LATIS: The Stellar Mass–Metallicity Relation of Star-forming Galaxies at z ∼ 2.5

We present the stellar mass–stellar metallicity relation for 3491 star-forming galaxies at 2 ≲ z ≲ 3 using rest-frame far-ultraviolet spectra from the Lyα Tomography IMACS Survey (LATIS). We fit stellar population synthesis models from the Binary Population And Spectral Synthesis code (v2.2.1) to medium-resolution (R ∼ 1000) and high signal-to-noise (>30 per 100 km s−1 over the wavelength range 1221–1800 Å) composite spectra of galaxies in bins of stellar mass to determine their stellar metallicity, primarily tracing Fe/H. We find a strong correlation between stellar mass and stellar metallicity, with stellar metallicity monotonically increasing with stellar mass at low masses and flattening at high masses (M * ≳ 1010.3 M ⊙). Additionally, we compare our stellar metallicity measurements with the gas-phase oxygen abundance of galaxies at similar redshift and estimate the average [α/Fe] ∼ 0.6. Such high α-enhancement indicates that high-redshift galaxies have not yet undergone significant iron enrichment through Type Ia supernovae. Moreover, we utilize an analytic chemical evolution model to constrain the mass loading parameter of galactic winds as a function of stellar mass. We find that as the stellar mass increases, the mass loading parameter decreases. The parameter then flattens or reaches a turning point at around M * ∼ 1010.5 M ⊙. Our findings may signal the onset of black-hole-driven outflows at z ∼ 2.5 for galaxies with M * ≳ 1010.5 M ⊙.


INTRODUCTION
Cosmic primordial gas is predominantly composed of hydrogen and helium, with small amounts of other light elements such as lithium.The gas accreted from the intergalactic medium (IGM) and circumgalactic medium (CGM) enables a galaxy to form stars, which produce heavy elements.Feedback processes such as stellar winds and supernova explosions expel some of these heavy elements into the interstellar medium (ISM) where new stars are born.Moreover, galactic outflows driven by supernovae and black hole feedback can transfer some enriched material to the IGM and CGM (e.g., Tremonti et al. 2004;Chisholm et al. 2018).Therefore, the metal content of galaxies is linked to their fundamental evolutionary processes (e.g., star formation and inflow/outflow), and determining its relationship to global properties, such as stellar mass (M * ), provides Corresponding author: Nima Chartab nchartab@carnegiescience.edu useful information for constraining galaxy evolution models (see review by Maiolino & Mannucci 2019).
The metal content of high redshift galaxies (2 ≲ z ≲ 3) is typically measured in the gas phase (metallicity of ISM) using strong rest-frame optical emission line diagnostics (e.g., Pettini & Pagel 2004) that are calibrated locally, suffering significant uncertainties.High-redshift galaxies have been observed to exhibit distinct physical conditions in their H ii regions compared to those at z = 0 (e.g., Erb et al. 2006;Steidel et al. 2014;Shapley et al. 2015), implying a potential evolution of metallicity calibrations with redshift.For direct estimates of gasphase metallicities, faint auroral lines (e.g., [Oiii] λ4363) need to be detected to determine the electron temperature of the ionized gas, which is now possible for a statistically significant number of high redshift galaxies thanks to the James Webb Space Telescope (e.g., Sanders et al. 2023, Curti et al. 2023, CECILIA survey;Strom et al. 2021, AURORA survey;Shapley et al. 2021).
Alternatively, stellar continuum emission can be used to measure the metal content of galaxies.Deep optical spectroscopy of z ∼ 2 − 3 galaxies allows us to access the rest-frame far-ultraviolet (FUV) part of the spectrum that contains important information about their underlying stellar population.Most of the emission in the FUV originates from short-lived O and B stars, and the inferred metallicities are expected to be similar to those derived for the ISM, out of which these young stars have recently formed.The photospheres of hot O and B stars, metal-dependent stellar winds, and interstellar lines all contribute to the FUV absorption features.These features usually have complex dependencies on age, metallicity, and initial mass function (IMF).A number of indices (e.g., 1425 Å and 1978 Å indices) have been identified that are optimized to depend only or mostly on metallicity (e.g., Leitherer et al. 2001;Rix et al. 2004;Sommariva et al. 2012).
In recent studies, full spectral fitting has been used to measure the stellar metallicity of high-redshift galaxies (e.g., Steidel et al. 2016;Cullen et al. 2019;Kriek et al. 2019;Topping et al. 2020;Kashino et al. 2022;Carnall et al. 2022).This method has the advantage of using all of the information in the spectra.Steidel et al. (2016) used a composite rest-frame UV and optical spectrum of 30 star-forming galaxies at z = 2.4 and found them to be α-enhanced relative to the solar abundances by a factor of 4-5.While Type Ia supernovae (SNe Ia) and core-collapse supernovae (CCSNe) both produce iron peak elements, α-elements are only produced by massive, short-lived stars that explode as CCSNe.Thus, the α/Fe abundance ratio is a powerful tool for constraining the relative contribution of SNe Ia and CCSNe.Based on [α/Fe] ∼ 0.6, they conclude that CCSNe dominate the enrichment of z ∼ 2 star-forming galaxies.SNe Ia provide enrichment over long timescales (1-3 Gyr) (Maoz & Mannucci 2012), while young high redshift galaxies with ages ≲ 1 Gyr are not sufficiently old for iron enrichment (Strom et al. 2017).
Although it is now well established that there is a strong positive correlation between the gas-phase metallicity and stellar mass of galaxies out to z ∼ 3.5 (e.g., Erb et al. 2006;Steidel et al. 2014;Sanders et al. 2020;Strom et al. 2022), there are only a few studies on the relationship between stellar metallicity and stellar mass, especially at high redshift (Cullen et al. 2019;Calabrò et al. 2021;Kashino et al. 2022).Cullen et al. (2019) used stacks of 681 star-forming galaxies at z = 2.5-5 from the VANDELS survey (Pentericci et al. 2018) with a spectral resolution of R ∼ 580, and Kashino et al. (2022) utilized 1336 star-forming galaxies at z = 1.6-3 drawn from zCOSMOS-deep survey (R ∼ 200) (Lilly et al. 2007) to study the stellar mass -stellar metallicity relation (hereafter stellar MZR) around cosmic noon.
Recently, the Lyα Tomography IMACS Survey (LATIS; Newman et al. 2020) has obtained deep optical spectroscopy of ∼ 3800 star-forming galaxies at 2 ≲ z ≲ 3 with a spectral resolution of R ∼ 1000.Compared to earlier studies of the stellar MZR, this data set is substantially larger and benefits from significantly higher resolution.
In this paper, we utilize the Binary Population And Spectral Synthesis code (BPASS v2.2.1;Eldridge et al. 2017;Stanway & Eldridge 2018) models to constrain the z ∼ 2.5 stellar MZR by fitting composite rest-frame FUV spectra of 3491 galaxies spanning a wide range of stellar masses, 10 9 M ⊙ ≤ M * ≤ 10 11.5 M ⊙ .We do not employ the latest release of BPASS models (v2.3;Byrne et al. 2022) with α-enhanced spectra.Although these αenhanced models are highly desired for applications at high redshift, they are limited to Main Sequence/Giant Branch stars while the spectra of OB and Wolf-Rayet stars remain unchanged from BPASS v2.2.1 and are not α-enhanced (private communication, C. Byrne), making them unsuitable for fitting FUV spectra dominated by OB stars.The paper is organized as follows.In Section 2, we present an overview of the LATIS survey and details of the sample used in this work.We then describe our spectral analysis to estimate the stellar metallicities in Sections 3 and 4. Our results are presented in Section 5. We discuss our results in Section 6 and summarize them in Section 7.
Throughout this work, we assume a flat ΛCDM cosmology with H 0 = 70 kms −1 Mpc −1 , Ω m0 = 0.3 and Ω Λ0 = 0.7.All magnitudes are expressed in the AB system, and the physical parameters are measured assuming a Chabrier (2003) IMF.We adopt the solar metallicity values of Z ⊙ = 0.0142 and 12 + log(O/H) ⊙ = 8.69 (Asplund et al. 2009).In this abundance scale, the solar oxygen and iron mass fractions are 0.00561 and 0.00128, respectively.

DATA
The LATIS survey is a five-year program (2017-2022) conducted using the Inamori-Magellan Areal Camera and Spectrograph (IMACS; Dressler et al. 2011) at the Magellan Baade telescope.The primary goal of LATIS is producing three-dimensional maps of the z ∼ 2.5 IGM at Mpc resolution, as traced by Lymanα absorption.LATIS densely sampled Lyman-break galaxies (LBGs) in three legacy fields, the Cosmic Evolution Survey (COSMOS; Scoville et al. 2007) and the Canada-France-Hawaii Telescope Legacy Survey (CFHTLS) D1/D4 fields.The final release catalog of CFHTLS (T0007)1 and the Ilbert et al. (2009) catalog of COSMOS were used for selection of LBGs in the CFHTLS D1/D4 and COSMOS fields, respectively.Galaxies were selected based on either their photometric redshifts or ugr optical colors with r−band magnitudes 23.0 < r < 24.8.The observations were performed using the custom grism and filter installed on IMACS resulting in a spectral resolving power of R = 990 and a spectral coverage of 3890-5830 Å.This spectral resolution is estimated at the midpoint of the wavelength coverage and is derived from the average size of targets, which are not point-like, in the two-dimensional spectra.For a full description of the survey strategy, observation and data reduction, we refer readers to Newman et al. (2020).
In the present paper, we use a sample of z > 1.5 galaxies with high-confidence spectroscopic redshifts based on multiple lines and a well-modeled spectrum (i.e., zqual=3 or zqual=4 defined in Newman et al. 2020).Out of the 7408 galaxies that were initially targeted with LATIS, 6568 received deep exposures as part of the main survey while 840 were observed on backup masks of bright candidates intended for poor conditions.A total of 5443 received a redshift quality flag (zqual) of 3 or 4 and had no severe data reduction problems.After further exclusion of 45 AGNs, 105 QSOs, and 44 spectra comprised of two heavily blended sources, we are left with 5249 galaxies, of which 3888 have a redshift z > 1.5.In the following sections we refine our sample further by applying stellar mass limits and requiring near-IR photometry, which narrows our sample to 3491 galaxies with 10 9 M ⊙ ≤ M * ≤ 10 11.5 M ⊙ .

Photometry
Multi-wavelength coverage from the UV to midinfrared is essential for fitting the spectral energy distributions (SEDs) of our sample and measuring physical parameters, such as stellar mass and star formation rate (SFR) (e.g., Chartab et al. 2023).The photometry of the LATIS galaxies in the COSMOS field is taken from the latest version of the COSMOS catalog (COS-MOS2020; Weaver et al. 2022).In the CFHTLS D1 and D4 fields, we construct new photometric catalogs by combining optical to near-infrared (NIR) photometry from multiple surveys.In both fields, we use ugriz images from the CFHTLS release T0007 and Spitzer/IRAC channel 1 and 2 mosaics from Annunziatella et al. (2018).In the D1 field, we use zY JHK s images from the VIDEO survey (Jarvis et al. 2013) while in the D4 field we use JHK s images from the WIRDS survey (Bielby et al. 2012).Images are registered and point spread function (PSF)-matched following the procedures described by Newman et al. (2012).
We use SWarp (Bertin 2010) to construct a χ 2 detection image from the NIR (VIDEO or WIRDS) images available in each field, and run SExtractor (Bertin & Arnouts 1996) in dual-image model to construct catalogs from the PSF-matched images.We measure colors in 2 ′′ diameter apertures and scale the fluxes to match FLUX AUTO in the r band.For the distinct treatment of the IRAC bands, see Newman et al. (2012).
When using the zeropoints provided by each survey, we find that the optical-NIR and NIR-IRAC colors of D1/D4 LATIS galaxies are inconsistent with those in the COSMOS field.This seems to arise, at least in part, from a flux-dependent PSF in the NIR images of the D1 and D4 fields; this effect means that matching the PSFs measured with stars does not match the PSFs applicable to faint galaxies.We therefore adjust the NIR and IRAC zeropoints to match the median optical-NIR and NIR-IRAC colors of faint galaxies to those measured in the COSMOS2015 catalog (Laigle et al. 2016, using the COSMOS2020 catalog instead results in insignificant differences).These zeropoint corrections are usually ≈ 0.1 mag and are fairly insensitive to the flux limits.After applying this procedure, we find that the average SEDs of the LATIS galaxies, as well as the derived distributions of stellar masses and SFRs (Section 2.2), are quite consistent across all fields.
We find that 93-97% of LATIS galaxies in the COS-MOS and D1 fields have clear matches in the COS-MOS2020 and our newly compiled (optical to NIR) D1 field catalogs.These matches are identified based on the criteria of having separations < 0. ′′ 5 and ∆r < 0.4 mag.The fraction is lower (60%) in the D4 field, because the WIRDS imaging does not cover the entire area.To ensure that the stellar mass measurements are robust, we require NIR coverage, which is particularly important in D1/D4 fields where the IRAC data are not as deep.Based on these criteria, we are left with approximately 3500 galaxies with 90% of the sample at 2.0 ≤ z ≤ 3.1 (Figure 1).

Stellar Mass & Star Formation Rates
We fit the UV to mid-infrared SED of galaxies to derive their physical parameters.We use the C++ version of the LePhare code (Arnouts et al. 1999;Ilbert et al. 2006) combined with a library of synthetic spectra generated by the Bruzual & Charlot (2003) population synthesis code.The redshifts are fixed to the spectroscopic redshifts from LATIS.Similar to the configuration employed in the COSMOS catalogs (Laigle et al. 2016;Weaver et al. 2022), the models incorporate exponentially declining star formation histories with nine efolding times in the range of 0.01 < τ < 30 Gyr and two delayed exponentially declining models (SFR ∝ te −t/τ ) with τ = 3 and 5 Gyrs.The delayed models are included since high-redshift galaxies likely exhibit star formation histories that differ substantially from simple exponential decays (e.g., Reddy et al. 2012;Kelson & Holden 2010).We adopt the Chabrier (2003) IMF, truncated at 0.1 and 100 M ⊙ , and the Calzetti et al. (2000) attenuation law to apply dust extinction (E(B − V) ≤ 1.1).The code also includes emission lines using the Kennicutt (1998) relation between SFR and UV luminosity, as described in Ilbert et al. (2009).Three different stellar metallicities are considered: Z * = 0.02, 0.008, and 0.004.
For each template, LePhare computes fluxes in all bands, then determines the template with a minimum χ 2 based on the model and observed fluxes.It also provides percentiles of the marginalized posterior probability distributions (probability ∝ e −χ 2 /2 ).In this work, we use the median posterior values for stellar mass, SFR, and sSFR (SFR/M * ).
Figure 2 shows the distribution of SED-derived stellar masses and SFRs.The main sequence (MS) relationship is also shown from Speagle et al. (2014) for comparison.Overall our sample is in agreement with the z ∼ 2.5 MS, although its average SFRs are slightly higher.We further compare to galaxies in the COSMOS2020 catalog, whose masses and star-formation rates were derived using similar techniques, thereby minimizing systematics.The green line in Figure 2 represents the median trend for the star-forming galaxies selected by At fixed stellar mass, our sample exhibits a ∼ 0.3 dex average enhancement in SFR compared to these galaxies.As a result, our sample predominantly represents star-forming galaxies located on and above the main sequence at cosmic noon.

Composite spectra
Most of the prominent absorption and emission features in the FUV are interstellar in origin.Typically, stellar features in the FUV spectrum are weak and require a high signal-to-noise ratio to be useful for measuring stellar metallicity.Gravitationally lensed systems often provide an adequate signal-to-noise ratio for individual galaxies.However, in our case, where individual galaxies lack sufficient signal-to-noise ratio, stacking techniques can be used to boost the signal strength.To construct a composite spectrum, the individual spectra are first corrected for Galactic extinction using the Schlafly & Finkbeiner (2011) dust map, followed by masking pixels with strong skyline contamination.The individual spectra are shifted to the rest frame using their spectroscopic redshifts and then normalized by their median fluxes within 1425-1500 Å to prevent the dominance of high-SFR galaxies in the composite spectrum.The normalized spectra are resampled onto a grid of wavelengths corresponding to ∆v = 100 km/s.This results in a wavelength spacing of ∼ 0.4 − 0.6 Å in a rest-frame FUV spectrum covering ∼ 1200 − 1800 Å.The composite spectrum is determined by the median of the normalized spectra, and the 1σ error is derived by bootstrapping. Figure 3 illustrates the high signal-to-noise composite spectrum of all 3491 galaxies with 10 9 M ⊙ ≤ M * ≤ 10 11.5 M ⊙ included in the present work (S/N=227 per 100 km/s at λ rest ∼ 1450 Å).

Stellar mass bins
To study the stellar MZR, we construct composite spectra within bins of stellar mass.The median uncertainty in stellar mass measurements is ∼ 0.15 dex, which provides a lower limit on the size of the stellar mass bins.Furthermore, we require that a composite spectrum has a median signal-to-noise ratio per 100 km/s of at least 30 over a wavelength range of 1221-1800 Å to ensure a reliable metallicity measurement.Over a similar wavelength range, we require that at least 25 galaxies contribute to every pixel of a composite spectrum in order to obtain reliable 1σ error estimates.These criteria leave us with nine stellar mass bins, which will be utilized to construct the stellar MZR in Section 5.2.

STELLAR POPULATION SYNTHESIS MODELS
To determine the stellar metallicity of galaxies, we compare the observed composite spectra with synthetic stellar population models.We use the models from the Binary Population And Spectral Synthesis code (BPASS v2.2.1;Eldridge et al. 2017;Stanway & Eldridge 2018) adopting the Chabrier ( 2003) IMF with a high-mass cutoff of 100 M ⊙ .The set of simple binary stellar population spectra (i.e., instantaneous starbursts) with an initially formed stellar mass of 10 6 M ⊙ at ages log(a ′ /years) = 6.0 − 11 in 0.1 dex increments is provided at a pixel resolution of 1 Å for 13 different stellar metallicities (10 −5 ≤ Z * ≤ 0.040).
We add nebular continuum to the stellar population models (BPASS+Nebular) using the photoionization code Cloudy v17.02 (Ferland et al. 2017).We model Hii regions with a spherical shell geometry of a fixed radius, adopting the BPASS spectrum as the incident spectrum.The metallicity of ionized gas is assumed to be Z neb = 0.5Z ⊙ , which is the average gas-phase metallicity of galaxies at the redshift interest of our study (e.g., Steidel et al. 2016;Strom et al. 2018;Sanders et al. 2021).Additionally, we use an ionization parameter of log U = −2.8 and electron density of n e = 250 cm −3 that are consistent with recent estimates for z ∼ 2 − 3 star-forming galaxies (e.g., Sanders et al. 2016;Strom et al. 2017;Topping et al. 2020).Due to the minimal contribution of the nebular continuum to the total flux (≲ 10%), our results are not sensitive to the assumptions of the Cloudy model.
We construct models assuming a constant star formation history by integrating burst models over their ages.The choice of constant star formation history is motivated by the fact that the FUV spectrum is mainly reliant on the past 100 Myr history of star formation, and the star formation history is expected to be relatively constant over this timescale for an ensemble of galaxies at z ∼ 2 − 3 (e.g., Reddy et al. 2012;Cullen et al. 2019).Therefore, model spectra, F(λ, a, Z * ), are constructed using the following equation: where BPASS+Nebular models (f BPASS ) are summed up over the age bins (a ′ ) such that N is the number of BPASS age bins within the age of a galaxy (a), ∆a ′ i are the widths of the age bins, and T (λ, a ′ i ) represents the transmission function due to dust attenuation.A population of stars younger than the lifetime of the stellar birth clouds suffers from greater attenuation than those outside the region (Charlot & Fall 2000).The typical lifetime of the birth clouds is 10 Myr (e.g., Blitz & Shu 1980).We therefore model the transmission function as follows (see also Carnall et al. 2022): where T ISM (λ) and T BC (λ) are the transmission functions of the ambient ISM and birth clouds, respectively.
In this paper, we do not intend to infer ISM dust attenuation parameters; they are merely nuisance parameters in our model and are coupled with possible flux calibration errors.To properly model any imperfections in spectroscopic calibration, a fifth-order polynomial function is used for the T ISM (λ) component; however, the T BC (λ) component is modeled using the Calzetti et al. (2000) dust attenuation, wherein the visual band attenuation A BC V is included as a free parameter.Figure 4 shows examples of the model spectra as the metallicity (Z * ) and A BC V are varied while the other parameters are fixed.It is evident that absorption features are much stronger as the stellar metallicity increases and that this has the largest effect on the spectrum.Moreover, the figure shows that although the FUV stellar absorption lines are relatively insensitive to the birth-cloud attenuation A BC V , the small variation in equivalent width of several features, including the 1718 Å and 1501 Å regions, can be used to capture any covariance between A BC V and Z * , given a high signal-to-noise spectrum.

FITTING SYNTHESIS MODEL SPECTRA
In this section, we present our Bayesian full spectrum fitting method, along with the associated parameters and priors.In order to accurately compare our    V at Z * = 0.005.For illustration purposes, we normalized all model spectra using the pseudo-continuum derived from the spline fit to the wavelength windows (black squares) suggested by Rix et al. (2004).model spectra with the observed composite spectra, it is important to ensure that the spectral resolution is consistent between the two.We first interpolate the BPASS+Nebular models to bins of 100 km/s, aligning them with the wavelength grid of the observed spectra.We then introduce a parameter, σ v , that accounts for the combined effects of galaxies' stellar velocity dispersion, the spectrograph resolution, and redshift uncertainties.We use estimates of these factors to place an informed prior on σ v while allowing for uncertainties in these estimates (In Appendix C we demonstrate the effect of adopting a fixed best estimate of σ v in our fitting).We convolve our model spectra with a Gaussian kernel in velocity space, resulting in model templates that we refer to as F(λ, Z * , a, σ v , A BC V ).We fit the BPASS+Nebular models ( F) to composite spectra within the framework of Bayesian inference.In the case of our problem, the parameter vector of the models is θ = (Z * , a, σ v , A BC V ), and assuming Gaussian and independent uncertainties in the composite spectra, the logarithm of the likelihood function can be constructed as, (3) where f (λ) and σ(λ) are the flux and 1σ uncertainty of the observed composite spectrum, respectively.Due to the fact that our models are provided for discrete values of the parameters, we linearly interpolate F models across different grids in order to define the likelihood at any given parameter value.Rather than treating T ISM as a free parameter in the model, we determine its value by fitting a fifth-order polynomial function to f (λ)/ F(λ, θ).This step is performed at every evaluation of the likelihood function during the Bayesian sampling process.This approach ensures that the model's continuum shape matches that of the observed spectra while keeping the computational complexity manageable.A BC V are given a uniform prior, but a logarithmic prior is applied for age and stellar metallicity.A Gaussian prior is adopted for σ v .In order to calculate this prior, we consider several factors: the spectral resolution of σ inst = 125 ± 10 km/s achieved for typical galaxies (representing the average and range over the rest-frame wavelength range of 1221-1800 Å, Section 2), random errors in redshifts of σ z = 93 ± 27 km/s (Newman et al. in prep), and a typical stellar velocity dispersion of σ * = 100±50 km/s (e.g., Barro et al. 2014).The BPASS models are distributed in wavelength bins of 1 Å, coarser than our pixels.Binning can be regarded as convolution by a top-hat kernel, which we approximate as a Gaussian with equal FWHM, leading to a model resolution of σ mod = 86 ± 10 km/s.Since σ v represents the kernel required to match the BPASS models to the observations, it can be computed as Combining our estimates of each term results in an effective resolution of σ v = 165 ± 35 km/s.Table 1 presents the list of the parameters and priors.In order to obtain the posterior distribution of the parameters, we use the MULTINEST nested sampling algorithm (Skilling 2006;Feroz & Hobson 2008;Feroz et al. 2009;Feroz & Skilling 2013;Feroz et al. 2019), which can be accessed via the PYMULTINEST interface (Buchner et al. 2014).
As shown in Figure 3, not every part of FUV spectra originates from stellar emission.Most of the prominent lines in a rest-UV spectrum are predominantly formed in the ISM and/or have a nebular origin.These lines are not included in the model spectra ( F).We therefore exclude the wavelength regions contaminated by interstellar absorption or nebular emission lines.The He ii λ1640 emission line is not excluded since this line originates primarily from hot stars that are well incorporated in BPASS models.Table 2 summarizes the masked regions ignored during the fitting.These masked wavelengths are also shown in Figure 3 as shaded regions.We note that our full spectral fitting procedure includes the 1370 and 1425 indices (Leitherer et al. 2001) which are attributed to Fe v λλ1360−1380 and O v λ1371, and to Si iii λ1417, C iii λ1427, and Fe v λ1430, respectively.The red lines in Figure 3 show some other stellar photosphere absorption lines within the wavelength range we considered for fitting (1221 − 1800 Å).
When conducting spectral fitting, we employ two iterations to address any potential discrepancies between the observed and model spectrum.During the first iteration, we exclusively mask out the regions where the ISM absorption lines are present, as specified in Table 2. Following the initial fit, we mask the pixels displaying 3σ deviations from the best-fit model and repeat the analysis.We find only a marginal change in our metallicity measurement during the second iteration.However, we find that this iterative process leads to more robust estimates of σ v that are closer to our prior expectation.As discussed in Section 2.3.1, we do not include λ > 1800 Å in our fitting procedure since we require the contribution of a minimum of 25 galaxies to every pixel of the composite spectrum in each stellar mass bin.Therefore, we do not use the 1935-2020 Å region known as the 1978 index (Rix et al. 2004), which is mainly sensitive to the iron abundance (Fe iii λλ1949 − 1966).We investigate the effect of excluding these regions, based on the composite spectrum of the entire sample, where there are ∼ 10× more galaxies than in the stellar mass bins, allowing us to perform spectral fitting over a wider range of wavelength, λ = 1221 − 2020 Å.We find that excluding the λ = 1800 − 2020 Å region from the full spectral fitting has a minimal effect on the inferred parameters.Full spectral fitting in the range of λ = 1221 − 2020 Å results in slightly higher Z * value, by ∼ 0.02 dex.

MEASUREMENTS OF STELLAR METALLICITY
We now use the framework built in Sections 3 and 4 to analyze composite LATIS spectra and infer their metallicities.We first consider the composite spectrum of the entire sample, and then turn to the composite spectra in stellar mass bins.5.1.Average metallicity of ∼ 10 10 M ⊙ star-forming galaxies at z ∼ 2.5 The signal-to-noise ratio of a composite spectrum can be maximized by stacking spectra of the whole sample, resulting in lower measurement uncertainties for the sample average metallicity.We fit the BPASS models to the composite spectrum of the full sample, which has a median stellar mass of M * ∼ 10 9.8 M ⊙ , and find that the average log(Z * /Z ⊙ ) is −0.87 +0.01 −0.01 .Figure 5 shows a corner plot demonstrating that all parameters are well constrained.Figure 6 shows that the models generally provide a good fit to the composite spectrum with root mean square (RMS) of residuals ∼ 2%.The RMS of the fractional noise for the composite spectrum is ∼ 1%.Given the very high signal-to-noise ratio of the full sample's composite spectrum, the RMS of the residuals is notably increased by model inaccuracies.Therefore, the error estimates of parameters derived by fitting the full composite spectrum may be underestimated, since they may not fully account for the model's limitations.Underestimation of errors is less of a concern for composite spectra in the stellar mass bins, where observational noise outweighs model inaccuracies.
The stellar photospheric absorption features in the rest-frame FUV spectrum are predominantly caused by transitions of highly ionized iron (e.g., Brandt et al. 1998).Therefore Z * can be translated to [Fe/H] by [Fe/H] ≈ log(Z * /Z ⊙ ) = −0.87+0.01 −0.01 .In this study, we measure the FUV-weighted stellar metallicity, which is expected to closely resemble the gas-phase metallicity measured in the ISM where recent star formation has occurred.In general, the FUV-weighted metallicity is expected to be marginally lower than the gas-phase metallicity by only ≲ 0.1 dex (Kashino et al. 2022).However, similar galaxies at this redshift have an average gas-phase oxygen abundance of [O/H] ∼ −0.3 +0.1 −0.1 , as derived from rest-frame optical strong lines (e.g., Erb et al. 2006;Steidel et al. 2014;Strom et al. 2018;Sanders et al. 2021).
To assess the gas-phase oxygen abundance within our sample, we measure [O/H] for the composite rest-frame optical spectra of 17 LATIS galaxies that were observed using Keck/MOSFIRE as part of the MOSDEF survey2 (Kriek et al. 2015).We ensure that these galaxies have Hα line detection with S/N ≥ 3 and are not flagged as AGNs based on X-ray emission or IRAC colors.Ad- ditionally, we require that log([Nii]λ6584/Hα) < −0.3 to exclude optical AGNs (e.g., Coil et al. 2015).The selected 17 galaxies have a median redshift of z ∼ 2.4 with a median stellar mass of 10 9.9 M ⊙ , which is comparable to the median mass of the entire sample.Additionally, the median of log(SFR[M ⊙ /yr]) for these 17 galaxies is 1.31, which is slightly lower than the 1.54 ob-served for the full sample.According to Sanders et al. (2018), this difference in SFR corresponds to a deviation of ∼ 0.03 dex in metallicity, which is within our measurement error, thereby ensuring the 17 galaxies are representative for purposes of metallicity comparison.We construct the composite spectrum and measure the flux of the emission lines following the method described in  This calibration relies on direct-method metallicities obtained from stacked spectra of local analogs of z ∼ 2 galaxies; as a result, there is no evaluation of its inherent scatter.We assume that the scatter is the same as Pettini & Pagel (2004) calibration, which is 0.18 dex.Since the composite spectrum is composed of 17 galaxies, the intrinsic error in the oxygen abundance of the composite spectrum reduces to 0.18/ √ 17 = 0.04 (e.g., Erb et al. 2006;Sanders et al. 2015).We add this in quadrature to the measurement uncertainty, yielding a final result of [O/H] = −0.30± 0.05.We emphasize that this error is purely statistical; the dominant systematic errors are discussed below.
There are substantial systematic uncertainties associated with gas-phase metallicity measurements.For instance, there is a systematic offset of ∼ 0.3 dex between "direct" electron-temperature-based and photoionization-model-based gas-phase metallicities (e.g., Kewley & Ellison 2008;Blanc et al. 2019).Furthermore, there is a 0.24 dex offset between gas-phase metallicities derived from collisionally excited lines and recombination lines (Steidel et al. 2016) in local H ii regions.These discrepancies are of particular importance when comparing gas-phase and stellar metallicity on an absolute scale.Consequently, the errors we have estimated for [O/Fe] represent measurement errors and do not encompass systematic errors.Taking into account the tendency of direct method gas-phase abundances to underestimate [O/H] (e.g., Cameron et al. 2023), it becomes plausible that the α-enhancement could be larger than our estimate.

Stellar Mass -Metallicity Relation
This section presents the relationship between stellar mass and stellar metallicity for our sample at z ∼ 2.5.The results of our full spectral fitting analysis on composite spectra in bins of stellar mass, as described in Section 2.3.1, are presented in Appendix A and summa- rized in Table 3.As shown in the left panel of Figure 8, the stellar M * − Z * relation exists at z ∼ 2.5, such that Z * increases with increasing stellar mass at low masses and flattens at M * ≳ 10 10.3 M ⊙ .This flattening has important implications for galactic winds in massive galaxies, which we will consider in Section 5.4.
Results from other studies at relatively similar redshifts are included in the left panel of Figure 8 for comparison.Our M * − Z * relation is in good agreement with that of zCOSMOS-deep survey from Kashino et al. (2022) at 1.6 < z < 3. Their sample shares a comparable range of redshifts with ours and was selected in a similar manner.Additionally, they employed stellar population models similar to ours.Therefore similar results are anticipated to be obtained.However, we find slightly lower metallicities compared to the VANDELS survey (Cullen et al. 2019).In our analysis, we utilize the BPASS v2.2.1 models, whereas, in Cullen et al. (2019), the WM-basic stellar population models were employed.They performed a comparative experiment using BPASS v2.1 models, which resulted in slightly lower metallicities by ∼ 0.09 dex.Therefore, it is recommended to adjust their measurements by ∼ 0.1 dex to ensure a meaningful comparison with our findings.After adjusting for this systematic offset, our result is in better agreement with Cullen et al. (2019).It should be highlighted that while our study focuses on z ∼ 2 − 3, the VANDELS survey explored a broader redshift range of z = 2.5 − 5.However, Cullen et al. (2019) did not observe significant redshift trends within this range.The lack of such trends suggests that the different redshift ranges between our studies may not be detrimental to our metallicity comparison.Moreover, our comparison supports their observation that stellar MZR does not evolve strongly over z ∼ 2 − 5.
Despite the generally satisfactory agreement with existing literature, we note that prior studies did not take into account the effect of extra attenuation in birth clouds, which may have an impact on the derived metallicity.The measurement of metallicity can be biased toward less obscured, lower metallicity environments if the extra attenuation present in the birth clouds of young stars is not taken into account (Chartab et al. 2022).To assess this effect, we repeat the analyses considering A BC V = 0 (Appendix B).We then determine the offset between the metallicity measurements and those obtained with additional dust attenuation considered for birth clouds.Our finding indicates that neglecting additional dust attenuation within birth clouds results in an underestimated stellar metallicity by ∼ 0.1 dex.
In the right panel of Figure 8, we compare our stellar MZR with the gas-phase M * − [O/H] relation from Sanders et al. (2021).We find that there is no clear trend of [O/Fe] with stellar mass, which is in agreement with previous studies (e.g., Strom et al. 2022;Kashino et al. 2022).Instead, the values remain relatively constant across all mass ranges, with an average value of [α/Fe] ∼ [O/Fe] ∼ 0.6 ± 0.1.This suggests that the galaxies in our sample are α-enhanced across the entire mass range.It is also worth noting that the average [α/Fe] value we find here is consistent with what we estimate for the composite spectrum of the entire sample in Section 5.1.
Previous studies of the gas-phase MZR at z ∼ 2 have not shown such a strong flattening at the massive end (right panel of Figure 8), although Sanders et al. (2021) reported signs of flattening for the most massive stellar mass bin at M * ∼ 10 10.5 M ⊙ .This difference could be partly explained by the smaller samples used in gas-phase studies, which are mostly restricted to M * < 10 10.5 M ⊙ , with only a handful of galaxies beyond this stellar mass.In contrast, our sample includes 219 galaxies in the last stellar mass bin of log(M * /M ⊙ ) = 10.4 − 11.5, and our overall larger sample also enables finer mass bins that better define the shape of the MZR.In stellar MZR studies, Cullen et al. (2019) analyzed a sample of M * ≲ 10 10 M ⊙ galaxies, the stellar mass range where we find no flattening consistent with their findings.The results of Kashino et al. (2022), although not explicitly stated in their paper, suggest signs of flattening at the massive end before a sharp rise (with large uncertainty) in their last stellar mass bin at M * ∼ 10 11 M ⊙ .
Moreover, the right panel of Figure 8 includes the stellar MZR at z ∼ 2 from Strom et al. (2022), who used a different method to constrain [Fe/H].They combined BPASS stellar population synthesis models with Cloudy photoionization models to infer [Fe/H] from the restoptical nebular spectra of galaxies in the Keck Baryonic  The number of galaxies in a stellar mass bin that is used to construct a composite spectrum.
b Median signal-to-noise ratio per pixel (100 km/s) over a wavelength range of 1221-1800 Å.
c The minimum number of galaxies contributing per pixel in the composite spectrum at the wavelength range of 1221-1800 Å.
Structure Survey (KBSS; Rudie et al. 2012;Steidel et al. 2014).We find slightly lower [Fe/H] along with a more complex mass dependence than the linear form used in the Strom et al. (2022) model, but the overall level agreement is still encouraging considering the very different data sets and methods underpinning these analyses.

UV Slope Dependence of the MZR
The flattening of the stellar MZR at high stellar masses observed in our study may have several possible explanations.One possibility is that in the highestmass bins, our sample is potentially biased against the most metal-rich galaxies.If such galaxies are also dusty and red, they may be underrepresented in our sample due to their faint FUV fluxes, which could lead to a flattening of the MZR at high masses.Another possibility is that the physical processes governing the stellar metallicity of galaxies change at high stellar masses.In this section, we explore the impact of dust on the MZR to assess the potential biases against red, dusty galaxies at high masses.In order to determine the UV continuum slope, denoted as β, we employ a power law (f λ ∝ λ β ) fit to the photometric data (outlined in Section 2.1) within the rest-frame range of λ = 1300 − 2600 Å.We require our galaxies to have at least three photometric measurements in this wavelength range.We do not use spectroscopic data to measure β since not all galaxies are covered in the desired wavelength range.Therefore, we rely on photometrically-measured β values to ensure consistent measurement across our sample.The median uncertainty in β is ∼ 0.14 (Figure 9).
To investigate the effect of dust content on the stellar MZR, we divide the sample into three subsamples corresponding to the tertiles of β. Figure 10 shows the stellar MZR for the subsamples representing the upper and lower tertiles of β.This comparison indicates no significant dependence of the MZR on β.If this independence holds over the full range of β present in a mass-limited sample, then the flattening of the MZR we observe cannot be attributed to color-related selection effects.However, we cannot guarantee that the independence of β and metallicity extends to extremely dusty massive galaxies.By examining Figure 9, we note that our sample closely follows the relation of the massselected star-forming sample in McLure et al. (2018) (green line) except for the very highest mass bin, where our sample is biased toward bluer galaxies.Nonetheless, at the turn-over stellar mass of the stellar MZR (M * ≈ 10 10.3 M ⊙ ), our sample does match with the trend reported by McLure et al. (2018).Therefore, we conclude that the observed flattening of the MZR is due to physical mechanisms rather than sample selection.

Mass-Dependent Trends in Galactic Winds
Analytic models provide valuable insight into the physical mechanisms governing the metallicity enrichment of galaxies.These models consider the interplay between key factors such as the inflow and outflow rates, SFR, the return of stellar material to the ISM, and changes in the total gas mass to unveil the complex mechanisms that shape the metallicity evolution of galaxies.When modeling the FUV-based stellar metallicity, which serves as a proxy for the iron abundance, it is important to incorporate the delay time distribution (DTD) of SNe Ia in order to track the evolution of iron-peak elements separately from the α-elements.The blue and red points represent the lower and upper tertiles of β subsamples, respectively.We find that there is no significant dependence of the MZR on β, suggesting that the flattening observed in the stellar MZR for high-mass galaxies is not heavily affected by selection effects related to color.
This distinction is crucial as iron-peak elements are produced by both CCSNe and SNe Ia, while α-elements are only produced by CCSNe.Incorporating the DTD for SNe Ia in the analysis allows for a more comprehensive understanding of the metallicity evolution in galaxies.
In this work, we adopt the one-zone model of Weinberg et al. (2017), which accounts for the DTD of SNe Ia.We can then constrain the mass loading factor (η = Ṁoutflow Ṁ * ) of z ∼ 2.5 galaxies based on the observed stellar MZR.
The Weinberg et al. (2017) model assumes that the star formation efficiency (SFE), which represents the ratio of the star formation to the gas mass available, remains constant over time ("linear Schmidt law").The DTD for enrichment from Type Ia SNe is modeled as an exponentially declining function, with a minimum delay time of 0.15 Gyr and an e-folding timescale of 1.5 Gyr.In this work, we adopt a simple constant star formation history and the fiducial values from Weinberg et al. (2017) for all free parameters of the model, except the mass loading parameter, SFE and CCSNe oxygen yield.We refer the reader to their paper for further details.
The adopted values for the iron and oxygen yields in Weinberg et al. (2017) result in a maximum [O/Fe] ∼ 0.4 for metal enrichment entirely from CCSNe.This value is smaller than that measured in z ∼ 2.5 galaxies (Section 5.2; Steidel et al. 2016), and it is also exceed by lowmetallicity stars in the Milky Way.For instance, Figure 11 shows that the abundances of APOGEE (DR17; Abdurro'uf et al. 2022) stars reach up to [O/Fe] ∼ 0.6 (see also Steidel et al. 2016), which is closer to the theoretical yield calculations of Nomoto et al. (2006).Therefore, we have adjusted the IMF-integrated oxygen yield to from 0.015 to 0.021 in order to achieve a plateau at [O/Fe] ≈ 0.6.It is important to highlight that this adjustment does not affect our main result, trends in the mass loading factor, which is based on iron abundances.
To estimate the SFE timescale, also referred to as the depletion timescale (M gas /SFR), we utilize the scaling relation introduced by Tacconi et al. (2020).They found that the depletion timescale depends mainly on the redshift and the offset from the main sequence.We compare our sample with a sample of star-forming galaxies (sSFR> 10 −10.1 yr −1 ; Pacifici et al. 2016) at z = 2 − 3 in COSMOS2020 catalog and find an average offset of ∼ 0.3 dex from the star-forming main sequence for our sample.By employing Equation 4 of Tacconi et al. (2020), we derive an average depletion timescale of ∼ 0.3 Gyr for our sample at z ∼ 2.5.
To explore the sensitivity of the analytic model to the mass-loading parameter η, we plot tracks of retain such high [α/Fe] values.This is consistent with the median age of our sample inferred from SED fitting in Section 2.2.
In the following, we aim to infer the mass loading parameter η and its dependence on stellar mass M * using the observed stellar MZR.We incorporate the M *dependent ages of galaxies into the chemical evolution modeling, estimated as the median of the characteristic star formation timescale, defined as τ * SF = M * /SFR, at a given stellar mass (Figure 12).Along with the other parameters described above, we can now use the Weinberg et al. (2017) model to compute [Fe/H] given η and M *3 .We utilize a cubic functional form for the mass loading parameter η(M * ) and determine the polynomial coefficients by fitting the model to the observed [Fe/H]-M * relation.The derived mass loading parameter, presented in Figure 13, indicates that it declines as the stellar mass increases, but flattens and possibly even increases for M * ≳ 10 10.5 M ⊙ .
We note that the absolute value of η is indeed sensitive to the various assumptions in the model, yet the overall trend remains consistent.Blue dashed and dotted lines in Figure 13 illustrate the trend for alternative mass-independent values of age, 200 Myr and 500 Myr, respectively.In all cases, a decline in η is observed as the stellar mass increases for low-mass galaxies; however, this trend plateaus/reverses at M * ∼ 10 10.5 M ⊙ .Sanders et al. (2021) derived the mass-dependence of η by applying a different chemical evolution modeling framework to observations of the gas-phase MZR at z = 2.3 and 3.3.They found that the metal mass-loading factor ζ(M * ), which differs from η by the ratio of the outflow and ISM metallicities, follows ζ(M * ) ∝ M −0.35±0.02* at low masses, and they showed that their data were compatible both with models in which ζ flattens at high masses M * ≳ 10 10.5 M ⊙ and models in which ζ continues to decline.We find that η ∝ M −0.35±0.01* at low masses M * ≲ 10 10.25 M ⊙ , in excellent agreement with Sanders et al. (2021), but our data clearly prefer a flattening at high masses.Our estimates of η are systematically higher than those of Sanders et al. (2021) by ∼ 0.3 dex, but we expect that the absolute value of η is significantly more sensitive to modeling assumptions than the slope of η(M * ).
To better understand the behavior of the mass loading parameter as a function of stellar mass, we compare our results with those obtained from the TNG50 cosmological simulation (orange dashed line in Figure 13).TNG50 is a cosmological hydrodynamic simulation within a 50 Mpc box sampled by 2160 3 gas cells, resulting in a baryon mass resolution of 8 × 10 4 M ⊙ and an average spatial resolution ∼ 100 − 200 pc in ISM gas.Nelson et al. (2019) measured the mass loading parameter at z ∼ 2 in TNG50 and found a similar nonmonotonic behavior as a function of galaxy stellar mass.Specifi-cally, they found that the mass loading parameter turns over and rises rapidly above 10 10.5 M ⊙ , consistent with our finding.In the simulation, this behavior is traced to winds driven by supermassive black holes, which dominate over stellar-driven winds in massive galaxies.

DISCUSSION
Our results show that the stellar MZR exists at z ∼ 2.5, such that stellar metallicity increases with increasing stellar mass at low masses.However, we also find evidence of flattening at the massive end of the MZR (M * ≳ 10 10.3 M ⊙ ).
In Figure 14, we compare stellar MZR measurements at z ∼ 0 with our result at z ∼ 2.5.Using full spectral fitting to the rest-frame optical spectra, Zahid et al. (2017) measured the stellar MZR for ∼ 200000 star-forming galaxies in the Sloan Digital Sky Survey (SDSS).Their result is in agreement with the relation derived by Kudritzki et al. (2016) based on analysis of individual supergiant stars.The results of Zahid et al. (2017) differ from those of Gallazzi et al. (2005), possibly due to the fact that the latter used Lick indices to measure metallicity and included data from both starforming and quiescent galaxies in their analysis.The slope of our stellar MZR at z ∼ 2.5 below the turnover mass (M * ∼ 10 10.3 M ⊙ ) is 0.40 ± 0.04, which agrees with the slope meausured at z ∼ 0 by Zahid et al. (2017).However our z ∼ 2.5 MZR shows a shift toward lower metallicities by ∼ 0.7 dex compared to the local stellar MZR.We note that the actual offset is greater, as z ∼ 0 metallicities in Zahid et al. (2017) are weighted by optical luminosity, and when converted to FUV-weighted for comparison with our measurements, they need to be shifted upwards by ∼ 0.1 dex to account for the differences in the stellar populations being traced (Kashino et al. 2022).Therefore, the average stellar metallicity of galaxies at fixed stellar mass increases by a factor of ∼ 5 from z ∼ 2.5 to 0 (see also Cullen et al. 2019).However, Sanders et al. (2021) reported that the gas-phase metallicity [O/H] only increased by a factor of 2 in the same time interval, leading to a discrepancy between the evolution of the FUV-weighted stellar metallicities and gas-phase oxygen abundances.While it is expected that these two metallicities evolve somewhat differently due to the different timescales that they trace, they should still closely follow each other with a difference of ∼ 0.1 − 0.2 dex (Kashino et al. 2022).Thus, the large observed difference in their evolution is likely due to delayed iron enrichment by SN Ia.
We note that our sample is primarily positioned on or above the star formation main sequence, with a median offset of ∼ 0.3 dex.At low redshifts, SFR is anticor- related with gas-phase metallicity at fixed stellar mass (e.g., Mannucci et al. 2010;Curti et al. 2020).Although the existence and form of the "fundamental metallicity relation" at higher redshifts is debated, Sanders et al. (2018) found that, at a fixed stellar mass, the relationship ∆ log(O/H) ∼ −0.15 × ∆ log(SFR) holds for starforming galaxies at z ∼ 2.3.If a similar relation applies to the iron abundance of young stars, then our estimates of the stellar metallicity would be ∼0.05dex lower than that of a mass-limited sample.This is a minor potential bias and is not likely to affect our main results on the shape of the stellar MZR.
Our MZR observations imply a transition in the characteristics of galactic winds around masses of M * ≈ 10 10.3 M ⊙ , but they cannot uniquely identify the physical origin of this transition.For further insights we turn to simulations and other observations.We find that our M * − η relation, derived using an analytic chemical evo-lution model, matches closely the one measured in the TNG50 simulation at z = 2.In the simulation, the mass loading attributed to star formation declines monotonically with stellar mass, and the upturn in η at high masses M * ≳ 10 10.5 M ⊙ is caused by winds driven by supermassive black holes (SMBHs; Nelson et al. 2019).In the TNG model, low-luminosity AGN in this population drive high-velocity winds that expel cool, dense gas from the nuclear ISM.These outflows eventually lead to the "inside-out" quenching of star formation.
Furthermore, independent observations support the widespread occurrence of AGN-driven outflows in a similar population of z ∼ 1-2 galaxies.Genzel et al. (2014) and Förster Schreiber et al. (2014) find a rapid increase at masses M * ≳ 10 10.9 M ⊙ of the incidence of broad, nuclear emission lines with line ratios indicative of shocks or AGN photoionization.This mass scale is higher than the mass at which we find a turnover of the MZR.How- ever, it is possible that AGN are driving winds in slightly lower-mass galaxies that are sufficient to affect the chemical evolution but are not as readily detectable via optical emission line diagnostics.This possibility might be tested by observations using the higher angular resolution afforded by extremely large telescopes.Taken together, these additional lines of evidence suggest that SMBH-driven winds in massive galaxies are the most likely origin of the flattening of the stellar MZR that we observe.

SUMMARY
In this paper, we present the stellar mass -stellar metallicity relation for 3491 star-forming galaxies at 2 ≲ z ≲ 3 using rest-frame FUV spectra from the LATIS survey.We utilize the BPASS (v2.2.1) synthesis models to fit high signal-to-noise composite spectra of galaxies.Our findings can be summarized as follows: • We find that the stellar metallicity increases monotonically with increasing stellar mass at lower masses but flattens at M * ≳ 10 10.3 M ⊙ .The slope of our z ∼ 2.5 MZR at the low mass end is similar to that of local stellar MZR.However, there is a significant offset of ∼ 0.7 dex toward lower metallicities compared to the local stellar MZR.
• We compare our measurements of stellar metallicity to the gas-phase oxygen abundance of galaxies at similar redshift.We find no clear trend in the [O/Fe] -M * relation.[O/Fe] remains relatively constant across the full mass rangewith an average value of [O/Fe] ∼ 0.6, suggesting that young galaxies at z ∼ 2.5 have yet to experience substantial iron enrichment resulting from SN Ia.
• Using an analytical model of chemical evolution, we constrain the mass loading parameter as a function of stellar mass.We find that the mass loading parameter decreases with stellar mass at low masses but plateaus or reverses at M * ∼ 10 10.5 M ⊙ .In combination with the TNG simulations and other observations, our results suggest that SMBH-driven outflows in massive galaxies at z ∼ 2.5 are responsible for the upturn in the M * −η relation, which in turn explains the flattening of the MZR at the massive end.Massive galaxies undergo strong SMBH-driven outflows, which remove the metal-rich gas from the ISM.
We thank the anonymous referee for providing insightful comments and suggestions that improved the quality of this work.This paper includes data gathered with the 6.5 meter Magellan Telescopes located at Las Campanas Observatory, Chile.We thank the staff at Las Campanas Observatory for their dedication and support.Software: NumPy (Harris et al. 2020), Matplotlib (Hunter2007),Astropy(AstropyCollaborationetal.2013, 2018, 2022), PyMultiNest (Buchner et al. 2014).

Figure 1 .
Figure 1.Spectroscopic redshift distribution for the entire LATIS sample.The vertical dashed lines indicate the median value along with the 5 th and 95 th percentiles.

Figure 2 .
Figure 2. SED-derived SFR as a function of stellar mass for the entire LATIS sample.For comparison, the mainsequence relation from Speagle et al. (2014) is also shown.The green line represents the median trend for the starforming galaxies at z = 2 − 3 in the COSMOS2020 catalog.sSFR> 10 −10.1 yr −1(Pacifici et al. 2016) and z = 2 − 3.At fixed stellar mass, our sample exhibits a ∼ 0.3 dex average enhancement in SFR compared to these galaxies.As a result, our sample predominantly represents star-forming galaxies located on and above the main sequence at cosmic noon.

Figure 3 .
Figure3.Composite rest-frame FUV spectrum of the entire LATIS sample.The top panel shows the number of galaxies contributing to each pixel in the composite spectrum, ranging from 536 to 3235.The composite spectrum is shown in the bottom panel; some prominent emission and absorption lines are labeled.Stellar and interstellar absorption features are colorcoded with red and green, respectively.Several nebular emission lines (blue), and fine structure emission lines (purple) are also included.The red shaded regions indicate wavelength ranges for 1370 and 1425 indices as defined inLeitherer et al. (2001); these indices cover several photospheric lines such as Fe v λλ1360 − 1380, O v λ1371, Si iii λ1417, C iii λ1427, and Fe v λ1430, which are not labeled for legibility.The gray line shows the 1σ error spectrum.

Figure 4 .
Figure 4.The BPASS v2.2.1 stellar population models at the average effective resolution of the LATIS composite spectra (σv = 250 km/s) assuming a constant SFR over 100 Myr.The top panel shows the variation of the model with metallicity at a fixed A BC V = 0. Various metallicities, ranging from 10 −4 to 0.01, are shown with different colors, as indicated in the legend.The bottom panel shows the sensitivity of the model to A BCV at Z * = 0.005.For illustration purposes, we normalized all model spectra using the pseudo-continuum derived from the spline fit to the wavelength windows (black squares) suggested byRix et al. (2004).

Figure 5 .
Figure 5. Corner plot showing the posterior distribution for fitted parameters on the composite spectrum of the entire LATIS sample.Contours correspond to 1σ, 2σ, 3σ and 4σ levels.The medians of the marginalized posteriors are indicated by red dashed lines.The black dashed lines on marginalized posteriors indicate the median values along with the 16 th and 84 th percentiles.

Figure 6 .
Figure 6.Similar to Figure 3 but the red line represents the best-fit model, while the orange line shows the unsmoothed BPASS model with its original resolution.The shaded regions in light blue, affected by interstellar absorption or nebular emission lines, are excluded from the fitting process.The bottom panel shows the fit residuals.Chartab et al. (2021).The resulting composite spectrum is shown in Figure7.We determine ⟨[Nii]λ6584/Hα⟩ = 0.13 ± 0.02, which corresponds to a gas-phase oxygen abundance of 12 + log(O/H) = 8.39 ± 0.03 ([O/H] = −0.30±0.03)using the calibrations ofBian et al. (2018).This calibration relies on direct-method metallicities obtained from stacked spectra of local analogs of z ∼ 2 galaxies; as a result, there is no evaluation of its inherent scatter.We assume that the scatter is the same asPettini & Pagel (2004) calibration, which is 0.18 dex.Since the composite spectrum is composed of 17 galaxies, the intrinsic error in the oxygen abundance of the composite spectrum reduces to 0.18/ √ 17 = 0.04 (e.g.,Erb et al. 2006;Sanders et al. 2015).We add this in quadrature to the measurement uncertainty, yielding a final result of [O/H] = −0.30± 0.05.We emphasize that this error is purely statistical; the dominant systematic errors are discussed below.Oxygen is the most abundant among the α elements.Therefore, a measurement of [O/Fe] can be taken as an approximation of [α/Fe].In this study, we find [Fe/H] = −0.87 ± 0.01 and [O/H] = −0.30± 0.05, resulting in [α/Fe] ∼ [O/Fe] = 0.57 ± 0.05 (statistical errors only).This result is consistent with the average α-enhancement reported for star-forming galaxies at z ∼ 2.5 (e.g.,Steidel et al. 2016;Strom et al. 2018;Cullen et al. 2019;Kashino et al. 2022), suggesting that

Figure 7 .
Figure 7. Composite spectrum of 17 LATIS galaxies in the rest-frame optical, obtained from Keck/MOSFIRE observations of the MOSDEF survey.Shaded regions indicate the corresponding errors obtained from bootstrapping.

Figure 8 .
Figure 8. Left: The stellar MZR for our sample of star-forming galaxies at z ∼ The data points represent the posterior median metallicity derived from composite spectra in bins of stellar mass, as described in Section 4. The error bars represent the 1σ uncertainty.For comparison, the MZR from previous studies at similar redshifts by Cullen et al. (2019) and Kashino et al. (2022) are shown in green and blue, respectively.Right: Comparison of our stellar [Fe/H] with the gas-phase [O/H] at z ∼ 2 measured by Sanders et al. (2021) and Strom et al. (2022).The blue dashed line indicates the stellar mass -stellar [Fe/H] relation inferred by Strom et al. 2022 based on photoionization modeling of rest-frame nebular spectra.

Figure 9 .Figure 10 .
Figure 9. UV spectral slope (β) plotted against stellar mass for the LATIS sample (blue points), with red diamonds indicating the median β in each stellar mass bin.The green line shows the M * − β relation from McLure et al. (2018) at 2 < z < 3 for comparison.The median uncertainty in the β measurements is displayed in the lower right corner.

Figure 11 .
Figure 11.Chemical evolution tracks in the [O/Fe] -[Fe/H] plane for different values of η.The different sizes of circles represent the ages of the galaxies along each track.The red square in the figure denotes our measurements for the entire LATIS sample at z ∼ 2.5.The gray density map illustrates the distribution of APOGEE (DR17) stars without a STAR BAD flag and with reliable abundance ratio values (O FE FLAG == 0 and FE H FLAG == 0).

Figure 12 .
Figure 12.The characteristic star formation timescale (τ * SF = M * /SFR) is shown as a function of stellar mass for our sample.The red diamonds represent the median values of τ * SF for each stellar mass bin, while the black line illustrates the best fit line to the median relationship.

Figure 13 .
Figure 13.The mass loading factor as a function of stellar mass at z ∼ 2.5, derived from the best-fit analytic model to our observed MZR.The relationship between the mass loading factor and the stellar mass is modeled as a cubic function.The shaded region in the main figure represents the 1σ error around the best fit.The blue dashed and dotted lines indicate the trends for constant, mass-independent age values of 200 and 500 Myr, respectively.The orange dashed line shows the trend from TNG50 simulation (Nelson et al. 2019) for outflows with a radial velocity greater than zero.The inset figure shows the best model MZR (blue line) along with the observed MZR (red diamonds).

Figure 14 .
Figure14.Comparison of our stellar MZR at z ∼ 2.5 with those ofZahid et al. (2017),Kudritzki et al. (2016) andGallazzi et al. (2005) measured at z ∼ 0. Our z ∼ 2.5 MZR exhibits a significant shift toward lower metallicities compared to the z ∼ 0 relation, while the slopes of the two relations are consistent.
N.C. and A.B.N. acknowledge support from the National Science Foundation under Grant No. 2108014.G.A.B. acknowledges the support from the ANID Basal project FB210003.

Figure C. 1 .
Figure C.1.Similar to the left panel of Figure 8, but with an orange curve representing the MZR when σv is fixed.

Table 1 .
Parameters and priors used for the model fitting

Table 3 .
Properties of composite spectra and measurements of stellar metallicity