The Stellar Mass versus Stellar Metallicity Relation of Star-forming Galaxies at 1.6 ≤ z ≤ 3.0 and Implications for the Evolution of the α-enhancement

The zCOSMOS-deep redshift survey has observed around 10⁴ galaxies in the central ∼0.8 deg² of the COSMOS field (Scoville et al. 2007). Here we provide a brief description and refer the reader to Lilly et al. (2007) and Kashino et al. (2021) for more details.

The observations were carried out using VLT/VIMOS (Le Fèvre et al. 2003) with the low-resolution blue grism with 10 slits, yielding a spectral resolution of R ∼ 200 and a spectral coverage of ≈3600–6700 Å. The selection of the targets was performed based on a then-current version of the COSMOS photometric catalog. All of the objects were color-selected through a BzK (Daddi et al. 2004) or ugr (Steidel et al. 2004) method with a blue magnitude cut B_AB < 25.25. These selection criteria isolate star-forming galaxies in a range 1.4 ≲ z ≲ 3.0 (Lilly et al. 2007). Redshifts were visually inspected in 2D and 1D reduced spectra by identifying multiple prominent spectral features in the rest-frame FUV window or Lyα emission line and break.

Figure 1 shows six representative examples of individual galaxy spectra spanning a range of redshifts and the confidence class. Note that the additional flux calibration described in Section 3.1 has been applied to these spectra. The wavelength regions without data are those impacted by strong sky lines or zeroth-order contamination. Those of class = 3 or 4 present robust redshift identification based on several strong absorption lines and/or a strong Lyα emission line. The identification of such spectral features becomes less secure for class = 2 spectra.

Zoom In Zoom Out Reset image size

We constructed the sample used in this paper from the full catalog of the zCOSMOS-deep survey (S. J. Lilly et al., in preparation). The sample is limited to those having a clear photometric counterpart in the COSMOS2015 photometric catalog (Laigle et al. 2016). Galaxies detected in X-rays (Civano et al. 2016) are also excluded to remove possible active galactic nuclei (AGNs) from the sample.

The redshift range for the current analysis is limited to 1.6 ≤ z ≤ 3.0 so that the VIMOS spectrum covers the range of λ_rest ≈ 1400–1700 Å for all the galaxies. We adopt all objects with a very secure zCOSMOS-deep redshift (Confidence Class = 3 or 4 ¹⁵ ) within the redshift range. For those with Class = 2, we use only those that are consistent to within ∣z_phot − z_spec∣/(1 + z_spec) ≤ 0.1 of the photometric redshift in the COSMOS2015 catalog. ¹⁶

The redshifts in these two categories are both estimated to be ≥99% reliable (Lilly et al., in preparation). We do not use any of the objects with less secure redshifts, nor any of those with broad emission lines (i.e., Class +10). Finally, we visually inspected individual sources and excluded 34 sources (2.5% of the remaining sample) for which the FUV continuum is barely detected (which are likely to be Lyα emitters with very faint continua), or that show obvious AGN signatures (e.g., strong Nv λ1240 emission lines), or that suffered from severe spectral contamination from companion sources that fall in the same slit.

The final sample consists of 1336 galaxies. Figure 2 shows the redshift distribution of the sample with the median redshift ${\left\langle z\right\rangle }_{\mathrm{med}}=2.22$. We note that the redshifts were all determined with strong absorption lines owing to mostly carbon and/or silicon ions in the ISM, but not with any kind of stellar iron lines. Here we ignore any possibility that our sample with secure redshifts may be biased toward those with strong iron features, although the strengths of the ISM absorption lines could be correlated at some level with the overall gas-phase metallicity and thus the ISM absorption strengths (Faisst et al. 2016).

Zoom In Zoom Out Reset image size

We derived stellar masses (M_*) for individual galaxies through spectral energy distribution (SED) fitting based on the photometry from the COSMOS2015 catalog, together with the precise spectroscopic redshifts. It should be noted that the stellar mass here denotes the mass of living stars at the time of observation, rather than the integral of the star formation rate (SFR). The SED fitting outlined in this section is done to the broadband photometry and is meant to derive stellar masses. We will independently perform another fitting to the composite FUV spectra using high-resolution model spectra to estimate the stellar metallicities (see Section 3.5).

Our SED fitting procedure uses the photometric fluxes measured within 32 broad-, intermediate-, and narrowband filters from GALEX near-UV to Spitzer/IRAC ch4, as listed in Table 3 of Laigle et al. (2016). For CFHT, Subaru, and UltraVISTA photometry, we used the fluxes measured in a 3''-diameter aperture and applied offsets provided in the catalog to convert them to the total fluxes. All the photometric bands whose rest-frame central wavelengths are within 1960 Å≤ λ_cen/(1 + z) ≤ 2440 Å were excluded in order to ensure that the SED fitting is not affected by the possible 2175 Å bump feature in the SED of the galaxies (see Kashino et al. 2021).

The stellar masses and SFRs are estimated using a Python code CIGALE (Burgarella et al. 2005; Noll et al. 2009; Boquien et al. 2019) for the SED fitting and adopted the stellar population synthesis models of Bruzual & Charlot (2003) with a Chabrier (2003) IMF. We considered delayed SFHs with an additional recent burst in order to model the long-term star formation that has formed the bulk of the stellar mass and the latest episode of star formation (e.g., Ciesla et al. 2016, 2017; Pearson et al. 2017),

$$\begin{eqnarray}&&\mathrm{SFR}(t)={\mathrm{SFR}}_{\mathrm{delayed}}(t)+{\mathrm{SFR}}_{\mathrm{burst}}(t),\end{eqnarray} \tag{ 1 }$$

where ${\mathrm{SFR}}_{\mathrm{delayed}}(t)\propto {{te}}^{-t/{\tau }_{0}}$ and ${\mathrm{SFR}}_{\mathrm{burst}}(t)\propto {e}^{-(t-{t}_{1})/{\tau }_{1}}$ if t > t₁ and SFR_burst(t) = 0 otherwise. The parameter t denotes the elapsed time since the onset of star formation, t₁ the galaxy age when the late episode of star formation onsets, and τ₀ and τ₁ the e-folding times of the main stellar population and the late starburst population.

In CIGALE, we also accounted for the effect on the photometry of nebular emission lines assuming a common ionization parameter $\mathrm{log}U=-2.8$ and dust emission based on the templates from Dale et al. (2014). The dust attenuation is accounted for using the prescription of Charlot & Fall (2000). The full list of the input parameters is presented in Table 1.

Table 1. Input Parameters of the SED Fitting with CIGALE

Parameter	Values
Delayed+Burst SFH
Age (main population) (Myr)	1000–4500 in steps of 250
e-folding time of the delayed SFH, τ0 (Myr)	1000, 2000, 3000
Age (starburst population) (Myr)	50, 100, 150, 200
e-folding time of the late starburst, τ1 (Myr)	10,000
Mass fraction of the late burst population	0.001, 0.003, 0.010, 0.020, 0.040, 0.100, 0.200, 0.400
Stellar Population: Bruzual & Charlot (2003)
Initial mass function	Chabrier (2003)
Metallicity, Z	0.008
Separation age between young and old populations (Myr)	10
Dust Attenuation: Charlot & Fall (2000)
AV in the diffuse ISM	0.0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8,
	0.9, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 3.5
δISM (power-law slope for the ISM)	−0.7
δBC (power-law slope for birth clouds)	−1.3
${A}_{V}^{\mathrm{ISM}}/({A}_{V}^{\mathrm{BC}}+{A}_{V}^{\mathrm{ISM}})$	0.44

Download table as:  ASCIITypeset image

It is known that SED fitting often leads to unrealistically young ages (a few × 100 Myr), and thus spuriously low stellar masses, for high-z star-forming galaxies when there is no limitation on the age. This is because the SED of such galaxies is dominated by the youngest stellar populations (e.g., Maraston et al. 2010). Moreover, it is not plausible that ∼ 10¹⁰ M_⊙ galaxies at z ∼ 2–3 have started to form stars just before we happen to observe them, given the fact that we can find star-forming galaxies continuously at higher redshifts. It has been shown that limiting the minimum age to a Gyr or more better recovers the stellar masses of high-z star-forming galaxies (Maraston et al. 2010; Pforr et al. 2012). To avoid artificially inferring unrealistically young ages, we thus limit the age of the main stellar population to be ≥ 1 Gyr. This corresponds to a formation redshift z ≳ 5 for galaxies at z ∼ 3 and should therefore not be far from reality.

The estimated stellar masses span the range $9.23\leqslant \mathrm{log}{M}_{* }/{M}_{\odot }\leqslant 10.85$ (the 2.5th–97.5th percentiles) with a median value of ${\left\langle \mathrm{log}({M}_{* }/{M}_{\odot })\right\rangle }_{\mathrm{med}}=10.0$. The reduced χ² values range mostly between 0.4 and 3.3 with the median value of 0.95. Figure 3 shows the distribution of the inferred M_* and instantaneous SFR. Note that there is a (spurious) sharp upper boundary of the SFR distribution in Figure 3. This is produced by the limited range of SFR probed by the adopted SFHs. This should be of no consequence, as we do not use the SFR values of the individual galaxies in the following analyses.

Zoom In Zoom Out Reset image size

For reference, the main-sequence M_*–SFR relations at similar redshifts are taken from Whitaker et al. (2014, z = 2.0–2.5). We also plot the local main sequence adapted from Renzini & Peng (2015). Overall, our sample is in good agreement with the main sequence at the earlier epoch, indicating that it should be a representative sample at these redshifts. We note that the magnitude limit, B_AB < 25.25, for the parent sample emerges as a cut in SFR at SFR ∼ 5 M_⊙ yr⁻¹ in Figure 3 (right panel). At $\mathrm{log}{M}_{* }/{M}_{\odot }\lesssim 9.6$, we are missing galaxies with lower SFRs owing to their low observed flux, and thus the results may be biased at the low-mass end.

We adopt a spectrophotometric correction to every single spectrum of our sample following the method described in Kashino et al. (2021). The spectra that are produced through the standard zCOSMOS-deep reduction pipeline were flux-calibrated based on standard-star observations. The nominal flux calibration, however, cannot perfectly correct for the effects of finite slit width, imperfect slit centering, and the effects of atmospheric dispersion.

We correct each spectrum with a smooth function of wavelength (see Equation (8) of Kashino et al. 2021) that is constructed using the differences between the actual photometric fluxes in four broadband and eight intermediate-band filters (shown in Figure 4) and the fluxes computed from the pipeline-processed spectrum in these same 12 filter bands. In doing so, we excluded any photometric bands that sampled the rest-frame wavelength of the Lyα emission line. This is because the flux of the strong Lyα line may be differently affected than the continuum flux owing to the possible extended shape of the emission and/or the overestimate of the sky level at these particular wavelengths.

Zoom In Zoom Out Reset image size

Examples of the flux-calibrated spectra are shown in Figure 1, demonstrating that the spectra are in agreement with the photometric fluxes. The left panel of Figure 4 shows how this spectrophotometric calibration works in the observed frame by comparing the spectral broad- and intermediate-band fluxes (${F}_{i}^{\mathrm{spec}}$, where i denotes a filter) recomputed in the flux-corrected spectra with the photometric fluxes (${F}_{i}^{\mathrm{phot}}$) from the COSMOS2015 catalog. The median values of the corrected ${F}_{i}^{\mathrm{spec}}$-to-${F}_{i}^{\mathrm{phot}}$ ratios in each band are all within ±0.015 dex, and the scatter is ∼ 0.02–0.05 dex, depending on the filters. This overall scatter (seen in the rest frame) is quite smaller than before correction ( ∼ 0.1–0.15 dex). Note that there is still scatter for a given filter because the correction used a smooth function of wavelength, so the effects of photometric noise in the different filters are still seen.

The right panel of Figure 4 shows the same data but now shifted to the rest-frame wavelengths. There is no systematic trend in the corrected flux ratios across the entire wavelength range for the stellar metallicity measurement (1274–2053 Å). The running medians with a window size of 100 Å are all within 0.005 dex over this wavelength range of interest. On the other hand, a significant systematic undercorrection is seen around 1216 Å most likely because of sky-subtraction issues associated with strong Lyα lines mentioned above.

This flux calibration is important to obtain composite spectra that correctly reflect the average shape of the galaxies' SEDs. We note, however, that the precision of this calibration is unlikely to be critical for our conclusions, because the stellar metallicities are measured based on the detailed shape of the spectra that results from the blending of numerous narrow stellar absorption lines, whereas the overall shape of the smooth continuum is fit with an arbitrary multiplicative λ-dependent function (see Section 3.5). In the remainder of the paper, the term "observed VIMOS spectra" will always refer to these accurately spectrophotometrically recalibrated spectra.

Precise determination of the spectroscopic redshifts is important to reduce the loss of the potential spectral resolution when spectra are stacked. The spectroscopic redshifts in the zCOSMOS-deep catalog have been determined by visually inspecting prominent emission and absorption features in each spectrum, but no systematic spectral fitting has been performed. As a consequence, the spectral redshifts may be more uncertain than the best estimates that can potentially be achieved from the existing spectra; thus, some improvements are possible.

We therefore made small adjustments to the spectroscopic redshifts by fitting a common template spectrum to each of the individual spectra. We constructed this template by stacking the observed spectra of the entire sample of 1336 galaxies used in this work by using the original spectroscopic redshifts from the catalog. In this template fitting, we applied an arbitrary normalization and a wavelength-dependent multiplicative factor that is intended to mimic the Calzetti et al. (2000) dust attenuation.

The differences between the revised and original spectroscopic redshifts, (z_orig − z_new)/(1 + z_new), have a Gaussian-like distribution with a standard deviation of 8.4 × 10⁻⁴ (252 km s⁻¹) and median of 6.4 × 10⁻⁶. We found that, by using the revised spectroscopic redshifts, the resulting stacked spectra are noticeably improved, showing sharper spectral features in both emission and absorption than those seen in the stacks based on the original redshifts. However, the adjustment of the spectroscopic redshifts has little effect on the metallicity measurements, and our conclusions do not change if the stacking is done using the original redshift values.

To infer the stellar metallicity as a function of stellar mass, we rely on stacked spectra from subsamples of galaxies defined by their stellar mass. The observed spectra are co-added as follows. We first transform all the individual spectra to the rest-frame wavelength based on their adjusted redshift (see Section 3.2) and rebin them to a common wavelength grid with a spacing of 1 Å. This is slightly oversampled relative to the pixel resolution of the original spectra (5.35 Å pixel^–1 in the observed frame) and equivalent to the wavelength grid spacing of the model spectra used for spectral fitting (see Section 3.4).

Each spectrum is then normalized by dividing by a fitted continuum of the form ${\lambda }^{\beta }\times {10}^{-0.4{A}_{\lambda }}$, where A_λ is the dust attenuation of the Charlot & Fall (2000) prescription. Here the overall dust attenuation is a free parameter for continuum fitting, independent of the result from the SED fitting in Section 2.3. By doing so, the effects of the variable overall shapes of the spectra are mitigated. We then take the mean value of the individual continuum-divided spectra at each wavelength grid while ignoring any spectral regions that are missing and/or contaminated, for example, by zeroth-order contamination or strong sky lines. Finally, to recover the global shape of the average SED, the stacked spectrum is multiplied by the mean of the fitted continua. We note that normalizing the individual spectra to their continuum level at λ_rest ∼ 1550 Å produces very similar results. Stacking without any normalization, however, yields more noisy composite spectra and more residuals with respect to the fitted models.

Figure 5 shows the composite spectrum (middle panel) of the entire sample. ¹⁷ The number of spectra used at each wavelength is indicated in the top panel. Note that, given the redshift range of our sample, the rest-frame wavelength range of 1400–1700 Å is covered by nearly all the input spectra, while shorter and longer wavelengths are less well sampled. Some prominent spectral features are clearly identified as marked by vertical lines.

Zoom In Zoom Out Reset image size

To derive stellar metallicities for the galaxies, we compare our observed composite spectra with model spectra, following the approaches described in Steidel et al. (2016) and Cullen et al. (2019).

We utilize the latest public data release of the population synthesis code "Binary Population and Spectral Synthesis" (BPASSv2.2.1; Eldridge et al. 2017; Stanway & Eldridge 2018). The package provides sets of single stellar population synthesis model spectra at a pixel resolution 1 Å pixel⁻¹ as a function of stellar ages for different IMFs and for discrete (overall) stellar metallicities (Z_* = 10⁻⁵, 10⁻⁴, 0.001, 0.002, 0.003, 0.004, 0.006, 0.008, 0.010, 0.014, 0.020, 0.030, and 0.040). It is important to note that although the BPASS models adopt a fixed abundance ratio based on the solar abundances, the model fit to the rest-frame FUV spectrum is mostly sensitive to the iron abundance of the young stellar component of the galaxies. This is because stellar photospheric absorption features in the wavelength range of interest are primarily due to transitions of highly ionized iron (Dean & Bruhweiler 1985; Brandt et al. 1998). We can therefore translate the inferred metallicity Z_* into [Fe/H].

We adopted the BPASSv2.2.1 model that used a Chabrier (2003) IMF with a high-mass cutoff of 100 M_⊙ and included binary star evolution, and we considered a continuous SFH with a duration of 100 Myr to construct a set of rest-frame UV template spectra. Note that the FUV spectrum is almost independent of the past SFH at ≳ 100 Myr previously, as it is dominated by massive young stars. Indeed, we checked that the results hardly change for different duration times between 10 and 1000 Myr. Figure 6 shows examples of the model spectra at different metallicities, demonstrating that the detailed shapes of the galaxy spectra across the rest FUV window are sensitive to the stellar metallicity. A feature of the BPASS binary models is that they model the broad He λ1640 emission line that originates in the winds from very massive stars (see Section 3.5 for a relevant description).

Zoom In Zoom Out Reset image size

We added nebular continuum and line emission to the templates, although it has only a minor contribution of ∼5% in the total continuum flux. The nebular continuum was computed using cloudy v17.00 (Ferland et al. 2017) adopting the BPASS spectrum itself as the incident spectrum. We assumed the electron density to be n_e = 300 cm⁻³ and the ionization parameter to be $\mathrm{log}U=-2.8$. These values are consistent with the recent estimates in z ∼ 2–3 star-forming galaxies (Sanders et al. 2016; Kashino et al. 2017; Strom et al. 2017, 2018) and local star-forming galaxies with specific SFR (sSFR) as high as our sample (∼5 Gyr⁻¹; Kashino & Inoue 2019). We change the gas-phase metallicity to calculate nebular emission according to the stellar metallicity, with an offset of $\mathrm{log}({Z}_{\mathrm{gas}}/{Z}_{* })=0.42$ as reported for z ∼ 2–3 galaxies by Strom et al. (2018). In other words, we considered an enhanced gas-phase metallicity, $12+\mathrm{log}({\rm{O}}/{\rm{H}})=8.69+\mathrm{log}({Z}_{* }/{Z}_{\odot })+0.42$, relative to each Z* of the BPASS template. We will indeed measure consistent $\mathrm{log}({Z}_{\mathrm{gas}}/{Z}_{* })$ or [O/Fe] values in Section 4.3 and discuss the evolution of [O/Fe] in Section 5.

Our goal is to estimate average stellar metallicities for binned subsamples of galaxies by using their composite spectra. To do so, we limited the wavelength range between 1274 Å ≤ λ_rest ≤ 2053 Å and excluded some narrow wavelength regions where the stellar continuum is impacted by interstellar absorption lines or nebular emission lines. The wavelength regions excluded from the fitting are summarized in Table 2. We note that the region of the He λ1640 line is not excluded because the BPASS models include the broad He λ1640 stellar emission line that originates from very hot stars. In fact, as shown in Section 4, the observed He λ1640 feature in the composite spectra could be attributed almost entirely to the contribution from stars that is predicted in the BPASS models.

Table 2. Rest-frame Wavelength Ranges Excluded from Fitting ^a

${\lambda }_{\min }$	${\lambda }_{\max }$	Interstellar b
(Å)	(Å)	Spectral Features
1294	1318	O i λ1302.17, Si ii λ1304.37, Si ii * λ1309.28
1326	1343	C ii λ1334.53
1353	1361	O i λ1355.60, O i λ1358.51
1385	1411	Si iv λλ1393.76, 1402.77
1518	1535	Si ii λ1526.71
1540	1559	C iv λλ1548.19, 1550.77
1600	1617	Fe ii λ1608.45
1658	1679	O iii] λλ1660.81, 1666.15, Al ii λ1670.79
1707	1712	Ni ii λ1709.60
1739	1744	Ni ii λ1741.55
1749	1754	Ni ii λ1751.91
1758	1763	C ii λ1760.5 (blended)
1806	1811	Si ii λ1808.01
1814	1820	Si ii λλ1816.93, 1817.45
1846	1871	Al iii λ1854.72, 1862.79
1880	1885	Si iii] λ1882.71
1890	1895	Si iii] λ1892.03
1898	1917	C iii] λ1906.68, 1908.73

Notes.

^a The wavelength range used for the fit is limited to λ_rest = 1274–2053 Å. ^b Wavelengths are all given in vacuum.

Download table as:  ASCIITypeset image

An added advantage of the full spectral fitting over using the traditional indices (e.g., the 1978 Å index) is that the former does not need to identify the pseudosmooth continuum. Contrarily, the latter requires it for measuring equivalent widths (EWs) of particular stellar absorption lines. The measurements of the EWs of faint features are quite sensitive to the assumed continuum level, while the determination of the pseudocontinuum has to be based on very limited wavelength regions free from narrow features, thus being easily affected by noise. The full spectral fitting is free from these difficulties.

We fit the constructed BPASS models to our composite spectra based on the maximum likelihood estimation. We assume that the associated errors are Gaussian and independent, so the logarithm of the likelihood ${ \mathcal L }$ is given by

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal L }=-\displaystyle \frac{1}{2}\sum _{i}\left[\displaystyle \frac{{({f}_{i}^{\mathrm{obs}}-{f}_{i}^{\mathrm{model}})}^{2}}{{\sigma }_{i}^{2}}+\mathrm{ln}(2\pi {\sigma }_{i}^{2})\right],\end{eqnarray} \tag{ 2 }$$

where ${f}_{i}^{\mathrm{obs}}$ is the observed composite spectrum at the ith wavelength grid, ${f}_{i}^{\mathrm{model}}$ is the model spectrum for a given set of parameters, and σ_i is the error on the observed flux. The summation is over all wavelength grids that are included in the fit.

For our models, we adopted eight free parameters: the logarithm of the overall stellar metallicity ($\mathrm{log}{Z}_{* }$), the Gaussian smoothing kernel σ_smooth applied to the model spectra, and six parameters for a fifth-order polynomial function that is multiplied by the template spectrum:

$$\begin{eqnarray}&&{f}_{i}^{\mathrm{model}}({Z}_{* })={f}_{i}^{\mathrm{BPASS}}({Z}_{* })\times \sum _{k=0}^{5}{p}_{k}{\left(\lambda [\mathring{\rm A} ]-1650\right)}^{k},\end{eqnarray} \tag{ 3 }$$

where ${f}_{i}^{\mathrm{BPASS}}({Z}_{* })$ is the template spectrum at given Z_* in arbitrary units, in which the nebular emission is accounted for. The polynomial term is ideally intended to reflect the change of the overall shape of the spectra due to dust attenuation. ¹⁸

To sample the posterior probability distribution of the model parameters, we employed a Markov Chain Monte Carlo (MCMC) technique using the emcee package for Python (Foreman-Mackey et al. 2013). In our fitting, we adopted a uniform prior probability function for each parameter. For stellar metallicity, a flat prior is adopted in log space (i.e., $\mathrm{log}({Z}_{* })$) between $-5\leqslant \mathrm{log}({Z}_{* })\leqslant \mathrm{log}(0.040)$ (which corresponds to $-3.15\leqslant \mathrm{log}({Z}_{* }/{Z}_{\odot })\leqslant 0.44$) according to the possible choices in the BPASS models. Since the models are only provided at 13 discrete metallicity values (see Section 3.4), we interpolated the flux values in $\mathrm{log}{Z}_{* }$–$\mathrm{log}(\mathrm{flux})$ space between the models. This enables us to generate a model at any metallicity value within the available range.

In the fitting, we need to match the spectral resolution of the model templates to that of the observed composite spectra. However, the spectral resolution of the composite spectra is not precisely known. It may vary from galaxy to galaxy within a subsample because of different velocity dispersions of the galaxies and different slit illumination profiles. We thus leave the smoothing scale σ_smooth, in units of Å, ¹⁹ as a free parameter rather than using a predetermined value. In the fitting, the model spectra are smoothed by a Gaussian kernel with σ_smooth, and the smoothing scale that gives the best fit to the data is constrained along with other parameters.

Before applying our procedure to the data, we evaluate the potential systematic uncertainties in the metallicity measurements that may come from the noise and the limited spectral resolution of the spectra.

As we will show below, the width of the smoothing kernel, σ_smooth, is estimated from the fits to be ∼ 3.3 Å. This is large enough to wash out individual narrow stellar absorption features in the model spectra. Still, the model spectra with different metallicities will retain their different characteristic shapes that result from blending of the narrow features after smoothing. Thus, the metallicity measurement is indeed possible. The possible systematic uncertainties, however, should be evaluated.

For this purpose, we constructed a set of artificial spectra using the BPASS models at different metallicities. We smoothed them with σ_smooth = 3.3 Å and modified the overall shape using the representative polynomial function as in Equation (3). We then added Gaussian noise with the moderate wavelength dependence that replicates the 1σ error spectra. We considered noise with four different levels: the best had the noise seen in the composite spectrum of the entire sample, i.e., the S/N at λ_rest ∼ 1450 Å, S/N₁₄₅₀ ≈ 190 per 1 Å. The other three are degraded to have twice (S/N₁₄₅₀ ≈ 95), five times (S/N₁₄₅₀ ≈ 38), or 10 times (S/N₁₄₅₀ ≈ 19) higher noise than the first one. The stacked spectra of the subsamples selected by stellar mass have noise levels between these values, depending on the number of galaxies in the bin and their average brightness.

We then attempted to measure the stellar metallicity in the same way as from the data in order to evaluate how accurately the input metallicities can be recovered. We repeated each setup 10 times with different random noise realizations. Here only a single fiducial duration of star formation (100 Myr) is considered.

Figure 7 shows the results of this exercise, separately, for four different noise levels. The difference between the inferred and input metallicities is shown as a function of the input value. The error bars indicate the 16th–84th percentiles of the inferred posterior for individual measurements. The results from the 10 trials with different random noises are shown separately with slight offsets in the x-axis. We also compile the posterior probability distributions of $\mathrm{log}{Z}_{* }$ from all 10 trials and derive the median and the 16th–84th percentiles that are shown by blue solid and dashed lines.

Zoom In Zoom Out Reset image size

Generally, the uncertainties reduce with increasing metallicity because the stellar absorption features in the spectrum become more prominent and sensitive to the metallicity value at higher metallicity. It can be seen that, with all these noise levels, the systematic biases in the metallicity measurements are negligible with no trend across the metallicity range of interest.

To investigate the relation between stellar mass and stellar metallicity, we divided the sample of 1336 star-forming galaxies into bins of stellar mass in two different binning schemes: first we equally split the sample into six mass bins, and then we divided the sample into bins with a constant width of ${\rm{\Delta }}\mathrm{log}{M}_{* }=0.2\,\mathrm{dex}$. After stacking, we performed the MCMC analysis to fit the BPASS template to the stacked spectra as described in Section 3.5.

Figure 9 shows the composite spectra and the corresponding best-fit models in the six equally populated subsamples. It is noticeable that the overall slope of the spectra becomes shallower for higher M_* owing to the increasing average dust attenuation. The overall shape and detailed features are both well reproduced across the entire wavelength range of interest by the model spectra.

Zoom In Zoom Out Reset image size

The results in different bins are summarized in Table 3. The stellar metallicities are measured to range over $-1.3\,\lesssim \mathrm{log}({Z}_{* }/{Z}_{\odot })\lesssim -0.4$, increasing with M_*. Note that the BPASS models cover this observed range of metallicity with a sufficiently large margin at either end, and thus our MCMC analysis should not be affected by the artificial limit of the explored metallicity range. We obtained consistent results from the fixed-width binning scheme.

Table 3. Stack Statistics and Metallicity Estimates ^a

$\mathrm{log}\left({M}_{* }/{M}_{\odot }\right)$		N	S/N1450 b	$\mathrm{log}\left({Z}_{* }/{Z}_{\odot }\right)$ c
Median	Min–Max			([Fe/H])
9.97	8.38–11.48	1336	189	$-{0.812}_{-0.008}^{+0.008}$
Binned Equally into Six Subsamples in M*
9.45	8.38–9.62	222	71	$-{1.045}_{-0.029}^{+0.029}$
9.74	9.63–9.83	222	72	$-{0.901}_{-0.030}^{+0.029}$
9.89	9.83–9.97	223	81	$-{0.848}_{-0.018}^{+0.017}$
10.06	9.97–10.14	223	86	$-{0.785}_{-0.018}^{+0.017}$
10.24	10.14–10.34	223	82	$-{0.717}_{-0.017}^{+0.017}$
10.52	10.34–11.48	223	73	$-{0.662}_{-0.012}^{+0.016}$
Binned into the 0.2 dex Width Intervals in M*
9.06	8.91–9.10	13	13	$-{1.281}_{-0.230}^{+0.159}$
9.24	9.10–9.30	35	28	$-{1.273}_{-0.180}^{+0.143}$
9.40	9.30–9.50	82	46	$-{1.087}_{-0.412}^{+0.095}$
9.62	9.50–9.70	169	67	$-{0.970}_{-0.031}^{+0.031}$
9.82	9.70–9.90	257	83	$-{0.837}_{-0.017}^{+0.017}$
10.00	9.90–10.10	287	98	$-{0.792}_{-0.016}^{+0.016}$
10.20	10.10–10.30	221	86	$-{0.711}_{-0.017}^{+0.017}$
10.37	10.30–10.50	147	62	$-{0.684}_{-0.017}^{+0.011}$
10.59	10.50–10.70	64	42	$-{0.672}_{-0.021}^{+0.025}$
10.78	10.70–10.90	29	27	$-{0.611}_{-0.037}^{+0.037}$
11.00	10.91–11.10	17	20	$-{0.393}_{-0.078}^{+0.072}$

Notes.

^a The top row is for the stack of the entire sample. The middle set of rows is for six equally populated bins of M_*, and the bottom set of rows is for stacks in 0.2 dex fixed-width M_* bins. ^b The S/N per unit pixel (1 Å) of the stacked spectra around 1450 Å, represented by the median S/N within λ_rest = 1430–1470 Å. ^c Metallicity estimates, represented by the median value of the posterior probability distribution function obtained from the MCMC analysis. The associated errors correspond to the 16th–84th percentiles. Because the fit is sensitive to the iron abundance, this can be translated into [Fe/H].

Download table as:  ASCIITypeset image

Figure 10 shows the relationship between stellar mass and stellar metallicity. The x-axis values of the points correspond to the median stellar masses within the bin, and the horizontal error bars indicate the minimum and maximum M_* values in each bin. The y-axis values and the error bars indicate the median and the 16th–84th percentiles of the posterior distribution of $\mathrm{log}({Z}_{* }/{Z}_{\odot })$. There is a tight correlation between these quantities, though the $\mathrm{log}({Z}_{* }/{Z}_{\odot })$ values in the lowest mass bins are relatively insecure. A linear fit to the results in log−log space, obtained from binning the sample into 0.2 dex width bins, yields

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}({Z}_{* }/{Z}_{\odot }) & = & -(0.81\pm 0.01)\\ & & +(0.32\pm 0.03)\mathrm{log}[{M}_{* }/({10}^{10}{M}_{\odot })],\end{array}\end{eqnarray} \tag{ 4 }$$

where the errors in the parentheses denote the nominal 1σ errors. This is shown as the blue line in Figure 10.

Zoom In Zoom Out Reset image size

This result is generally in good agreement with that of Cullen et al. (2019), who used stacks of in total 681 star-forming galaxies at z = 2.5–5 ($\left\langle z\right\rangle \approx 3.5$). There is a difference at the low-mass end, where the Cullen et al. (2019) data show a tentative sign of flattening at $8.8\lesssim \mathrm{log}{M}_{* }/{M}_{\odot }\lesssim 9.5$, which, however, is not seen in our data (it should be noted that their lowest mass point is only an upper limit). Overall, our metallicities are slightly lower at given M_*, which appears to be counter to the expected redshift evolution since the Cullen et al. (2019) sample is generally at somewhat higher redshifts. This apparent discrepancy may be accounted for by a systematic effect from the methodology that we now explain.

The method in Cullen et al. (2019) is very similar to ours and was based on stacks of rest-frame FUV spectra and model templates with various metallicities. The authors, however, used the Starburst99 high-resolution WM-basic stellar population models for their fiducial results. They compared their fiducial estimates to the ones derived using the BPASSv2.1 models and found that the use of the BPASS models leads to systematically lower metallicities by ∼ 0.1 dex. Accounting for this systematic offset will then bring the two results into better agreement. Note that the lower-metallicity limit of the Starburst99 model is $\mathrm{log}({Z}_{* }/{Z}_{\odot })=-1.15$, which is not sufficiently low for unbiased fitting for our sample, and thus we adopted the BPASS model in this work.

Our results also appear to be lower by ∼ 0.4 dex than the average measurement of $\mathrm{log}({Z}_{* }/{Z}_{\odot })=-0.425\pm 0.159$ ²³ at M_* ∼ 10¹⁰ M_⊙ obtained by Halliday et al. (2008) using the Rix et al. (2004) 1978 Å index. This offset suggests that there possibly remain substantial systematic uncertainties in the stellar metallicity measurements for different samples and methodologies.

It should also be noted that the stellar metallicities may be underestimated because the integrated stellar emission would be biased toward stars formed in the less obscured and thus lower-metallicity environments. Such systematic biases in the stellar metallicity measurements will be explored in future papers.

We now compare our result with other work at low redshifts to see how the stellar MZR evolves through cosmic time. Zahid et al. (2017) established the stellar MZR in the range 10^8.5 ≲ M_*/M_⊙ ≲ 10¹¹ using ∼2 × 10⁵ star-forming galaxies at z < 0.25 (median $\left\langle z\right\rangle =0.08$) from SDSS. The authors employed a full spectral fitting approach using the observed composite spectra in the rest-frame optical. Their result is in agreement with the measurements obtained from spectroscopy of individual stars in nearby galaxies compiled by Kudritzki et al. (2016). Furthermore, comparing to the stellar metallicities of dwarf galaxies, measured by Kirby et al. (2013), revealed that the M_*–Z_* correlation holds smoothly down to M_* ∼ 10⁴ M_⊙.

We note that the metallicities from Kirby et al. (2013) are intended to purely reflect the iron abundance (Fe/H). We assume that other local metallicity measurements also reflect the iron abundance. Particularly, this approximation can be justified for the Zahid et al. (2017) measurements that are obtained from comparison between the observed rest optical spectra and the models assuming solar α/Fe. This is because the inferred total metallicity (Z_*) is known to be related to [Fe/H] via $[\mathrm{Fe}/{\rm{H}}]=\mathrm{log}({Z}_{* }/{Z}_{\odot })-A[\alpha /\mathrm{Fe}]$ with A ≈ 0.8–0.9 (Trager et al. 2000). Indeed, we will find [α/Fe] ≈ [O/Fe] ∼ 0 for those galaxies.

Figure 11 compares these low-redshift measurements with our results at high redshift. Our high-redshift stellar MZR is clearly offset below the local relation (at our sampled masses) with ${\rm{\Delta }}\mathrm{log}Z\approx 0.8\,\mathrm{dex}$ at a given mass. We should note that the local MZR is not yet corrected for the underestimation that arises from the measurements based on the rest-frame optical integrated light of the galaxies. We will take account of this correction, which is about +0.1 dex (see Appendix), in the subsequent sections, where we will compare high- and low-redshift MZRs more precisely.

Zoom In Zoom Out Reset image size

We employed the empirical parameterization introduced by Curti et al. (2020, Equation (2)) to express the local stellar MZR:

$$\begin{eqnarray}&&\mathrm{log}\left(\displaystyle \frac{Z}{{Z}_{\odot }}\right)={{ \mathcal Z }}_{0}-\displaystyle \frac{\gamma }{\beta }\mathrm{log}\left(1+{\left(\displaystyle \frac{{M}_{* }}{{M}_{0}}\right)}^{-\beta }\right),\end{eqnarray} \tag{ 5 }$$

where Z₀ is the asymptotic metallicity at the massive end, M₀ is a characteristic mass where the relation begins to flatten, γ is the power-law slope of the relation at M_* ≪ M₀, and β determines the width of the transition region between the two extremes. Fitting the measurements from Zahid et al. (2017) yields ${({{ \mathcal Z }}_{0},\,{M}_{0},\,\gamma ,\,\beta )}_{z\sim 0}$ = (0.049, 10.24, 0.40, 6.36). In Figure 11, this best-fit relation clearly represents the local measurements.

Comparing with the local relation, our result at z = 1.6–3.0 is found to have a slightly shallower slope, 0.32 ± 0.03 (see above Equation (4)). As noted above, the high-redshift relation is offset to lower metallicities and, perhaps for this reason, does not show evidence of a saturation at the high-mass end.

As already noted, the gas-phase and stellar "metallicities" can be more or less translated, at both high and low redshifts, into oxygen and iron abundances, respectively. Comparison of the two may therefore give some insight into the dependence and evolution of the relative abundance ratio, O/Fe, often used as a proxy of the α-enhancement.

The gas-phase metallicity reflects the instantaneous oxygen abundance at the time of observation. Measuring the gas-phase iron abundance is, however, quite challenging because of the faintness of iron emission lines, as well as the uncertainties in dust depletion factors; we therefore need to rely on the stellar metallicities. The stellar metallicities of our high-z galaxies estimated from the rest-frame FUV spectra, which are dominated by short-lived, recently formed stars, offer a more "instantaneous" measurement, similar to the gas-phase O/H. As we discuss in Appendix, we find that the FUV-weighted [Fe/H] is consistent within ∼0.02 dex with the gas-phase values. We can thus adopt the ratio of the gas-phase oxygen abundance to the stellar iron abundance of z ∼ 2.2 galaxies, [O/Fe], as a proxy of the instantaneous α-enhancement of the gas in the galaxies for the remainder of the paper.

In contrast, the iron abundance based on the rest-frame optical spectra could be lower than the instantaneous value because of the substantial contribution to the integrated light from older stars. This is the case of the stellar MZR at z ∼ 0 from Zahid et al. (2017). Therefore, the comparison between the two [Fe/H] measurements is not straightforward. As shown in Appendix, the offset with respect to the instantaneous [Fe/H] is probably around ≈ − 0.1 dex for low-redshift galaxies. In the subsequent sections, we thus shift the local stellar MZR by ${\rm{\Delta }}\mathrm{log}({Z}_{* })=+0.1\,\mathrm{dex}$ so that it better reflects the instantaneous iron metallicity. This correction, however, does not have any significant impact on our conclusions.

We take the gas-phase [O/H] MZRs at z ∼ 0 and z ∼ 2.2, respectively, from Curti et al. (2020, Equation (2)) and Sanders et al. (2020). The former uses the latest accurate calibration between the optical strong-line ratios and the metallicity determined from the direct method. The latter is purely based on the direct method O/H measure in individual 18 galaxies at z = 1.7–3.6 (median $\left\langle z\right\rangle =2.17$). An important caveat is that the sample of Sanders et al. (2020) is not representative of our zCOSMOS-deep sample, but rather a compilation from the literature of different surveys. The following analysis is thus based on an assumption that the result of Sanders et al. (2020) and ours both independently represent the same, typical star-forming galaxy population at these epochs.

Sanders et al. (2020) corrected their O/H measurements for the residuals around the epoch's main sequence, ${{\rm{\Delta }}}_{\mathrm{MS}}\mathrm{log}(\mathrm{SFR})$, of their sample galaxies to obtain more representative metallicities at given redshift. The corrections achieve ∼0.2 dex on average. These corrections are, however, questionable at some level owing to systematic uncertainties in determining the shape of the main sequence and the ${{\rm{\Delta }}}_{\mathrm{MS}}\mathrm{log}(\mathrm{SFR})$ dependence of the metallicity. In particular, the authors adopted a locally derived relation, ${\rm{\Delta }}({\rm{O}}/{\rm{H}})=-0.29\times {{\rm{\Delta }}}_{\mathrm{MS}}\mathrm{log}(\mathrm{SFR})$. However, this would be a function of redshift, and indeed a weaker dependence, or none at all, has been claimed for z ∼ 2–3 galaxies in other papers (Onodera et al. 2016; Sanders et al. 2021). We therefore adopt here their direct measurements of O/H (Equation (7) of Sanders et al. 2020) with a moderate constant correction of ${\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})=+0.1$. Although this choice is more or less arbitrary, it is indeed within the statistical error of the direct measurement.

A similar bias might also be expected in our own sample. However, we do not find a significant offset from the main sequence at this epoch given by Whitaker et al. (2014) in the M_*–SFR diagram of Figure 3. We therefore assume that the observed galaxies and their stellar metallicities are reasonably representative of galaxies of the given M_* at these epochs.

In the following analysis, we ignore the possible depletion of oxygen onto dust grains in the ISM, which would cause some underestimate of the gas-phase O/H. The depletion of O is estimated to be ∼ − 0.09 ± 0.05 dex in z ∼ 2 star-forming galaxies (Steidel et al. 2016). This is within the observed uncertainties.

In Figure 12, we compare the empirical fits to the four MZRs (stellar and gas phase; z ∼ 0 and z ∼ 2.2). All the metallicities are here normalized to the solar values. It can be seen that the evolution in [Fe/H] between z ∼ 0 and z ∼ 2.2 is larger than that seen in [O/H]. At M_* ∼ 10¹⁰ M_⊙, the [Fe/H] metallicities at z ∼ 2.2 are about 0.9 dex lower than locally (including the +0.1 dex shift in the local stellar MZR; see above), whereas the [O/H] metallicities are only about 0.4 dex lower than that at z ∼ 0. In other words, the offset between the instantaneous [O/H] and [Fe/H] metallicities evidently increases with redshift; put another way, at M_* = 10¹⁰ M_⊙, [O/Fe] ≈ 0 is indicated locally, but [O/Fe] = 0.47 ± 0.12 at z ∼ 2.2. Given the similar slopes of the Fe and O MZRs at z ∼ 2.2, the average [O/Fe] is consistent with being independent of M_*. The implied value of [O/Fe] at z ∼ 2.2 is close to the average [O/Fe] ≈ 0.4 reported for star-forming galaxies at z ∼ 2.3 (Strom et al. 2018; Cullen et al. 2021). It also approaches the predicted maximum value (∼0.6) that is predicted for metal enrichment from CCSNe (Nomoto et al. 2006; Kobayashi et al. 2020) alone, suggesting that the production of iron by SNe Ia had not progressed far at these redshifts.

Zoom In Zoom Out Reset image size

We can in principle eliminate M_* from the gas and stellar MZRs (taken from Figure 12) to yield the relations between [O/Fe] and either [O/H] or [Fe/H]. This is done in Figure 13, which shows the [O/Fe] as a function of [Fe/H] (left panel) and [O/H] (right panel). The mass ranges of the z ∼ 0 and z ∼ 2.2 relations are limited to $\mathrm{log}({M}_{* }/{M}_{\odot })=8.5\mbox{--}11.0$ and 8.8–10.7, respectively.

Zoom In Zoom Out Reset image size

Of course, this procedure of eliminating M_* from the MZR fits may not produce the same [O/Fe] versus metallicity relations as would be obtained by considering individual galaxies, because it does not consider the scatter in the observed quantities at a given M_*. In other words, this procedure is tantamount to assuming that galactic mass is the primary driver of the variations in the [Fe/H] and [O/H] within the observed sample. In the current study, we have no alternative to this procedure, because we relied on stacked spectra (in stellar mass bins) for our metallicity measurements.

The results show that these relations have negative slopes in both panels and at both redshifts. Interestingly, the low-redshift [O/Fe]–[Fe/H] (and [O/Fe]–[O/H]) relations in Figure 13 appear to broadly line up with the overall evolutionary vectors from z ∼ 2.2 to z ∼ 0. The high-redshift [O/Fe]–[Fe/H] (and [O/Fe]–[O/H]) relation is also consistent with having the same slope, but this slope is quite uncertain.

Uncertainties in the slopes of the input MZR can propagate to have a large effect on the resultant slopes of the inferred [O/Fe] versus metallicity relations. We constructed confidence intervals on these slopes by considering the range of possible linear relations that are obtained using random combinations of fits (within the uncertainties) to the input M_*–[O/H] and M_*–[Fe/H] relations. These are shown in Figure 13. Flat or even positive correlations between [O/Fe] and overall metallicities are evidently allowed by our available high-redshift data and analysis methods. We do not try to evaluate the range of possible slopes at low redshift because, while the statistical errors in the fits are very small, it is difficult for us to assess any systematic uncertainties in those taken from the independent studies.

Given the large uncertainties in the slopes of these relations, we focus instead on the [O/Fe] and [Fe/H] (or [O/H]) values at a single fiducial mass of M_* = 10¹⁰ M_⊙ (close to the middle of our high-z mass range) as being representative of each redshift. These representative values for M_* = 10¹⁰ M_⊙ at high and low redshifts are shown by the two blue and black squares in each of the panels of Figure 13.

In what follows, we adopt the gas-regulator model of Lilly et al. (2013) in which the SFR is instantaneously regulated by the mass of gas (M_gas) present in some reservoir, via the star formation efficiency (SFE = SFR/M_gas), and with a wind-driven mass loss that scales with the SFR via a "mass-loading" factor η. Mass conservation then straightforwardly gives (see Equation (9) of Lilly et al. 2013)

$$\begin{eqnarray}&&{\rm{\Phi }}=(1-r+\eta )\cdot \mathrm{SFR}+{\dot{M}}_{\mathrm{gas}},\end{eqnarray} \tag{ 6 }$$

where Φ is the mass inflow rate and r is the recycling factor (or called the return fraction), i.e., the fraction of mass that is formed into stars and then at later times returned to the ISM.

Obviously, once the two parameters, SFE and mass loading η, are specified (possibly as a function of mass and/or redshift), the gas accretion history (and gas content history M_gas(t)) of a given system follows completely from its SFH, since ${\dot{M}}_{\mathrm{gas}}$ will also be given by the change in SFR(t). Note that it is the changing gas reservoir that distinguishes this "gas-regulator" model from the "bathtub" models of Bouché et al. (2010) and Davé et al. (2012), in which ${\dot{M}}_{\mathrm{gas}}$ is set to be zero.

This means that, within the context of the gas-regulator model, the instantaneous metallicity of the gas reservoir will also be completely determined by the SFH once the (possibly mass- and/or epoch-dependent) SFE and η are specified, along with an assumption of the metallicity of the inflowing gas; we will here assume for simplicity that this is zero.

In the following, we consider the instantaneous gas-phase metallicity as the ratio of metal mass in the gas phase and gas mass at given time. We do not consider metal depletion onto dust grains, and thus the gas-phase metal mass represents all metals except those locked in surviving stars and those metals gone in the wind.

We can assume that the oxygen is produced only by CCSNe that occur "promptly" (with zero time delay) after the birth of the progenitor stars. The change of oxygen mass M_O in the gas phase is therefore expressed as

$$\begin{eqnarray}&&\displaystyle \frac{{{dM}}_{{\rm{O}}}}{{dt}}={y}_{{\rm{O}}}^{\mathrm{CC}}\cdot \mathrm{SFR}-{Z}_{{\rm{O}}}(1-r+\eta )\cdot \mathrm{SFR}+{\rm{\Phi }}{Z}_{{\rm{O}}}^{\inf },\end{eqnarray} \tag{ 7 }$$

where ${y}_{{\rm{O}}}^{\mathrm{CC}}$ is the IMF-weighted oxygen yield defined as the oxygen mass synthesized and then returned into the ISM per unit mass formed, ²⁴ Z_O( = M_O/M_gas) is the oxygen abundance, and ${Z}_{{\rm{O}}}^{\inf }$ is the oxygen abundance of the infalling gas.

We use Equation (7) to numerically track the chemical evolution. However, it is also useful to formalize the metallicity too. The change in Z_O is then obtained by eliminating Φ using Equation (6) as

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{{dZ}}_{{\rm{O}}}}{{dt}} & = & {y}_{{\rm{O}}}^{\mathrm{CC}}-({Z}_{{\rm{O}}}-{Z}_{{\rm{O}}}^{\inf })(1-r+\eta )\mathrm{SFE}\\ & & -\,({Z}_{{\rm{O}}}-{Z}_{{\rm{O}}}^{\inf })\displaystyle \frac{{\dot{M}}_{\mathrm{gas}}}{{M}_{\mathrm{gas}}}.\end{array}\end{eqnarray} \tag{ 8 }$$

The metallicity in an equilibrium condition is thus derived by setting dZ_O/dt to zero, i.e.,

$$\begin{eqnarray}&&\begin{array}{l}{Z}_{{\rm{O}}}^{\mathrm{eq}}={Z}_{{\rm{O}}}^{\inf }\\ \ +\,\displaystyle \frac{{y}_{{\rm{O}}}^{\mathrm{CC}}}{1-r+\eta +{\mathrm{SFE}}^{-1}\left((1-r)\mathrm{sSFR}+{\dot{M}}_{\mathrm{gas}}/{M}_{\mathrm{gas}}\right)}.\end{array}\end{eqnarray} \tag{ 9 }$$

Lilly et al. (2013) showed that the timescale for driving Z_O toward ${Z}_{{\rm{O}}}^{\mathrm{eq}}$ is shorter than the timescale on which the equilibrium conditions are varying. Assuming equilibrium is therefore a good approximation.

For iron, we must consider the substantial amount of iron that is produced by SNe Ia, which occur with some considerable delay after the birth of their stellar progenitors. In this analysis, we adopt the expression for the distribution of the SN Ia delay time, f_Ia(t), that was formalized by Greggio (2005) for a single stellar population (see also Greggio et al. 2008; Greggio 2010). We consider contributions from the single-degenerate and double-degenerate channels. Note that f_Ia(t) is normalized so that the time integration equals 1. Figure 14 shows the adopted delay time distribution (DTD), compared with a simple exponential parameterization used in Andrews et al. (2017) and Weinberg et al. (2017). The adopted one is more sensitive to SFR at earlier times in the past.

Zoom In Zoom Out Reset image size

Using f_Ia, the change of the gas-phase iron mass M_Fe is written as

$$\begin{eqnarray}\begin{array}{rcl}\frac{{{dM}}_{\mathrm{Fe}}}{{dt}} & = & {y}_{\mathrm{Fe}}^{\mathrm{CC}}\cdot \mathrm{SFR}-{Z}_{\mathrm{Fe}}(1-r+\eta )\cdot \mathrm{SFR}+{\rm{\Phi }}{Z}_{\mathrm{Fe}}^{\inf }\\ & & +{y}_{\mathrm{Fe}}^{\mathrm{Ia}}{\int }_{0}^{t}\mathrm{SFR}({t}^{{\prime} }){f}_{\mathrm{Ia}}(t-{t}^{{\prime} }){{dt}}^{{\prime} },\end{array}\end{eqnarray} \tag{ 10 }$$

where ${y}_{\mathrm{Fe}}^{\mathrm{CC}}$ is the CCSN iron yield and ${Z}_{\mathrm{Fe}}^{\inf }$ is the iron abundance of the infalling gas. The last term denotes the contribution from SNe Ia, where ${y}_{\mathrm{Fe}}^{\mathrm{Ia}}$ is the time-integrated SN Ia yield of iron for unit mass formed. ²⁵

The steady-state iron metallicity is then written as

$$\begin{eqnarray}\begin{array}{rcl}{Z}_{\mathrm{Fe}}^{\mathrm{eq}} & = & {Z}_{\mathrm{Fe}}^{\inf }\\ & & +\,\displaystyle \frac{{y}_{\mathrm{Fe}}^{\mathrm{CC}}+({y}_{\mathrm{Fe}}^{\mathrm{Ia}}/\mathrm{SFR}){\displaystyle \int }_{0}^{t}\mathrm{SFR}({t}^{{\prime} }){f}_{\mathrm{Ia}}(t-{t}^{{\prime} }){{dt}}^{{\prime} }}{1-r+\eta +{\mathrm{SFE}}^{-1}\left((1-r)\mathrm{sSFR}+{\dot{M}}_{\mathrm{gas}}/{M}_{\mathrm{gas}}\right)}.\end{array}\end{eqnarray} \tag{ 11 }$$

Assuming ${Z}_{{\rm{O}}}^{\inf }={Z}_{\mathrm{Fe}}^{\inf }=0$, Equations (9) and (11) yield

$$\begin{eqnarray}&&\displaystyle \frac{{Z}_{{\rm{O}}}^{\mathrm{eq}}}{{Z}_{\mathrm{Fe}}^{\mathrm{eq}}}=\displaystyle \frac{{y}_{{\rm{O}}}^{\mathrm{CC}}}{{y}_{\mathrm{Fe}}^{\mathrm{CC}}+({y}_{\mathrm{Fe}}^{\mathrm{Ia}}/\mathrm{SFR}){\displaystyle \int }_{0}^{t}\mathrm{SFR}({t}^{{\prime} }){f}_{\mathrm{Ia}}(t-{t}^{{\prime} }){{dt}}^{{\prime} }}\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\approx \,\displaystyle \frac{{y}_{{\rm{O}}}^{\mathrm{CC}}\mathrm{SFR}}{{y}_{\mathrm{Fe}}^{\mathrm{Ia}}{\displaystyle \int }_{0}^{t}\mathrm{SFR}({t}^{{\prime} }){f}_{\mathrm{Ia}}(t-{t}^{{\prime} }){{dt}}^{{\prime} }},\end{eqnarray} \tag{ 13 }$$

where the approximation holds when the SNe Ia dominate the iron production. This indicates that the [O/Fe] is approximately proportional to the number ratio of the CCSNe and SNe Ia at any time. In other words, as the denominator is some kind of average SFR in the past, the [O/Fe] would be tightly correlated with the sSFR. We will see this in Section 5.3.

The parameters in these equations are not very well constrained from observations. We therefore basically follow the "fiducial" choice of Andrews et al. (2017) and Weinberg et al. (2017): r = 0.4 and ${y}_{{\rm{O}}}^{\mathrm{CC}}=0.017$ (see also Vincenzo et al. 2016 for the IMF-weighted yield and return mass fraction). As discussed below (Section 5.2.2), we adjusted the values of the CCSN and SN Ia iron yields, taking ${y}_{\mathrm{Fe}}^{\mathrm{CC}}=0.00081$ (instead of 0.0012) and ${y}_{\mathrm{Fe}}^{\mathrm{Ia}}=0.0022$ (instead of 0.0017), together with the parameters determining the SFE and η. As noted above, for simplicity we set both ${Z}_{{\rm{O}}}^{\inf }$ and ${Z}_{\mathrm{Fe}}^{\inf }$ to zero.

5.2.1. Choice of Star Formation Histories

Again for simplicity, we construct a representative set of SFHs by integrating the evolving main sequence of star-forming galaxies across cosmic time. We adopted the local M_*–SFR relation (z ≈ 0.08) derived by Renzini & Peng (2015) and the redshift evolution at fixed M_* as follows:

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{logSFR}\,({M}_{\odot }\,{\mathrm{yr}}^{-1}) & = & 0.76\times \mathrm{log}{M}_{* }({M}_{\odot })-7.55\\ & & +\,3.0\times \mathrm{log}(1+z)\end{array}\end{eqnarray} \tag{ 14 }$$

for z ≤ 2.4 and

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{logSFR}\,({M}_{\odot }\,{\mathrm{yr}}^{-1}) & = & \mathrm{logSFR}(z=2.4)\\ & & +\,1.2\times \mathrm{log}\left(\displaystyle \frac{1+z}{1+2.4}\right)\end{array}\end{eqnarray} \tag{ 15 }$$

for z > 2.4 (see, e.g., Lehnert et al. 2015; Tasca et al. 2015).

The inferred evolutionary tracks in the M_*–SFR plane and SFHs are shown in Figure 15. In the former, we mark the values of our representative galaxies at z = 0.08 and z = 2.22. These simulated galaxies at z = 2.22 are in broad agreement with our sample galaxies, as shown by gray dots.

Zoom In Zoom Out Reset image size

5.2.2. Parameters of the Regulator Systems and Iron Yields

In order to derive the time-dependent gas content of each model galaxy from their individual SFHs and thereby compute the corresponding chemical evolution, we need to define the two parameters of the gas-regulator system, the mass-loading factor η and SFE. We assume that these values scale with the instantaneous M_* of the system, but that they are not redshift dependent, i.e.,

$$\begin{eqnarray}&&\eta =\max \left\{{\eta }_{10}\times {({M}_{* }/({10}^{10}\,{M}_{\odot }))}^{a},{\eta }_{\min }\right\},\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\mathrm{SFE}={\mathrm{SFE}}_{10}\times {({M}_{* }/({10}^{10}\,{M}_{\odot }))}^{b}.\end{eqnarray} \tag{ 17 }$$

Here we consider the minimum value for η at high masses (assuming a negative a) for better representation of the saturation feature of the local MZR at the high-mass end. Although this is an arbitrary treatment, the bending in the average η toward high masses is seen owing to effects of AGNs in simulations (Nelson et al. 2019).

Now Equations (7) and (10) can be used to compute the evolution with time of [O/H] and [Fe/H] for any arbitrary SFR(t), i.e., for each of the representative SFH identified above. However, using the fiducial iron yields given by Andrews et al. (2017) and Weinberg et al. (2017) gives an [Fe/H] MZR that is slightly higher at z ∼ 2.2 but lower locally than those observed at each redshift shown in Figure 12. We also notice that the value of ${y}_{\mathrm{Fe}}^{\mathrm{CC}}=0.0012$ adopted in the above papers gives [O/Fe] (=0.5 with ${y}_{{\rm{O}}}^{\mathrm{CC}}=0.017$) from pure CCSNe that is smaller than the values seen in the literature (∼0.6; e.g., Nomoto et al. 2006). We therefore allowed the CCSN and SN Ia iron yields to vary in order to better reproduce the observed MZRs.

We used an MCMC algorithm to determine values of the six parameters to be $({\eta }_{10},a,{\eta }_{\min },{\mathrm{SFE}}_{10}({\mathrm{Gyr}}^{-1}),b,{y}_{\mathrm{Fe}}^{\mathrm{CC}},{y}_{\mathrm{Fe}}^{\mathrm{Ia}})$ =(2.58, −0.267, 2.34, 0.310, 0.043, 0.00083, 0.00215) that, using our representative SFHs, well reproduce the M_*–[O/H] and M_*–[Fe/H] relations at z ∼ 0 and z ∼ 2.2 shown in Figure 11.

We note that these values are not far from those suggested from observations and/or simulations. The mass-loading factor and its scaling are similar to what is found in Muratov et al. (2015, see their Equation (8)). The typical (total gas) depletion timescales (t_dep = 1/SFE) of several Gyr have been observed (Bigiel et al. 2008) and reproduced in simulations (Semenov et al. 2017). From Equations (14) and (17), we obtain ${M}_{\mathrm{gas}}\propto {M}_{* }^{0.72}$, which is in broad agreement with the scaling relation between the molecular gas mass and stellar mass ($\sim {M}_{* }^{0.59};$ Tacconi et al. 2020). The corrections of the iron yields are also small ( ∼ 30%) with respect to the fiducial values.

5.2.3. Calculated Chemical Evolutionary Tracks

In Figure 16, we show the chemical evolution tracks for the input SFHs in the iron and oxygen MZR diagrams. The individual lines correspond to each of the SFHs that were shown in Figure 15. As there, the positions of the galaxies at z = 2.22 and at z = 0.08 are marked (again, black diamonds and red circles, respectively). At z = 2.22 we limit those to having M_* higher than ≈ 10^8.8 M_⊙, which is the lower mass limit of our high-z sample.

Zoom In Zoom Out Reset image size

This figure illustrates how the models well reproduce the observed MZRs in both oxygen and iron, including their slopes, at both z ∼ 0 and z ∼ 2.2. This was, of course, to be expected since both the SFE(M_*) and η(M_*) functions of the regulator model and the iron yields have been adjusted so as to match these overall MZRs at both z ∼ 2 and z ∼ 0, although we stress that the adopted parameters are, in our view, completely reasonable.

Rather, our interest in constructing these models is in examining and understanding the expected relations between the overall metallicity and the α-enhancement at these two redshifts, especially in comparison with the evolutionary offset between them. We address this in the following sections of the paper.

In Figure 17, we show the α-enhancement versus overall metallicity of these same models, i.e., [O/Fe] as a function of either [Fe/H] (left panel) or [O/H] (right panel). As in Figure 16, the model galaxies with M_* ≳ 10⁹ M_⊙ are marked at z = 2.2 and at z = 0 in order to compare with the observations (as in Figure 13). Both the evolutionary tracks in this diagram and the variations within the population at fixed epoch (locus of black and red points) may be read from this diagram.

Zoom In Zoom Out Reset image size

It can be seen that the modeled evolutionary tracks are in good agreement with the overall change in both metallicity and [O/Fe] between z ∼ 2.2 and the present epoch. This is not surprising because the model parameters were tuned (within reasonable ranges) to fit the oxygen and iron MZRs, and thus implicitly the [O/Fe], at both redshifts.

The models present variations across mass that show a negative slope between these quantities. The observed slopes are similar to what are observed at these redshifts. However, as mentioned in Section 4.3, large uncertainties in the slopes do not enable us to make any robust statement. We thus focus on the representative [O/Fe] and [Fe/H] (or [O/H]) values at our single fiducial mass of M_* = 10¹⁰ M_⊙ at each redshift (squares in Figure 17).

It is noticeable that, at least for the range of representative SFHs considered, the model tracks all follow a relatively narrow path in the [O/Fe]–[Fe/H] plane. The tightness of this path is enhanced because, at a given epoch, the loci of (model) galaxies of different masses are close to being parallel to the individual evolutionary tracks.

We now turn to investigating whether the observed offset between z ∼ 0 and z ∼ 2.2 (i.e., the "evolution vector") in this diagram can be explained by the change with redshift in some other parameter, and whether that parameter, if present, might also be responsible for producing the variation within the population at a single epoch. While the [O/Fe] (or generally α/Fe) is often loosely considered as being indicative of the "age" of a stellar system, because of the time delay in producing much of the iron via SNe Ia, it should really reflect, especially in a continuous "flow-through" scenario, the sSFR. The sSFR reflects the ratio of the current SFR to some average SFR in the past. This ratio will thus determine the relative number of CCSNe and SNe Ia at any point in time, which will then determine the "instantaneous" α-enhancement of the gas (see Equation (13)).

To explore this, we show in Figure 18 the same evolutionary tracks (as in Figure 17) but now color-coded by the instantaneous sSFR (top panels) and M_* (bottom panels) of the model galaxy. Now the different dependences of the evolutionary tracks on these quantities are quite obvious: it is clear that the [O/Fe] is very tightly correlated with the sSFR across the whole range of these evolutionary tracks and is largely unaffected by either the metallicity or the mass of the system. This is indicated by the striking horizontal banding of the sSFR-coded colors in the top panels. In contrast, the inclined M_*-coded color banding in the bottom panels shows that the metallicity (whether [O/H] or [Fe/H]) depends on both the stellar mass and [O/Fe]. In other words, using the tight relation between [O/Fe] and sSFR, the metallicity depends on the mass and sSFR together.

Zoom In Zoom Out Reset image size

Figure 19 shows the same data in a different way, demonstrating that, at least in the gas-regulated models and representative SFHs considered in this paper, there is remarkably little scatter between [O/Fe] and the instantaneous sSFR of the model galaxies. The [O/Fe] scales as ≈ sSFR^0.34 at $\mathrm{log}(\mathrm{sSFR}\,[{\mathrm{Gyr}}^{-1}])\lt 0.5$.

Zoom In Zoom Out Reset image size

These results support the idea that the tight relation between [O/Fe] and sSFR is quite fundamental, being almost independent of the epoch, the shape of the SFHs, and thus, in our modeling, the present-day mass and thus the values of SFE and η.

Lastly, we compare these results to the concept of the so-called "fundamental metallicity relation" (FMR; e.g., Mannucci et al. 2010; Lara-López et al. 2010; see also Andrews & Martini 2013; Lara-López et al. 2013). The FMR established that the SFR appears to be a second parameter in the gas-phase [O/H] MZR at low redshift, in the sense that higher-SFR galaxies at a given mass have lower oxygen metallicities. Further, it was then shown that the observed evolutionary change in [O/H] to high redshift (at a given mass) was the same as that obtained by simply extrapolating the trend with SFR established at z ∼ 0 to the much higher SFR seen at high redshift.

An important insight into the FMR came from the introduction of the gas-regulated model of galaxies (Lilly et al. 2013). Assuming a quasi-equilibrium state, the (gas-phase) metallicity of the galaxy is set "instantaneously" by two considerations. The first is the specific rate at which the system is being fed by gas (dM_gas/dt/M_gas), which itself can be inferred from the sSFR. The second is the values of the two regulator parameters: the SFE and mass loading η. The latter parameters determine the gas content of the system that is necessary to achieve the required sSFR (see Lilly et al. 2013 for details and discussion). The evolving sSFR of main-sequence galaxies can be viewed in this framework as the primary driver of the observed redshift evolution of the gas metallicity, but the positive MZR at a given epoch is the consequence of the change in sSFR along the main sequence, together with the mass dependences of the SFE and mass loading η.

It is important to clarify that the dependence of [O/H] on the mass-dependent SFE and η (see Lilly et al. 2013) contrasts with the fact that [O/Fe] is determined almost entirely by the sSFR alone, acting as a proxy of the number ratio of CCSNe and SNe Ia at a given time. The variations in SFE and η (with mass) cause the variations in O/H at fixed sSFR (thus at fixed O/Fe) and thus the range in the evolutionary tracks in the [O/Fe]–metallicity diagrams (Figure 17). Note that, in the context of the gas-regulator model, if SFE and η are constant, the gas-phase metallicity is also determined by sSFR alone. The [O/Fe]–metallicity tracks in Figure 17 will thus be completely independent of SFHs, having no scatter, while, at a single epoch, the values of [O/Fe] and metallicities vary along this single path if the sSFR changes with mass along the main sequence, or a population has a scatter in SFR at fixed mass as expected in the real universe.

Figure 20 shows data for solar neighbor FGK stars in the Milky Way (MW) adapted from Ramírez et al. (2013) and Amarsi et al. (2019). These stars are separated into the thin-disk, thick-disk, and halo populations. It is known that stars in the so-called "thick disk" are typically older, have lower metallicity, and are α-enhanced (i.e., have higher [O/Fe]) compared to stars in the "thin disk." It is indeed clear that thin- and thick-disk stars occupy on average different regions in the [O/Fe]–metallicity plane. Halo stars present large scatters in both metallicity and α-enhancement but mostly have lower Z and higher α than the thin-disk stars.

Zoom In Zoom Out Reset image size

Galactic stars, especially the thick-disk stars, clearly lie along the locus extending between the two M_* ∼ 10¹⁰ M_⊙ squares in Figure 20. In particular, it can be seen that our best-estimate locations of M_* = 10¹⁰ M_⊙ galaxies at z ∼ 2 in the two [O/Fe]–metallicity diagrams lie at the upper tail of the sequence of thick-disk stars in our own Galaxy.

Figure 20 also shows that the simple models constructed in Section 5 successfully reproduce the locus of these Galactic stars. This gives added confidence that they are reasonable.

Martig et al. (2016) found the median ages of the thick-disk population of stars to be ≈9–5 Gyr, decreasing from the inner to the outer disk. Kilic et al. (2017) also estimated the average age of the thick-disk stars to be ≈ 9 Gyr. It is thus quite plausible that those Galactic thick-disk stars located around our high-z measurement (i.e., the relatively metal-poor α-enhanced stars) formed around z ∼ 2–3, i.e., some 10 Gyr ago and some 3.5 Gyr after the big bang. We have already discussed how, in terms of the simple flow-through scenario of chemical evolution discussed in Section 5, it is natural to get these high α-enhancements ∼10 Gyr ago. Note that, in our modeling, present-day galaxies with the same mass as the Galaxy (6.4 × 10¹⁰ M_⊙; McMillan 2011) had a mass around this fiducial value of M_* = 10¹⁰ M_⊙ at z ∼ 2.

The evident agreement of [O/Fe] and [Fe/H] (or of [O/Fe] and [O/H]) between the old stars in the Galactic thick disk and the high-redshift galaxies with masses comparable to those expected for MW progenitors provides, in our opinion, a beautifully direct link between the results of Galactic archeology and observations of galaxies in the high-z universe that has long been assumed but rarely established directly (see Cullen et al. 2021 for a similar conclusion).

In the subsection above, we found a general agreement between the chemical abundances seen in high-z (z ∼ 2) galaxies and the MW thick-disk population, plus the estimated ages of the latter, and it supports the idea that, when observing z ∼ 2 galaxies, we are witnessing the formation of stars that will constitute the thick-disk population in MW-like galaxies. However, as we will see below, some caveats and questions are posed when a detailed comparison is made with bulge stars in the Galaxy.

Recently, Queiroz et al. (2020) have presented a detailed study of chemical abundances from APOGEE data (Abolfathi et al. 2018) combined with stellar distances from Gaia (Gaia Collaboration et al. 2018). They constructed detailed distributions in [α/Fe] versus [Fe/H] plots at several radial distances from the center of the Galaxy and at three vertical distances from the Galactic plane (see their Figure 6). Note that they measure [α/Fe] using different multiple α-elements, and the comparison to [O/Fe] probably needs some offset (order 0.1 dex; McWilliam et al. 2008; Jönsson et al. 2020). They showed that there are basically two disjoint populations in the [α/Fe]–[Fe/H] plot producing, at each spatial location, a bimodal distribution of stars (the so-called α-dichotomy). Broadly speaking, these two distinct populations comprise a metal-rich (solar or more), low-α (close to solar) component and a lower-metallicity, high-α component. Globally, ∼50% of bulge stars belong to the metal-rich component (Zoccali et al. 2017).

In the bulge region (R < 2 kpc), the supersolar metal-rich component with [α/Fe] spread around zero dominates. Moving radially outward in the plane of the disk, the metal-rich component shifts slightly to lower (i.e., subsolar) metallicities without significant change in [α/Fe], while the metal-poor, high-α component progressively vanishes. Instead, moving vertically away from the plane, it is the metal-rich component that tends to vanish, leaving the metal-poor component dominant, without much change in its abundances. This variation in the spatial manifestation of the bimodality may therefore be related to the "thin" and "thick" disk components, at least at intermediate radii. However, the dominant metal-rich component in the central bulge region challenges this simple picture and poses questions for the beautiful concordance between z ∼ 2 galaxies and the thick-disk stars shown in Section 6.1.

It is clear that we find no evidence in our zCOSMOS-deep sample at z ∼ 2.2 for star formation with abundances comparable to the metal-rich low-α component. Some studies, however, have indicated that this metal-rich component is dominated by stars of age ∼10 Gyr or older, as demonstrated by multiband Hubble Space Telescope photometry extending from the near-UV to the near-IR (Brown et al. 2009; Renzini et al. 2018). ²⁶ Based on these indications, this component seems to have formed at z ∼ 2. The obvious question is, why are such abundances not seen in our z ∼ 2 zCOSMOS-deep galaxies? If we suppose that MW progenitor(s) should indeed have been found in our sample of main-sequence star-forming galaxies at z ∼ 2.2, then there are some possible explanations for this observational discrepancy.

The first possibility is that the central metal-rich component in the MW had not in fact formed by z ∼ 2, i.e., that the ages of these stars have been overestimated. This possibility would be consistent with our observational results and also with our models, since these suggest that galaxies with SFHs that follow the evolution of the main sequence cannot have [α/Fe] ∼ 0 at z ∼ 2. The main problem with this first possibility is the observational evidence for old ages for these metal-rich low-α stars. One way out could come from assuming that the central regions of the MW are unusual and not representative of the general galaxy population seen in the high-redshift universe.

A second possibility is that metal-rich stars are indeed forming at z ∼ 2, as indicated by their ages, but that these metal-rich stars will not have contributed to the rest-frame FUV spectra on account of being heavily obscured by dust, which may be common in a supersolar-metallicity environment.

A second possibility is that metal-rich stars are indeed forming at z ∼ 2, as indicated by their ages, but that these metal-rich stars have contributed to the rest-frame FUV spectra on account of being heavily obscured by dust, which may be common in a supersolar-metallicity environment. The global dust obscuration in the FUV estimated from SED fitting is on average A_FUV ≈ 1.4 at M_* = 10¹⁰ M_⊙, increasing with M_* from ∼0 ( ∼ 10⁹ M_⊙) to ∼3 ( ∼ 10¹¹ M_⊙). This means that on average 70%–80% of the light from SFR in the galaxies at M_* = 10¹⁰ M_⊙ has been absorbed. If there is a large variation in dust obscuration across the galaxy, increasing toward regions that are more metal-rich and solar-α, the light from this SFR may have been preferentially, and possibly completely, removed from the observed FUV spectra.

We can, however, test the hypothesis of substantial heavily obscured star formation making almost no contribution to the FUV light by examining the far-infrared emission of our sample using the public super-deblended catalog (Jin et al. 2018). Only 54 sources out of the entire sample of 1336 galaxies have FIR luminosity detected at > 3σ. Even for these FIR-selected galaxies the nominal FIR-based SFRs (median ∼ 140 M_⊙ yr⁻¹; estimated in the same way as in Kashino et al. 2021) are only 50% higher than the SED-based SFR estimates ( ∼ 100 M_⊙ yr⁻¹). There is therefore no compelling evidence with the available data of substantial "hidden" SFR across the entire sample.

The third possibility is that the formation of the metal-rich component was completed well before z ∼ 2, and that such stars are present in the zCOSMOS-deep galaxies but already old enough so as not to contribute to the ultraviolet light. This may be partly supported by observations that have revealed metal-rich (≳solar) stars in quiescent galaxies at z ∼ 2 (Estrada-Carpenter et al. 2019; Morishita et al. 2019). However, such a scenario seems implausible to us for several reasons. Pushing star formation to earlier epochs will clearly tend to increase the α-enhancement, not decrease it, exacerbating the problems identified in the previous paragraph. Not least, the stars would have had to form in a still shorter time at even higher sSFR. Indeed, Onodera et al. (2015) found an enhanced average [α/Fe] ( ≈ 0.3), along with supersolar [Z/H] ( ≈ 0.24), for quenched galaxies at z ∼ 1.6.

Lastly, we should also consider a bias in the selection of our sample, namely, that galaxies containing metal-rich star formation may be not represented in our blue-selected (B_AB ≲ 25) zCOSMOS-deep sample, which excludes objects with low FUV continuum. It is naturally expected that galaxies either dominated by metal-rich star formation or having experienced starbursts are dust rich, and thus their rest-frame UV emission may well disappear below the selection limit. A highly complete study, including UV faint sources, is desired to gain a more robust conclusion.

We return to emphasize the basic difficulties (required in both the second and third possibilities discussed above) of achieving low (solar) α-enhancement at z ∼ 2, at least using our fiducial SN Ia DTD that provides the satisfactory match presented in Section 5. We showed in Section 5.3 that the SFHs that are obtained by integrating the main sequence simply do not yield gas abundances with low α-enhancements at these redshifts. This is ultimately because of the high sSFR of the main sequence at z ∼ 2. Except in extremely contrived scenarios, low α-enhancement stars must be formed at low sSFR. If a significant mass of stars is to be formed, then this star formation must be maintained over an extended period of time, i.e., over timescales of order sSFR⁻¹. This is naturally achieved at late epochs in the universe, but it is hard to see how this can happen at much earlier times. Quite independent of our own observational results, this conundrum represents a basic puzzle about the origin of the low α-enhancement metal-rich central population identified in the central kiloparsec of our Galaxy by Queiroz et al. (2020) if these stars are truly over 10 Gyr old rather than ≲ 4 Gyr old (i.e., formed after z ∼ 0.5) as indicated by our model (Figure 17).

It is clear that our models, which employ representative smooth SFHs, produce one continuous sequence of stars in the α–metallicity plane, rather than the two distinct components as observed (Queiroz et al. 2020). This suggests that the construction of the MW was a more complex phenomenon. Recently, an SFH with an extended gap ( ∼ 6 Gyr) between the early formation of the metal-poor high-α and the metal-rich low-α components has been proposed to reproduce the α-dichotomy of the MW stars (Lian et al. 2020; Queiroz et al. 2021). As star formation of the metal-poor component quenched, the ISM was enriched only in Fe by ongoing production of SNe Ia on a ∼ 1 Gyr timescale, and the metal-rich, low-α populations were produced from the low-α gas when star formation resumed after the gap. In this scenario, the MW would have been in the gap at z ∼ 2–3. But this is manifestly not the case for the majority of our sample galaxies, as they actually have high sSFRs. This would suggest that the α-dichotomy may be a characteristic of the MW rather than a common feature, though some cases of "rejuvenation" after a hiatus in star formation have been proposed (Mancini et al. 2019). We should also mention another work by Clarke et al. (2019), who proposed clumpy star formation led by disc fragmentation as the origin of high-α populations, and claimed that the α-dichotomy is common in massive galaxies as high-z galaxies commonly show clumpy morphologies.

Lastly, there is another possibility to be mentioned, namely, that the SN Ia DTD may shift to shorter delays with increasing metallicity, in particular in the supersolar regime. For example, in the double-degenerate scenario for SN Ia progenitors, the event results from the merging of two white dwarfs (WD) as they spiral in owing to gravitational wave radiation. Thus, the DTD is controlled by the distribution of the binary WD separations as they emerge from their last common envelope event. As the delay time scales as the fourth power of such separation (e.g., Greggio 2005), a ∼30% reduction in the distribution of WD separations would result in a shift by a factor of ∼5 in the DTD, e.g., moving the median delay time from, say, ∼1 Gyr to ∼200 Myr. The amount of orbital shrinkage during common envelope events is notoriously hard to predict; nevertheless, we note that stellar sizes are larger at high metallicity, and hence stars fill their Roche lobe at an earlier evolutionary phase compared to the case at lower metallicity. An effect of metallicity on the WD orbital separations and on the resulting DTD is therefore to be expected, though hard to quantitatively predict, and the few attempts in this direction remain inconclusive (Meng et al. 2011; Meng & Yang 2012; Kistler et al. 2013). The old ages, high metallicity, and low [α/Fe] of bulge stars could then be more easily reconciled if the DTD were to move to substantially shorter delays in the supersolar regime.

It is clear to us that integrated light spectroscopy of high-redshift galaxies, such as presented in this paper, can highlight the problem, but it is unlikely to enable us to understand the complexities revealed by individual star abundances across the body of the MW. It is possible that high-resolution imaging and integrated field spectroscopy with the forthcoming JWST instruments will finally give clues to the solution of this puzzle.