CO Emission, Molecular Gas, and Metallicity in Main-Sequence Star-Forming Galaxies at $z\sim2.3$

We present observations of CO(3-2) in 13 main-sequence $z=2.0-2.5$ star-forming galaxies at $\log(M_*/M_{\odot})=10.2-10.6$ that span a wide range in metallicity (O/H) based on rest-optical spectroscopy. We find that CO(3-2)/SFR decreases with decreasing metallicity, implying that the CO luminosity per unit gas mass is lower in low-metallicity galaxies at $z\sim2$. We constrain the CO-to-H$_2$ conversion factor ($\alpha_{\text{CO}}$) and find that $\alpha_{\text{CO}}$ inversely correlates with metallicity at $z\sim2$. We derive molecular gas masses ($M_{\text{mol}}$) and characterize the relations among $M_*$, SFR, $M_{\text{mol}}$, and metallicity. At $z\sim2$, $M_{\text{mol}}$ increases and molecular gas fraction ($M_{\text{mol}}$/$M_*$) decrease with increasing $M_*$, with a significant secondary dependence on SFR. Galaxies at $z\sim2$ lie on a near-linear molecular KS law that is well-described by a constant depletion time of 700 Myr. We find that the scatter about the mean SFR-$M_*$, O/H-$M_*$, and $M_{\text{mol}}$-$M_*$ relations is correlated such that, at fixed $M_*$, $z\sim2$ galaxies with larger $M_{\text{mol}}$ have higher SFR and lower O/H. We thus confirm the existence of a fundamental metallicity relation at $z\sim2$ where O/H is inversely correlated with both SFR and $M_{\text{mol}}$ at fixed $M_*$. These results suggest that the scatter of the $z\sim2$ star-forming main sequence, mass-metallicity relation, and $M_{\text{mol}}$-$M_*$ relation are primarily driven by stochastic variations in gas inflow rates. We place constraints on the mass loading of galactic outflows and perform a metal budget analysis, finding that massive $z\sim2$ star-forming galaxies retain only 30% of metals produced, implying that a large mass of metals resides in the circumgalactic medium.


INTRODUCTION
The cold gas mass and abundance of heavy elements (i.e., metallicity) are fundamental properties of galaxies that are key to understanding galaxy formation and evolution. Molecular gas clouds are the sites of star formation such that galaxy star-formation rates (SFRs) email: rlsand@ucdavis.edu * The data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. † NHFP Hubble Fellow depend on the mass of molecular gas (M mol ) available in the interstellar medium (ISM). This close tie is manifested in a tight correlation between the surface densities of SFR (Σ SFR ) and molecular gas mass (Σ mol ) known as the molecular Kennicutt-Schmidt relation (.e.g, Kennicutt 1998;Bigiel et al. 2008;Leroy et al. 2008;de los Reyes & Kennicutt 2019;Kennicutt & De Los Reyes 2021). Likewise, metallicity is connected to star formation as new metals are produced via nucleosynthesis and returned to the ISM through the processes of stellar evolution and death, such that the trend of increasing metallicity over time traces the buildup of galaxy stellar mass (M * ).
Galaxy-integrated SFR, metallicity, and M mol have been found to depend strongly on M * in star-forming populations, resulting in a number of scaling relations: the star-forming main sequence (MS) in which SFR increases with increasing M * (e.g., Brinchmann et al. 2004;Noeske et al. 2007); the mass-metallicity relation (MZR) in which ISM metallicity, traced by the gas-phase oxygen abundance, increases with increasing M * (e.g., Tremonti et al. 2004;Andrews & Martini 2013;Curti et al. 2020); and the M mol -M * relation in which more massive galaxies have larger molecular gas reservoirs (e.g., Bothwell et al. 2014;Saintonge et al. 2017). At z ∼ 0, these scaling relations display secondary dependences that connect these properties in 3-dimensional parameter spaces. At fixed M * , O/H decreases with increasing SFR, forming the SFR−O/H−M * "Fundamental Metallicity Relation" (SFR-FMR; e.g., Mannucci et al. 2010;Cresci et al. 2019;Curti et al. 2020). Similarly, O/H and M mol are inversely related at fixed M * in the M mol −O/H−M * relation, or "Gas-FMR" (Bothwell et al. 2016b,a). A Gas-FMR has also been found at z ∼ 0 using atomic hydrogen gas masses (Bothwell et al. 2013;Hughes et al. 2013;Lara-Lopez et al. 2013;Brown et al. 2018). Since SFR is determined by the amount of cold gas present, it is likely that the Gas-FMR is a more fundamental relation from which the SFR-FMR emerges. These higherorder relations are realized as correlated scatter around pairs of scaling relations: the SFR-FMR is represented by an anti-correlation between residuals around the MS and residuals around the MZR, while a signature of the Gas-FMR is that residuals around the M mol -M * relation are anti-correlated with those about the MZR.
Characterizing the interrelation between M * , SFR, M mol , and metallicity is of central importance to understanding the cycle of baryons that governs galaxy growth. The class of gas-regulator, equilibrium, or bathtub models of galaxy formation indicate that the gas fraction (M gas /M * ) and ISM metallicity of a galaxy is governed by the rates of gas accretion (Ṁ in ) and outflow (Ṁ out ) relative to the SFR (e.g., Peeples & Shankar 2011;Davé et al. 2012;Lilly et al. 2013;Peng & Maiolino 2014), such that the mass-loading factor of outflows (η out =Ṁ in /SFR) can be constrained using measurements of both metallicity and gas fraction. Such models suggest that the SFR-FMR and Gas-FMR arise as a direct result of baryon cycling. Freshly accreted unenriched gas will expand the gas reservoir leading to a larger SFR, while simultaneously diluting metals in the ISM leading to lower O/H. The SFR-FMR and Gas-FMR are also ubiquitous features of numerical simulations of galaxy formation in a cosmological context that include feedback (e.g., Davé et al. 2017Davé et al. , 2019De Rossi et al. 2017;Torrey et al. 2018Torrey et al. , 2019. The existence of both a SFR-FMR and Gas-FMR, and the associated cor-related residuals around scaling relations, is thus a signature of a self-regulated baryon cycle governing galaxy growth. Searching for the SFR-FMR and Gas-FMR at high redshift presents an opportunity to understand baryon cycling during an epoch when galaxy formation was proceeding rapidly, when gas inflows and outflows are expected to be more intense than locally on average. There has been great progress in characterizing metallicity scaling relations at high redshift in the past several years thanks to extensive rest-optical spectroscopic surveys of representative star-forming galaxies at z ∼ 1 − 3 (e.g., Steidel et al. 2014;Kriek et al. 2015;Momcheva et al. 2016;Kashino et al. 2019). The MZR has been found to exist out to z ∼ 3.5, evolving toward lower O/H at fixed M * with increasing redshift (e.g., Erb et al. 2006a;Maiolino et al. 2008;Troncoso et al. 2014;Cullen et al. 2014Cullen et al. , 2021Sanders et al. 2015Sanders et al. , 2021Topping et al. 2021). The existence of a SFR-FMR has also been confirmed at z = 1.5 − 2.5 (Zahid et al. 2014;Sanders et al. 2018Sanders et al. , 2021Henry et al. 2021). M mol is commonly inferred from indirect tracers including CO line emission and dust continuum emission. The former requires a conversion factor between CO luminosity and M mol (α CO ; e.g., Wolfire et al. 2010;Schruba et al. 2012;Bolatto et al. 2013;Accurso et al. 2017), while the latter requires an assumed dust-to-gas ratio (e.g., Sandstrom et al. 2013;De Vis et al. 2019) or an empirical calibration between Rayleigh-Jeans tail dust emission and M mol (Scoville et al. 2016). These conversion factors have been found to depend strongly on ISM metallicity in the local universe, such that knowledge of the metallicity is key to accurate inferences of M mol . Progress in understanding the cold gas properties of high-redshift galaxies has been challenging due to the difficulty of detecting these tracers at cosmological redshifts.
CO line emission is the most common tracer of molecular gas in the nearby universe, but has proven difficult to measure in high-redshift main-sequence galaxies. Early work at z > 1 was limited to extreme sources including submillimeter galaxies (SMGs) and ultraluminous infrared galaxies (ULIRGs) (e.g., Greve et al. 2005;Tacconi et al. 2008), and more typical galaxies could only be reached with strong gravitational lensing Coppin et al. 2007;Danielson et al. 2011;Saintonge et al. 2013). The sensitive IRAM NOEMA interferometer enabled targeted surveys to reach mainsequence galaxies at z ∼ 1 − 2 Magnelli et al. 2012), culminating with the PHIBSS survey that measured CO for ∼ 50 main-sequence galaxies at z = 1 − 2.5 .
More recently, the Atacama Large Millimeter Array (ALMA) has provided further improvement in sensitivity at millimeter wavelengths. Using ALMA, the AS-PECS blind spectral scan survey has detected CO(3−2) and CO(2−1) for dozens of main-sequence galaxies at z = 1 − 3 in the Hubble Ultra Deep Field Boogaard et al. 2019), while targeted campaigns measured CO for a number of other z > 1 galaxies (e.g., Silverman et al. 2015Silverman et al. , 2018. Main-sequence CO(1−0) detections at z = 2 − 3 have also been obtained with deep VLA observations (Pavesi et al. 2018;Riechers et al. 2020). However, these existing samples of main-sequence CO-detected galaxies at high redshift still only probe the most massive mainsequence galaxies, with nearly all targets at z > 2 having stellar masses above 10 10.7 M . This fact highlights an important disconnect between high-redshift galaxies with cold gas information and the large samples at z ∼ 1 − 3 with detailed rest-optical spectroscopy spanning log(M * /M ) = 9.0 − 10.5, for which robust metallicity constraints are available. It is necessary to obtain CO observations of lower-mass galaxies that are more typical of high-redshift star-forming populations and critically have metallicity measurements in order to search for a Gas-FMR at high redshift, improve constraints on baryon cycling in early galaxies, and understand how to reliably estimate M mol in main-sequence high-redshift galaxies when CO and dust tracers are not available.
In this work, we analyze CO(3−2) observations obtained with ALMA for 13 near main-sequence galaxies at z ∼ 2.3 with a mean stellar mass of 10 10.4 M . These targets also have deep rest-optical spectroscopy from the MOSDEF survey , providing robust determinations of gas-phase metallicity and SFR. This unique combination of measurements allows us to carry out a systematic analysis of the relations among M * , SFR, M mol , and metallicity for typical starforming galaxies at z ∼ 2, and place constraints on gas flows and baryon cycling in these early galaxies.
This paper is organized as follows. In Section 2, we present the observations and describe how derived properties are calculated. The results are presented in Section 3. The sample properties are described in Sec. 3.1. Empirical relations between CO(3−2) luminosity and galaxy properties are explored in Sec. 3.2. We place constraints on the relation between α CO and O/H at z ∼ 2 in Sec. 3.3. In Sec. 3.4, we derive molecular gas masses and characterize the relationships between M * , SFR, M mol , and metallicity. We discuss these results in Section 4, and summarize our conclusions in Section 5.
Throughout, we assume a standard ΛCDM cosmology with H 0 =70 km s −1 Mpc −1 , Ω m =0.3, and Ω Λ =0.7. Magnitudes are in the AB system (Oke & Gunn 1983) and wavelengths are given in air. Stellar masses and SFRs are on the Chabrier (2003) initial mass function (IMF) scale. The term metallicity refers to the gas-phase oxygen abundance unless otherwise stated, where solar metallicity is 12+log(O/H) = 8.69 (Asplund et al. 2021). Molecular gas masses include a 36% contribution from helium and metals.

Sample Selection and Rest-Optical Spectroscopic Observations
Our sample was selected from the MOSFIRE Deep Evolution Field (MOSDEF) survey , a deep near-infrared (rest-frame optical) spectroscopic survey of ∼1500 galaxies at 1.4 ≤ z ≤ 3.8 using the MOSFIRE instrument (McLean et al. 2012) on the 10 m Keck I telescope. Survey targets were selected in the five CANDELS extragalactic legacy fields (Grogin et al. 2011;Koekemoer et al. 2011) from the photometric catalogs of the 3D-HST survey (Brammer et al. 2012;Skelton et al. 2014;Momcheva et al. 2016) down to fixed H-band magnitude measured from HST /WFC3 F160W imaging in three redshift intervals, with limits of H AB ≤ 24.0, 24.5, and 25.0 for targets at 1.37 ≤ z ≤ 1.70, 2.09 ≤ z ≤ 2.61, and 2.95 ≤ z ≤ 3.80, respectively. These redshift intervals are chosen such that the strong rest-frame optical emission lines [O ii]λλ3726,3729,Hβ,[O iii]λλ4959,5007,Hα,[N ii]λ6584, and [S ii]λλ6716,6731 fall in windows of atmospheric transmission in the near-infrared. Targets were selected based on pre-existing spectroscopic or HST grism redshifts when available, and photometric redshifts otherwise. The completed MOSDEF survey measured robust redshifts for ∼ 1300 targeted galaxies. See Kriek et al. (2015) for a full description of the survey design, observations, and data reduction.
Using the MOSDEF catalogs, we selected a sample of 14 star-forming galaxies in the middle redshift interval at z ∼ 2.3 where the CO(3−2) emission line can be observed in Band 3 with ALMA. Targets were required to be in the COSMOS field for accessibility with ALMA, have spectroscopic redshifts in the range 2.0 ≤ z ≤ 2.6, have detections at S/N≥3 for the Hβ, [O iii]λ5007,Hα,and [N ii]λ6584 lines to ensure robust metallicity and SFR constraints, and have stellar masses within ±0.1 dex of log(M * /M ) = 10.5 1 . Active galactic nuclei (AGN) were removed from the sample, having been identified based on their X-ray and IR properties (Coil et al. 2015;Azadi et al. 2017Azadi et al. , 2018Leung et al. 2019) and when log([N ii]/Hα) > −0.3. This selection yielded a sample of 14 star-forming galaxies at z = 2.08 − 2.47 for followup ALMA observations, the properties of which are presented in Table 1.

ALMA Observations and CO(3−2) Measurements
These 14 targets were observed in ALMA Cycle 6 Program 2018.1.01128.S (PI: R. Sanders) with the Band 3 receiver in 21 scheduling blocks over 7 Nov. to 29 Dec. 2018. The observations were designed such that the CO(3−2) line (ν rest = 345.796 GHz) fell within one spectral window in Band 3, between 99 and 113 GHz at the redshifts of our targets. The requested spectral configuration provided a bandpass of 1.875 GHz per spectral window with a native resolution of 7.8125 MHz, corresponding to ≈ 22 km s −1 at z ∼ 2.3. Observations were carried out with 43 12 m antennas in the C43-3 to C43-5 configurations, providing beam major axes of 0.7 −2.5 with a mean beam size of 1.5 . On-source integration times were 10 to 170 minutes per target reaching sensitivies of 86 to 456 µJy/beam integrated over 50 km s −1 .
The ALMA data were calibrated and imaged using the Common Astronomy Software Applications (CASA) package (McMullin et al. 2007). We utlized flux, phase, and bandpass calibrated measurement sets provided by the North American ALMA Regional Center at NRAO. The CASA task tclean was used to produce clean imaged datacubes with a velocity resolution of 50 km s −1 (≈ 1/3 of the rest-optical line FWHM; Price et al. 2016) using natural weighting. We searched for detections by collapsing the channels within ±200 km s −1 of the systemic redshift determined from the rest-optical lines and searching for positive peaks in these moment-0 maps exceeding 3σ at the spatial location of each target galaxy. Figure 1 shows CO(3−2) contours superposed on HST false-color images for our sample (top panel of each subplot). We found CO(3−2) emission is robustly detected in 6 targets with peak S/N>5 in the moment-0 maps, tentatively detected in two targets with peak S/N=3−4, and not detected in the remaining 6 targets.
1 The stellar masses used to select the sample observed with ALMA were not corrected for the contribution of strong emission lines to the broadband photometry. The stellar masses used in this paper have been derived from emission-line corrected photometry . Accordingly, the final stellar mass range differs slightly from that originally used in target selection.
For all detected sources, we perform optimal extractions to produce one-dimensional science spectra (Horne 1986) after converting the datacube intensity units from Jy/beam to Jy/pixel using the beam size in each frequency channel. All but one source (ID 3324) have spatial profiles that are consistent with being unresolved. We extract these sources using a 2-dimensional Gaussian profile set by the beam shape. ID 3324 presents an unsual spatial profile with two peaks separated by ≈1 beam width (1.5 ). Spectra extracted separately at each peak location yield 1D lines with roughly comparable fluxes and velocity centroids that are consistent with one another and the rest-optical redshift at less than the 1σ level. Given the relatively small area probed by our observations, the probability of a S/N>5 noise fluctuation located spatially and spectrally close to our targets is small. However, we cannot rule out the possibility that the second peak of ID 3324 is in fact spurious. We proceed by extracting the total line flux for ID 3324 using a spatial profile of two point sources aligned with the two peak locations, but note that our results in Section 3 do not significantly change if we exclude ID 3324 entirely or only include the CO(3−2) flux from the peak that is more spatially coincident with the HST imaging centroid. Spectra of non-detections are extracted using a point source profile at the target centroid location from HST /WFC3 F160W imaging (rest-optical). Employing boxcar extractions yields consistent results with slightly larger uncertainties. The extracted 1D CO(3−2) spectra are displayed in the bottom panels of each subplot in Figure 1.
Integrated CO(3−2) line fluxes are measured by fitting single Gaussian profiles to the extracted 1D spectra, with the uncertainty on the flux estimated from the covariance matrix. For non-detections, we fit a single Gaussian where the velocity width is fixed to be the same as that of the strong rest-optical lines (Hα and [O iii]λ5007) and the centroid is allowed to vary ±50 km s −1 from the systemic redshift. The resulting uncertainties are used to infer 3σ upper limits on the CO(3−2) flux. Employing bandpass integrations or double Gaussian profiles for sources that may be spectrally double-peaked (IDs 19753 and 2672) yield line fluxes that are consistent within the uncertainties with the single Gaussian values. The 6 robust detections have integrated CO(3−2) S/N=3.8−5.4, while the two tentatively detected lines have significances of 2.8σ and 2.3σ.
One source, ID 19753, is robustly detected in CO(3−2) but is not included in the analysis because of ambiguity in whether the rest-optical line emission is associated with the region dominating the CO(3−2) emission. ID 19753 displays a clumpy morphology with  two dominant components displaying a large color difference, with a red clump to the east and a blue western clump that dominates the UV. It is probable that the rest-optical line emission originates from the UV-bright blue component since the Hα/Hβ ratio implies relatively low nebular reddening (E(B-V) gas =0.16). Dust continuum emission from ALMA Band 6 observations are offset from the blue component and are more spatially coincident with the red component (Shivaei et al. 2022). Due to the large beam size of the Band 3 observations, we cannot determine whether the CO(3−2) emission is co-spatial with the blue or red component, though it is likely that the CO emission aligns with the dust emission. Given the likelihood of the CO(3−2) and restoptical lines of ID 19753 arising from distinct spatial components, we remove this source and proceed with a final sample of 13 galaxies.

Redshifts, Line Ratios, and SFRs
Rest-optical emission line fluxes are measured by fitting Gaussian line profiles to the slit-loss corrected 1D science spectra from the MOSDEF survey  Gaussian where the width of the three lines are tied but the [N ii] centroids can vary slightly from their expected position relative to Hα. In all cases, the continuum is set to the best-fit model from SED fitting such that stellar Balmer absorption is accounted for in the best-fit Hα and Hβ fluxes. Systemic redshifts are inferred from the best-fit centroids of the highest-S/N line (Hα or [O iii]λ5007).
Emission-line fluxes are corrected for the effects of dust by using the Hα/Hβ ratio to calculate E(B-V) gas , assuming an intrinsic ratio of 2.86 (Osterbrock & Ferland 2006) and the Milky Way extinction curve of Cardelli et al. (1989). Reddy et al. (2020) derived a nebular attenuation curve directly from z ∼ 2 MOS-DEF measurements that is consistent with the Milky Way curve. All emission-line ratios utilized for metallicity estimates are calculated using dust-corrected line fluxes. Star-formation rates (SFRs) are derived from dust-corrected Hα luminosities using the Hα conversion factor of Hao et al. (2011) adjusted to a Chabrier (2003) IMF. The star-formation rate surface density is calculated as Σ SFR = SFR/2πR 2 eff , where R eff is the restoptical half-light elliptical semimajor axis derived from HST /WFC3 F160W imaging as cataloged by van der Wel et al. (2014). We define the offset from the starforming main sequence at fixed M * as a function of redshift as ∆log(SFR) MS ≡log(SFR/SFR MS (M * , z)), where we adopt the SFR MS (M * , z) parameterization of Speagle et al. (2014), which matches the mean SFR-M * relation of MOSDEF star-forming galaxies at z ∼ 2.3 and z ∼ 3.3 ).

Stellar Masses
Stellar masses are derived from spectral energy distribution (SED) fitting of the extensive broadband photometry in the COSMOS field from the 3D-HST survey catalogs (Skelton et al. 2014;Momcheva et al. 2016). Photometry was fit using the SED fitting code FAST ) in combination with the flexible stellar population synthesis models of Conroy et al. (2009), assuming constant star-formation histories, solar stellar metallicity, a Chabrier (2003) IMF, and a Calzetti et al. (2000) attenuation curve. Because these models do not include a nebular emission component, photometric measurements were corrected for the contribution of strong rest-optical emission lines prior to fitting as described in Sanders et al. (2021), a necessary step due to the large emission-line equivalent widths common at z > 2 ( Reddy et al. 2018). Correcting the broadband photometry results in lower stellar masses by 0.12 dex on average in our sample, The two galaxies with the highest rest-optical equivalent widths have mass estimates reduced by 0.35 dex. The resulting best-fit stellar continuum model is used in the emission-line fitting, as described above.

Gas-Phase Metallicities
Gas-phase metallicities, given as 12+log(O/H), are derived from reddening-corrected rest-frame optical line ratios. We employ metallicity calibrations derived from composite spectra of extreme local galaxies that are analogs of z ∼ 2 galaxies from Bian et al. (2018). The functional form of the calibrations we use are given in Appendix A and Table A1. These calibrations reflect the significant evolution of H ii region ionization conditions between z ∼ 0 and z > 1 (e.g., Steidel et al. 2014Steidel et al. , 2016Shapley et al. 2015Shapley et al. , 2019Sanders et al. 2016Sanders et al. , 2020aStrom et al. 2017Strom et al. , 2018Runco et al. 2020;Topping et al. 2020a,b) and provide a good match to the existing sample of z > 1 galaxies with direct-method metallicities (Sanders et al. 2020b). Since we are interested in measuring O/H, we consider metallicities based on line ratios involving α elements (i.e., O, Ne) to be more robust than those involving N since N/O has a large scatter at fixed O/H in H ii regions and galaxies (∼0.2 dex; e.g., Pilyugin & Thuan 2011;Strom et al. 2017) due to the differing nucleosynthetic production channels of N and α elements. Accordingly, we establish a hierarchy of preferred metallicity indicators based on the emission lines available for a given source.
The  [NII] such that O3N2 and N2 metallicities are available. In the following analysis, the oxygen abundance of each source is taken to be the best available estimate based on the preference of O2O3Ne3, then O3N2, then N2 unless otherwise specified. We obtain consistent results if we instead use O3N2 metallicities uniformly for the full sample, with the main quantitative difference being that the slopes of the anti-correlations presented in Sec (1) where log(M 0 (z)/M ) = 9.90 + 2.06 × log(1 + z).

Composite CO(3−2) Spectra
Given that CO(3−2) emission is not robustly detected for approximately half of our sample, we employ spectral stacking techniques to include information from all of the targets in our analysis. We create composite spectra by converting the extracted 1D spectra from flux density to luminosity density and shifting the frequency axis to velocity offset based on the systemic redshift measured from the rest-optical lines. We then combine the individual spectra by taking the unweighted mean luminosity density at intervals of 50 km s −1 velocity offset. We do not use any weighting because the on-source integration times and sensitivities varied widely from source-tosource, such that the few sources with long integration times dominate the stacks if inverse-variance weighting is used. The composite CO(3−2) spectrum is then converted back to flux density using the mean redshift of the included individual galaxies. The composite CO(3−2) line flux is measured by fitting a Gaussian profile to the composite 1D spectrum, and is presented in Table 1.
We create composite CO(3−2) spectra of all of the galaxies in our sample (13 sources; labeled "stack-all") and the subset that are not robustly detected (8 sources; "stack-nondet"). We also stack in two bins of various galaxy properties divided such that there are roughly equal numbers of galaxies per bin. The galaxy properties according to which we bin include offset from the MS (∆log(SFR) MS ; "stack-delsfms"), and metallicity (12+log(O/H); "stack-oh"). Splitting the sample in Σ SFR results in identical bins as stack-delsfms, while binning according to offset from the mass-metallicity relation (∆log(O/H) MZR ) yields the same results as stackoh. The mass, SFR, effective radius, and metallicity associated with each composite CO(3−2) spectrum and given in Table 1 are the mean log(M * ), log(SFR), restoptical R eff , and 12+log(O/H) of the individual galaxies included in each stack.
Example composite CO(3−2) 1D spectra are shown in the bottom row of Figure 2 for the stack-all and stacknondet samples. The integrated CO(3−2) significance is ≥ 4.8σ for all stacked spectra used in this paper (stacknondet has the lowest significance). The top row displays 2D stacks constructed by producing datacubes using tclean with a uv taper of 3 and a spatial sampling of 0.5 /pixel to homogenize the beam sizes, integrating the channels within ±300 km s −1 of the systemic redshift, and taking the mean flux value at each spatial position relative to the HST imaging centroid. While we utilize composite CO(3−2) fluxes from the stacked 1D spectra in the following analysis, line fluxes derived by spatially integrating the 2D composites agree with the 1D fluxes within 20%.  . Example 2D (top) and 1D (bottom) composite CO(3−2) spectra for the stack-all (left; 13 galaxies) and stack-nondet (right; 8 galaxies) samples. The white crosshair is aligned with the expected spatial location of emission based on HST imaging centroids. The integrated CO(3−2) line S/N is given in the upper right corner of the bottom panels. The similarity of the noise patterns in the left and right columns arises because there are 8 galaxies in common between the two stacks.

Literature CO Sample
We draw from the literature a supplementary sample of star-forming galaxies at z = 2 − 3 with both CO and metallicity constraints available. We selected galaxies at z = 2 − 3 with CO(3−2) measurements from the PHIBSS ) and ASPECS  surveys. This literature CO sample was then cross-matched with rest-optical spectroscopic surveys to identify objects for which metallicity constraints could be derived. Objects with [N ii]/Hα > 0.5 were excluded as probable AGN. One additional object at z = 1.99 with rest-optical spectroscopy from Shapley et al. (2020) was included with a CO(1−0) measurement from Riechers et al. (2020). This literature sample comprises 23 galaxies with CO measurements (19 CO detections and 4 upper limits), 16 (12 CO detections and 4 upper limits) of which have gas-phase metallicity constraints from rest-optical emission lines (1 from O3O2Ne3, 6 from O3N2, and 9 from N2). The properties of the literature sample are presented in Table C2. Stellar masses and SFRs are taken from the literature sources and converted to a Chabrier (2003) IMF when necessary. Line ratios are calculated using the rest-optical line fluxes  Note-ID numbers are from the 3D-HST v4.1 photometric catalog ( from the given spectroscopic O/H references, which we then use to derive metallicities as described in Sec. 2.5. This sample of z = 2 − 3 literature sources complements the MOSDEF-ALMA sample in our analysis.

Sample Properties
The redshift, M * , and SFR properties of the MOSDEF-ALMA and z = 2 − 3 literature CO samples are shown in Figure 3. The two samples are well matched in redshift (top panel) with mean redshifts of 2.25 and 2.30, respectively. The combined sample spans z = 2.0 − 2.7 with a mean redshift of 2.28±0.17, where the uncertainty represents the sample standard deviation.
The middle panel of Fig Collectively, the combined sample is representative in both SFR and metallicity of typical z ∼ 2.3 star-forming galaxies in the mass range log(M * /M ) = 10.3 − 11.0, falling on the mean MS and MZR at this redshift. The addition of the MOSDEF-ALMA observations expands the sample of z ∼ 2 galaxies with CO measurements to lower masses and metallicities than were available to date. The wide range of offsets both above and below the mean MS and MZ relations makes MOSDEF-ALMA  an ideal sample with which to search for correlated residuals around mass-scaling relations that are signposts of self-regulated baryon cycling.

CO(3−2) Luminosity and Galaxy Properties
Before deriving molecular gas masses, we first explore empirical relations between observed CO(3−2) luminosity and global galaxy properties with the goal of assessing whether the CO production efficiency, and thus the CO-to-H 2 conversion factor, varies systematically. Integrated CO(3−2) line fluxes, S CO , are converted to total line luminosities using the relation of Solomon et al. (1997): where ν obs is the observed line frequency in GHz and D L is the luminosity distance in Mpc. The MOSDEF-ALMA CO(3−2) luminosities are given in Table 1.
The relation between L CO(3−2) and stellar mass is shown in the top panel of Figure 4. The points are colorcoded by offset from the star-forming main sequence, ∆log(SFR) MS . We find a significant positive correlation between L CO(3−2) and M * , with an evident secondary dependence on ∆log(SFR) MS such that galaxies at fixed M * with higher L CO(3−2) have higher SFR. Such a secondary dependence is expected if galaxies above (below) the MS have larger (smaller) molecular gas reservoirs than those on the MS, as has been observed in the local universe (e.g., Saintonge et al. 2016Saintonge et al. , 2017Saintonge & Catinella 2022) and at high redshifts with dust-and CO-based measurements (e.g., Tacconi et al. 2013Tacconi et al. , 2018Tacconi et al. , 2020Genzel et al. 2015;Scoville et al. 2017;Aravena et al. 2020). Using an orthogonal distance regression, we fit a function that is linear in log(M * ) and ∆log(SFR) MS to the combined sample of individual galaxies (excluding limits), and find a best-fit relation of: (3) Since this relation is a function of ∆log(SFR) MS , the fact that the MOSDEF-ALMA sample lies above the MS on average does not introduce bias. The stack of the MOSDEF-ALMA non-detections is 1σ consistent with this best-fit relation, suggesting that the exclusion of limits when fitting has not significantly biased this result. After accounting for measurement uncertainties, the remaining intrinsic scatter in L CO(3−2) around the best-fit function is σ int = 0.17 dex. This equation can thus be used to predict the CO(3−2) luminosity of z ∼ 2  galaxies to within a factor of ≈1.5 using M * and SFR alone. If the SFR term is neglected such that L CO(3−2) is a function of M * alone (gray line, top panel of Fig. 4), the intrinsic scatter is 0.31 dex.
The ratio of L CO(3−2) to SFR is a proxy of the CO production per unit M mol since SFR is tightly coupled to M mol via the molecular KS law (e.g., Kennicutt 1998 Fig. 4 implies that CO(3−2) emission is more efficiently produced per unit molecular gas mass in highmass galaxies. Figure 5 presents L CO(3−2) /SFR as a function of metallicity. There is a positive correlation between these two quantities, with significantly smaller scatter in L CO(3−2) /SFR at fixed O/H compared to that at fixed M * , suggesting that the relation between L CO(3−2) /SFR and O/H is more fundamental such that the correlation with M * emerges because of the mass-metallicity relation. Applying a Spearman correlation test to the individually detected sources yields ρ S = 0.45 and a p-value = 0.05, indicating a 2σ correlation. If the one significant outlier is excluded (discussed below), the correlation is stronger with ρ S = 0.65 and p-value = 0.003, indicating 3σ significance. Fitting a linear function to the individual galaxies, excluding limits, yields: There is no significant secondary dependence on ∆log(SFR) MS , suggesting that the best-fit relation is robust even though the MOSDEF-ALMA sample is offset above the MS on average. The composite spectra agree with the best-fit relation at the 1σ level, implying that our exclusion of limits when fitting has not significantly biased the result. The intrinsic scatter about this best-fit line is 0.12 dex after accounting for measurement uncertainties.
The single significant outlier in Fig. 5, lying ∼ 4σ below the best-fit line, is ID 13701 in the MOSDEF-ALMA sample, a potential merger with tidal tail-like features and two distinct brightness peaks in HST imaging (Fig. 1). It is possible that the low L CO(3−2) /SFR of ID 13701 is a result of an enhanced star-formation efficiency (SFR/M mol ) during a major merger (Di Matteo et al. 2007;Sargent et al. 2014; Reyes 2021) since L CO(3−2) traces M mol , though it is unclear how significant such enhancement is expected to be in high-redshift mergers (Fensch et al. 2017).
The positive correlation between L CO(3−2) /SFR and O/H demonstrates that the production efficiency of CO(3−2) is a function of metallicity, with lowmetallicity galaxies outputting less CO(3−2) per unit SFR than metal-rich galaxies. If the star-formation efficiency (or, equivalently, depletion timescale) is similar across our sample, then the correlation in Fig. 5 implies a decreasing CO luminosity per unit M mol with decreasing metallicity in z ∼ 2 star-forming galaxies. This same trend has been observed in local galaxy samples, manifesting as a metallicity-dependent CO-to-H 2 conversion factor (α CO ) with α CO increasing with decreasing metallicity (Wilson 1995;Arimoto et al. 1996;Wolfire et al. 2010;Schruba et al. 2012;Bolatto et al. 2013;Accurso et al. 2017). Higher α CO at lower O/H, corresponding to lower CO luminosity per unit gas mass, arises because CO is dissociated at greater depths into H 2 clouds in low-metallicity environments. CO molecules are photodissociated by far-UV photons and require dust-shielding to prevent dissociation (e.g., Wolfire et al. 2010;Glover & Mac Low 2011). Metal-poor galaxies have lower dust-to-gas ratios (e.g., Sandstrom et al. 2013;De Vis et al. 2019) and low-metallicity massive stars produce harder and more intense UV radiation fields, such that there is an increasing fraction of COdark H 2 gas with decreasing metallicity. The trend in Fig. 5 thus suggests that a metallicity-dependent α CO is required to accurately translate observed CO luminosities of z ∼ 2 galaxies into molecular gas masses.
3.3. The CO-to-H 2 Conversion Factor at z ∼ 2 Given the evidence for a metallicity-dependent α CO at z ∼ 2 presented above, we utilize a combined sample of MOSDEF-ALMA and literature targets at z > 1 with CO and rest-optical spectra to constrain the form of the α CO −O/H relation at high redshifts. Genzel et al. (2012) presented evidence for a α CO -O/H relation at high redshift using a sample of ∼ 40 z > 1 main-sequence galaxies, but the majority of their sample lacked spectroscopic metallicity constraints and instead had metallicities inferred from the MZR, with only a subset of 9 CO-detected galaxies having metallicities based on [N ii]/Hα ratios. The analysis below represents a significant improvement on the foundational results of Genzel et al. (2012) by using a sample composed entirely of galaxies with spectroscopic metallicities.
To maximize the sample size for the α CO analysis, we supplement the combined z = 2 − 3 MOSDEF-ALMA and literature CO+O/H sample with 14 additional galaxies at z > 1 with the necessary observations that were not included in our primary sample because they either fall outside of the target redshift range (i.e., at z = 1 − 2 or z > 3), are extreme starbursts with ∆log(SFR) MS ≥10, or are strongly gravitationally lensed. Sources in this supplementary z > 1 sample are given in Table C2. Thus, this expanded sample is only used for the α CO analysis, while the remainder of the paper employs only the more homogeneous z = 2 − 3 MOSDEF-ALMA and literature samples described in Sec. 2.7. The α CO sample includes 43 galaxies spanning z = 1.08 − 3.22 with a mean redshift of 2.12. Of these 43 sources, all have metallicity constraints from rest-optical line ratios, 33 have CO detections, and 10 have CO upper limits.
The CO-to-H 2 conversion factor α CO is defined as where M mol is the total molecular gas mass including a 36% contribution from Helium, and r J1 is the excitation correction factor to convert higher-J CO luminosities to that of the ground transition. We adopt r 31 = 0.55 ) and r 21 = 0.76 . Our α CO sample contains 33 galaxies with CO(3−2) measurements, 7 with CO(2−1) and 3 with CO(1−0). Our results are thus most sensitive to the assumed r 31 value, whereas the adopted r 21 will have only a minor effect. Constraining α CO requires an estimate of M mol that is independent of CO emission. We employ two methods of estimating M mol : from dynamical masses and the molecular KS law. Both techniques require the galaxy size in the calculations, so we required a half-light radius measurement in our selection criteria. We operate under the assumption that the effective radii of rest-optical continuum emission, Hα emission, and molecular gas are similar for z ∼ 2 galaxies, as has been found at z ∼ 1 − 2 by comparing resolved Hα and CO sizes to HST imaging Förster Schreiber et al. 2019). We describe the calculations of M mol using each method below.

Dynamical mass method
The dynamical mass M dyn provides a measure of the total mass in a galaxy, including stars, gas, and dark matter. We assume that z ∼ 2 galaxies are dominated by baryons within the effective radius such that dark matter is negligible (Genzel et al. 2017Price et al. 2021), and that the total gas mass in z > 1 galaxy disks is dominated by the molecular component . Under these assumptions, M mol can be inferred from M dyn using In the α CO sample, 13 objects have M dyn constraints based on measured rotation, either from Hα velocity mapping (7 galaxies; Förster , CO velocity mapping (1 galaxy; Swinbank et al. 2011), or forward-modeling of tilted emission lines in slit spectra (5 galaxies; Price et al. 2016Price et al. , 2020. For targets without measured rotation, M dyn can be estimated via the virial mass equation where G is the gravitational constant, σ v is the velocity dispersion measured from emission line widths (corrected for instrumental resolution), k is the virial coefficient, and R eff is the half-light elliptical semimajor axis. The value of the virial coefficient depends on the mass and velocity distribution of stars and gas in the source. Rather than assuming a value for k from the literature, we calibrate k such that dispersion-based M dyn estimated using equation 7 matches rotation-based M dyn on average for the 13 targets with measured rotation. This process yields a best-fit value of k = 8.6, very close to the virial coefficient for an exponential profile (Sérsic index n = 1) from Cappellari et al. (2006). Our bestfit virial coefficient is somewhat higher than has been used for disks and spheroids in past work (k ≈ 3 − 6; e.g., Pettini et al. 2001;Cappellari et al. 2006;Erb et al. 2006b;Price et al. 2016), though this difference is at least partially explained by the fact that we do not perform inclination corrections such that our best-fit k encompasses both the virial coefficient and an average inclincation correction factor. Excluding the 13 objects with measured rotation, we find that 26 of the remaining galaxies have published line widths, for which we calculate M dyn using equation 7. We thus have M dyn constraints for a sample of 39 galaxies (13 from rotation, 26 from line widths), for which we estimate M mol with equation 6 and α dyn CO using equation 5. Seven objects have M * >M dyn such that the inferred M mol and α dyn CO are unphysically negative, which we attribute to the large systematic scatter in deriving M dyn from spatially unresolved measurements. For these sources, we instead calculate 3σ upper limits on M mol and α dyn CO .
The molecular gas mass can also be estimated by applying the molecular KS relation between the surface densities of SFR (Σ SFR ) and molecular gas mass (Σ mol ) (e.g., Kennicutt 1998 converted Σ SFR to Σ mol based on the adopted KS law with Σ SFR in units of M yr −1 kpc −2 and Σ mol in units of M pc −2 , and finally derived M mol using where the same value of r eff is used in the first and last steps. We adopt the molecular KS law of Tacconi et al. (2013) adjusted for the different r 31 assumed in that work, with log(c) = −2.96 and n = 1.05. This M mol estimate is then applied in equation 5 to infer α KS CO .
3.3.2. The αCO−O/H relation at z ∼ 2 Figure 6 shows α CO as a function of metallicity using the dynamical mass method (left panel) and the KS law method (right panel). The black diamonds show mean values of the individual objects (excluding limits) in 3 bins of O/H. With both methods, we find a trend of decreasing α CO with increasing O/H. In each panel, only one source at 12+log(O/H) > 8.7 has inferred α CO > 10, while several galaxies exceed this value at lower metallicities. Likewise, α CO lower than the Milky Way value of 4.36 is common at 12+log(O/H) > 8.6, while at lower metallicities there are only two galaxies that have α dyn CO <4.36 and none with α KS CO below the Galactic value. A Spearman correlation test indicates ρ S = −0.55 and p-value = 0.0098 for α dyn CO , and ρ S = −0.62 and p-value = 1.1 × 10 −4 for α KS CO , rejecting the null hypothesis of no correlation at 2.6σ and 3.3σ, respectively.
The best-fit relations to the individual galaxies (excluding limits; thick gray lines in Fig. 6) where x = 12 + log(O/H) − 8.7 and α CO is in units of M (K km s −1 pc 2 ) −1 . The MOSDEF-ALMA stacks are consistent with the best-fit relations, suggesting that the exclusion of limits has not introduced a significant bias. These best-fit relations are fully consistent with one another despite the different methods employed to estimate α CO . They are also generally consistent with existing α CO −O/H calibrations (Genzel et al. 2012(Genzel et al. , 2015Bolatto et al. 2013;Accurso et al. 2017), though the z ∼ 2 data favor a steeper slope than that of Bolatto et al. (2013) and are in particularly good agreement with the Accurso et al. (2017)  ). If we instead use the z ∼ 0 calibrations of Pettini & Pagel (2004) or Curti et al. (2017Curti et al. ( , 2020, the metallicities decrease by ≈ 0.15 dex on average. Given that the majority of the sample has CO(3−2) measurements, the value of the r 31 excitation correction factor, which is not well constrained at high redshift, can shift galaxies vertically. If we instead assume r 31 that is 1.5 times larger (i.e., r 31 = 0.84 as found by Riechers et al. 2020 for five ∼ 10 11 M galaxies at z ≈ 2.6), the inferred α CO would decrease by 0.18 dex. Finally, the assumed molecular KS law and virial coefficient affect α KS CO and α dyn CO , respectively. While these systematic uncertainties are a significant concern for the absolute normalization of the α CO −O/H relation, they each shift the majority of galaxies in our sample in the same direction unless there is a strong systematic dependence of r 31 or the virial coefficient on metallicity. It is plausible that r 31 increases with decreasing metallicity since the combination of reduced dust shielding and harder and more intense UV radiation fields at low metallicity is expected to more highly excite CO. In this scenario, the true α CO (O/H) slope would be steeper than our derivation that assumes a constant r 31 , though more high-redshift galaxies with detections of both CO(3−2) and CO(1−0) are needed to robustly assess the impact of excitation Riechers et al. 2020). With the currently available constraints, we conclude that α CO depends on O/H at z ∼ 2 with similar dependence to that observed at z ∼ 0. It is crucial to take this metallicity dependence into account to derive accurate gas masses for subsolar metallicity galaxies that are common in the high-redshift universe.

The Molecular Gas Content of z ∼ 2
Star-Forming Galaxies In this section, we derive molecular gas masses for the MOSDEF-ALMA and z = 2 − 3 literature samples, and explore relationships between molecular gas content, stellar mass, SFR, and metallicity. Based on the good agreement with the z ∼ 2 galaxies in Fig. 6, we adopt the Accurso et al. (2017) relation between α CO and metallicity 2 with a floor at the Milky Way value of α MW CO = 4.36 M (K km s −1 pc 2 ) −1 such that super-solar metallicity galaxies do not have significantly smaller values of α CO : where x = 12 + log(O/H). Molecular gas masses are calculated assuming this α CO (O/H) relation and r 31 = 0.55, where the latter is adopted for consistency with Tacconi et al. (2018). For the 7 literature sources with-out spectroscopic metallicity constraints, O/H from the MZR(M * , z) of equation 1 is used for M mol calculations, and these sources are excluded from any plots and analyses that directly include metallicity. The derived molecular gas properties of the MOSDEF-ALMA sample and stacks are presented in Table 1, while those of the z = 2 − 3 literature sample are given in Table C2.

Molecular gas content and stellar mass
The molecular gas masses, M mol , and molecular gas fractions, µ mol =M mol /M * , are shown as a function of M * in the top and middle panels of Figure 7, with points color-coded by offset from the star-forming main sequence. On average, the MOSDEF-ALMA sample has lower M mol and higher µ mol than the more-massive literature sample, indicating a positive correlation between M mol and M * and an anti-correlation between µ mol and M * . M mol is typically larger than M * below 10 10.8 M , with an average µ mol ≈2 at 10 10.4 M . Consequently, gas mass is expected to dominate the baryonic mass of z ∼ 2 galaxies at 10 10.4 M .
There is a significant additional dependence on SFR such that galaxies with higher SFR at fixed M * have larger M mol , as has been observed in previous works (e.g., Tacconi et al. 2013Tacconi et al. , 2018Tacconi et al. , 2020Scoville et al. 2017;Liu et al. 2019). The black line shows the Tacconi The bottom panel of Figure 7 presents the molecular gas depletion timescale, t depl =M mol /SFR. We do not find evidence for any significant dependence of t depl on M * across log(M * /M ) = 10.2 − 11.2. In the sample averages, we find slightly lower t depl at lower M * , though the MOSDEF-ALMA and literature mean values are consistent with one another at the ≈1.5σ level. A dependence on offset from the MS is apparent, with higher-SFR galaxies at fixed M * having shorter t depl . The gray points once again show the MOSDEF-ALMA stacks corrected for their offset from the MS according to the dependence of the Tacconi et al. (2018) scaling relation (t depl ∝ ∆SFR −0.44 MS ). The SFR-corrected stack of all MOSDEF-ALMA targets has log(t depl /Gyr) = −0.20 ± 0.10 at 10 10.4 M , while the z = 2 − 3 literature mean is log(t depl /Gyr) = −0.09 ± 0.07 at 10 10.85 M . We thus find that the molecular gas depletion timescale of z ∼ 2 main-sequence star-forming galaxies at log(M * /M ) ∼ 10 − 11.5 is t depl = 600 − 800 Myr with no significant M * dependence, shorter by a factor of ∼ 2 than that of main-sequence z ∼ 0 galaxies at similar M * with t depl ≈ 1.4 Gyr .
3.4.2. The z ∼ 2 molecular KS law Figure 8 shows Σ SFR vs. Σ mol for z = 2 − 3 galaxies. 4 We find a tight correlation between these quantities that is well fit by a linear relation with constant t depl = 700 Myr (orange dashed line), in agreement with the findings of Tacconi et al. (2013) based on the PHIBSS z = 1.0 − 1.5 and z = 2.0 − 2.5 samples, the latter of which is included in our literature sample. The best-fit relation from fitting the individual galaxies (ex- Our best-fit relation is fully consistent with the z = 1 − 2.5 relation from PHIBSS ). The best-fit slope at z = 2−3 is consistent with the molecular KS law in the local universe that is generally found to be near-linear with either integrated or spatially-resolved observations (e.g., Kennicutt 1998;Bigiel et al. 2008;Leroy et al. 2008; de los Reyes & Kennicutt 2019; Kennicutt & De Los Reyes 2021), although the normalization is offset from what is appropriate for typical integrated z ∼ 0 main-sequence galaxies with t depl ≈ 1 − 1.5 Gyr (gray dotted line).
We find that the z ∼ 2 molecular KS law is very tight over nearly 3 orders of magnitude in Σ mol and Σ SFR , with an intrisic scatter of 0.07 dex (18%) in Σ SFR around the best-fit line after accounting for measurement uncertainties. The MOSDEF-ALMA results demonstrate that z ∼ 2 galaxies down to 10 10.2 M fall on a well-defined KS relation, implying that equation 14 can be used to accurately infer M mol (and M dyn , assuming negligible dark matter) for existing large samples of z ∼ 2 galaxies with SFR and size measurements in deep extragalactic legacy fields. This implication is especially beneficial for existing large rest-optical spectroscopic surveys, which have mean sample stellar masses of ≈ 10 10 M (e.g., MOSDEF, Kriek et al. 2015;KBSS-MOSFIRE, Steidel et al. 2014).

Correlated residuals around mean scaling relations
We now search for correlations among the residuals around three mean scaling relations: the SFR-M * relation (i.e., star-forming main sequence; MS), the massmetallicity relation (MZR), and the M mol -M * relation. Models of galaxy growth governed by a self-regulating baryon cycle predict that the scatter around these scaling relations will be (anti-)correlated (e.g.  Figure 7. We find a positive correlation between these quantities, indicating that galaxies with higher SFR have larger M gas at fixed M * . The correlation is highly significant: a Spearman correlation test yields ρ S = 0.79 with a p-value=5.6 × 10 −5 . The trends of stacked MOSDEF-ALMA data in two bins of ∆log(SFR) MS (red squares) and the individual objects are consistent with the SFR dependence of the Tacconi  in which galaxies at fixed M * with higher SFR have lower O/H. The inverse correlation of ∆log(O/H) MZR on ∆log(SFR) MS has been observed with much higher significance using both individual galaxies and stacked spectra of a considerably larger (N ∼ 250) MOSDEF z ∼ 2.3 star-forming sample (Sanders et al. 2018. Here, we confirm the existence of a SFR-FMR among the 13 MOSDEF-ALMA targets. There is not a significant correlation present for the literature CO sample in the right panel of Figure 9. This apparent lack of a SFR-FMR may be due to the small sample size or the impact of systematic uncertainties on O/H, or could be the true physical scenario for this sample. The small sample size of only 16 galaxies can make it difficult to recover the secondary dependence of O/H on SFR in addition to the primary dependence on M * , considering that there is both intrinsic scatter in the SFR-FMR and several literature sources have fairly large statistical uncertainties on O/H. The metallicities of the literature sample are derived from a heterogeneous set of strong-line ratios: 9 based on N2, 6 from O3N2, and 1 from O3O2Ne3. While we have made efforts to control systematics between O/H derived from each of these indicators (see Appendix A), metallicites based on N2 (and O3N2 to a lesser extent) are strongly affected by object-to-object variations in N/O that are commonly a factor of 2 at fixed O/H (Pilyugin & Thuan 2011;Strom et al. 2017, Fig . A1). In contrast, the MOSDEF-ALMA sample is more uniform in metallicity indicator, with 11 based on O3O2Ne3 (unaffected by N/O) and 2 derived using O3N2 (less strongly affected by N/O than N2).
It is also possible that there truly is no dependence of O/H on SFR at fixed M * in the literature sample.
In the local universe, the SFR-dependence in the SFR-FMR significantly weakens and may disappear entirely or even invert at high masses (> 10 10.5 M ; e.g., Mannucci et al. 2010;Yates et al. 2012). That the massive literature sample (mean M * = 10 10.85 M ) displays no SFR dependence while the less-massive MOSDEF-ALMA sample (mean M * = 10 10.4 M ) does could thus indicate a behavior similar to that of the z ∼ 0 SFR-FMR at high masses. Uniform full-coverage of the restoptical lines for a larger sample of such massive z ∼ 2 star-forming galaxies is required to discern whether the lack of a correlation is the result of systematic effects or is real.
The top panel of Figure (Fig. 8), such that the SFR-FMR likely emerges because this more fundamental Gas-FMR exists among M * , M mol , and O/H. We thus expect that if a SFR-FMR is present, a Gas-FMR will be present as well, as found in the MOSDEF-ALMA sample. We do not find evidence for a Gas-FMR among the more-massive z = 2 − 3 literature sample, consistent with the lack of a SFR-FMR in the right panel of Fig. 9.
Interpreting the top panel of Figure 10  where O/H is inferred from the MZR(M * , z) given in equation 1, and finally multiplying by r 31 . When evaluated at z = 2.3, the resulting L CO(3−2) (M * , z) relation is similar to the best-fit L CO(3−2) vs. M * relation in the top panel of Figure 4, such that we obtain consistent results if that relation is instead assumed.
The bottom panel of Figure 10 displays ∆log(L CO(3−2) ) T18,MS vs. ∆log(O/H) MZR . The distribution of MOSDEF-ALMA targets and stacks is flat, such that L CO(3−2) has no dependence on O/H at fixed M * . Consequently, the anti-correlation between ∆ log(M mol ) T18,MS and ∆log(O/H) MZR is significantly affected by how α CO depends on metallicity. If α CO has no dependence on O/H, then no Gas-FMR is found in the z ∼ 2 MOSDEF-ALMA sample. However, we have presented strong evidence for an O/H-dependent α CO in Fig. 6. Furthermore, the inverse relation between α CO and O/H is well-established in the local universe (e.g., Bolatto et al. 2013;Accurso et al. 2017) and is thought to be driven by physics that appear to behave similarly at higher redshifts, namely a correlation between metallicity and dust content (Reddy et al. 2010;Shapley et al. 2020Shapley et al. , 2021Shivaei et al. 2022, Popping et al., in prep.), and an anti-correlation between metallicity and UV radiation field intensity. We thus conclude that the evidence for an anti-correlation between M mol and O/H at fixed M * in the MOSDEF-ALMA sample is robust, though the exact slope of this inverse relation will depend on the α CO (O/H) relation assumed. Our analysis of the relations among M * , SFR, metallicity, and gas content provides strong evidence for the existence of baryon cycling processes that govern the growth of z ∼ 2 massive galaxies (log(M * /M ) ∼ 10.0 − 11.5). The typical depletion times of 700 Myr are only ≈25% of the age of the universe at z = 2.3. Since a star-forming main sequence exists at this redshift, accretion of gas from the IGM and/or CGM is required to sustain star formation in these galaxies, as has been pointed out by previous studies (e.g., Tacconi et al. 2013). However, the gas fractions are not large enough to explain the ISM metallicities at z ∼ 2. For the typical µ mol ≈ 1.5 and assuming neutral atomic gas is negligible, a simple accreting box model with an accretion rate equal to the SFR predicts 12+log(O/H)=9.1, assuming the stellar O yield of core-collapse SNe from a Chabrier (2003) IMF. This value is considerably higher than the average observed 12+log(O/H) = 8.6 − 8.7.
Considering the large amount of evidence for ubiquitous high-velocity outflowing gas around high-redshift star-forming galaxies (e.g., Steidel et al. 2010;Förster Schreiber et al. 2019;Weldon et al. 2022), strong gas outflows present a probable candidate mechanism to further decrease the effective yield and drive model ISM metallicities lower toward observed values.

Implications for outflow mass loading
The mass outflow rate (Ṁ out ) is often parameterized in the outflow mass loading factor (η out =Ṁ out /SFR) that encapsulates the efficiency with which feedbackdriven galactic winds remove mass (and metals) from galaxies. We place quantitative constraints on η out by employing the ideal gas-regulator model of Lilly et al. (2013). In this formalism, the ISM metallicity Z ISM is a function of η out and the gas fraction µ gas =M gas /M * : where y is the stellar metal yield, R is the fraction of stellar mass returned to the ISM through stellar evolutionary processes, and Z 0 is the metallicity of accreting gas. We assume Z 0 Z ISM such that accreting metals are negligible. We adopt a stellar oxygen yield of y O = 0.033 as a mass fraction (equivalent to 12+log(O/H) y = 9.45 by number density) and R = 0.45, both values appropriate for core-collapse SNe enrichment with a Chabrier (2003) IMF (Vincenzo et al. 2016). Finally, we assume that the total gas mass is dominated by the molecular component in high-redshift star-forming galaxies such that M gas ≈M mol , as found in studies of the atomic-to-molecular hydrogen fraction using semi-analytic models (Lagos et al. 2011;Popping et al. 2014Popping et al. , 2015. The top panel of Figure 11 displays metallicity as a function of µ mol . The black lines show gas-regulator models of equation 15 with η out varying between 0.5 and 8. We find that the z ∼ 2 CO sample is generally described by models with η out ranging from 0.5 − 4, with an average value of η out ≈ 2. Only 5 galaxies lie in the unphysical region of parameter space that would mathematically require η out < 0, though 4 of these are ≈ 1σ consistent with the physical region. The data are incompatible with significantly lower stellar yields (i.e., y O = 0.015 for a Salpeter (1955) IMF), which would place a large fraction of the z ∼ 2 sample in the unphysical regime. Such low yields could only be accommodated if Z 0 is comparable to Z ISM . Suzuki et al. (2021) found that z ∼ 3.3 galaxies at ∼ 10 10.5 M have metallicities significantly lower than z ∼ 0 and z ∼ 1.6 galaxies at fixed µ gas (where µ gas includes both molecular and atomic components for local galaxies). These authors concluded that there is a sharp increase in η out at fixed µ gas between z ∼ 1.6 and z ∼ 3. We find that our z ∼ 2 sample, with a mean value of 12+log(O/H) = 8.6 − 8.7 at µ gas =1.6, falls near the distribution of the local and z ∼ 1.6 galaxies at fixed µ gas , implying no strong evolution in η out at fixed µ gas over z = 0 − 2.5. The apparent evolution in η out at fixed µ gas at z > 3 found by Suzuki et al. (2021) may be due to systematics associated with converting dust mass to M mol and the metallicity calibrations these authors adopted. Obtaining joint CO and metallicity constraints for z > 3 main-sequence galaxies would aid in understanding whether η out is different at z > 3. Using our measured values of O/H and µ mol , we solve equation 15 for η out and present the resulting values as a function of M * in the bottom panel of Figure 11. The majority of the combined z ∼ 2 sample has η out = 1 − 4, with an average value of ≈2 and no apparent dependence on M * over log(M * /M ) ∼ 10.2 − 11.2. The blue line displays an inference of the average η out as a function of M * for lower-mass main-sequence galaxies, obtained by inputting the best-fit z ∼ 2.3 MZR of Sanders et al. (2021) and the Tacconi et al. (2018) µ mol scaling relation into equation 15, assuming the same y O and R as above. As found by Sanders et al. (2021), η out decreases with increasing M * up to 10 10.5 M as η out ∝ M −1/3 * , at which point the blue line aligns with the average value of the higher-mass CO sample. The CO-based results suggest that η out does not continue falling with increasing M * above 10 10.5 M , but instead flattens out. A similar behavior is inferred on average for z ∼ 0 galaxies (dashed black line). Since η out > 1 even at the highest masses, the flattening of η out with respect to M * is responsible for the flattening of the MZR at high masses (Fig. 3). As Sanders et al. (2021) noted, the lower-mass inferences based on the MZR suggest that, at fixed M * η out is larger at z ∼ 2 than at z ∼ 0. This does not appear to be the case at > 10 10.5 M , where both lowand high-redshift galaxies are inferred to have similar η out on average.
In the EAGLE simulations, Mitchell et al. (2020) found that η out falls roughly as a power law with M * even out to very high masses when AGN feedback is not included, while including AGN feedback causes η out to flatten out above halo masses of ∼ 10 12 M (M * ∼ 10 10 M ). Nelson et al. (2019) see a flattening of η out toward high masses and eventual steep increase of η out with M * at very high masses that are attributed to AGN feedback in IllustrisTNG. Other analyses of cosmological simulations have found that the flattening of the MZR at high masses requires the inclusion of AGN feedback (e.g., De Rossi et al. 2015Ma et al. 2016).
While we have excluded AGN from our sample, the tracers used to identify AGN (X-ray luminosity, near-IR colors tracing hot dust, and rest-optical line ratios) can only identify galaxies with contemporaneous high rates of black hole accretion. Galaxies with AGN activity in the recent past but low levels of black hole accre-tion in the present would correctly be identified as starforming, though AGN feedback may have a lasting effect on ISM gas via material removal or gas heating that makes it easier for SNe feedback to drive outflows. This scenario is consistent with the idea that super-massive black hole accretion rates trace SFRs (in turn tracing gas inflow rates) on average (Madau & Dickinson 2014), while maintaining significant variability such that galaxies alternate between an "on" AGN phase and "off" starformation dominated phase (Hickox et al. 2014;Heckman & Best 2014). We speculate that the transition from a power-law η out (M * ) to a flat η out occurs at the mass above which the AGN duty cycle becomes high enough that the imprint of AGN feedback on ISM gas does not have time to be erased between consecutive AGN "on" phases, providing a physical driver for observed MZR flattening at high masses. The occurrence rate of strong AGN rises steeply above the mass at which the MZR begins to flatten at z ∼ 0 (e.g., Kauffmann et al. 2004;Aird et al. 2019). Using the KMOS 3D sample at z = 0.6−2.7, Förster  found that the occurrence rates of both AGN and AGN-driven outflows grow steeply with increasing M * and exceed 10% at M * > 10 10.5 M , similar to the mass where we find the onset of MZR and η out flattening at z ∼ 2. While these results are generally consistent with intermittent AGN feedback as a driver of high-mass MZR flattening at z ∼ 2, robustly evaluting this scenario will require better quantitative constraints on the turnover mass of the z ∼ 2 MZR to directly compare to measures of AGN frequency as a function of M * . 4.1.2. The role of accretion rate in scaling relations at The (anti-)correlated residuals around the MS, MZR, and M mol -M * relation displayed in Figures 9 and 10 demonstrate the existence of a multi-dimensional relation among M * , SFR, M mol , and metallicity in the z ∼ 2 MOSDEF-ALMA sample. This relation is such that, at fixed M * , galaxies with higher gas masses have higher SFR and lower metallicity. Such trends are found in a range of cosmological simulations (e.g., Ma et al. 2016;Davé et al. 2017Davé et al. , 2019De Rossi et al. 2017;Torrey et al. 2018Torrey et al. , 2019, and are a staple of models of galaxy evolution based on a self-regulating baryon cycle (Davé et al. 2012;Lilly et al. 2013;Peng & Maiolino 2014). These trends have also been observed in z ∼ 0 star-forming galaxy populations as the SFR-FMR and Gas-FMR (Mannucci et al. 2010;Lara-López et al. 2010;Bothwell et al. 2013Bothwell et al. , 2016bHughes et al. 2013;Lara-Lopez et al. 2013;Brown et al. 2018). Evidence for a SFR-FMR at z > 1 has been found previously (Zahid et al. 2014;Sanders et al. 2018Sanders et al. , 2021Henry et al. 2021), while hints of a Gas-FMR have been found at z ∼ 1.4 (Seko et al. 2016). Here, we confirm the interdependent nature of these four properties and the simultaneous existence of the SFR-FMR and Gas-FMR at high redshift for the first time.
The correlated residuals in M mol , SFR, and metallicity can be explained if deviations of M mol from the mean reflect variations in the gas inflow rate,Ṁ in . The observed trends then indicate that higherṀ in drives larger SFR due to the availability of more cold gas while simultaneously diluting metals in the ISM as metal-poor gas is mixed in, lowering O/H. This scenario implies that the scatter of the MS, MZR, and M mol -M * relations are all primarily driven by variations in accretion rates, a conclusion reached by previous works modeling the SFR-FMR and Gas-FMR (Dayal et al. 2013;Ford et al. 2014;Brown et al. 2018;Torrey et al. 2019). The simultaneous existence of a correlation between ∆ log(M mol ) T18,MS and ∆log(SFR) MS and an anti-correlation between ∆ log(M mol ) T18,MS and ∆log(O/H) MZR requires that SFR and metallicity respond to the accretion of fresh gas quickly, faster than the timescale of significant variations in the inflow rate, as is seen in cosmological simulations (Muratov et al. 2015;Torrey et al. 2018Torrey et al. , 2019. Indeed, the dynamical timescales of z ∼ 2 galaxy disks (on which galaxyaveraged SFR and metallicity can vary) are a few tens of Myr (e.g., Förster Reddy et al. 2012), while halo dynamical times (on whichṀ in from smooth gas accretion is expected to vary) are roughly an order of magnitude larger (e.g., Torrey et al. 2018). The overall tightness of these three scaling elations (1σ scatter=0.1 − 0.3 dex at fixed M * ) then suggests that large variations inṀ in averaged over the SFR and O/H response timescale (i.e., disk dynamical time, ∼ 30 − 50 Myr) are rare.
Comparing the observed quantitative slopes of these correlated residual relations to those produced in numerical simulations and semi-analytic models of galaxy formation presents an avenue to make high-order tests of accretion and feedback models (e.g., Davé et al. 2017;Torrey et al. 2019;Pandya et al. 2020). The slopes are not well-constrained with the current z ∼ 2 MOSDEF-ALMA sample due to the small sample size and data quality. This work motivates a larger and deeper survey of CO in z ∼ 2 galaxies with accompanying metallicity measurements to quantitatively constrain the secondary dependence on gas mass.

ISM metal mass and retention
The total metal mass in the ISM can be estimated from the combination of metallicity and gas mass. We convert from the number density 12+log(O/H) to a mass fraction Z O = M O /M gas assuming He makes up 36% of the gas mass, and then use the ratio of the oxygen and total metal stellar yields y O /y Z = 0.6 5 to infer the total gas-phase ISM metal mass fraction Z ISM = M Z /M gas . We then multiply Z ISM by M mol , assuming M mol ≈M gas for the z ∼ 2 sources, to obtain the mass of metals in the ISM, M Z,ISM . The top panel of Figure 12 shows that M Z,ISM increases with increasing M * on average (i.e., for the composites), as expected since both metallicity and M mol increase with increasing M * . The typical ISM metal mass of z ∼ 2 star-forming galaxies is M Z,ISM ≈ 10 8.7 M at M * = 10 10.5 M . The ratio of gas-phase ISM metal mass to galaxy stellar mass has been used as a measure of a galaxy's ability to retain metals since M * is approximately proportional to the total metals produced over a galaxy lifetime (Ma et al. 2016;Torrey et al. 2019). Torrey et al. (2019) defined this quantity as the ISM metal retention efficiency, γ Z,ISM = M Z,ISM /M * , and argued that evolving gas fractions, not outflow efficiencies, control the evolution of the MZR because higher redshift galaxies in Il-lustrisTNG have higher γ Z,ISM at fixed M * and are thus inferred to be more efficient at retaining gas-phase metals than their low-redshift counterparts. We plot γ Z,ISM vs. M * in the bottom panel of Figure 12, and find that the combined z ∼ 2 sample indeed displays higher ISM metal retention efficiencies at fixed M * by a factor of ∼ 3 than the mean value at z ∼ 0 (black dashed line).
This result implies that z ∼ 2 galaxies at log(M * /M ) ∼ 10.5 − 11 are less efficient at ejecting metals, seemingly in conflict with the result based on gas-regulator models that η out is the same on average at z ∼ 2 and z ∼ 0 over this mass range (Fig. 11, bottom  panel). In their analysis of the MZR in the IllustrisTNG simulations, Torrey et al. (2019) pointed out the tension between γ Z,ISM that implies that gas fraction primarily governs the MZR shape and its evolution (see also Ma et al. 2016), and analyses based on gas-regulator or equilibrium models in which outflow efficiency is inferred to play the larger role (e.g., Peeples & Shankar 2011;Davé et al. 2012;Lilly et al. 2013;Sanders et al. 2021). Torrey et al. (2019) argued that the importance of outflows is overestimated when using current gas-regulator models because they are not properly treating potentially important physical processes known to occur in numerical simulations, notably the recycling of enriched previously-outflowing material (Finlator & Davé 2008;Muratov et al. 2015Muratov et al. , 2017Anglés-Alcázar et al. 2017;Pandya et al. 2020). Theoretical work has shown close connections between η out and µ gas or Σ gas (Hayward & Hopkins 2017;Kim et al. 2020), such that the action of outflows vs. gas fraction in setting the ISM metallicity may not be cleanly separable. We will revisit the apparent tension between Figures 11 and 12 at the end of Sec. 4.3.

A budget of metals in massive z ∼ 2 star-forming galaxies
We now combine the ISM metal mass calculated above in Sec. 4.2 with estimates of the metal mass in stars and dust to perform a simple metal budget analysis. We carry out this analysis on the sample averages of the MOSDEF-ALMA and z = 2 − 3 literature samples. We estimate the mass of metals locked in stars, M Z,stars , by multiplying the stellar metallicity by M * , assuming the stellar metallicity is equal to the current gas-phase metallicity. This approach likely overestimates the true M Z,stars since ISM metallicity increases with time on average such that only recently formed stars will have metallicity equal to that of the current ISM, while previous generations of stars should have lower metallicity. The mass of metals in dust, M Z,dust , is taken to be the total dust mass since dust grains are essentially entirely composed of metals. To estimate M Z,dust , we combine the gas mass with the dust-to-gas ratio as a function of O/H from De Vis et al. (2019) 6 , which is consistent with current dust-to-gas constraints at z ∼ 2 (Shapley et al. 2020, Popping et al., in prep.). As before, we assume that M gas ≈M mol in high-redshift star-forming galaxies (Lagos et al. 2011;Popping et al. 2014Popping et al. , 2015Tacconi et al. 2018). The lifetime metal production is calculated by first dividing M * by (1-R) to infer the lifetime total stellar mass produced, and multiplying the resulting value by the stellar metal yield. We adopt y Z = 0.051 and R = 0.45 for a Chabrier (2003) IMF (Vincenzo et al. 2016).
The left panel of Figure 13 displays the mass of metals in the ISM, stars, dust, and the sum total of metals accounted for inside galaxies compared to the lifetime metal mass produced for the z ∼ 2 samples. We find that the majority of metals inside z ∼ 2 galaxies reside in the gas-phase ISM, even at very high M * . This is unlike what is found at z ∼ 0, where only low-mass galaxies at M * < 10 10 M have most of their metals in the ISM, while the majority of metals in more massive local galaxies resides in stars (Peeples et al. 2014;Oppenheimer et al. 2016;Muratov et al. 2017). This difference can be attributed to the steep increase of µ gas with redshift: z ∼ 2 galaxies have roughly an order of magnitude larger M gas than z ∼ 0 galaxies at fixed M * . Dust appears to be a subdominant destination for metals at z ∼ 2, similarly to z ∼ 0. Our dust mass estimates, based on a dust-to-gas ratio, are ∼ 0.2 dex lower than those derived from ALMA dust-continuum measurements of z ∼ 2.3 galaxies at similar M * by Shivaei et al. (2022). Even if M Z,dust was larger by that amount, it would remain smaller than M Z,ISM and M Z,stars .
The fraction of metals accounted for relative to the lifetime metal mass produced is presented in the right panel of Figure 13 Metal mass has been estimated for the gas-phase ISM (blue circle), stars (yellow star), interstellar dust (red "x"), and the sum of these three components representing the total metal mass accounted for within the galaxy (orange square). The black line shows the estimated total metal mass produced. Right: Fraction of total produced metals accounted for in the ISM, stars, and dust, obtained by dividing the metal mass in each phase by the lifetime metal mass produced. The orange error bar shows the cumulative uncertainty on the sum. The approx70% of produced metals that are unaccounted for are inferred to be lost to the CGM and IGM.
the ISM, stars, and dust makes up only 30 ± 4% of the total metals produced, implying that the majority of produced metals have been ejected from the galaxies in our samples and reside in either the CGM or IGM. This value is very similar to the 25% average metal retention fraction derived for z ∼ 0 galaxies by Peeples et al. (2014) using similar metal budgeting techniques, and is also close to that of ∼ 10 10 M galaxies at z = 2 in the FIRE simulations (Muratov et al. 2017). The metal loss fraction of 70% does not show variation across the mass range spanned by the z ∼ 2 samples, again similar to the behavior at z ∼ 0 (Peeples et al. 2014;Oppenheimer et al. 2016) despite large differences in global properties such as µ gas , SFR, and Σ SFR that are thought to be related to outflow strength. This analysis predicts a large mass of metals in the CGM of z ∼ 2 star-forming galaxies, with M Z,CGM ∼ 10 9 M at M * = 10 10.5 M . Significant metal absorption from low-and high-ionization species has been observed in the CGM of z ∼ 2 galaxies (e.g., Adelberger et al. 2005;Steidel et al. 2010;Rudie et al. 2019). Recently, Rudie et al. (2019) found that the mass of metals in the CGM of M * ∼ 10 10 M z ∼ 2 galaxies is > 25% of M Z,ISM alongside evidence that metal ejection from the halo into the IGM may be common. Improved observational constraints on the CGM metal mass of high-redshift galaxies are needed to see if the missing metals in Figure 13 can be accounted for. The fact that the total metal retention fraction is found to be ≈30% for galaxies at log(M * /M ) ∼ 10.5 at both z ∼ 0 and z ∼ 2 implies that they have a similar efficiency of metal ejection, consistent with our inferences of η out at both redshifts in the bottom panel of Fig. 11. We can now understand the apparent tension in the evolution of the ISM metal retention efficiency, γ Z,ISM , at fixed M * in Fig. 12 toward higher metal retention with increasing redshift. As a consequence of their large gas fractions (i.e., M gas > M * ), z ∼ 2 galaxies store the largest fraction of their metals in the gas-phase ISM, with smaller contributions from stars and dust. In contrast, gas fractions are low at z ∼ 0 (M gas ≈ 0.2M * at ∼ 10 10.5 M ), such that the vast majority of metals are stored in stars and only a small part is in the ISM for local galaxies in this mass range (Peeples et al. 2014). Thus, the evolution of γ Z,ISM at fixed M * in Fig. 12 reflects an increasing fraction of total metals stored in the gas-phase ISM with increasing redshift as a consequence of higher µ gas , rather than a higher overall metal retention (and lower metal ejection efficiency) at higher redshift.
We therefore do not find any tension between analyses based on gas-regulator models and metal retention fractions, but stress that galaxy metal retention (and its evolution with redshift) cannot be properly evaluated using metals in the gas-phase ISM alone because the relative contribution of this phase to the total metal mass changes as a function of both M * and redshift. This finding suggests that the missing physics (e.g., outflow recycling) in current gas-regulator models may not be required to obtain a reasonably accurate understanding of the origin and evolution of the MZR, though this problem ultimately requires the development of analytic models including these high-order processes to evaluate the impact of their inclusion or exclusion. In Fig. 11, inferences at lower masses from the MZR suggest that η out increases with increasing redshift below 10 10 M . Accordingly, we expect the total metal retention fraction to be lower at z ∼ 2 than at z ∼ 0 at fixed M * for low-mass systems. This motivates an analysis of the metal budget at z ∼ 2 across a wider mass range, which we will address in future work.

Cold gas content scaling relations at z ∼ 2
Recent studies have calibrated scaling relations of M mol , µ mol , and t depl as a function of redshift, M * , and SFR (or, equivalently, offset from the MS) by combining samples with Rayleigh-Jeans dust continuum and/or CO line emission observations spanning z = 0 − 4 (Genzel et al. 2015;Scoville et al. 2017;Tacconi et al. 2018Tacconi et al. , 2020Liu et al. 2019;Wang et al. 2022). These scaling relations potentially have great utility by providing estimates of the cold gas content of galaxies based on galaxy properties that are easier to directly constrain. However, at z ≥ 2, the samples used to calibrate such scaling relations are composed of massive galaxies (M * 10 10.5 ) such that using these relations to estimate gas masses of lower-mass high-redshift galaxies requires extrapolation. Furthermore, the vast majority (and in some works entirety) of galaxies in the z > 1 calibration samples have gas mass measurements based on dust continuum emission, while only a small fraction have CObased M mol . Using the combination of the lower-mass MOSDEF-ALMA and more-massive literature CO samples in this work, we can evaluate which cold gas scaling relations are reliable at z ∼ 2.
The scaling relations of Scoville et al. (2017), Tacconi et al. (2018), andLiu et al. (2019) are displayed in each panel of Figure 7, all evaluated at z = 2.3 for galaxies on the MS. We find that the scaling relation of Tacconi et al. (2018) reliably reproduces CO-based M mol and µ mol of main-sequence galaxies on average (gray square and colored diamond) down to 10 10.4 M , in particular having a slope that matches the observed trend with M * . In contrast, the Scoville et al. (2017) and Liu et al. (2019) relations have slopes as a function of M * that are too shallow for M mol and too steep for µ mol , suggesting that these relations will overestimate gas mass when extrapolating to M * 10 10 M . The normalization of the Scoville et al. (2017) relation is too high across the full range of masses probed by the z ∼ 2 CO sample, while the Liu et al. (2019) relation underestimates M mol and µ mol at ∼ 10 11 M . The good match between the relation of Tacconi et al. (2018) and the CO-based mean values suggests that the Tacconi et al. (2018) parameterization is more reliable when extrapolating to lower masses.
We found little dependence of t depl on M * in Figure 7, with lower-mass galaxies having marginally smaller t depl . Both the Scoville et al. (2017) relation that has no mass dependence and the relation of Tacconi et al. (2018)  ) is clearly too steep to match the CO-based t depl and would severely overestimate t depl at lower M * . The Scoville et al. (2017) and Tacconi et al. (2018) relations bracket the observations in normalization, with the former slightly higher and the latter slightly lower.
The SFR dependence of M mol and µ mol at fixed M * is nearly identical for Tacconi et al. (2018) and Liu et al. (2019) (µ mol ∝ (SFR/SFR MS ) 0.53 ), matching the observed trend in Figure 9, left panel. Scoville et al. (2017) find µ mol ∝ (SFR/SFR MS ) 0.32 , significantly shallower than the CO-based constraints at z ∼ 2. In summary, we find that the Tacconi et al. (2018) scaling relations provide the best overall match to the CO-based observations for estimates of M mol , µ mol , and t depl as a function of M * and SFR at z ∼ 2, and provide a particularly good match in mass-dependence suggesting extrapolation to lower masses is reliable. 7 4.5. Prospects for detecting CO of low-mass high-redshift galaxies with ALMA With the MOSDEF-ALMA sample, we have pushed the stellar mass of z ∼ 2 main-sequence galaxies with CO-based M mol down by a factor of 2, from 10 10.7 M to 10 10.4 M . It is of interest to investigate the prospects of detecting CO emission in even lower-mass z ∼ 2 galaxies to expand the dynamic range of parameter space probed by CO in the M * -SFR plane. We estimated the integrated CO(3−2) line flux, S CO(3−2) , for z = 2.3 galaxies on the MS by combining equations 4 and 1 to estimate L CO(3−2) /SFR as a function of M * , using the z = 2.3 MS of Speagle et al. (2014) to convert to L CO(3−2) , and finally obtaining S CO(3−2) with equation 2. We then inferred S CO(3−2) for offsets above and below the MS using the best-fit SFR dependence of L CO(3−2) at fixed M * : L CO(3−2) ∝ (SFR/SFR MS ) 0.71 (equation 3).
The predicted S CO(3−2) at z = 2.3 across the M * -SFR plane is displayed in Figure 14. We converted S  to the peak line flux density assuming a velocity FWHM of 150 km s −1 , the average value for z ∼ 2.3 galaxies in the MOSDEF survey spanning log(M * /M ) = 9.0−10.5 (Price et al. 2016). We then used this peak flux density to estimate the required integration time with ALMA to detect CO(3−2) across the M * -SFR plane. We employed the ALMA Sensitivity Calculator 8 to calculate the RMS sensitivity for different on-source integration times. We assumed an observed frequency of 104.8 GHz (CO(3−2) at z = 2.3) that falls in ALMA Receiver Band 3, an equatorial target (e.g., the COSMOS field), and a bandwidth of 50 km s −1 equal to 1/3 of the line FWHM for the RMS calculations. Using Monte Carlo simulations, we found that integrated line S/N=4 is obtained if the RMS value integrated over 1/3 of the line FWHM is 3.5 times smaller than the line peak. For a given integration time, we identified the contour in the M * -SFR plane for which the line peak is 3.5 times the RMS value, above which galaxies will have integrated CO(3−2) S/N>4. The blue, green, orange, and red lines show contours where the integrated line flux has S/N=4 in on-source integration times of 1, 2.5, 5, and 10 hours, respectively. Integrations of ∼ 1 h on-source reach the current mass limit of 10 10.4 M on the MS. At z ∼ 2.3, a survey with 5 h on-source will reach 10 10 M on the MS, 10 9.7 M for moderate starbursts (3×MS), and 10 9.4 M for extreme starbursts (10×MS). Such a survey could nearly double the parameter space containing galaxies with CO detections at z ∼ 2 relative to the current sample, probing a regime that significantly overlaps with samples from large spectroscopic surveys for which detailed metallicity, SFR, dust reddening, and ionization constraints are available. Detecting CO(3−2) in main-sequence galaxies at masses significantly below 10 10 M would require prohibitively long integrations with ALMA and is thus not a feasible route to constrain the gas masses of these targets. More promising strategies to reach lower masses include measuring CO in gravitationally-lensed galaxies, where even modest magnification factors of 3 − 5 would lead to significant gains in efficiency, or deriving M mol from the Rayleigh-Jeans dust-continuum, which has proven to be detectable with ALMA in shorter integration times than are required for CO (e.g., Kaasinen et al. 2019;Aravena et al. 2020;Suzuki et al. 2021;Shivaei et al. 2022).

SUMMARY AND CONCLUSIONS
We have analyzed ALMA observations of CO(3−2) for a sample of 13 moderate-mass main-sequence galaxies at z ∼ 2.3 that uniquely have near-infrared spectra from the MOSDEF survey covering the full set of strong rest-optical emission lines, from which we have derived gas-phase metallicities. We supplemented the MOSDEF-ALMA sample with a sample of more massive z = 2−3 main-sequence galaxies from the literature with existing CO and metallicity constraints. The combination of cold gas content and metallicity information provides a powerful tool to constrain baryon cycling in the high-redshift universe. Our main conclusions are as follows: (i) We characterized the dependence of L CO(3−2) on M * , SFR, and O/H at z ∼ 2, finding that L CO(3−2) /SFR increases with O/H in a tight relation (Fig. 5). This result implies that CO luminosity per unit gas mass is lower in low-metallicity galaxies, carrying important implications for the conversion factor between CO luminosity and molecular gas mass.
(ii) We estimated the CO-to-H 2 conversion factor α CO for a sample of z = 1.1 − 3.2 galaxies with spectroscopic metallicity constraints using two techniques based on dynamical masses and the molecular KS law. With both methods, we found a significant dependence of α CO on O/H such that α CO increases with decreasing O/H (Fig. 6), a trend that is physically driven by less dust shielding and more intense UV radiation fields in lowmetallicity systems. The z ∼ 2 relation is consistent with z ∼ 0 α CO (O/H) relations within the uncertainties.
(iii) We found that M mol increases and µ mol decreases with increasing M * , while both properties display a strong secondary positive correlation with SFR at fixed M * (Fig. 7). The scaling relation of Tacconi et al. (2018Tacconi et al. ( , 2020 reproduces the observed M * dependence better than other published µ mol scaling relations (Scoville et al. 2017;Liu et al. 2019), suggesting it provides a more reliable extrapolation if applied to lower-mass z ∼ 2 samples. The molecular depletion timescales are 700 Myr on average and do not display any dependence on M * across log(M * /M ) = 10.2 − 11.4.
(iv) A tight near-linear molecular KS law exists at z ∼ 2 (Fig. 8, equation 14), providing a reliable means of estimating M mol indirectly from SFR and size measurements for large high-redshift star-forming samples.
(v) In the z ∼ 2 MOSDEF-ALMA sample, we found that residuals around the SFR-M * (star-forming main sequence), O/H-M * (mass-metallicity relation), and M mol -M * relations are correlated (Figs. 9 and 10). These correlations are such that, at fixed M * , galaxies with larger M mol have higher SFR and lower O/H. This result confirms the simultaneous existence of a SFR-FMR and Gas-FMR at z ∼ 2, both of which have been observed at z ∼ 0. These results suggest that the scatter of both the star-forming main sequence and massmetallicity relation are driven by stochastic variations in gas inflow rates that are traced by variations in gas fraction at fixed M * . Better gas mass constraints spanning a wider range of M * and SFR are required to obtain the quantitative form of the high-redshift Gas-FMR and its evolution across cosmic history, which could provide high-order tests of gas accretion and feedback models in numerical simulations.
(vi) We used gas-regulator models to infer the outflow mass-loading factors of the z ∼ 2 sample, finding η out = 1 − 4 with a typical value of 2.5 with no significant dependence on M * above 10 10.2 M (Fig. 11). A similar flattening of η out at high M * is observed at z ∼ 0. We conclude that the high-mass flattening of the MZR is driven by a flattening in η out , which in turn may be caused by the onset of cyclical AGN feedback in high-mass galaxies.
(vii) We performed a metal budgeting analysis, estimating the mass of metals found inside massive z ∼ 2 main-sequence galaxies in the gas-phase ISM, dust, and stars (Fig. 13). Comparing the sum of these phases to the total metal mass produced over a galaxy lifetime showed that massive z ∼ 2 star-forming galaxies retain only 30% of produced metals, implying that two thirds of produced metals are ejected into the CGM and/or IGM and demonstrating the important role of feedbackdriven outflows in removing metals from galaxies. Future constraints on CGM metal masses at z ∼ 2 can determine whether the missing metals can be accounted for.
(viii) We find that the M mol , µ gas , and t depl scaling relations of Tacconi et al. (2018) provide the best overall match to the z ∼ 2 CO samples in terms of M * and SFR dependence, suggesting that these relations can be extrapolated to lower masses without incurring large systematic biases. The µ gas scaling relations of Scoville et al. (2016) and Liu et al. (2019) are too steep with M * such that they likely overestimate the molecular gas content of low-mass systems at z ∼ 2.
(ix) We explored the possibility of detecting CO(3−2) in low-mass z ∼ 2 galaxies with ALMA (Fig. 14). Mainsequence galaxies at M * = 10 10 M can be detected in 5 hours on source, while starbursts a factor of 5 times above the main-sequence can be reached with the same depth down to 10 9.5 M . Detecting CO at lower masses is inefficient with ALMA due to the combined effect of decreasing M mol with decreasing M * and large α CO at low metallicity. Targeting lensed galaxies for CO and deriving M mol from dust continuum observations present promising possibilites to reach low masses at high redshift.
These results demonstrate the power of combining measurements of metallicity and gas mass to constrain baryon cycling for high-redshift galaxies. Existing samples at z > 2 with these measurements are small and largely have low-S/N gas measurements (CO or dust continuum) that often require stacking to produce significant detections. Combining samples from the literature potentially introduces systematic effects due to the variety of selection criteria and rest-optical line ratios available for metallicity determinations. The quest for a detailed understanding of baryon cycling at high redshifts would greatly benefit from a deeper survey of CO emission in a large and systematically selected sample of galaxies with accompanying rest-optical spectroscopy. ALMA currently has the capabilities to obtain the required measurements for a large sample at z ∼ 2.
This paper makes use of the following ALMA data: ADS/JAO.ALMA#2018.1.01128.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. Support for this work was provided by NASA through the NASA Hubble Fellowship grant #HST-HF2-51469.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astron-omy, Incorporated, under NASA contract NAS5-26555. We acknowledge support from NSF AAG grants AST-1312780, 1312547, 1312764, and 1313171, archival grant AR-13907 provided by NASA through the Space Telescope Science Institute, and grant NNX16AF54G from the NASA ADAP program. We also acknowledge a NASA contract supporting the "WFIRST Extragalactic Potential Observations (EXPO) Science Investigation Team" (15-WFIRST15-0004), administered by GSFC. We additionally acknowledge the 3D-HST collaboration for providing spectroscopic and photometric catalogs used in the MOSDEF survey. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.

A. METALLICITY CALIBRATIONS
The set of calibrations used to convert strong-line ratios to metallicities deserves careful consideration in analyses where metallicity plays a key role, especially when a uniform set of emission lines is not available for the entire sample. We adopt calibrations derived from the local sample of z ∼ 2 analogs from Bian et al. (2018) that Sanders et al. (2020b) found to be a good match to actual z ∼ 2 galaxies with direct-method abundances for line ratios of [O iii], [O ii], [Ne iii], and Hβ. Here, we describe some important differences in our use of the Bian et al. (2018) calibration set relative to what is reported in that work. Table A1 provides the metallicity calibrations adopted in this work as polynomial coefficients for use in the following equation: where R is the line ratio under consideration, x=12+log(O/H), and c n are the coefficients. Differences from the values in Bian et al. (2018) (Storey & Zeippen 2000).
While strong-line ratios involving the [N ii]λ6584 line are commonly used to derive the oxygen abundance (e.g., N2 and O3N2; Pettini & Pagel 2004;Marino et al. 2013), N-based calibrations are sensitive to N/H such that differences in N/O at fixed O/H between the calibrating and target samples will systematically bias the O/H estimates. Unfortunately, Sanders et al. (2020b) were unable to test the N2 and O3N2 analog calibrations from Bian et al. (2018) due to the small number of galaxies at z > 1 with both direct-method metallicities and [N ii] detections. Figure  , and Hβ (O3O2Ne3). Gray squares show composite spectra of z ∼ 2.3 star-forming galaxies from Sanders et al. (2021). Small circles display individual z ∼ 2.3 galaxies, color-coded by the difference in [N ii]/[O ii] relative to the average value at fixed M * . The solid black line denotes a one-to-one relation. The dashed pink line presents the best-fit offset between each pair of metallicities, as fit to the composite spectra. The N2 and O3N2 calibrations given in Table A1 have been renormalized such that N2, O3N2, and O3O2Ne3 yield consistent metallicities on average for the z ∼ 2.3 sample.
For reference, the leading constant factors of the original Bian et al. (2018) calibrations are 8.82 and 8.97, respectively. The coefficients in Table A1 match these renormalized relations and yield consistent metallicities on average from N2, O3N2, or O3O2Ne3 for z ∼ 2.3 galaxies.  Figure A1 can be explained if the Bian et al. (2018) analog calibrating sample has higher N/O at fixed O/H than the z ∼ 2.3 sample. Figure  3 star-forming galaxies and highlights the care that is needed when using local analogs in place of actual high-redshift samples. These analogs were selected to lie in the same region of the [N ii] BPT diagram as z ∼ 2 samples. This criterion has clearly failed to select galaxies with identical ISM conditions to those at z ∼ 2, consistent with recent work that found local and high-redshift galaxies occupying the same region of the BPT diagram maintain some distinct differences in their ISM conditions (Runco et al. 2020). Improving metallicity calibrations for use at high redshifts will require more careful selection of local analogs and, ultimately, an expanded sample of high-redshift sources with directmethod metallicities that will soon be enabled by the James Webb Space Telescope.

B. PARAMETERIZATION OF THE EVOVLING MASS-METALLICITY RELATION
We utilize recent observational determinations of the mass-metallicity relation at different redshifts to obtain a functional form for the evolving mass-metallicity relation parameterized by M * and z. We use the highredshift measurements of Sanders et al. (2021) based on composite spectra of z ∼ 2.3 and z ∼ 3.3 star-forming galaxies from the MOSDEF survey. Metallicities of the high-redshift samples are based on O3O2Ne3 and de- rived using the Bian et al. (2018) z ∼ 2 analog calibrations. We add the binned data of Curti et al. (2020) for a sample of ∼150,000 star-forming galaxies at z med = 0.08 from the Sloan Digital Sky Survey (SDSS) to set the local baseline. Metallicities of the local sample were determined using calibrations based on stacked spectra of normal z ∼ 0 galaxies derived in Curti et al. (2020). Both calibration sets are on a T e -based metallicity scale. We fit the z ∼ 0, z ∼ 2.3, and z ∼ 3.3 data with a function of the form MZR(M * , z) = Z 0 −γ/β ×log(1+[M * /M 0 (z)] −β ) (B4) where Z 0 is the high-mass asymptotic metallicity, γ is the low-mass power law slope, and β controls the sharpness of the turnover that occurs at the mass log(M 0 (z)/M ) = m 0 + m 1 * log(1 + z). This is the same functional form used to fit the z ∼ 0 SFR-FMR by Curti et al. (2020), but with the secondary parameter being (1 + z) instead of SFR. We find the following best-fit parameters: Z 0 = 8.80, γ = 0.30, β = 1.08, m 0 = 9.90, and m 1 = 2.06. This fitting function yields a nearly identical relation to the Curti et al. (2020) bestfit z ∼ 0 MZR when evaluated at z = 0.08. Figure B3 shows the best-fit MZR(M * , z) alongsize the z ∼ 2.3 and z ∼ 3.3 data (green and blue, respectively; z ∼ 0 data is not displayed for clarity), displaying a good fit across the entire mass range. To test the applicability of this fitting function across a wider range of redshifts, we compare to samples at z ∼ 0.8 (Zahid 8.5 9  Figure B3. The mass-metallicity relation for samples of star-forming galaxies at a range of redshifts: z ∼ 0.8 (Zahid et al. 2011), z ∼ 1.5 (Topping et al. 2021), and z ∼ 2.3 and 3.3 . Solid black lines show the best-fit MZR(M * , z) to data at z = 0.08 (Curti et al. 2020), z ∼ 2.3, and z ∼ 3.3, assuming the functional form given in equation B4. Additional black lines display the best-fit MZR(M * , z) evaluated at z = 0.8 and z = 1.5 in comparison to observational MZR data at these redshifts that were not used in the fitting process, demonstrating that the best-fit function is robust across the intermediate redshift range as well. Dashed gray lines display the MZR(M * , z) function of Genzel et al. (2015), evaluated at z = [0.8, 1.5, 2.3, 3.3].
et al. 2011, DEEP2) and z ∼ 1.5 (Topping et al. 2021, MOSDEF). The metallicities of the z ∼ 0.8 sample have been reevaluated using the R 23 and O 32 z ∼ 0 calibrations of Curti et al. (2020), while those of the z ∼ 1.5 sample have been rederived with the renormalized highredshift N2 calibration given in Table A1. Our best-fit MZR(M * , z) function provides a good match at intermediate redshifts as well, demonstrating that it is widely applicable over z = 0 − 3.3 for galaxies above 10 9 M . The gray dashed lines show the MZR(M * , z) derived by Genzel et al. (2015) and adopted by Tacconi et al. (2018Tacconi et al. ( , 2020, evaluated at z = [0.8, 1.5, 2.3, 3.3]. Our new relation yields higher metallicity at fixed M * and redshift across the full mass range (0.15 dex in O/H at 10 10.5 M and z = 2.3). The difference primarily stems from the use of larger and more representative highredshift datasets in this work (see Sanders et al. 2021 for a discussion) and different choices of metallicity calibrations, in particular the use of different calibrations at low and high redshifts in this work to reflect evolving ISM ionization conditions.