The Gas-phase Mass–Metallicity Relation for Massive Galaxies at z ∼ 0.7 with the LEGA-C Survey

The massive end of the gas-phase mass–metallicity relation (MZR) is a sensitive probe of active galactic nuclei (AGN) feedback that is a crucial but highly uncertain component of galaxy evolution models. In this paper, we extend the z ∼ 0.7 MZR by ∼0.5 dex up to log (M ⋆/M ⊙) ∼ 11.1. We use extremely deep VLT VIMOS spectra from the Large Early Galaxy Astrophysics Census (LEGA-C) survey to measure metallicities for 145 galaxies. The LEGA-C MZR matches the normalization of the z ∼ 0.8 DEEP2 MZR where they overlap, so we combine the two to create an MZR spanning from 9.3 to 11.1 log (M ⋆/M ⊙). The LEGA-C+DEEP2 MZR at z ∼ 0.7 is offset to slightly lower metallicities (0.05–0.13 dex) than the z ∼ 0 MZR, but it otherwise mirrors the established power-law rise at low/intermediate stellar masses and asymptotic flattening at high stellar masses. We compare the LEGA-C+DEEP2 MZR to the MZR from two cosmological simulations (IllustrisTNG and SIMBA), which predict qualitatively different metallicity trends for high-mass galaxies. This comparison highlights that our extended MZR provides a crucial observational constraint for galaxy evolution models in a mass regime where the MZR is very sensitive to choices about the implementation of AGN feedback.


INTRODUCTION
The mass-metallicity relation (MZR), the tight correlation between galaxy stellar mass and gas-phase oxygen abundance, is sensitive to the physical processes that govern the baryon cycle and shape galaxy evolution.Infalling gas fuels star formation and triggers the production of metals (Pagel 1997).Feedback from star formation can then drive large scale outflows that eject gas and metals from a galaxy into the intergalactic or circumgalactic media (see Veilleux et al. 2005 andTumlinson et al. 2017 andreferences therein), where it may eventually reaccrete onto the galaxy (Oppenheimer et al. 2010).
At high masses, the steep fall off of the stellar mass function strongly implies a disruption of the baryon cycle and the quenching of star formation (Di Matteo et al. 2005).Many galaxy evolution models rely on some form of AGN feedback to shut off star formation (e.g., Springel et al. 2005;Somerville & Davé 2015), but the dominant mechanism-such as thermal feedback that heats and ejects gas from the interstellar medium (Schaye et al. 2015) or radio-mode feedback that prevents cooling from the circumgalactic medium (Kereš et al. 2005)-remains difficult to pin down with current observations.However, differences in AGN feedback prescriptions in cosmological simulations result in significantly different shapes of the high mass end of the MZR at intermediate redshifts (Davé et al. 2017;Torrey et al. 2019), making it a key observational constraint for differentiating between AGN feedback mechanisms.
The MZR has been widely studied in the local universe (Tremonti et al. 2004;Andrews & Martini 2013;Curti et al. 2020), at intermediate redshifts (z = 0.5 − 1.5; e.g., Pérez-Montero et al. 2009;Zahid et al. 2011;Topping et al. 2021), beyond cosmic noon with ground-based spectroscopic samples (z = 1.5 − 3.5, e.g., Maiolino et al. 2008;Shapley et al. 2017;Sanders et al. 2018Sanders et al. , 2020;;Topping et al. 2021;Sanders et al. 2021), and recently extending up to z ∼ 8 with the advent of JWST (Langeroodi et al. 2022;Curti et al. 2023;Nakajima et al. 2023).In the local universe, metallicity increases with stellar mass according to a power law as 12 + log(O/H) ∝ M 1/3 below a characteristic turnover mass (∼ 10 10.25 M ), and then the MZR flattens as metallicities of the most massive galaxies asymptote to a constant value.Generally, metallicity decreases at fixed stellar mass with increasing redshift.The shape and normalization of the MZR evolve with redshift (e.g., Lilly et al. 2003, Moustakas et al. 2011, Sanders et al. 2020, Sanders et al. 2021), with the metallicities of lower mass galaxies evolving more dramatically across cosmic time (Zahid et al. 2013).However, above the turnover mass (log M /M 10.25), the MZR has only been reliably measured locally due to the difficulty of detecting the required emission lines in metal-rich galaxies at higher redshifts.
In this work, we take advantage of the ultra-deep spectra of a nearly mass-complete sample at 0.6 < z < 1.0 obtained by the Large Early Galaxy Astrophysics Census (LEGA-C van der Wel et al. 2016;Straatman et al. 2018) survey.The primary science driver of the LEGA-C survey was the study of stellar populations and dynamics of massive galaxies at z ∼ 0.8, so the spectra have a much higher signal-to-noise ratio (S/N) per pixel than typical surveys of emission line galaxies.LEGA-C's high S/N enables detection of even weak emission lines, significantly mitigating the observational bias towards objects with bright emission lines and systematically lower metallicities.Despite the extremely long integration time required to achieve the high S/N per pixel, the LEGA-C survey is large enough to construct a statistically robust MZR up to 10 11.1 M -high enough to constrain the effects of AGN feedback on the MZR.
The paper is structured as follows: §2 describes the LEGA-C data set, sample selection methods, and Bayesian metallicity estimation method; §3 presents the MZR at z ∼ 0.7, and §4 compares the LEGA-C MZR with cosmological hydrodynamic simulations.Through-out this paper we assume the following cosmology: H 0 = 70 km s −1 Mpc −1 , Ω M = 0.3, and Ω Λ = 0.7.

The LEGA-C Survey
The LEGA-C survey used a K -band selection to target a nearly mass-selected sample of 3,741 high-mass (M 10 10 M ) galaxies with existing Ultra-VISTA photometry at z = 0.6 − 1.The spectra were obtained with 20 hours of integration time each on the Very Large Telescope (VLT) with the VIsible MultiObject Spectrograph (VIMOS Le Fèvre et al. 2003) and resulted in spectra over the wavelength range 5800−9500 Å with an average S/N of 20 Å−1 .Details of the survey design and sample selection can be found in van der Wel et al. ( 2016), Straatman et al. (2018), andvan der Wel et al. (2021).LEGA-C's K -band selection and high S/N is well-suited for measuring metallicity even for objects with weak emission lines and significant stellar Balmer absorption lines that would otherwise render the measurement of weak Balmer lines challenging.This opens up novel parameter space including the high mass, high metallicity, and low specific star formation rate (sSFR) regime.
We adopt the stellar masses and line flux measurements from the LEGA-C catalog (DR3; van der Wel et al. 2022).Stellar masses were determined from Ultra-VISTA photometry (Muzzin et al. 2013a, Muzzin et al. 2013b) using the FAST code (Kriek et al. 2009) and Bruzual & Charlot (2003) stellar population libraries, assuming a Chabrier (2003) initial mass function and Calzetti et al. (2000) dust extinction law, with an exponentially declining star formation rate (SFR).If we had chosen to adopt the MAGPHYS stellar masses (da Cunha et al. 2008, as measured for LEGA-C galaxies by de Graaff et al. 2021), this would systematically increase stellar masses by ∼ 0.2 dex.
In order to properly account for absorption underlying the Balmer emission lines used in this analysis, the stellar continuum must be modeled and subtracted.The stellar continua were fit using pPXF (Cappellari & Emsellem 2004, Cappellari 2012, Cappellari 2017) with MILES empirical stellar population templates (Vazdekis et al. 2010).For further details see van der Wel et al. (2016) and Bezanson et al. (2018).We use emission line fluxes from continuum-subtracted spectra (v3.11) fit by Maseda et al. (2021), which used Platefit (Tremonti et al. 2004;Brinchmann et al. 2004;Gunawardhana et al. 2020) to simultaneously fit emission lines as Gaussians.These Gaussians were fixed to the LEGA-C spectroscopic redshifts but allowed a ±300 km s −1 velocity offset between the Balmer lines and forbidden lines such Gray points are those with S/N < 3, which are excluded from the sample.The remaining points follow the same color scheme as the right panel, described below.We do not find a trend of Hγ S/N with stellar mass.Only a small fraction of these objects have Hγ S/N < 3, which is required for dereddening.Right panel: The same sample on the mass-excitation (MEx) plot of Juneau et al. (2011), which plots R3 against stellar mass and separates star-forming galaxies and AGN based on their line ratios.The red points show objects that are determined to be MEx AGN by the cut of Juneau et al. (2011).Orange points show objects determined to be other forms of AGN, such as broad-line, radio, IR, and through the "blue diagram" of Lamareille et al. (2009).Green crosses show objects with line fluxes dominated by DIG, and blue points show star-forming objects with line fluxes dominated by HII regions.
as [O ii]λ3727.Uncertainties on the line fluxes were determined by fitting to fluxes perturbed by wavelengthdependent noise vectors.The standard deviation of these fits was combined with the uncertainty in the overall flux calibration to give a line flux uncertainty.

Sample Selection
We require [O ii]λ3727 and [O iii]λ5007 for our analysis, which are only simultaneously accessible for a small subset of the spectroscopic sample.The high spectral resolution of LEGA-C results in rather limited wavelength coverage.Of the ∼ 3500 galaxies in the primary LEGA-C sample, only 425 galaxies include both spectral features.This selection limits the redshift range of the sample analyzed to 0.568 < z < 0.894.We further limit our sample to star-forming galaxies using UVJ restframe colors, adopting the Muzzin et al. (2013a) colorcolor cuts designed for the Ultra-VISTA data set, leaving 245 objects.Finally, because we utilize Hβ and Hγ to deredden line fluxes, we require a S/N 3σ in each line flux.Since our Bayesian metallicity estimation method naturally handles upper limits on [O ii]λ3727 and [O iii]λ5007 line fluxes, we do not apply a S/N cut on them to avoid biasing our sample.However, we visually inspect the spectra that fall into this category (18 objects), rejecting half of them for incomplete line coverage.This leaves 211 galaxies.
The left panel of Fig. 1 shows the Hγ S/N as a function of stellar mass, with a horizontal line at S/N = 3 showing our selection.The gray points are those with Hγ S/N < 3; the remaining points follow the same coloring convention as the right panel of the figure, described below.After we select star-forming galaxies, we are limited by our ability to measure the Balmer lines, because we do not place an S/N cut on [O ii]λ3727 or [O iii]λ5007.We do not find a trend of Hγ S/N with mass, and we therefore reach a high level of completeness.We measure each emission line, on average, 92% of the time.Since the LEGA-C survey was not designed to only target galaxies with the brightest emission lines, its extreme depth enables analysis of a much more complete set of star-forming galaxies than previous surveys.The high Hγ detection fraction allows us to use the Balmer decrement to deredden spectra, instead of needing to assume an A V (see Savaglio et al. 2005).For sample selection purposes only, we dereddened the line fluxes using the Hγ/Hβ Balmer decrement (e.g., equation A14 from Momcheva et al. 2013).However, when we measure our metallicity, we simultaneously solve for 12 + log(O/H), A V , and dereddened Hβ flux, in both cases adopting a Cardelli et al. (1989) extinction law with R V = 3.1.
Gas-phase metallicities cannot be derived from the line ratios of AGN as the emission lines are not thermally excited, so they must be excluded.We are unable to use the Baldwin-Phillips-Terlevich diagnostic diagram (BPT; Baldwin et al. 1981) to remove AGN from our sample as we do not have coverage of [N ii]λ6583 or Hα.Instead, we identify AGN using the mass-excitation diagram (Juneau et al. 2011), shown in the right panel of Fig. 1.There are four populations on this panel.The red points denote objects that are shown to be mass-excitation AGN by Juneau et al. (2011).The orange points denote objects that are determined to be AGN by other means (mid-IR diagnostics, X-ray fluxes, radio emission, and broad emission lines; S. Vervalcke priv.communication).We also used cuts from the R3 vs. R2 diagram (see Section 2.4 of Lamareille 2010) that only requires lines at the bluer end of the optical range to remove additional likely AGN (the LINER, Seyfert 2, and star-forming/Seyfert 2 classes).Finally, we remove galaxies whose emission is dominated by diffuse ionized gas (DIG) via the Σ Hα cut suggested by Zhang et al. (2017), and retain galaxies with log(Σ Hα [erg s −1 kpc −2 ]) < 39.These are the objects shown as green crosses in Fig. 1.For the purposes here, we calculate Hα via the intrinsic ratio between dereddened Hβ and Hα.The removal of these objects has a minimal effect on the best-fit MZR.The blue points show galaxies that are star-forming and make up the final sample used in the rest of this paper.
Of the 211 galaxies with coverage of the emission lines of interest, 52 are identified as AGN (orange and red points in Fig. 1), and 159 are identified as star-forming.Of these 159 star-forming galaxies, 145 are found to have line fluxes dominated by HII regions.This latter group constituting the sample used in the rest of the analysis.The median redshift of our final sample is z = 0.70, with stellar masses from log(M /M ) = 9.81 to 11.1.The median log(SFR [M yr −1 ]) = 1.1, ranging from −0.14 to 1.89.The median log(sSFR [yr −1 ]) = −9.3,ranging from −10.4 to −8.6.and SDSS (gray contours) samples.There is a slight offset in the LEGA-C data to higher line ratios for R2, R3, and R23, but not for O32.We note that the median line ratio uncertainties are smaller than the symbol size if reddening uncertainties are not included.Note the smooth trend in line ratios and the lack of sharp features at log(M /M ) > 10.5 that would be expected if AGN feedback dramatically alters the gas-phase metallicity at high mass.

Emission Line Ratios
It is useful to compare emission line ratios between studies as a first step, given that they are not affected by the choice of metallicity calibration.Fig. 2 shows line ratios versus stellar mass for four common metallicitysensitive diagnostics (R2, R3, R23, and O32) for LEGA-C galaxies whose line fluxes are dominated by HII regions.The LEGA-C line ratios are shown in blue circles with the blue line and shaded region indicating the running median and one-sigma scatter, respectively.Line flux ratios from the Sloan Digital Sky Survey (SDSS, York et al. 2000) MPA-JHU catalog (Kauffmann et al. 2003, Brinchmann et al. 2004, Tremonti et al. 2004) from the local universe are shown as gray contours.The LEGA-C trends generally agree with the SDSS data with a slight offset to higher R2, R3, and R23 in agreement with the more detailed analysis of LEGA-C emission line ratios of Helton et al. (2022).The similarity between the line ratios of LEGA-C and SDSS galaxies motivates our adoption of the Curti et al. (2017) strong line metallicity calibration that is calibrated to direct method metallicities of stacked SDSS spectra.Furthermore, smoothness of the trends and the absence of sharp features at log(M /M ) > 10.5 implies that the MZR should be smooth in this regime regardless of the metallicity calibration-in contrast to some simulation pre- dictions with certain AGN feedback prescriptions (see Section 4 for more discussion).

Measuring Gas-Phase Metallicities
We estimate metallicities using a Bayesian framework (see Wang et al. 2017, Wang et al. 2019) to account for covariance between the dust correction and metallicity.We use the emcee package (Foreman-Mackey et al. 2013) to run a Markov Chain Monte Carlo (MCMC) that samples combinations of metallicity, nebular dust extinction, and dereddened Hβ flux (f dered Hβ ), which can be used to generate [O ii]λ3727, Hγ, Hβ, and [O iii]λ5007 line fluxes given a metallicity calibration (Curti et al. (2017) in our case) and assuming case B recombination for the Balmer line ratio.We adopt a flat prior on metallicity (12 + log(O/H) = [7.6,9.3]) 1 and on nebular extinction (A V = [0, 4]).We use the Jeffreys' prior distribution for f dered Hβ and limit its range to [0 to 10 5 ] in units of 10 −19 erg s −1 cm −2 .The χ 2 statistic for our MCMC is given by where f dered i is the dereddened line flux of emission line i, σ f dered i is its measured uncertainty, R i is the ratio of the dereddened flux of emission line i to the dereddened Hβ flux, and σ Ri is the intrinsic scatter in the R i -metallicity relation from Curti et al. (2017) (where σ Ri = 0 for the Balmer lines because we assume no uncertainty in the case B recombination Hβ/Hγ line ratio).
With our four emission lines ([O ii]λ3727, Hγ, Hβ, and [O iii]λ5007), it is possible to use R2 and R3 as the line ratios of interest (R i in Equation 2) or R23 and O32.We 1 Formally, our metallicity range extends beyond the calibrated maximum metallicity of the Curti et al. ( 2017) calibrations (12 + log(O/H) = 8.85), so as to not impose an artificial ceiling on metallicity.We treat objects with metallicities well-beyond this value with caution, which happens to be only a single object in our sample, whose inclusion does not affect the shape of the MZR or any of our main conclusions.
found similar results for both approaches, but the R2-R3 combination produced slightly smaller uncertainties, so we adopt those results as our fiducial measurements going forward.The dereddened line fluxes, A V , and metallicities for individual galaxies are provided in Table 1.
As a check on the Bayesian approach, we also calculated metallicities using the more common approach of first determining the reddening from the Balmer decrement and then estimating metallicity from the dereddened line flux ratios.We found similar median estimated metallicities (average discrepancy of 0.03 dex), but underestimated uncertainties (average discrepancy of 0.05 dex) using the Kobulnicky & Kewley (2004) and Curti et al. (2017) R23 calibrations (assuming the upper branch of R23 and after adjusting to the same calibration scale using the transformations from Teimoorinia et al. 2021) compared to our Bayesian approach.
3. MASS-METALLICITY RELATION AT z ∼ 0.7 Fig. 3 shows the LEGA-C MZR in comparison to other z ∼ 0.7 MZRs and the z ∼ 0 SDSS MZR.The left panel of Fig. 3 shows individual LEGA-C galaxies as blue circles, the running median (with a window size of 0.5 dex) as a blue line, and the 16-84th percentile region as a blue band.It also shows the Zahid et al. (2011) MZR at z ∼ 0.8 (red) using data from DEEP2 (Davis et al. 2003, Newman et al. 2013) and the Curti et al. (2020) z ∼ 0 MZR for SDSS galaxies (gray). 2 The LEGA-C MZR matches the DEEP2 MZR almost perfectly in the mass range probed by both surveys (log(M /M ) = 10.0-10.6)and smoothly extends an additional 0.5 dex in stellar mass.We find that  2019) MZR from stacks of eBOSS galaxies (orange) at z ∼ 0.7 is shown in orange with the most massive stellar mass bin indicated with a dotted line due to its low S/N (as noted by Huang et al. 2019).Overall, the eBOSS MZR agrees well with the LEGA-C+DEEP2 MZR where they overlap.
the LEGA-C MZR has slightly smaller scatter (0.123 dex) than the DEEP2 MZR (0.225 dex).The LEGA-C and DEEP2 objects show significant overlap in their sSFR distributions, suggesting that these two samples are probing similar populations of star-forming galaxies though at somewhat different mass ranges.Given these similarities, we do a combined fit to the LEGA-C and DEEP2 MZRs to form a single z ∼ 0.7 MZR that spans 1.7 dex in stellar mass (blue line in the right panel of Fig. 3).
The LEGA-C, DEEP2, and combined LEGA-C+DEEP2 MZRs agree well with the Huang et al. ( 2019) z ∼ 0.7 MZR3 (orange line in the right panel of Fig. 3) that was constructed using stacked spectra from the eBOSS survey (Dawson et al. 2016).The slight turnover at the highest mass point could be due to AGN contamination in the stack due to the difficulty of identifying AGN in individual eBOSS spectra.With the LEGA-C data, we found that including AGN would introduce a downturn in the MZR at the highest masses, hence the need for removing them.
To determine the best fit MZR across the full mass range by combining LEGA-C and DEEP2 MZR, we use the LMFIT package (Newville et al. 2014) and the functional form given in Equation 2 of Curti et al. (2020): where M 0 is the characteristic turnover mass, Z 0 is the saturation metallicity that the MZR asymptotically approaches, γ is the power law index of the MZR below M 0 , and β characterizes how quickly the power law begins to flatten.We use a bin size of 0.16 dex to smooth over small-scale variations while retaining the large-scale trend.
The uncertainties on the combined fit were found by resampling the median metallicities according to their uncertainties and refitting, but the uncertainties on the fit are not shown for clarity because they are on the order of a few hundredths of a dex at all masses.The parameters of the best fit LEGA-C+DEEP2 z ∼ 0.7 MZR are Z 0 = 8.74 ± 0.01, log(M 0 /M ) = 10.13 ± 0.05, γ = 0.30 ± 0.02, and β = 1.99 ± 0.05.The LEGA-C+DEEP2 MZR lies at an average −0.07 dex vertical offset from the SDSS MZR, with a −0.05 dex offset at the high mass end and a −0.13 dex offset at the low mass end.The values of log(M 0 /M ) and γ are within 1σ of the SDSS MZR.
The offset of the LEGA-C+DEEP2 MZR from the z ∼ 0 MZR (∼ −0.1 dex below the turnover and ∼ −0.05 dex above it) is significantly less than offsets found by previous studies at similar redshifts (e.g., Pérez-Montero et al. 2009;Lamareille et al. 2009).The offsets in these studies of ∼ −0.2 dex could be due to a bias towards galaxies with stronger emission lines, their smaller sample sizes, or lower masses.Both of the aforementioned studies used data from VVDS (Le Fèvre et al. 2005), which also used the VIMOS spectrograph but before it was upgraded to have better red sensitivity.
Finally, we checked to see if accounting for the secondorder dependence of metallicity on SFR (in addition to stellar mass) would further reduce the scatter in metallicity, but we did not find evidence for this so-called fundamental metallicity relation (FMR; Mannucci et al. 2010;Lara-López et al. 2010) in our LEGA-C sample.The FMR has been found to hold out to z ∼ 3.3 (Sanders et al. 2021), but all MZRs at z ∼ 0.7 and beyond are at lower masses than LEGA-C.The FMR is difficult to constrain at the high-mass end because there is little dynamic range in SFR, and the trend may reverse depending on the adopted metallicity calibration (Curti et al. 2020), so perhaps our finding is not too surprising.Further complicating the issue is that LEGA-C does not cover Hα, which would provide the most robust short timescale SFR measurements.Given that the measured LEGA-C SFRs probe both longer timescales and different regions of the galaxy than metallicity indicators, we expect that this added scatter could erase additional correlations.We tested whether the scatter in the LEGA-C MZR correlates with other SFR indicators (e.g., UV+IR, spectral energy distribution fitting) but did not find evidence for a clear fundamental metallicity relation.

DISCUSSION AND CONCLUSIONS
The LEGA-C MZR probes a very interesting stellar mass regime at z ∼ 0.7.The high-mass end of the MZR is important to constrain because of its sensitivity to AGN feedback in cosmological hydrodynamical simulations (De Rossi et al. 2017, Davé et al. 2019, Torrey et al. 2019), which are typically tuned to match the evolving stellar mass function (Weinberger et al. 2018 et al. 2020et al. , Zinger et al. 2020)).Because the exact implementation and strength of AGN feedback is hotly debated (Davé et al. 2017, Torrey et al. 2019), the observed MZR provides an important constraint and can be used to calibrate simulations.
In Fig. 4 we compare the LEGA-C+DEEP2 MZR (blue) to the MZRs from the IllustrisTNG (orange; Torrey et al. 2019) and SIMBA (green; Davé et al. 2019) simulations at similar redshifts.We focus mainly on the difference in the shape of the MZRs due to uncertainties in the overall normalization of stellar masses, metallicities, and stellar yields, so we normalized the metallicities of the IllustrisTNG and SIMBA MZRs to match the LEGA-C+DEEP2 MZR at log(M /M ) = 10.4.The IllustrisTNG MZR is shifted upwards 0.02 dex, the SIMBA MZR is shifted upwards 0.15 dex.The LEGA-C+DEEP2 MZR flattens at log(O/H) ∼ 8.75 and qual-itatively disagrees with both simulations at the highmass end-the range of the MZR where LEGA-C becomes crucial.The SIMBA MZR uses an sSFR cut of sSF R > −0.38Gyr −1 at z = 0.7.The IllustrisTNG MZR uses mass-weighted metallicity as opposed to SFRweighted metallicity to avoid highlighting exclusively star-forming gas (Torrey et al. 2019).
The IllustrisTNG simulation predicts a downturn in metallicity at log(M /M ) = 10.6, in sharp contrast to the observed LEGA-C MZR, which flattens but does not decrease.This downturn, found at z = 1 and z = 2 but not at z = 0, is driven by a transition from thermal to kinetic AGN feedback in the models (Weinberger et al. 2018) that removes large quantities of metal-rich gas from the centers of galaxies.Given the disagreement with the measured MZR, perhaps the ejective feedback within IllustrisTNG is removing too much metal-rich gas.
On the other hand, the SIMBA simulation predicts that the metallicity continues to increase up to log(M /M ) = 11 before flattening.Davé et al. (2017) speculate that the precursor of SIMBA simulation (i.e., MUFASA) may underestimate the metalloading of winds or overestimate wind recycling in highmass galaxies, and this problem may also be present in the SIMBA simulations.The median SIMBA MZR at z ∼ 0.7 (shown here) is robust to the SFR cut used to select star-forming galaxies.The SIMBA MZR is in closer agreement with the LEGA-C+DEEP2 MZR at the high-mass end than the MUFASA z ∼ 1 MZR due to SIMBA's on-the-fly dust model, suggesting that the treatment of dust in simulations is important for re-producing the observed MZR (Davé et al. 2019).It is also possible that SIMBA's implementation of quenching by heating halo gas to stop accretion (mimicking radio-mode AGN feedback as opposed to the ejective AGN feedback of IllustrisTNG) may result in effectively closed-box evolution to high metallicities.
Measuring the massive end of the MZR at z ∼ 0.7 is a significant step towards understanding how AGN feedback affects galaxies.Further observation of LEGA-C galaxies with PFS (Takada et al. 2014;Greene et al. 2022) to measure Hα and [N ii]λ6583 would allow AGN removal through the classic BPT diagram, as well as calculations of metallicity through other calibrations.Additionally, constraining the high-mass end of the MZR at z ∼ 1 − 2 with PFS and MOONRISE (Maiolino et al. 2020) is a natural extension of this work and will help to more concretely define the evolution of the MZR at high masses.At even higher redshifts, JWST can probe the metallicities of the galaxies that will evolve into high mass objects at low redshift.Finally, more detailed comparisons to cosmological simulations and theory will help improve our understanding of AGN feedback.

Figure 1 .
Figure1.Left panel: Hγ S/N for non-quiescent LEGA-C objects as a function of stellar mass.Gray points are those with S/N < 3, which are excluded from the sample.The remaining points follow the same color scheme as the right panel, described below.We do not find a trend of Hγ S/N with stellar mass.Only a small fraction of these objects have Hγ S/N < 3, which is required for dereddening.Right panel: The same sample on the mass-excitation (MEx) plot ofJuneau et al. (2011), which plots R3 against stellar mass and separates star-forming galaxies and AGN based on their line ratios.The red points show objects that are determined to be MEx AGN by the cut ofJuneau et al. (2011).Orange points show objects determined to be other forms of AGN, such as broad-line, radio, IR, and through the "blue diagram" ofLamareille et al. (2009).Green crosses show objects with line fluxes dominated by DIG, and blue points show star-forming objects with line fluxes dominated by HII regions.

Figure 2 .
Figure2.R2, R3, R23, and O32 versus stellar mass for the LEGA-C sample (blue circles with the running median and 16-84th percentiles as blue lines and bands, respectively) and SDSS (gray contours) samples.There is a slight offset in the LEGA-C data to higher line ratios for R2, R3, and R23, but not for O32.We note that the median line ratio uncertainties are smaller than the symbol size if reddening uncertainties are not included.Note the smooth trend in line ratios and the lack of sharp features at log(M /M ) > 10.5 that would be expected if AGN feedback dramatically alters the gas-phase metallicity at high mass.

Figure 3 .
Figure 3. Left: The LEGA-C mass-metallicity relation with individual measurements shown as blue circles with error bars indicating the inner 68 percentile highest probability density interval; the median and the 16-84th percentile range are shown as a blue line and band, respectively.The z ∼ 0 SDSS MZR (Curti et al. 2020) is shown in gray for comparison.The LEGA-C MZR shows a −0.05 dex offset from the SDSS MZR at the high mass end, and a −0.07 dex offset at the low mass end.The z ∼ 0.8 DEEP2 MZR (Zahid et al. 2011) (red line with red error band showing the 16-84th percentile region) agrees well with the LEGA-C MZR where they overlap in mass.Right: the LEGA-C and DEEP2 running median MZRs are shown as blue and red circles with error bars indicating the uncertainty on the median.The combined LEGA-C+DEEP2 MZR is shown as a blue line.The Huang et al. (2019) MZR from stacks of eBOSS galaxies (orange) at z ∼ 0.7 is shown in orange with the most massive stellar mass bin indicated with a dotted line due to its low S/N (as noted byHuang et al. 2019).Overall, the eBOSS MZR agrees well with the LEGA-C+DEEP2 MZR where they overlap.

Figure 4 .
Figure4.A comparison between our observed LEGA-C+DEEP2 MZR (blue) and MZRs from the IllustrisTNG simulation (orange;Torrey et al. 2019) and the SIMBA simulation (green;Davé et al. 2019) that shows the extreme sensitivity of the high mass end of the MZR to the impleof AGN feedback in simulations.Given the large uncertainty in the overall normalization of the metallicity scale, we renormalized the simulated MZRs to match the LEGA-C+DEEP2 MZR at log(M /M ) = 10.4.The Il-lustrisTNG MZR shows a significant decrease in metallicity above log(M /M ) = 10.6 in contrast to the LEGA-C MZR, likely due to AGN feedback ejecting too much metal-rich gas from central regions.The SIMBA MZR, on the other hand, continues to increase to log(M /M ) = 11, which may be caused by insufficient metal-loading of winds, excessive recycling of winds, or radio-mode AGN feedback that results in closed-box-like metallicity evolution.

Table 1 .
A snippet of the results from our Bayesian metallicity estimation method.Columns, from left to right, are the galaxy's LEGA-C ID and slit mask, spectroscopic redshift, stellar mass, dereddened [O ii]λ3727 and [O iii]λ5007 flux relative to Hβ, nebular extinction AV , and gas-phase metallicity, the latter four columns provide the inner 68 percentile highest probability density interval from the MCMC posterior distribution.The full table is available in machine-readable format.