The VANDELS survey: the stellar metallicities of star-forming galaxies at 2.5 < z < 5.0
F. Cullen R. J. McLure J. S. Dunlop S. Khochfar R. Davé R. Amorin M. Bolzonella A. C. Carnall M. Castellano A. Cimatti M. Cirasuolo G. Cresci J. P. U. Fynbo F. Fontanot A. Gargiulo B. Garilli L. Guaita N. Hathi P. Hibon F. Mannucci F. Marchi D. J. McLeod L. Pentericci L. Pozzetti A. E. Shapley M. Talia G. Zamorani




firstpage–lastpage 2016


Pseudo-electromagnetic fields in topological semimetals
    Dmitry I. Pikulin
    Accepted – . Received 2017 May 30
========================================================







We present the results of a study utilising ultra-deep, rest-frame UV, spectroscopy to quantify the relationship between stellar mass and stellar metallicity for 681 star-forming galaxies at 2.5<z<5.0 (⟨ z ⟩ = 3.5 ± 0.6) drawn from the VANDELS survey.
Via a comparison with high-resolution stellar population synthesis models, we determine stellar metallicities (Z_∗, here a proxy for the iron abundance) for a set of high signal-to-noise ratio composite spectra formed from subsamples selected by mass and redshift.
Across the stellar mass range 8.5 < log(⟨ M_∗⟩/M_⊙) < 10.2 we find a strong correlation between stellar metallicity (Z_∗/Z_⊙) and stellar mass, with stellar metallicity monotonically increasing from Z_∗/Z_⊙ < 0.09 at ⟨ M_∗⟩ = 3.2 × 10^8M_⊙ to  Z_∗/Z_⊙ = 0.27 at ⟨ M_∗⟩ = 1.7 × 10^10M_⊙.
In contrast, at a given stellar mass, we find no evidence for significant 
metallicity evolution across the redshift range of our sample.
However, comparing our results to the z=0 stellar mass-metallicity relation for star-forming galaxies, we find that the ⟨ z ⟩ = 3.5 relation is consistent with being shifted to lower metallicities by ≃ 0.6 dex at all stellar masses.
Contrasting our derived stellar metallicities with estimates of the gas-phase metallicities of galaxies at similar redshifts and stellar masses, we find evidence for enhanced O/Fe ratios in z ≳ 2.5 star-forming galaxies of the order (O/Fe) ≳ 1.8 × (O/Fe)_⊙.
Finally, by comparing our results to the predictions of three cosmological simulations, we find that the ⟨ z ⟩ = 3.5 stellar mass-metallicity relation is consistent with current predictions for how outflow strength scales with galaxy stellar mass.
This conclusion is supported by an analysis of one-zone analytic chemical evolution models, and suggests that the mass loading parameter (η=Ṁ_outflow/M_∗) scales as η∝ M_∗^β with β≃ -0.4.
 


galaxies: metallicity - galaxies: high redshift - galaxies: evolution - 
galaxies: star-forming






§ INTRODUCTION


The relationship between the stellar mass and metallicity of galaxies as a function of cosmic epoch provides a fundamental constraint on models of galaxy formation <cit.>.
The metallicity of a galaxy is determined by the integrated past history of star formation, the fraction of metal-enriched gas lost from the interstellar medium (ISM) through outflows, and the dilution of enriched ISM gas by pristine inflows from the intergalactic medium (IGM).
Metallicity scaling relations therefore provide rich information regarding some of the key physical processes governing the evolution of galaxies.

In recent years, significant efforts have been made towards determining the gas-phase metal content of galaxies at 2 ≲ z ≲ 4 from the ratios of strong optical nebular emission lines <cit.>.
Despite this, the substantial uncertainties inherent in converting nebular line ratios to abundances have frustrated efforts to reach a consensus.
In addition to the known inconsistency between different metallicity calibrations in the local Universe <cit.>, the problem is compounded at higher redshifts by the evolution in the physical conditions in  regions, potentially rendering local calibrations unusable <cit.>.
Disconcertingly, it is still unclear whether the correlation between nebular line ratios and galaxy stellar mass in high-redshift galaxies is dictated primarily by the gas-phase metallicity or is purely driven by the ionizing properties of the massive stellar population <cit.>.
To definitively address these concerns, direct estimates of gas-phase metallicity are needed at high-redshift, requiring the the detection of faint oxygen auroral lines such as  and , which are extremely challenging observations that have only been achieved a handful of times at z≳2 <cit.>.

An alternative method for probing the metal content of galaxies utilizes their stellar continuum emission directly.
Abundances derived in this way are referred to as stellar (Z_∗) as opposed to gas-phase (Z_g) metallicities.
When available, stellar metallicities are a useful independent probe of the metal content of galaxies. 
In the local Universe, several authors have investigated the stellar mass-metallicity relationship for large statistical galaxy samples, primarily from the extensive Sloan Digital Sky Survey (SDSS) dataset <cit.>, integral field spectroscopic (IFS) surveys <cit.> and, at lower stellar masses, from stellar spectroscopy in local dwarf galaxies <cit.>.
In the majority of these cases, metallicities have been derived from rest-frame optical continuum observations.
Recently, <cit.> presented a compilation of Z_∗ estimates in the local Universe spanning ≃ 7 orders of magnitude in stellar mass (ranging from M_∗≃ 10^4 M_⊙ to M_∗≃ 10^11M_⊙) finding evidence for a continuous relation which rises from  ≃ -2.5 at the lowest stellar masses, up to  ≃ 0.0, flattening at stellar masses above M_∗≃ 10^10M_⊙.
Determining whether or not a similar relation was already in place at earlier cosmic epochs is clearly of significant interest.

Unfortunately, estimates of Z_∗ for galaxies at z≳2 are rare, primarily because the measurement requires a high signal-to-noise ratio (S/N) detection of the stellar continuum, an expensive and challenging observation for faint high-redshift sources.
Nevertheless, estimates have been made for small samples of lensed galaxies at z ∼ 2-3 <cit.>, a number of unlensed sources at z > 3 <cit.> and a composite spectrum of star-forming galaxies at z≃ 2 <cit.>.
In all of these cases, the observations have been taken with ground-based optical spectrographs and therefore the stellar metallicities are based on rest-frame far-ultraviolet (FUV) spectra.
The distinction is important, because FUV-based metallicities are considered to be, to a first approximation, a measurement of the iron abundance in the photospheres of the young, massive, O- and B-type stars in the galaxy <cit.>.
This is not necessarily the case for optical-based stellar metallicities, which trace more evolved stars and sample longer star-formation timescales, and is certainly not the case for estimates of Z_g derived from optical nebular emission lines, which trace the young stellar population, but mainly the oxygen and nitrogen abundance.

Most commonly, stellar metallicities at high-redshift have been determined using a calibration of FUV photospheric absorption indices developed by <cit.>, and later extended in <cit.>.
However, the results from these studies have generally been inconclusive, especially with respect to the scaling relation between metallicity and stellar mass.
This has been, in part, a consequence the of small sample sizes and subsequent lack of dynamic range, but is also a consequence of the fact that these indices are also prone to significant contamination by interstellar medium (ISM) absorption <cit.>.
More recently, an alternative method for measuring Z_∗, by fitting the full FUV spectrum at λ≤ 2000 Å, has been demonstrated by <cit.>  using the BPASSv2.0 stellar population synthesis models including massive stellar binaries <cit.>.
The method is a natural extension of the previous methods, making more complete and consistent use of the available data and should, in principle, place more robust constraints on Z_∗ than the standard indices approach <cit.>.

Applying the full FUV spectral fitting method to a composite spectrum of 30 star-forming galaxies at z=2.4, <cit.> found a low FUV-based stellar metallicity of Z_∗/Z_⊙≃0.1, roughly a factor of five lower than the gas-phase metallicity derived from rest-frame optical nebular emission lines (Z_g/Z_⊙≃0.5), which they interpreted as evidence for super-solar oxygen-to-iron ratios (O/Fe) in high-redshift star-forming galaxies.
This result was subsequently confirmed using a larger statistical sample of 150 star-forming galaxies at z∼2-3 in <cit.> albeit via a slightly different method in which Z_∗ and Z_g were simultaneously estimated from fitting optical nebular emission lines.
However, although <cit.> established that star-forming galaxies at z≳2 have significantly sub-solar stellar metallicities (iron abundances), they did not explore the existence of a scaling relation of Z_∗ with stellar mass.
Therefore, determining the stellar mass-metallicity relation at z≳2 remains a key open question in the study of galaxy evolution.

Fortunately, further progress can be made in this area using the deep spectroscopic data taken as part of the VANDELS survey <cit.>.
VANDELS is an ESO public spectroscopic survey providing exceptionally deep integration times (of up to 80 hours) for ∼ 2000 galaxies at z ≳ 1.
At redshifts 2.5 < z < 5.0, the VANDELS spectra cover the rest-frame FUV, allowing stellar metallicities to be estimated using the methods described above.
Moreover, since all of the VANDELS targets are drawn from the CDFS and UDS survey fields, accurate stellar masses can be determined for all sources using the deep multi-wavelength photometric catalogues available in these regions.
In this paper, we present a study of the stellar mass-metallicity relationship at 2.5 < z < 5.0 using deep rest-frame FUV spectra from the VANDELS survey.

The structure of the paper is as follows.
In Section <ref> we describe the selection of our VANDELS star-forming galaxy sample along with the other relevant datasets used in this work.
In Section <ref> we begin by reviewing the metallicity information contained within rest-frame FUV galaxy spectra, and outline our method for determining FUV-based stellar metallicities.
In Section <ref> we present our stellar mass-metallicity relation and compare to the predictions from three independent cosmological simulations.
In Section <ref> we discuss some of the implications of our results, including an investigation into which physical parameters are potentially driving the observed stellar metallicities in our sample, and a discussion of O/Fe ratios in high-redshift star-forming galaxies.
Finally, in Section <ref> we summarize our results and conclusions.
Throughout this paper metallicities are quoted relative to the solar abundance taken from <cit.> which has a bulk composition by mass of Z_∗=0.0142, an iron mass fraction of Z_∗,Fe=0.0013 and an oxygen mass fraction of Z_∗,O=0.0058, and we assume the following cosmology: Ω_M =0.3, Ω_Λ =0.7, H_0 =70 km s^-1 Mpc^-1.
All magnitudes are quoted in the AB magnitude system <cit.>.

	







§ DATA


In this section we present the various observed and simulated datasets used in this work.
The dataset of primary interest, described in Section <ref>, is the rest-frame FUV spectra of star-forming galaxies at  2.5 < z < 5.0 provided by the VANDELS survey.
These high-redshift data are supplemented by a selection of FUV spectra of starburst galaxies and star-forming regions in the local Universe which are described in Section <ref>.
Finally, in Section <ref>, we discuss the data extracted from two state-of-the-art cosmological simulations <cit.> which are used to both test our method and to aid in the analysis and interpretation of our results.




 §.§ VANDELS



The primary spectroscopic data were obtained as part of the VANDELS ESO public spectroscopic survey <cit.>.
The VANDELS survey is a deep optical spectroscopic survey of the CANDELS CDFS and UDS fields with the VIMOS spectrograph on ESO's Very Large Telescope (VLT) targeting massive passive galaxies at 1.0 ≤ z ≤ 2.5, bright star-forming galaxies at 2.4 ≤ z ≤ 5.5 and fainter star-forming galaxies at 3.0 ≤ z ≤ 7.0.
In this work, we focus exclusively on the star-forming galaxies at z ≤ 5.

All galaxies were drawn from four independent H-band-selected catalogues.
The CDFS and UDS regions are covered by the CANDELS survey, and benefit from extensive WFC3/IR imaging <cit.>.
In these fields, photometry catalogues produced by the CANDELS team were used <cit.>.
Within the wider-field regions, two bespoke PSF-homogenized catalogues were produced, primarily from publicly available ground-based imaging <cit.>.
From these catalogues objects were pre-selected as potential spectroscopic targets based on robust photometric redshift estimates.
The bright star-forming galaxies were chosen to satisfy 2.4 ≤ z_phot≤ 5.0 with an i-band magnitude of i ≤ 25.
This i-band constraint was chosen to ensure the final spectra had sufficient signal-to-noise (S/N) to allow detailed analyses of individual objects.
The faint star-forming galaxies relevant to our study (i.e. in the redshift interval 3.0 ≤ z_phot≤ 5.5) were selected to have 25 ≤ H ≤ 27  i ≤ 27.5 in the HST regions [In the wider-field regions the limiting i-band magnitude is slightly brighter at i ≤ 26.0.].  
As discussed in <cit.> (illustrated in their Fig. 3), this selection ensured that the sample of star-forming galaxies is consistent with being drawn from the main sequence of star-forming galaxies at all redshifts.
An additional constraint, enforced by the observing strategy, was placed on the i-band magnitudes so that slits were efficiently allocated to objects requiring 20, 40 and 80hr of integration time in a roughly 1:2:1 ratio <cit.>.
The final spectroscopic sample of faint star-forming galaxies was an unbiased random (approximately 1 in 4) subsample of the pre-selection catalogue.
    
The observations and reduction of the VIMOS spectra are described in detail in the first data release paper <cit.>.
Briefly, observations were obtained with the ESO-VLT VIMOS spectrograph using the medium resolution grism which covers the wavelength range 4800 < λ_obs < 10000 Å with a resolution of R=580 and a dispersion of 2.5 Å per pixel.
The median seeing across all observations was FWHM=0.7.
The spectra were reduced using the Easy-Life pipeline <cit.> which is an update of the algorithms originally described in <cit.>.
Easy-Life performs all standard data reduction procedures including the removal of bad CCD pixels and sky subtraction at the individual exposure level.
The 2D spectrum for each individual exposure is extracted and resampled onto a common linear wavelength scale; these individual-exposure 2D spectra are then combined to create the final 2D spectrum for each object.
The 1D spectra are extracted, flux calibrated using the spectrophotometric standard stars and corrected for telluric absorption.
A final correction for atmospheric and galactic extinction is applied and the 1D spectra are normalized to the i-band photometry.
     
     

We draw our sample from the second VANDELS data release[Data available through the consortium database at the link http://vandels.inaf.it/db] which comprises roughly 65% (1362/2106) of the full survey.
Redshifts for all of the spectra have been determined by members of the VANDELS team and assigned a redshift quality flag as described in <cit.>.
We only consider galaxies at 2.3 ≤ z_spec≤ 5.0 and with a redshift quality flag 3 or 4 (corresponding to a ≥95% probability of being correct), leaving a total of 681 galaxies.
Fig. <ref> shows the distribution of our final sample in the mass-redshift plane compared to the underlying distribution of all 2.4 ≤ z_phot≤ 5.0 star-forming galaxies in the pre-selection catalogue.
We have extended our selection to z_spec≥ 2.3 so that we include galaxies whose photometric redshift was slightly overestimated at the lower redshift boundary of the star-forming galaxy selection (z_phot=2.4).
The upper redshift limit of z=5 is enforced primarily because above this redshift the FUV wavelength coverage becomes prohibitively small.
The average redshift of our final sample is ⟨ z ⟩ = 3.5 ± 0.6.
Stellar masses for the full VANDELS sample were derived from SED fitting using <cit.> templates with solar metallicity and assuming exponentially declining star-formation histories as described in <cit.>.
We also derived stellar masses for our 2.5 < z < 5.0 sample using fast++, a rewrite of fast <cit.> described in <cit.>.
We adopted the same set of parameters as <cit.>, but allowed for rising star-formation histories using a delayed exponentially declining model (M_∗∝ te^t/τ) and a range of metallicities values (0.3-2.5 ×Z_⊙).
The stellar masses returned from these two procedures were fully consistent, but for this work we adopt the fast++ values.

As discussed in <cit.>, the pre-selection catalogue contains star-forming galaxies consistent with being drawn from the star-forming main sequence across the full redshift range.
Therefore, since the final selection is a random subsample of galaxies from the pre-selection catalogue, the galaxies in our final sample can be considered typical star-forming galaxies for their stellar mass.
The stellar mass range of the final sample is 8.2 < <11.0 with a median specific star-formation rate of sSFR=4.4 ^-1.




  §.§.§ Composite FUV spectra


The typical continuum S/N is too low to extract metallicity information, therefore we focus our analysis on stacked spectra in bins of stellar mass and redshift.
The mass and redshift bins are illustrated in Fig. <ref> and detailed in Table <ref>.
In total we chose four redshift bins which were selected to encompass roughly equivalent intervals of cosmic time (≈ 400 Myr), with three stellar-mass bins per redshift.
We also analyse seven composite spectra stacked in bins of stellar mass across the full redshift range, with the stellar mass ranges given in Table <ref>.
The mass bins were chosen manually with the aim of keeping the bin widths as narrow as possible whilst also ensuring that they contained enough galaxies to have an acceptable S/N.

The final stacked composite spectra were formed by first shifting each contributing spectrum into the rest-frame using z_spec. 
To correct for redshift differences, the flux of each spectrum was scaled to the flux that would be observed at the mean redshift of the stack.
For the stellar-mass redshift stacks these flux-correction factors are all ≲ 5 %; for the stellar-mass only stacks the maximum average correction is 15±7% which occurs in the highest redshift bin.
We have confirmed that our results do not change significantly if, instead of this flux-correction method, the spectra are simply normalized and then combined.
The primary benefit of preserving the absolute flux values is that an error spectrum can be robustly estimated.
The individual spectra were then resampled onto a common wavelength grid, which varied slightly depending on the redshift range (see Table <ref>), with a dispersion of 1 Å per pixel. 
The final flux at each dispersion point was taken as the median of all the individual flux values after rejecting 3σ outliers, and the 1σ error was calculated from bootstrap re-sampling of the individual flux values.
An example of one composite spectrum is shown in Fig. <ref>.
The effective spectral resolution element of the composites is 3.0 Å.




 §.§ FOS-GHRS Local Sample


The main high-redshift sample was supplemented by a sample of local galaxies for which published rest-frame FUV spectra are available.
The motivation for including a local reference sample was that the majority of published stellar metallicities in the local Universe are derived from rest-frame optical spectra <cit.>, and it is not clear that optical and FUV light will originate from the same population of stars or trace the same abundance type.
Even in the cases where rest-frame UV spectra are analysed, the focus is primarily on deriving gas-phase metallicities (not stellar metallicities) from the depths of ISM absorption features <cit.>.
To our knowledge, there are no local examples of stellar metallicities derived solely from global FUV continuum fitting.
For this reason, we have constructed a comparison sample of local starbursts and star-forming galaxies from the FUV spectroscopic atlas presented in <cit.> which is comprised of spectra taken with the Faint Object Spectrograph (FOS) and Goddard High Resolution Spectrograph (GHRS).

Our sample is drawn from the 46 rest-frame UV spectra observed with the Faint Object Spectrograph (FOS) and the Goddard High Resolution Spectrograph (GHRS) on-board the HST presented in <cit.>.
These 46 spectra are drawn from 28 individual galaxies (i.e. some spectra are just different regions of the same galaxy).
The spectral resolution of the individual spectra spans the range 0.5-3Å depending on the instrument, aperture size and physical extent of the object being observed.
For consistency with the the VANDELS data, we smooth all spectra to a common 3Å resolution.
We also require that the spectra have wavelength coverage in the interval 1410 ≤λ≤ 1450 Å to enable an analysis of normalized composite spectra, and finally that the corresponding galaxy has a measurement of absolute K-band magnitude that we can use to estimate the stellar mass using the relation from <cit.>.
This selection leaves 26/46 of the original spectra from 18/28 of the original local starbursts and star-forming galaxies presented in <cit.>.
A list of the individual spectra used in this work is given in the Appendix (Table <ref>). 
Composite spectra were constructed in three bins of stellar mass as outlined in Table <ref>.
To create the composites, all individual spectra were median combined on the same wavelength grid used for the VANDELS galaxies.
In this case, an error spectrum was estimated by propagating the error spectra of each of the individual spectra (bootstrapping could not be applied in this case since the individual spectra were normalized).



 §.§ Simulation Data


To compare the observations to simulation predictions we have extracted data from two state-of-the-art cosmological hydrodynamical simulations:  <cit.> and  <cit.>.
In Section <ref> we use the simulation data to assess the accuracy our adopted method for measuring stellar metallicities and in Section <ref> we compare the observed stellar mass-metallicity relationship to the simulation predictions.
Having two independent simulations allows us to account for systematic variations between their methodologies and theoretical predictions.
Here, we give a brief overview of the simulation details but refer the reader to the references above for further information.



  §.§.§ FiBY


The  simulation suite is a set of high-resolution cosmological hydrodynamical simulations using a modified version of the GADGET code used in the Overwhelmingly Large Simulations (OWLS) project <cit.>. 
The code tracks metal pollution for 11 elements: H, He, C, N, O, Ne, Mg, Si, S, Ca and Fe and calculates the cooling of gas based on line-cooling in photoionization equilibrium for these elements <cit.> using tables pre-calculated with CLOUDY v07.02 <cit.>. 
Furthermore, the simulation incorporates full non-equilibrium primordial chemistry networks <cit.> including molecular cooling functions for H_2 and HD. 
Star formation is modelled using the pressure law implementation of <cit.>, which yields results consistent with the Schmidt-Kennicutt law <cit.>. 
The simulations include feedback from stars by injecting thermal energy into the neighboring particles <cit.>. 
Element yields from Type Ia and core-collapse supernova (CCSNe) are implemented following the prescription of <cit.>.
In this work we will focus on the _XL simulation which has individual gas and star particle masses of  = 5.68 and covers a co-moving volume of (32  Mpc)^3.
The lowest redshift simulated in the  simulation is z=4, and this is the closest redshift to the mean redshift of our sample (⟨ z ⟩=3.5).
Therefore, for our  comparison sample, we extracted the 591 galaxies in _XL with  > 8.0 at z=4.
The maximum galaxy stellar mass is M_∗=3.8 × 10^10 M_⊙, with a median value of 4.0 × 10^8 M_⊙.



  §.§.§ Simba


The  simulation suite is based on the GIZMO cosmological gravity plus hydrodynamics solver <cit.> and is described in detail in <cit.>.
The code implements a H_2-based star formation rate, with H_2 fractions calculated using the sub-grid model of <cit.>.
The chemical enrichment model tracks the same 11 elements as in , with radiative cooling and photoionization heating modeled using the Grackle-3.1 library <cit.>, and element yields calculated for CCSNe, Type-Ia SNe and Asymptotic Giant Branch stars following the prescriptions of <cit.>, <cit.> and <cit.>, respectively.
 employs a mass outflow rate scaling with galaxy stellar mass motivated by the results derived from the extremely high resolution FIRE zoom simulations <cit.>.
In this work we focus on the data from the m50n1024 simulation which has a co-moving box length of 50 h^-1Mpc and individual gas/star element resolution of 2.28× 10^6 M_⊙.
For consistency with , we focus on the z=4 snapshot in the simulation.
From this, we extract a sample of 1749 galaxies down to a minimum stellar mass of 4.5 × 10^8 M_⊙, with a maximum stellar mass of 1.7 × 10^11 M_⊙ and a median value of 8.2 × 10^8 M_⊙.

    




§ MEASURING STELLAR METALLICITIES


We now describe our method for estimating stellar metallicities from rest-frame FUV spectra.
In Section <ref> we describe the metallicity information content of the rest-frame FUV spectra of star-forming galaxies.
The stellar population synthesis models we use to compare to the data are discussed in Section <ref>.
Then, in Section <ref>, we describe the statistical method we employ to robustly constrain the metallicity.
Finally, in Section <ref> we discuss various tests of our method. 



 §.§ Stellar Metallicity Indicators in the UV


To begin with, it is worthwhile reviewing what metallicity information is contained within the FUV spectrum of a star-forming galaxy.
Figs. <ref> and <ref> serve as a useful references for this discussion.
From simple photon energy arguments, we know that the rest-frame FUV is dominated by the light from young, massive, O- and B-type stars.
These stars emit a continuum which contains absorption (and potentially emission) features due to elements within the stellar photosphere and the expanding stellar wind. 
The strength of these features is naturally a strong function of the total photospheric metallicity <cit.>.

The most prominent absorption features are typically the strong stellar wind lines produced from the stellar atmospheres, namely ,  and .
These lines are often extremely broad, blueshifted by ∼ 1000  kms^-1 and exhibit P-Cygni line profiles. 
They are useful metallicity diagnostics since the strength of the stellar winds, and hence the line profiles, have a strong metallicity dependence which is linked to the metallicity-dependent mass-loss rates <cit.>.
Naturally, however, these lines are also sensitive to the IMF and star-formation history.

Outside these strong wind features, the rest-frame FUV shortward of λ=2000Å contains an abundance of stellar photospheric absorption features, primarily due to transitions of highly ionized iron <cit.>.
The individual absorption features have low equivalent widths which require high spectral resolution and S/N to be seen individually, but the resulting strong photospheric line blanketing can be seen in unresolved O- and B-star populations in local starburst galaxies, star-forming regions and even high-redshift galaxies <cit.>.

    
    

These photospheric absorption features have been utilized in the past as metallicity indicators, for example <cit.> developed the `1978 index' at 1935-2050Å, a wavelength regime which is believed to be dominated by  transitions in early B-type stars. 
The 1978 index, which is thought to be a direct probe of the iron abundance ([Fe/H]), was used to derive the stellar metallicity of a composite spectrum of 75 star-forming galaxies at z≈2 by <cit.>.
Additional wavelength regions (e.g. , , , ) have been investigated by other authors and shown to be similarly sensitive to the photospheric abundance <cit.>.
These individual absorption features are analogous to the Lick index system developed for constraining abundances from rest-frame optical spectra <cit.>.
However, <cit.> have recently demonstrated that this approach is limited by the fact that the many of these indices are significantly contaminated by ISM absorption lines.

An alternative approach, demonstrated most recently by <cit.>, is to simply fit the full FUV spectra directly, thereby constraining the FUV abundance by using all of the metallicity-sensitive regions simultaneously.
Indeed, similar approaches have proven successful for constraining stellar metallicities at optical wavelengths <cit.>.
The obvious merit of this approach is that potentially relevant information is not being thrown away <cit.>.
On the other hand, the fact that the absorption features across the full FUV spectrum are a result of a number of different element species means there is some ambiguity as to what abundance is actually being measured.
<cit.> (among others) have pointed out that, because the vast majority of the FUV line-blanketing is a result of transitions in highly-ionized iron, a metallicity derived by fitting the FUV spectrum is, to a first approximation, a measure of the iron abundance in the stellar photosphere.
The point is demonstrated in Table 1 of <cit.>, which provides a comprehensive list of the photospheric spectral lines identified in star-forming galaxies both locally and from lensed high-redshift sources <cit.>.
In this work we will adopt a full spectra fitting approach similar to <cit.>, which is outlined in detail in Section <ref>.
Throughout this paper we refer to the stellar metallicity derived from rest-frame FUV observations as Z_∗/Z_⊙ (or ), but emphasize that this should be understood as a proxy for the iron abundance. 

A final important point to note is that, while these FUV absorption features are strongly dependent on metallicity, the low-order shape of the FUV continuum is not (in contrast to the case at optical wavelengths).
By far the dominant factor in determining the FUV continuum shape is the wavelength-dependent dust attenuation law.
Interestingly, the attenuation law derived from local starbursts by <cit.> still appears to be a good approximation for the average shape of the attenuation law for star-forming galaxies at high-redshift with  ≳ 9.5 <cit.> , although this is still a matter of debate <cit.>, and object by object variation, and even some stellar mass dependence, is expected <cit.>.
Crucially, however, the lack of metallicity dependence on the global continuum shape (i.e the overall curvature of the FUV spectrum) negates any strong degeneracies between the metallicity and the wavelength-dependence/normalisation of the dust attenuation.



 §.§ Stellar Population Synthesis Models


To derive metallicities from the observed FUV spectra we require models to compare to the data.
A number of different stellar population synthesis models exist <cit.>, but in this work  we have opted to use the Starburst99 (SB99) high-resolution WM-Basic theoretical stellar library described in <cit.>.
Our choice was motivated partly for clarity, and also by the fact that the WM-Basic SB99 spectra are the highest resolution FUV spectral models available, and have been extensively tested against individual spectra of hot stars in the Galaxy and Magellanic Clouds.
We argue that the WM-Basic SB99 models still represent the most robust stellar population models available for comparing to spectroscopic data in the FUV.
However, the SB99 models do not account for some physical effects that are known to be important in massive star evolution, such as the prevalence of binary stars and the phenomenon of `quasi-homogeneous evolution' at low metallicities <cit.>.
In Appendix <ref> we have demonstrated that using stellar populations models which do account for these phenomena <cit.> would not affect our main results.

We considered constant star-formation rate models from the latest version of SB99 <cit.>, assuming the weaker-wind Geneva tracks with stellar rotation and single-star evolution.
We considered the Geneva evolution tracks at Z_∗=(0.001, 0.002, 0.008, 0.014, 0.040), and an η=2.3 IMF slope in the mass range  corresponding to the default <cit.> IMF.
The high-resolution WM-basic spectra provided by Starbust99 cover the wavelength region 900-3000 Å at a dispersion of 0.4 Å/pixel, so we smoothed and resampled the models to match the resolution (≈ 3 Å) and sampling (1 Å/pixel) of the VANDELS composites (Fig. <ref>).
We considered constant star-formation models over timescales of 100, 300 and 500  but adopt the 100  models as our fiducial set (we discuss the motivation for this choice below).

Nebular continuum emission was added to all models using cloudy v17.00 <cit.> following the method described in <cit.>, assuming  region parameters consistent with the current best estimates at z≃2-3 <cit.>.
When fitting the observed data we considered models with both maximal nebular contribution (f_esc=0%) and zero nebular contribution (f_esc=100%).
However, since the average escape fraction is measured to be relatively low in high-redshift star-forming galaxies <cit.>, we adopt the models including maximal nebular continuum as our fiducial set.
We also note that, since the nebular continuum only acts to alter the shape of the continuum it is degenerate with the dust prescription, not with the stellar metallicity.
Assuming f_esc=100% does not change the results presented here.

Finally, in Section <ref>, we will discuss how we tested our FUV-fitting method using synthetic FUV spectra derived from  and  simulation data.
For generating these synthetic spectra we also used WM-Basic SB99 models, but in this case considered the instantaneous burst models generated at 1000 time-steps split logarithmically between 1  and 1.5 .
These instantaneous burst models are required so that the synthetic spectra can be constructed to match the known star-formation histories and chemical-abundance histories of the simulations.
The instantaneous burst models assume a total stellar mass within the burst of 10^6 M_⊙ with all other parameters (e.g. IMF, metallicity range) being identical to the constant star-formation rate models.
As will be described in more detail below, the synthetic spectra were constructed from these instantaneous burst models using the average star-formation and chemical-abundance histories of the simulated galaxies following the method outlined in <cit.>.



 §.§ A statistical estimate of the stellar metallicity



To fit the SB99 models to our data we adopted a Bayesian forward-modeling approach.
The method follows the familiar approach of model fitting using Bayes' theorem with,

    P(θ|D)∝ P(D|θ)P(θ),

where P(θ|D) is the posterior probability on the input parameters θ given the data (D), P(D|θ) is the likelihood function (L), and P(θ) is the prior. 
Assuming the error bars are Gaussian and independent, the logarithm of the likelihood function is given by,

    ln(L) = -1/2∑_i[(f_i-f(θ)_i)^2/σ_i^2 + ln(2πσ_i^2)]

where f is the observed flux, f(θ) is the model flux for a given set of parameters θ, and σ is the error on the observed flux.
The summation is over all wavelength pixels included in the fit.
For our model we adopted four free parameters: the logarithm of the stellar metallicity (log(Z_∗/Z_⊙)) and three parameters used to fit the overall continuum shape, based on a physically motivated parameterization of the dust attenuation law.

    

The dust law parameterization, taken from <cit.>, is a modification of the <cit.> attenuation law for starburst galaxies which allows for the slope of the curve to be modified and includes a prescription for the UV bump at 2175Å. 
It is given by:

    A_λ,mod=A_V/R_V,mod[k_λ,CalzR_V,mod/R_V,Calz(λ/5500Å)^δ+D_λ],

where A_V is the absolute V-band attenuation (i.e. the normalization of the attenuation law), k_λ,Calz is the total-to-selective attenuation curve for the <cit.> law and δ is the power law exponent used to modify the slope of the <cit.> law.
R_V,mod is the modified total-to-selective attenuation ratio which is simply a function of δ <cit.> and R_V,Calz=4.05.
Finally, D(λ) is the Drude profile which is commonly used as the functional form of the UV bump at 2175Å and is given by:

    D_λ(B)=Bλ^2Δ^2/[λ^2-λ_c^2]^2-λ^2Δ^2

where λ_c is the central wavelength of the feature (2175Å), Δ is the width (350Å) which are both held fixed, and B is the amplitude, which we allowed to vary[Although the peak of the UV bump at 2175Å. is outside the range of our fitting it can still potentially affect the shape of the UV continuum due to its width (350Å)].
This dust attenuation prescription therefore has three free parameters: A_V (normalisation), δ (slope) and B (UV bump strength).
When δ=0 and B=0 the parameterization simply becomes the <cit.> law.
We note that the role of the three dust parameters is to fit the shape of the continuum, which, as discussed above, has little dependence on the stellar metallicity (Fig. <ref>). 

To perform the fitting, we used the nested sampling code multinest <cit.>[We accessed multinest via the python interface pymultinest <cit.>.], which is an implementation of the nested sampling algorithm described in <cit.>.
Nested sampling is an alternative to the traditional Markov Chain Monte Carlo (MCMC) method of sampling the posterior distribution for a given Bayesian inference problem, which enables the extraction of 1D posterior distributions for any given parameter in the model by marginalizing over all other free parameters.
We adopted simple flat priors on all the model parameters.
The prior in  is imposed by the SB99 models to be -1.15 <  < 1.45. 
As the models are only provided at five fixed metallicity values, we linearly interpolated the logarithmic flux values between the models to generate a model at any metallicity within the prescribed range.
We experimented with other interpolation schemes (e.g. interpolating in flux linear flux values) but found that this did not strongly affect our results, and performed simple tests (e.g. recovering the known metallicity of the Starburst99 templates) to ensure the interpolation scheme was not strongly biasing the recovered metallicities.
For A_V and B we considered values in the range 0 - 5 and for δ we considered -1.0 < δ < 1.0.

A final point to note is that not all wavelength pixels were included in the fitting.
As discussed above (and illustrated in Fig. <ref>) a significant portion of the rest-frame FUV spectra is contaminated by features not related to stellar emission (e.g. interstellar absorption lines and nebular emission lines).
We excluded these wavelength regions, using only the spectral windows sensitive to photospheric absorption and stellar-wind features outlined in Table 3 of <cit.>.
Finally, we also only considered wavelengths redward of  at 1216Å to avoid the additional uncertainties related to the intergalactic medium  absorption and dust correction at shorter wavelengths.
An example of a fit to one spectrum (VANDELS-m5) is shown in Fig. <ref>. 
We list the best-fitting log(Z_∗) values and errors in the last column of Table 1.
These quoted errors are the statistical uncertainties only, a discussion of the various systematic uncertainties, which we find to be at the ≃ 10 % level, is given in the appendix.




 §.§ Accuracy of the estimate and systematic uncertainties


Before discussing the stellar metallicities derived from the VANDELS data, it is worth considering how reliable we expect our method to be.
We performed two tests to explore any potential systematic uncertainties related to our adoption of a simplified star-formation and chemical-abundance history.
Firstly, we used data from the  and  simulations to construct synthetic FUV spectra based on realistic star formation and chemical-abundance histories, and compared our recovered  values to the true FUV-weighted stellar metallicities.
Secondly, we derived stellar metallicities from the FUV spectra of local star-forming regions and starburst galaxies, which we compared to the published gas-phase metallicities of the same objects.
The results of both of these independent tests are discussed below.

    



  §.§.§ Synthetic spectra based on simulation data


In previous analyses of the FUV spectra of star-forming galaxies at high redshift, it has been common to assume models with constant star-formation histories <cit.>.
Generally, a star-formation timescale of 100  Myr is assumed.
It is further assumed that the metallicity of the stars does not vary strongly across the star-formation history.
These assumptions make modeling the FUV spectra of high-redshift galaxies relatively simple compared to longer-wavelength regimes and/or lower-redshift galaxies, where the timescales involved are longer, so that more complex star formation and chemical evolution histories must be considered <cit.>.

Nevertheless, it is worth testing the validity of these assumptions against predictions from simulations.
To do this we have used the  and  simulation data described in Section <ref>.
Because our composite spectra are averages across a population, we focused on the average star formation and chemical abundances histories of galaxies in the  and  simulations. 
We considered two stellar mass bins which encompass the observed data, a low mass bin with 8.5 <  < 9.5 and a high mass bin with  > 9.5.
As a result of their different volumes the maximum mass is different in the two simulations (10^10.5M_⊙ in  and 10^11.2M_⊙ in ).
The average star formation and chemical evolution histories for these simulated galaxies are shown in Fig. <ref>.

     

In accordance with common assumptions, it can be seen that, although the average star formation rates are smoothly rising over the full formation history, all are relatively constant within the past 100 .
Even in the case of the most extreme variation (the high-mass bin of  galaxies) the 100  averaged star formation rate is 29±4 M_⊙yr^-1 (ranging from 23 - 35 M_⊙yr^-1).
In fact, for the star formation rates, one could argue that there is relatively little variation over roughly 300-500  timescales in all cases.
This is in good agreement with previous analyses of the average star formation histories of simulated galaxies at similar redshifts <cit.>.

An invariance on timescales of ≃ 100  is also true for the chemical abundances, which overall follow a very similar trend to the star formation rates.
There is, in general, a slightly steeper rise over the longer timescales, although even in the most extreme case (again the high-mass bin of  galaxies) the increase is only a factor ≃ 2 over 500 .
Interestingly, despite the exact details of the global histories varying quite significantly, all share a common theme: rapid evolution over the first  (i.e. from formation up to z≃6) followed by a more gradual evolution over the next 500  (i.e. z=4-6).
Based on the simulation data it would seem that assuming a simple constant star formation rate and chemical abundance model over 100  timescales is a reasonable approximation for the composite galaxy spectra in our sample, since, across 100 , both the average star formation rates and chemical abundances of the galaxy population are relatively invariant in the two stellar mass bins.

However, below we will also demonstrate explicitly that our results should not be strongly affected by assuming longer timescales.
This is mainly due to the fact that, assuming constant star formation, the model FUV spectra reach an equilibrium (i.e. become time invariant) after ≃ 50  <cit.>, and therefore any timescale ≳ 50  should yield similar results.
We also note that while these assumptions hold for galaxies on average, the same might not be true for individual galaxies.
On an individual basis, the star formation histories will not necessarily be smooth, and bursts of star formation will likely play a more significant role <cit.>.

As an explicit test, we constructed four synthetic FUV spectra based on the four star-formation and chemical-abundance histories illustrated in Fig. <ref>, using the SB99 WM-Basic instantaneous burst models scaled appropriately to match the integrated stellar mass.
To account for the effect of dust, we also attenuated each individual burst model according to its metallicity. 
Under our simple prescription, the absolute attenuation at 1500 Å is related to the metallicity via A_1500 = -2.24 ×log(Z_∗) - 2.16, with a random scatter of σ_A_1500=0.5 magnitudes (in this case Z_∗ is the total metallicity not the iron abundance).
This scaling relation is motivated by results from the FiBY simulation (Khochfar et al. in preparation), and implies ≃ 2 magnitudes of FUV attenuation for stars formed within solar metallicity environments.
We note that, despite its simplicity, the relation yields global (i.e. galaxy averaged) A_1500 values in reasonable agreement with the the attenuation versus stellar mass relation presented in <cit.>.
Finally, for simplicity, we used a Calzetti attenuation curve to convert A_1500 into A_λ across the full FUV spectral range.
The synthetic spectra were smoothed to the resolution of the VANDELS data and gaussian noise was added assuming a S/N per pixel of 15 (comparable to the observed stacks).

One unfortunate aspect of this procedure is that, for a portion of all the star-formation histories, the metallicity is below the minimum imposed by the Starburst99 models (Z_∗/Z_⊙ = 0.07; log(Z_∗/Z_⊙) = -1.15).
In these cases we are forced to simply use the lowest metallicity SB99 spectrum.
This is particularly a problem for the low metallicities predicted in the low-mass bin of the  simulation (the black solid curve in the right-hand panel of Fig. <ref>).
For the other three cases the issue is somewhat mitigated, since the FUV light is dominated by stars with metallicities within the SB99 parameter space.
Nevertheless, we note that this effect could still introduce some small bias into the recovered  values.

The results of the test are shown in Fig. <ref>, and generally confirm the conclusions drawn from Fig. <ref>.
For the three cases in which the FUV-weighted metallicity is within the model parameter space the true value is well recovered, with derived (true)  values of -0.78 ± 0.06 (0.78), -0.63 ± 0.07 (0.62) and -0.80 ± 0.06 (0.83). 
As mentioned above, the low mass galaxies in the  simulation have an average FUV-weighted metallicity of log(Z_∗/Z_⊙)=-1.88, below the lower limit of the models.
In this case, obviously somewhat by design, the solution bumps up against the edge of the parameter space.
For a similar case within the observed sample we would consider this an upper limit on .
It can also be seen that adopting constant star-formation rate models with longer timescales (300 and 500 ) does not have a strong systematic affect on the recovered metallicity although, in general, the 100  Myr models perform best.

Finally, we note that, although it is encouraging, this test cannot provide a definitive justification of our method, since we are building and fitting the FUV spectra with, fundamentally, the same set of SB99 models. 
However, it does serve to demonstrate that our simplified assumptions regarding the star formation and chemical-abundance histories should not significantly bias the recovered metallicity values.



  §.§.§ The FOS-GHRS local Universe sample


    

       


Given that the test described above, using synthetic spectra, is not truly independent, it is desirable to perform an independent test using real data.
To this end, we have used the above method to fit both the individual and composite spectra from the local FOS-GHRS galaxy sample described in Section <ref>.
The idea here is to compare the measured UV stellar metallicities to the published gas-phase metallicities of the same galaxies, under the assumption that the two parameters should correlate.
This should be a reasonable assumption, since the nebular emission lines used to determine the gas-phase metallicity are expected to originate from the hot, ionized, gas surrounding the same young, massive O- and B-type stars which dominate the FUV stellar emission.
Indeed, such a correlation is commonly observed in the  regions of Local Group galaxies, and is consistent with a 1:1 relation <cit.>.
Moreover, for the integrated emission of local galaxies (out to z≃0.3), a 1:1 correlation is observed for the young (< 2  old) stellar population <cit.>.

As discussed in Section <ref>, we produced three composite FOS-GHRS spectra in bins of absolute K-band magnitude (M_K, a proxy for the stellar mass) as outlined in Table <ref>.
Spectral fits to the M_K composite spectra are shown in Fig. <ref>.
We find that the stellar metallicities of the two brightest composites (M_K= -23.1 and -25.1; or log(M_∗/M_⊙)=10.33 and 11.13) are similar with  ≃ 0.2 (or Z/Z_⊙≃1.6) in excellent agreement with their median gas-phase metallicities.
For the faint, low-mass, composite with M_K= -18.5 ( =8.5) we find a stellar metallicity of  ≃ -0.58 (or Z/Z_⊙≃0.26), roughly a factor six lower) and, again, the stellar metallicity is in good agreement with the median gas-phase metallicity.
We note that the fits to the low-mass and high-mass composites, while generally good, do not appear to match the observed spectra in the range ≃ 1200-1400 Å despite having statistically acceptable reduced χ-squared values (χ^2_r ≃ 0.7).
These fits would likely be improved with a more complex dust model, as the current discrepancy could plausibly be a result of strong variation in the dust law at shorter wavelengths due to the relatively low number of objects in each composite. 
However, as the current fits are statistically acceptable, we did not consider a more complex dust model in this analysis.

We find a similar result for the individual galaxies, which are fitted in the same way (although often using only a restricted portion of the full wavelength range), albeit with a much larger scatter.
In fact, the results in Fig. <ref> are consistent with the stellar and gas-phase metallicities scattering about a 1:1 relation.
The Spearman rank correlation coefficient between the two quantities is 0.70, with the probability of obtaining this value by chance formally zero (p=3.1 × 10^-4).
We note that although a correlation between stellar and gas-phase abundances is to be expected, a direct 1:1 correspondence is only expected if the O/Fe ratios (i.e. chemical abundance patterns) in these galaxies are close to the solar value, we will return to this issue in the case of high-redshift galaxies in Section <ref>.

In summary, we have performed two independent tests of our adopted method.
Firstly, we have shown that the assumption of a constant star-formation history with a timescale of 100 , at a single metallicity, is sufficient to recover the FUV-weighted stellar metallicities of synthetic galaxy spectra built from more realistic star formation and chemical-abundance histories.
Furthermore, by applying our method to composite FUV spectra of local galaxies we have shown that we are capable of reproducing the expected trend between stellar and gas-phase abundances.
One further potential source of bias which we do not discuss in detail here, relates to the effect of redshift uncertainties on the metallicity-sensitive continuum features when constructing composite spectra. 
However, we do not expect this effect to significantly affect our derived metallicities, and a brief discussion is provided in Appendix <ref>.

 	




§ THE STELLAR METALLICITIES OF STAR-FORMING GALAXIES AT 2.5 < 𝐳 < 5.0


In this section we discuss stellar mass-metallicity relationship at 2.5 < z < 5.0 derived from the VANDELS composite spectra.
Before discussing the observational results however, it is important to consider that fact that the stellar metallicities derived using our method are FUV-weighted, and it is therefore worth briefly discussing how FUV-based metallicities (tracing the recent star-formation history) might be biased with respect to the mass-weighted metallicities (which trace the global star-formation history) at these redshifts.  

    



 §.§ FUV weighted versus stellar-mass weighted metallicities


From the  and  simulation data, we were able to derive both the mass-weighted and FUV-weighted stellar metallicities for each simulated galaxy and we present the results of the comparison in Fig. <ref>, which shows the median mass-weighted and FUV-weighted stellar mass-metallicity relations for both simulations.
The mass-weighted metallicities of each simulated galaxy were derived simply using the mass and metallicity information of each star particle.
For the FUV-weighted metallicities we used the FUV luminosity of each star particle at 1500 Å taken from the Starburst99 instantaneous burst models.
In both cases, the FUV-weighted relation is offset to lower metallicity, although crucially the shape of the relation is unaffected.
The offset is due to the fact that the youngest stars, which dominate the FUV flux, are generally more metal-enriched and therefore, according to our attenuation prescription, are more heavily attenuated and do not contribute as much to the FUV.
For example, it can be seen from Fig. <ref> that the intrinsic FUV-weighted relation (i.e. assuming no dust attenuation) for the  galaxies  (red-dotted line), closely follows the mass-weighted relation.
Therefore, any preferential attenuation of the youngest, most metal-enriched, stars will result in a bias to lower metallicities when using FUV spectra.
There is clearly some systematic uncertainty here related to our method of attenuating the individual star-particle spectra, however, based on reasonable dust assumptions, we expect the FUV-weighted metallicities of high-redshift star-forming galaxies to be biased low with respect to the mass-weighted values by ≃ 0.05-0.1 dex, but predict that the shape of the mass-metallicity relation to be similar in either case.
Finally, we note that throughout this section we will also compare our results to the relations extracted from the FIRE simulation presented in <cit.>[<cit.> assume a solar metallicity of 0.02 and therefore we have converted their relations to our assumed value of 0.0142.].
Since these FIRE relations are for mass-weighted metallicities, we have scaled them down by 0.05 dex.



 §.§ The 2.5 < 𝐳 < 5.0 stellar mass-metallicity relation


The 2.5 < z < 5.0 stellar mass-metallicity relation derived from the VANDELS sample is shown in Fig. <ref>.
The vertical error bars show the 16th and 84th percentiles of the posterior distribution on  and the horizontal error bars show the range of stellar mass in each bin.
The first thing to note is that a clear monotonic increase in  with  is observed for our sample.
Focusing on the purely stellar-mass stacks,  evolves from a lower limit of < -1.04 (68%) at  ≃ 8.5 to -0.57 ± 0.06 at  ≃ 10.2, an increase by a factor ≥3 in the average metallicity across roughly two decades of stellar mass. 
In Fig. <ref> we show how this metallicity evolution affects the form of the FUV spectrum by showing the best-fitting SB99 models overlaid on the composite spectra (the effect is subtle, but as the average mass of the composite spectra increases, the continuum becomes less smooth and the undulations resulting from metal line blanketing become more pronounced).
At all masses, the stellar metallicities (to first order a proxy for the iron abundance) are significantly sub-solar (≲ 25%). 
From an analysis of the redshift and mass tacks, we do not find strong evidence for a significant trend with redshift at a given stellar mass in the VANDELS sample. 
It can be seen from Fig. <ref> that the redshift-mass stacks are consistent with scattering around the stellar mass-only relation.

With respect to the normalization of the stellar mass-metallicity relation, our results are consistent with recent results from the independent KBSS-MOSFIRE survey <cit.>.
The orange data point in Fig. <ref> shows the result derived from a stack of 30 star-forming galaxies at z=2.40±0.11 presented in <cit.> who estimated Z_∗/Z_⊙≃ 0.14 using a similar method to our own, by fitting Starburst99 WM-Basic models to the composite FUV spectrum, with the additional requirement that the best-fitting stellar model should also account for the observed nebular optical emission lines (a constraint we currently do not have for our sample).
<cit.> find a lower best-fitting metallicity model using the alternative BPASSv2 stellar models, which is consistent with what we report in Appendix <ref>. 
Despite the slight difference in method, average redshift, and the stellar models employed, the consistency is encouraging, and provides further evidence that the stellar metallicities (or iron abundances) of galaxies at high redshift are significantly sub-solar (Z_∗/Z_⊙≲ 0.25).



 §.§ Redshift evolution of the stellar mass-metallicity relation



Extending the comparison to lower redshifts, we show the evolution of the mass-metallicity relation between 2.3 < z < 5.0 (⟨ z ⟩ = 3.5) and z=0 in Fig. <ref>.
For a our main local Universe comparison sample we have used a recent determination of the stellar mass-metallicity relation for ∼ 200,000 star-forming galaxies in the Sloan Digital Sky survey by <cit.> (Z17).
The solid line (connecting the triangular data points) in Fig. <ref> shows the z=0 relation and the dashed line shows an arbitrary downward shift of this relation by 0.6 dex.
This simple shift matches the VANDELS data remarkably well, implying an increase in metallicity of a factor ∼4 from z ∼ 3.5 to the present day without any obvious dependence on stellar mass.

However, it is important to note that there will be some biases associated with this simple comparison. 
For example, whereas our metallicities are determined from the rest-frame FUV, tracing the massive OB stellar population, and are sensitive to the iron abundance,  derive metallicities by fitting to rest-frame optical spectra in the wavelength range ≃ 4000 - 7000 Å.
It is not immediately obvious that the methods trace either the same stellar population or abundance type.
We can partially address this issue via a comparison with the metallicities derived from the composite FOS-GHRS spectra (blue data points in Fig. <ref>).
In this case we know that the method, at least, is not biasing the comparison since the metallicities have been derived in identical ways.
Although the mass sampling is not ideal, from the three stellar-mass bins available we find the derived stellar metallicities are in approximate agreement with the z=0 relation, especially when taking into account that, when measured for individual galaxies, the scatter in the local stellar mass-metallicity relation is of the order ≳ 0.1 dex <cit.>.
There are again a number of caveats, primarily the fact that, for the high-mass galaxies, the FOS-GHRS observations are often taken for individual  regions within the galaxy rather than the galaxy globally.
However, one can make the argument that these regions are likely to dominate the global FUV emission of these galaxies.
Taken as a whole, the data presented in Fig. <ref> provide strong evidence for an evolution in stellar metallicities of roughly a factor 4 between 2.3 < z < 5.0 and z=0.

Another potential complication worth acknowledging is the fact that, although comparing the stellar metallicities of star-forming galaxies at these redshifts is clearly of interest, it is likely that the two datasets do not form an evolutionary sequence since, depending on their stellar mass, star-forming galaxies at high redshift will likely be the progenitors of present-day quiescent galaxies <cit.>.
Nevertheless, the stellar metallicities of early-type galaxies with M_∗>10^10M_⊙ estimated from rest-frame optical absorption features also fall predominantly within the range -0.2<  <0.2 <cit.>, similar to the values for the star-forming population, and therefore a similar conclusion would be drawn from this comparison.

In summary, we find clear evidence for a monotonic increase in stellar metallicity with galaxy stellar mass across the redshift range 2.3 < z < 5.0.
Within this narrow redshift window we find no strong evidence for redshift evolution at fixed stellar mass.
However, a direct comparison to the stellar metallicities of star-forming galaxies in the local Universe implies that metallicities increase by a factor of 4, independent of mass, between z ∼ 3.5 and z=0.



§ DISCUSSION


Here we discuss some of the implications of our results.
Firstly, we compare the observed data to the predictions of three cosmological simulations.
Then, using a simple one-zone chemical evolution model, we attempt to identify the key physical parameters driving both the significantly sub-solar stellar metallicities, as well as the observed mass dependence in our sample.
Finally, we compare the stellar metallicities derived here with published gas-phase metallicities at similar redshifts in order to investigate suggestions of enhanced O/Fe ratios in star-forming galaxies at z ≳ 2.5.



 §.§ Comparison to cosmological simulations


The left-hand panel of Fig. <ref> compares the absolute  values to simulation predictions for the FUV-weighted stellar mass-metallicity relation. 
In general, the absolute metallicity values predicted in the simulations are either systematically too large () or too small ( and FIRE).
Given these discrepancies are generally a factor ≲ 2, they could be explained by uncertainties in the element yields from core-collapse supernovae (CCSNe), SNe Ia and asymptotic giant branch (AGB) star winds used in the simulations <cit.>.
Indeed, yields in the cosmological simulation Mufasa <cit.> were arbitrarily reduced by a factor of two to match published gas-phase oxygen abundance data.
The correction became unnecessary in  due to a more detailed model for how elements become locked in dust grains <cit.>, however this comparison would suggest that the iron-peak element yields may still not consistent with current data.
Another explanation, illustrated in Section <ref> below, could be related to the strength of galactic outflows in the simulations.

The shape of the relation, on the other hand, appears to be relatively well reproduced by all three simulations.
The right-hand panel of Fig. <ref> shows the data and simulated relations normalized at  = 9.56 (the stellar mass of the VANDELS-m4 composite).
The exact slope of the relation based on the observed data is somewhat uncertain given the upper limit at the lowest stellar mass bin, but the simulations are clearly in relatively good agreement.
The  data, in particular, match the shape of the relation extremely well.

Determining what governs the overall shape of the stellar mass-metallicity relation in simulations is not straightforward due to the complex interplay of star formation, inflows and outflows which determine the metallicity evolution.
However, as we discuss in Section <ref> below, one of the key physical parameters governing the shape and normalization of the relation is the strength and mass-scaling of galactic outflows, parameterized by the mass-loading parameter η = Ṁ_outflow/Ṁ_∗. 
The prescription for how the mass-loading parameter scales with stellar mass (the η-M_∗ relation) in  is, in fact, based on the results of the high-resolution FIRE zoom simulations, which include a detailed stellar feedback model accounting for the effects of radiation pressure, supernovae, stellar winds and photoionization and photoelectric heating <cit.>.
The fact that  and FIRE use a similar outflow model could account for the good agreement in the shape of the MZR between the two simulations.
The feedback model used in  is based on an earlier prescription for supernova feedback <cit.> which results in generally weaker stellar winds (i.e. smaller mass-loading parameters), although crucially the scaling of η with stellar mass is similar to the  and FIRE simulations (S. Khochfar private communication).
The similarity of the slope of the η-M_∗ relation in each simulation could therefore account for the consistency in their predictions for the shape of the MZR. 
On the other hand, the differences in the absolute  values are likely to be a result of both the different normalizations of the η-M_∗ relations, as well as differences in the assumed stellar yields.

Finally, we note that the lack of observed redshift evolution between 2.5 < z < 5.0 is also in reasonable agreement with predictions from simulations. 
For example, the FIRE simulation predicts ≃ 0.16 dex of evolution in  at fixed  across 2.3 < z < 5.0 <cit.>, a value roughly comparable to the 1-2σ error bars on the derived metallicities for the mass-redshift composites.
This is shown by the grey dotted lines in Fig. <ref>, where we have scaled the FIRE metallicities upwards by a factor of 2 such that the z=2.3 and z=5.0 relations bracket our data.
However, given the fact that we are only utilizing ∼ 65 % of the full VANDELS dataset, evidence for a redshift evolution will be an area that can be explored in more detail in future work. 
Fig. <ref> also shows the prediction for the redshift evolution of the stellar mass-metallicity relation to z=0 from the FIRE simulation[Unfortunately, we cannot produce similar predictions from the  and  simulations since they terminate at redshifts z=4 and z=1 respectively.].
Again, it can be seen that the predicted evolution of the stellar MZR by ≃ 0.6 dex is in remarkable agreement with the observations.
We also note that an evolution of this magnitude is predicted from the EAGLE simulations <cit.>, although we do not show this relation in Fig. <ref> for clarity.
In summary, the current predictions from cosmological simulations seem to be in broad agreement with our data.




 §.§ Analytic one-zone chemical evolution models


Although simple analytic models lack the rigor of detailed cosmological simulations, they can still provide useful physical insights into galactic chemical evolution. 
The primary advantage of using analytic approximations is, perhaps, the ability to quickly determine the sensitivity of chemical evolution to different physical parameters, a procedure which is computationally expensive in the case of detailed simulations.
Such analytic models are typically specified by a gas-accretion history, a star-formation law, and prescriptions for gas inflows/outflows and nucleosynthetic yields.
Numerous analytic models exists <cit.> however,
in this section, we compare our results with the analytic one-zone chemical evolution model presented in <cit.> (hereafter ).
The  model is particularly suited to our purpose since, as discussed, the FUV stellar metallicities we derive are a proxy for the iron abundance, and  incorporate a realistic delay time distribution (DTD) for Type-Ia SNe, allowing the evolution of iron peak elements to be tracked separately to that of the α elements.

    


Briefly, the key elements of the  model are as follows.
The star-formation rate (Ṁ_∗) is assumed to be proportional to the gas mass (M_gas) with the star formation efficiency (SFE =Ṁ_∗/M_gas) assumed to be constant (the so-called `linear Schmidt law').
For the purposes of this discussion we will refer to the inverse of the star-formation efficiency as the gas depletion timescale (t_dep=M_gas/Ṁ_∗) as this is a quantity that has been extensively studied out to z≃6 <cit.>.
Gas outflows are assumed to be a constant multiple of the star-formation rate, with the constant of proportionality commonly referred to as the mass-loading parameter (η = Ṁ_outflow/Ṁ_∗).
Therefore, for a specified star-formation history, the gas depletion timescale and mass-loading parameter define the total gas-accretion history.
The gas accreted from the IGM is assumed to be pristine (Z=0) and the outflowing gas is assumed to be at the same metallicity as the ambient star-forming ISM[We note that this is one important simplifying assumption of the model, since there is strong evidence that outflows are in fact metal-enriched with respect to the average ISM metallicity <cit.>.].
Three star-formation histories are explicitly solved for in  but we only consider the linear-exponential case here[In the linear-exponential case, the star-formation history is modeled as the product of a linear rise and exponential decline: Ṁ_∗(t) ∝ te^-tτ_ sfh.].
For our purposes, this essentially equates to a linearly rising star-formation history.
Elements produced by core-collapse supernova (CCSNe) are instantaneously recycled into the star-forming ISM, whereas enrichment from Type Ia SNe is delayed as specified by the DTD.
The form of the DTD is exponentially decreasing in time with a minimum delay time of 0.15 Gyr and an e-folding timescale of 1.5 Gyr[Again, it should be noted that this DTD is not necessarily consistent with current observational data <cit.> and represents another source of systematic uncertainty when considering the predictions of the  model.].
The models assumes a <cit.> IMF and the CCSNe yields of <cit.> and <cit.> which predict that 1.5M_⊙ of oxygen and 0.12 M_⊙ of iron are returned to the ISM for every 100 M_⊙ of star formation.
The net return of iron from Type Ia SNe, assuming the fiducial model parameters, is 0.17 M_⊙ for every 100 M_⊙ of star formation.
These nucleosynthetic yields are assumed to be independent of stellar metallicity.
In general yields will be metallicity dependent, however, for iron, the metallicity dependence is not expected to be strong <cit.>.
A fixed fraction of stellar mass formed into stars is assumed to be recycled from the envelopes of CCSNe progenitors at its original metallicity and is referred to as the recycling parameter (r=0.4).
All free parameters excluding the star formation rate are assumed to be fixed with time[The fiducial model parameters are taken from a numerical implementation of the  model presented in <cit.>, and are adopted here unless explicitly stated.
They are calibrated to reproduce the [O/Fe] versus [Fe/H] distribution of local thin disk, thick disk and halo stars.].
Using these various approximations, it is possible to derive analytic equations for the evolution of chemical abundances. 
We refer interested readers to the original  paper for detailed derivations.

As discussed in , the main governing parameters in their model are the mass-loading parameter (or outflow efficiency) η and the star-formation efficiency, parameterized by the familiar gas-depletion timescale (t_dep=M_gas/Ṁ_∗).
The left-hand panel of Fig. <ref> shows the predicted time evolution of  for different values of η and t_dep. 
Here,  refers to the iron abundance ([Fe/H]) in the  models.
It can clearly be seen that the assumed value of the mass-loading parameter has a pronounced effect, strongly influencing the typical  value reached after ≃ 2.5  of evolution (corresponding to z ≃ 2.5). 
The reason is simply that by increasing the mass-loading parameter (from η=3-10 in our example) ISM enrichment is suppressed by the removal of more enriched gas from the galaxy.
In contrast, it is difficult to ascribe the low absolute abundances to a combination of relatively weak stellar winds (η=3) with long depletion timescales  (i.e. low star-formation efficiency, which will lead to lower metallicities simply due to the fact that fewer metals are being formed in stars).
As can be seen from Fig. <ref>, even depletion timescales of 2  are not sufficient to reach agreement with the data. 
In reality, timescales of ≃ 2  are more typical of local star-forming galaxies <cit.>, and it has been found from observations of cold gas in high-redshift galaxies that the depletion timescale appears to decrease with increasing redshift, with typical values in the range 200 - 700  across z=1-6 <cit.>.
Therefore, within the context of the  model, the mass of metals produced assuming realistic values of the depletion timescale is too large to be consistent with the observations unless relatively large mass-loading parameters (i.e. high outflow efficiencies, η≃ 10 on average across our sample) are assumed.
We note that the importance of strong outflows in shaping both the stellar and gas-phase metallicity is also found in many semi-analytic galaxy evolution models <cit.>.

Further insights can be gained by investigating the mass dependence of the stellar metallicity. 
However, since the  model is only parameterized in terms of η, to simulate mass dependence a mapping between η and mass is required.
Motivated by the outflow scaling relations from the FIRE simulations presented in <cit.>[We note that these scaling relation have since been updated by <cit.> and <cit.> but we consider this simple parameterization sufficient for our purposes here.], we assume η is related to stellar mass via

    η = α(M_∗/10^10M_⊙)^β,

where M_∗ is the galaxy stellar mass, α is the mass-loading parameter at M_∗=10^10M_⊙ and β is the power-law exponent.
<cit.> find a redshift-independent relation with α=3.55 and β=-0.351.
Equation <ref> can be easily incorporated into the  model in order evaluate the values of α and β (and hence the η(M_∗) relation) most consistent with our data.

The right-hand panel of Fig. <ref> shows the difference in metallicity (Δ) as a function of stellar mass for different parameterizations of the stellar mass versus η relationship assuming both the <cit.> normalization (α=3.6) and a slightly higher normalization implied by the  models (a value of α≃6 is consistent with the highest mass - i.e. highest metallicity - data point in the left-hand panel of Fig. <ref>).
The Δ values are calculated with respect to  ≃ 9.56.
The data favor a η-M_∗ scaling relation with a similar power-law slope to the <cit.> relation, irrespective of the absolute normalization.
A simple χ^2 analysis returns a best-fitting power-law slope of β=(-0.45, -0.40) for the α=(3.6, 6.0) cases respectively.
In either case the mass-loading parameter decreases by roughly a factor of 5 across the stellar mass range 10^8 < M_∗/M_⊙ < 10^10 with a median value of η=8 for α=3.6 and η=15 for α=6.0.
We note there is some degeneracy here between α and the assumed yields, however the corrections to the  would have to be quite large to make η=3.0 consistent with the observed data (e.g. both the CCSNe and Type Ia SNe yields would have to be decreased by a factor 3).

In summary, whilst acknowledging the caveats associated with simple analytic models, we find that the  model supports a scenario in which the significantly sub-solar iron abundances of the galaxies in our sample, and the dependence of iron abundance on stellar mass, are a consequence of strong outflows which scale with  to the power ≃ -0.4.
In this scenario, the mass of gas outflowing from the galaxy, at M_∗ = 10^9 M_⊙, is predicted to be roughly an order-of-magnitude greater than the gas mass being incorporated into new stars.
Finally we note that the inferences drawn from the  model are broadly similar to the those drawn from hydrodynamical simulations.
Although these methods are not fully independent, since simulations still rely some analytic sub-grid recipes, the agreement is encouraging and highlights the utility of simple analytic prescriptions for chemical evolution.



 §.§ Enhanced O/Fe ratios at high redshift


    

Important insights can be gained by comparing the stellar metallicities of our sample to published gas-phase metallicities at similar redshifts.
The comparison is of interest because, while the stellar metallicities derived here are sensitive to elements dominating photospheric absorption in massive stars (namely iron), the gas-phase metallicities are determined from the nebular emission lines emitted from the surrounding  regions and are therefore sensitive to the dominant nebular coolants (namely oxygen).
Under the assumption that the chemical abundances of massive stars should be similar to the surrounding  region gas-out of which they presumably formed-then comparing stellar and gas-phase metallicities should be a good proxy for the O/Fe ratio in these galaxies <cit.>.

Non-solar O/Fe ratios should be expected for galaxies with constant or rising star-formation histories and relatively young ages (<500  - 1 ). 
This results from the fact that the release of Fe into the ISM is determined by both the rate of Type-Ia SNe and CCSNe, while the release of O (and other α elements) is determined only by the rate of CCSNe.
Crucially, the Type-Ia SNe rate is sensitive to the star-formation rate of a galaxy over a large range of epoch before the time of observations <cit.>, while CCSNe enrich the ISM almost instantaneously (within a few  of a given star-formation episode).
Therefore, for constant or rising SFHs, the abundance of the ISM will always be dominated by CCSNe products at young ages.
The yields of <cit.> predict, for a Salpeter IMF <cit.>, that CCSNe yields produce (O/Fe) ≃ 4 - 6 × (O/Fe)_⊙ for initial stellar metallicities ≃ 0.1 - 2 × Z_⊙.
The yields of the <cit.> and <cit.> fiducial model predict that CCSNe produce 1.5M_⊙ of O and 0.12M_⊙ of Fe for every 100M_⊙ of star formation assuming a Kroupa IMF <cit.>, equivalent to (O/Fe) ≃ 3 × (O/Fe)_⊙.
Given the growing observational and theoretical evidence for constant or rising star-formation histories at the redshift of our sample <cit.> we might therefore expect to see (O/H) ≃ 3-6 × (Fe/H).

Indeed, this result has already been reported for star-forming galaxies at ⟨ z ⟩=2.3 from the KBSS-MOSFIRE sample.
<cit.> present an analysis of the optical spectra of 148 galaxies allowing them to constrain [Fe/H] and [O/H] for individual galaxies, finding an average [O/Fe] =0.42 (i.e. (O/H) ≃ 2.6 × (Fe/H)), somewhat lower than the initial [O/Fe] ≃ 0.7 reported for a stack of 30 star-forming galaxies at the same redshift derived from a combined analysis of FUV and optical spectra presented in <cit.>.
Nevertheless, both values are comparable to the expected range of [O/Fe] from the CCSNe yield models described above. 

Although we are unable to directly derive [O/H] for the galaxies in our sample, it is still instructive to compare the stellar metallicities reported here to [O/H] values at similar redshifts from the literature.
To this end we show, in Fig. <ref>, the gas-phase mass metallicity relationship for a sample of 34 galaxies at 3 < z < 5 from the AMAZE/LSD survey <cit.> compared to the stellar metallicities derived for our VANDELS sample.
The AMAZE/LSD sample spans the stellar mass range 8.5 <  < 11.0 and is composed of Lyman-break selected galaxies with a bias to the most highly star-forming galaxies at a given mass <cit.>.
The metallicities were derived from a combination of the ,  and  nebular emission lines using the semi-empirical metallicity calibration of <cit.> (M08).
Fig. <ref> shows that the stellar metallicities (Fe/H) of star-forming galaxies at z≃3.5 are systematically lower than the gas-phase metallicities (O/H) by ≃ 0.25 dex at all stellar masses compared to the mean relation found by <cit.>.
This value is comparable to the prediction of the  chemical evolution model described above, which, assuming the best-fitting η-M_∗ relation, predicts [O/Fe] ≃0.2 dex at z=3.5 (i.e. after ≃1.8  of evolution) across all stellar masses[This value is lower than that predicted directly from the yields because it incorporates the assumed star formation history and Type-Ia SNe DTD.].

However, there are a number of caveats that should be acknowledged with respect to this comparison. 
Firstly, we note that the galaxies in the <cit.> sample have, on average, higher star-formation rates at fixed stellar mass than the galaxies in our sample (the median sSFR of the AMAZE/LSD galaxies is 5.4 ^-1 compared to 4.4 ^-1 for our sample).
Assuming a stellar mass-SFR-metallicity relation <cit.> is in place at high redshifts (as claimed recently by <cit.> at z=2.3) then, due to their higher star-formation rates, the AMAZE/LSD gas-phase metallicities will be biased low with respect to the true gas-phase metallicities of our sample.
In this case the offset of [O/Fe] =0.25 dex would be a lower limit.
Secondly, the comparison is affected by zero-point offsets in both the stellar and gas-phase metallicity values.
For example, in Appendix <ref> we show how using an alternative stellar population model <cit.> can result in a systematic shift of ≃ - 0.1 dex in .
Additionally, the zero-point offset in gas-phase metallicity calibrations, which has been well-studied in the local Universe <cit.>, implies that the calibration adopted by <cit.> (M08) returns metallicities intermediate between other calibrations with offsets of up to ≃± 0.2 dex across the stellar mass range shown in Fig. <ref> <cit.>.

Clearly, this comparison is currently limited to qualitative statements, and is hampered by the different selection biases affecting the VANDELS and AMAZE/LSD samples. 
A more definitive answer will have to wait until rest-frame optical observations of VANDELS galaxies are available, from which [O/H] can be directly estimated.
Nevertheless, adopting currently published values suggests that the O/Fe ratio is enhanced in high-redshift galaxies relative to the solar values, with (O/Fe) ≳ 1.8 × (O/Fe)_⊙.
One major practical implication of this result is that stellar population models which assume solar chemical abundance ratios are likely not valid for modeling the full spectra energy distribution of galaxies at high-redshift <cit.>[Although the Starburst99 models used here assume a solar abundance pattern, the fact that the FUV spectral features are primarily determined by Fe abundance means our recovered  values should not be strongly dependent on the detailed abundance pattern.
Nevertheless, in future it would clearly be desirable to test this assumption.].



§ CONCLUSIONS


In this paper we have presented the results of a study of the stellar mass-metallicity relationship for a sample of 681 star-forming galaxies at 2.3 < z < 5.0 spanning the stellar masses range 8.2 <  < 11.0, based on deep rest-frame FUV spectroscopic data from the VANDELS survey <cit.>.
Stellar metallicities (to first order a proxy for the iron abundance) have been derived by fitting high-resolution theoretical Starburst99 WM-Basic stellar population synthesis models <cit.> to high signal-to-noise composite spectra in bins of stellar mass and redshift.
Our method is tested using both synthetic spectra generated from simulation data and FUV spectra of galaxies in the local Universe.
Finally, we have investigated the normalization and shape of the stellar mass-metallicity relationship at 2.3 < z < 5.0, and compared our results to predictions from state-of-the-art cosmological simulations and simple one-zone analytic models for chemical evolution.
Our main results can be summarised as follows:



    
  * We find a strong correlation between stellar metallicity and stellar mass at 2.3 < z < 5.0 (⟨ z ⟩ = 3.5) with  monotonically increasing from < -1.04 for galaxies with ⟨ M_∗⟩ = 3.2 × 10^8M_⊙ to -0.57 at ⟨ M_∗⟩ = 1.7 × 10^10M_⊙.
    Across the full mass range  the metallicities (iron abundances) we derive are significantly sub-solar with Z_∗≲ 0.25 Z_⊙.

    
  * We do not observe a strong relationship between stellar metallicity and redshift within our current sample.
    Within the error bars, the metallicities derived for composite spectra binned in redshift and stellar mass are consistent with scattering around the stellar-mass only relationship.
    
    
  * The stellar mass-metallicity relation we derive at ⟨ z ⟩ = 3.5 is offset to lower metallicities by ≃ 0.6 dex compated to the z=0 relation for star-forming galaxies derived from optical spectroscopic data <cit.>.
    There are some systematic uncertainties related to this comparison (e.g. probing different stellar populations and/or abundance types), but similar results are obtained by comparing to the FUV-based metallicities at z=0 from FOS-GHRS data. 

    
  * The normalization of the ⟨ z ⟩ = 3.5 stellar mass-metallicity relation is not well reproduced by simulations.
    Absolute metallicity values tend to fall either systematically above (e.g ) or below (e.g. , FIRE) our data.
    However, given the relatively large uncertainties in stellar yields, it is hard to draw strong conclusions from this observation.
    On the other hand, the shape of the relation is generally well recovered.
    We postulate that this is likely due to the fact that the scaling between outflow efficiency (parameterized by the mass loading parameter η = Ṁ_outflow/Ṁ_∗) and stellar mass is similar across all of the simulations.

    
  * We investigate the parameters affecting the shape and normalization of the stellar mass-metallicity relation using detailed one-zone analytic chemical evolution models <cit.>. 
    We find that these models can reproduce the observed shape the relation by assuming a power-law scaling relation between η and M_∗  that is similar to that derived from fully cosmological simulations <cit.> and predicted from analytic models of momentum-driven winds <cit.>.
    A comparison of the models and simulations to our data suggest that η∝ M_∗^β with β≃ -0.4.
    
    
  * Furthermore, simulations which include state-of-the-art prescriptions for stellar feedback (FIRE and ), as well as the  model, suggest that an average mass-loading parameter of ⟨η⟩≃ 10 is required for consistency with our data, ranging from η≃ 6 at M_∗ = 10^10M_⊙, to η≃ 40 at M_∗ = 10^8M_⊙.
    Although the absolute normalization of the η-M_∗ relation is subject to yield uncertainties, the reasonable agreement with predictions from cosmological simulations is encouraging.

    
  * Finally, by comparing the stellar metallicities derived here (which trace Fe/H)  to published gas-phase metallicities at similar redshifts <cit.>, we find that the gas-phase abundances (which trace O/H) are likely enhanced by at least a factor of ≃ 2.
    This comparison provides further evidence that star-forming galaxies at these epochs are α-enhanced systems, although it will require direct [O/H] estimates for the galaxies in our sample to provide a more robust measurement.





§ ACKNOWLEDGMENTS


FC, RJM, JSD, SK, AC and DJM acknowledge the support of the UK Science and Technology Facilities Council.
This work is based on data products from observations made with ESO Telescopes at La Silla Paranal Observatory under ESO programme ID 194.A-2003(E-Q).
AC acknowledges the grants PRIN MIUR 2015, ASI n.I/023/12/0 and ASI n.2018-23-HH.0
GC has been supported by the INAF PRIN-SKA 2017 program 1.05.01.88.04.
This research made use of Astropy, a community-developed core Python package for Astronomy <cit.>, NumPy and SciPy <cit.>, Matplotlib <cit.>, IPython <cit.> and NASA's Astrophysics Data System Bibliographic Services. 



mnras                      
       





§ FOS-GHRS LOCAL SAMPLE


 

    

The local FUV spectra are taken from a complication of 28 local starbursts and star-forming galaxies presented in <cit.> .
 present 46 rest-frame UV spectra observed with the Faint Object Spectrograph (FOS) and the Goddard High Resolution Spectrograph (GHRS) of the HST.
The spectral resolution of the individual spectra spans the range 0.5-3Å depending on the instrument, aperture size and physical extent of the object being observed.
For consistency with the the VANDELS data, we smooth all spectra to a common 3Å resolution.
We also require that the spectra have wavelength coverage in the interval 1410 ≤λ≤ 1450 Å to enable an analysis of normalized composite spectra, and finally that the corresponding galaxy has a measurement of absolute K-band magnitude which we can use to estimate the stellar mass.
This selection leaves 26/46 of the original spectra from 18/28 of the original local starbursts and star-forming galaxies presented in .
A list of the individual spectra used in this work is presented in Table <ref>. 



§ EFFECT OF REDSHIFT UNCERTAINTIES


One potential source of bias in our analysis relates to the effect of redshift uncertainties when stacking.
The redshifts for the individual VANDELS spectra are constrained primarily by the strong interstellar absorption features and/or  emission line.
These features originate from outflowing gas and are therefore often redshifted or blueshifted by up to hundreds of kms^-1 with respect to the systemic redshift of the galaxy.
For example, from a sample of 89 galaxies at z∼2-3, <cit.> find an average velocity offset for redshifts derived from interstellar absorption lines of ⟨Δ v_IS⟩ = -164 ± 16 kms^-1, and of ⟨Δ v_Lyα⟩ = +445 ± 27 kms^-1 for redshifts derived from  emission.
Each galaxy in our sample will therefore have some unknown redshift offset which could potentially result in a systematic bias when deriving metallicities.
In particular, the weak continuum features which are used to constrain the metallicity via continuum fitting could be `washed out' to make the continuum look increasingly featureless, and therefore bias measurements towards lower metallicity solutions.

We tested the effect of redshift uncertainties using simulated galaxy spectra.
For each of the five default Starburst99 WM-Basic templates (Z_∗=0.001, 0.002, 0.008, 0.014, 0.040) we produced 100 versions (i.e. roughly the number of galaxies in our mass stacks) each with a random velocity offset distributed uniformly within the range -150 < Δ v  < 500 kms^-1.
The 100 spectra were then median combined to mimic stacking spectra with random, unknown, redshift uncertainties.
These spectra were convolved to the resolution of the VANDELS data and gaussian noise was added assuming a S/N per pixel of 15. 
The metallicity was then derived using the method applied to the observed data as described in Section <ref>.
For each template the process was repeated 100 times and the distribution of recovered metallicities is shown in Fig. <ref>.
It can be seen that although there is a systematic offset to lower metallicities for all but the lowest metallicity template (for which it is impossible to measure a lower metallicity), the effect is relatively small (≤ 10 %).
We therefore conclude that the low metallicities derived for our VANDELS sample are robust against the effect of redshift uncertainties.



§ BPASSV2.1 STELLAR POPULATION SYNTHESIS MODELS


In our main analysis we focused on using the Starburst99 WM-Basic stellar population synthesis models to fit the VANDELS spectra for the reasons outlined in Section <ref>.
However, we have also performed the same analysis using the BPASSv2.1
 stellar population synthesis models described in <cit.>.
It is well known that using different models will cause systematic offsets in derived parameters <cit.> and it is useful to have an estimate of how significant this offset my might be in the case of FUV-derived metallicities. 

The unique feature of the BPASSv2.1 models is the inclusion of massive binary star evolution, which can have a strong effect on the predicted UV spectrum, particularly at low metallicities.
For example, at Z_∗/Z_⊙ < 0.35, a phenomenon known as `quasi-homogeneous evolution' can occur <cit.>, which results in stars living longer on the main-sequence and becoming hotter than single stars at similar mass and metallicity.
QHE can have a profound effect on the estimated ages and ionizing photon output of stellar populations.
Stellar metallicities should be less affected, particularly since the O-star models in BPASSv2.1 are also generated with WM-Basic code, albeit with a different parameter set to Starburst99.
Nevertheless, it worth testing this assumption.

	

We considered a 100  constant star formation BPASSv2.1 model which included binary star evolution and had an IMF cutoff of 100 M_⊙, with IMF index of -1.3 between 0.1 - 0.5 M_⊙ and -2.35 above 0.5 M_⊙.
The model is referred to as BPASSv2.1-100bin.
We considered the following set of metallicity values Z_∗ =
(0.001, 0.002, 0.003, 0.004, 0.006, 0.008, 0.010, 0.014, 0.020,
0.030, 0.040), interpolating between the models in the same way as for Starburst99.
Each of the seven mass stacks were fitted with the BPASSv2.1-100bin model in the same way as was done for the Starburst99 models (see Section <ref> for details).
Fig. <ref> compares the metallicities derived from the BPASSv2.1-100bin model and the metallicities derived from Starburst99.
It can be seen that the metallicity estimates are clearly correlated, with the BPASSv2.1 values offset to slightly lower metallicity by ∼ 0.09 dex (roughly a factor of 1.2).
However, it is clear from this figure that adopting the BPASSv2.1 model estimates would not affect our main results.



^1SUPAScottish Universities Physics Alliance, Institute for Astronomy, University of Edinburgh, Royal Observatory, Edinburgh EH9 3HJ

^2University of the Western Cape, Bellville, Cape Town 7535, South Africa

^3South African Astronomical Observatories, Observatory, Cape Town 7925, South Africa

^4Instituto de Investigación Multidisciplinar en Ciencia y Tecnología, Universidad de La Serena, Raúl Bitrán 1305, La Serena, Chile

^5Departamento de Física y Astronomía, Universidad de La Serena, Av. Juan Cisternas 1200 Norte, La Serena, Chile

^6INAF - Osservatorio Astronomico di Bologna, via P. Gobetti 93/3,
I-40129, Bologna, Italy

^7INAF-Osservatorio Astronomico di Roma, Via Frascati 33, I-00040 Monte Porzio Catone (RM), Italy

^8University of Bologna, Department of Physics and Astronomy (DIFA) 
Via Gobetti 93/2- 40129, Bologna, Italy

^9INAF - Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125, Firenze, Italy

^10European Southern Observatory, Karl-Schwarzschild-Str. 2, 86748 Garching b. München, Germany

^11The Cosmic Dawn Center, Niels Bohr Institute, University of Copenhagen, Juliane Maries Vej 30, DK-2100 Copenhagen Ø, Denmark

^12INAF-Astronomical Observatory of Trieste, via G.B. Tiepolo 11, 34143 Trieste, Italy

^13INAF-IASF Milano, via Bassini 15, I-20133, Milano, Italy

^14 Núcleo de Astronomía, Facultad de Ingeniería, Universidad Diego Portales, Av. Ejército 441, Santiago, Chile

^15Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

^16European Southern Observatory (Chile)

^17Department of Physics and Astronomy, University of California, Los Angeles, 430 Portola Plaza, Los Angeles, CA 90095, USA



