The Heavy Metal Survey: The Evolution of Stellar Metallicities, Abundance Ratios, and Ages of Massive Quiescent Galaxies since z ∼ 2

In the present-day Universe, nearly all massive galaxies have quiescent stellar populations. Archaeological studies have revealed that they formed the bulk of their stars rapidly at z > 2 (e.g., Gallazzi et al. 2005; Thomas et al. 2005; McDermid et al. 2015). Thus, these studies imply that quiescent galaxies should already exist when the Universe was only a fraction of its current age.

Quiescent galaxies have indeed been identified and studied out to z ∼ 5 (e.g., Franx et al. 2003; Cimatti et al. 2004, 2008; Daddi et al. 2005; Kriek et al. 2006; Kriek et al. 2009; van de Sande et al. 2013; Belli et al. 2014; Glazebrook et al. 2017; Schreiber et al. 2018; Valentino et al. 2020; Antwi-Danso et al. 2023; Carnall et al. 2023; Nanayakkara et al. 2024) and are thought to dominate the massive galaxy population out to z = 2.5 (e.g., Muzzin et al. 2013; Tomczak et al. 2014; McLeod et al. 2021). Interestingly, quiescent galaxies at higher redshifts are significantly more compact than the population at z ∼ 0, with their half-mass radii growing by a factor of 3 since z ∼ 2 (Daddi et al. 2005; Dokkum et al. 2008; Damjanov et al. 2009; van der Wel et al. 2014; Suess et al. 2019a, 2019b; Mowla et al. 2019; Miller et al. 2023). Distant quiescent galaxies, however, have central mass densities similar to their nearby analogs, and the size difference is primarily attributed to the larger envelopes of the lower-redshift galaxies (e.g., Bezanson et al. 2009; Hopkins et al. 2009; van Dokkum et al. 2010; van de Sande et al. 2013; Barro et al. 2017; Ji & Giavalisco 2022). Furthermore, the color gradients of quiescent galaxies become stronger with time, with the galaxy outskirts becoming progressively bluer than the centers (Suess et al. 2019a, 2019b, 2020; Miller et al. 2023).

One way to understand the structural evolution of massive quiescent galaxies is through the "inside-out" growth scenario, in which compact galaxies still grow via gas-poor minor mergers after they have become quiescent (e.g., Bezanson et al. 2009; Naab et al. 2009; van Dokkum et al. 2010; Oser et al. 2012). An alternative explanation is suggested by the fact that the stellar mass function shows a tenfold increase in the number of massive quiescent galaxies since z ∼ 2 (McLeod et al. 2021). If galaxies that quench at lower redshift have larger sizes and stronger color gradients, the observed evolution of average properties could instead be explained by population growth (i.e., progenitor bias; e.g., Khochfar & Silk 2009; van Dokkum et al. 2010; Carollo et al. 2013; Poggianti et al. 2013). The relative importance of minor mergers versus progenitor bias remains a central question in massive galaxy evolution.

The evolution in elemental abundances and stellar population properties over cosmic time provides a unique insight into the assembly and star formation histories of massive quiescent galaxies (e.g., Matteucci 1994; Trager et al. 2000; Choi et al. 2014; Conroy et al. 2014; Peng et al. 2015; Spitoni et al. 2017; Maiolino & Mannucci 2019; Trussler et al. 2020; Beverage et al. 2021). Unfortunately, measuring elemental abundances beyond the low-redshift Universe poses a considerable challenge. At higher redshifts, key absorption features are faint and shifted to near-infrared (NIR) wavelengths. Recently, ultradeep spectroscopic surveys have started pushing these measurements beyond the low-redshift Universe and finding that the chemical properties of massive quiescent galaxies at z ∼ 0.7 closely resemble today's massive early-type galaxies (Choi et al. 2014; Leethochawalit et al. 2018; Bevacqua et al. 2023; Beverage et al. 2023). These studies show that the most massive galaxies were already in place 6 billion yr ago, while the low-mass quiescent population has continued to grow since z ∼ 0.7 (see Leethochawalit et al. 2019; Beverage et al. 2023).

While minimal evolution is found out to z ∼ 0.7, the picture at higher redshift is less clear. The few existing stellar metallicity measurements of massive quiescent galaxies at z ≳ 1.4 show systematically lower [Fe/H] compared to z ∼ 0, but there is large scatter depending on the methods used to derive the metallicities (Lonoce et al. 2015; Onodera et al. 2015; Kriek et al. 2016, 2019; Carnall et al. 2022; Saracco et al. 2023; Zhuang et al. 2023). At z > 2, the only two existing measurements show strong α-enhancement, with [Mg/Fe] ∼ 0.5–0.6 dex, indicative of extremely short star formation timescales of <100 Myr (Kriek et al. 2016; Jafariyazani et al. 2020). To further clarify the picture at z > 1 and assess the assembly of the quiescent galaxy population over the past 11 billion yr, larger galaxy samples are needed.

To that end, we have executed the Keck Heavy Metal Survey, an ultradeep spectroscopic survey with LRIS (Oke et al. 1995; Rockosi et al. 2010) and MOSFIRE (McLean et al. 2010, 2012) on the Keck I Telescope (see Kriek et al. 2024). This survey collected high signal-to-noise ratio (S/N) spectra of 21 massive quiescent galaxies at 1.4 ≲ z ≲ 2.2, covering multiple Balmer and metal absorption lines.

The Heavy Metal Survey has more than tripled the number of individual galaxies at z ≳ 1.4 with measurements of stellar metallicities, ages, and elemental abundance ratios. Preliminary results of five Heavy Metal galaxies were presented in Kriek et al. (2019) and are updated in the current work. The survey design is presented in the companion survey paper (Kriek et al. 2024). In Section 2, we describe our full-spectrum modeling technique. In Section 3, we present the stellar population results. In Section 4, we assess how the metal content of quiescent galaxies evolves between z ∼ 2.2 and z ∼ 0 by comparing our results with the LEGA-C survey (z ∼ 0.7) and the Sloan Digital Sky Survey (SDSS; z ∼ 0). In Section 5, we consider the implications of our results in the buildup of the massive quiescent galaxy population. In Section 6, we summarize our results. Throughout this work, we assume a flat ΛCDM cosmology with Ω_m = 0.3 and H₀ = 70 km s⁻¹ Mpc⁻¹.

2.1. The Heavy Metal Survey

2.2. Stellar Ages and Elemental Abundances

In the fitting presented here, we assume a model with a single burst of star formation, approximated as a simple stellar population (SSP), and a fixed Kroupa (2001) initial mass function. We note that alf also allows for a young second SSP component to account for late-time star formation. In the end, the full-spectrum model has 26 free parameters: velocity offset, velocity dispersion, SSP age, isochrone metallicity, the abundances of 19 individual elements (Fe, O, C, N, Na, Mg, Si, K, Ca, Ti, V, Cr, Mn, Co, Ni, Cu, Sr, Ba, and Eu), and the emission line strengths of H, [O iii], and [N ii].

The fitting is done in the rest frame, and the instrumental resolution of the data is accounted for by convolving the grid of models to match the resolution of the spectrum. We measure this instrumental resolution for each spectrum by finding the average full width at half-maximum (FWHM) of 20–25 skylines in the variance spectrum. The average instrumental dispersions are σ_inst = 62, 36, and 33 km s ⁻¹ for the LRIS-Red, MOSFIRE-J, and MOSFIRE-H spectra, respectively. We note that the instrumental dispersion is lower than the template resolution in alf. This template resolution can be an issue if the intrinsic velocity dispersion is lower than 100 km s ⁻¹. However, as we will find, the intrinsic velocity dispersions of the Heavy Metal galaxies are significantly higher than 100 km s ⁻¹. The spectral continuum is also removed from the observations by fitting a high-order polynomial to the ratio of the data to the model (where $n\equiv ({\lambda }_{{rmmax}}{\lambda }_{\min })/100\,\mathring{\rm A}$). We test the impact of the polynomial on the results by assuming different orders than the default (n = 8 and n = 6). We find that the polynomial order does not make a meaningful impact on the results, to well within the 1σ uncertainties. Furthermore, Conroy et al. (2018) test the impact of the normalizing polynomial on the fitting results and find that if fluxing issues are approximately linear on <50 Å scales, the default polynomial assumption is sufficient. Once normalized, the spectrum is fit using a Markov Chain Monte Carlo (MCMC; Foreman-Mackey et al. 2013) with 1024 walkers and a burn-in of 20,000 steps. The priors of the MCMC fit are uniform, and the walkers are initialized around the solar abundance pattern.

After fitting all 19 galaxies, we visually inspect the best-fit models and the corresponding posterior distribution functions (PDFs) and remove any galaxy with poorly constrained parameters. We removed three galaxies because their PDFs do not follow normal distributions and/or they have ill-defined peaks. Two of these galaxies lack spectral coverage of key absorption features (D4000 break), leading to unreliable abundance measurements, and the third galaxy has strong [O ii] and hydrogen emission. We also remove two galaxies with best-fit stellar ages <1 Gyr, as the stellar population models are not optimized for these young ages. As a sanity check, we also inspect correlations between all fitted parameters to ensure the relevant parameters are not being driven by other poorly constrained parameters. After removing these five galaxies, we are left with elemental abundances and stellar age measurements for 10 galaxies at z ∼ 1.4 and four galaxies at z ∼ 2.1.

Five of the original 19 galaxies in this study were originally presented in Kriek et al. (2019). In this study, we have since reextracted the 1D spectra using an improved optimal weighting procedure. We also assume a single SSP when fitting, whereas Kriek et al. (2019) assume a two-burst model. Nonetheless, the results presented here are consistent to within 1σ of what is found in Kriek et al. (2019).

Example spectra with corresponding best-fit models are shown in Figure 1, while the remaining spectra and fits can be found in Appendix A. Corner plots corresponding to the galaxies in Figure 1 are included in Appendix B. The best-fit parameters are listed in Table 1 for all 19 galaxies. Objects with poorly constrained abundances are noted and are removed from the following analysis. In Appendix D, we conduct a mock recovery test to check for systematic uncertainties in the alf fitting as a function of S/N. Systematic uncertainties remain <0.1 dex even at the lowest S/N (8 Å⁻¹). We refer to Kriek et al. (2024) for tabulated details of the observations and spectral energy distribution fitting results.

Zoom In Zoom Out Reset image size

Table 1. Heavy Metal Quiescent Galaxy Parameters

ID a	z	σ	Age	[Fe/H]	[Mg/Fe]
		(km s−1)	(Gyr)
23621	1.358	${199}_{-41}^{+38}$	${2.1}_{-0.7}^{+1.0}$	$-{0.41}_{-0.41}^{+0.36}$	${0.31}_{-0.37}^{+0.34}$
25407	1.359	${201}_{-28}^{+27}$	${2.5}_{-0.6}^{+0.9}$	$-{0.47}_{-0.24}^{+0.26}$	${0.64}_{-0.26}^{+0.24}$
26888	1.359	${166}_{-27}^{+33}$	${3.7}_{-1.1}^{+2.3}$	$-{0.10}_{-0.26}^{+0.20}$	${0.27}_{-0.23}^{+0.22}$
217249 b , c	1.377	${185}_{-24}^{+21}$	${0.8}_{-0.1}^{+0.2}$	$-{0.38}_{-0.26}^{+0.21}$	${0.34}_{-0.25}^{+0.25}$
214695 b	1.396	${210}_{-25}^{+24}$	${3.4}_{-0.9}^{+2.2}$	$-{0.29}_{-0.27}^{+0.24}$	${0.34}_{-0.19}^{+0.18}$
213947 b	1.397	${216}_{-13}^{+16}$	${0.9}_{-0.1}^{+0.1}$	$-{0.17}_{-0.18}^{+0.17}$	${0.38}_{-0.21}^{+0.20}$
213931 b	1.399	${371}_{-13}^{+13}$	${2.6}_{-0.4}^{+0.6}$	$-{0.12}_{-0.12}^{+0.13}$	${0.27}_{-0.07}^{+0.07}$
214340 b	1.418	${149}_{-23}^{+19}$	${1.4}_{-0.3}^{+0.4}$	$-{0.34}_{-0.20}^{+0.20}$	${0.24}_{-0.23}^{+0.23}$
23482	1.419	${289}_{-35}^{+38}$	${2.6}_{-0.7}^{+1.3}$	$-{0.09}_{-0.26}^{+0.24}$	${0.29}_{-0.22}^{+0.23}$
25702	1.419	${243}_{-20}^{+21}$	${1.8}_{-0.4}^{+0.7}$	$-{0.24}_{-0.17}^{+0.18}$	${0.32}_{-0.16}^{+0.15}$
23351	1.420	${257}_{-13}^{+14}$	${1.3}_{-0.2}^{+0.3}$	$-{0.07}_{-0.17}^{+0.15}$	${0.12}_{-0.13}^{+0.14}$
59375 c	1.550	${333}_{-43}^{+42}$	${9.5}_{-4.9}^{+3.4}$	${0.21}_{-0.30}^{+0.30}$	${0.36}_{-0.29}^{+0.28}$
60736 c	1.861	${270}_{-27}^{+30}$	${2.1}_{-0.6}^{+1.7}$	$-{0.49}_{-0.25}^{+0.23}$	${0.60}_{-0.30}^{+0.26}$
104779	2.123	${186}_{-31}^{+30}$	${1.5}_{-0.4}^{+0.9}$	$-{0.15}_{-0.28}^{+0.24}$	${0.55}_{-0.23}^{+0.23}$
56163 c	2.160	${213}_{-48}^{+54}$	${2.9}_{-1.1}^{+1.1}$	$-{0.43}_{-0.32}^{+0.32}$	${0.29}_{-0.31}^{+0.32}$
106812	2.230	${187}_{-45}^{+41}$	${1.5}_{-0.4}^{+1.1}$	$-{0.44}_{-0.28}^{+0.28}$	${0.41}_{-0.32}^{+0.36}$
108899	2.233	${306}_{-23}^{+24}$	${1.2}_{-0.1}^{+0.2}$	$-{0.38}_{-0.30}^{+0.22}$	${0.48}_{-0.29}^{+0.25}$
107590 c	2.234	${173}_{-23}^{+25}$	${0.8}_{-0.1}^{+0.1}$	$-{0.42}_{-0.22}^{+0.18}$	${0.44}_{-0.30}^{+0.26}$
103236	2.243	${197}_{-35}^{+37}$	${1.5}_{-0.4}^{+1.0}$	$-{0.23}_{-0.23}^{+0.22}$	${0.50}_{-0.23}^{+0.25}$

Notes.

^a See Kriek et al. (2024) for coordinates, magnitudes, and stellar population and structural properties. ^b Appears in Kriek et al. (2019). ^c Poorly constrained, removed from analysis.

Download table as:  ASCIITypeset image

In addition to elemental abundances and stellar ages, the full-spectrum fitting provides measurements of the stellar velocity dispersions. We find that the velocity dispersions of the Heavy Metal galaxies agree with previous measurements of the dynamical mass–stellar mass relation at similar redshifts (e.g., Belli et al. 2016), and we explore the implications of the Heavy Metal dynamical masses in the companion survey paper (Kriek et al. 2024).

Figure 2 shows the variation in [Fe/H], [Mg/Fe], and formation time (t_form) as a function of velocity dispersion σ_v for the Heavy Metal galaxies at z ∼ 1.4 (blue triangles) and z ∼ 2.1 (green stars). For comparison, we also show trends of quiescent galaxies at z ∼ 0 and z ∼ 0.7 from Beverage et al. (2023) in Figure 2. The z ∼ 0 trends are derived by refitting the stacked spectra of a sample of 5234 quiescent galaxies in SDSS from Conroy et al. (2014), while the z ∼ 0.7 trends are based on 135 quiescent galaxies in the LEGA-C survey. Initially, Beverage et al. (2021) reported a −0.2 dex offset in [Fe/H] between LEGA-C and SDSS galaxies based on the spectra from the second LEGA-C data release (DR2). However, with the subsequent release, DR3, which doubled the sample size and improved data reduction (see van der Wel et al. 2021), Beverage et al. (2023) no longer found an offset in [Fe/H] between z ∼ 0.7 and z ∼ 0. Figure 2 also includes results from other studies of quiescent galaxies at z ∼ 1.5 (Carnall et al. 2022; Zhuang et al. 2023) and z ∼ 2.1 (Kriek et al. 2016; Jafariyazani et al. 2020). ¹¹ We note that several other works have derived stellar metallicities of distant quiescent galaxies by fitting low-resolution grism spectra (e.g., Morishita et al. 2018; Estrada-Carpenter et al. 2019). However, such studies rely on the shape of the stellar continuum to constrain the metallicities, leading to high degeneracies with stellar age; therefore, we do not include these results in our comparison. To ensure unbiased comparisons, we restrict the analysis to studies that utilize alf.

Zoom In Zoom Out Reset image size

In the left panel of Figure 2, we show that the Heavy Metal galaxies have lower [Fe/H] compared to the z ∼ 0 and z ∼ 0.7 galaxies with similar velocity dispersions. Both the z ∼ 1.4 and z ∼ 2.1 galaxy populations have Fe abundances between −0.5 and −0.1 dex, with typical uncertainties of 0.2 dex. There is a hint that the z ∼ 2.1 galaxies are slightly offset to lower [Fe/H] compared to z ∼ 1.4. We fit a σ–[Fe/H] relation to the z ∼ 1.4 population using a simple linear regression, with confidence intervals determined by perturbing the [Fe/H] of each data point according to their uncertainties and refitting 1000 times. The resulting relation is shown in the left panel of Figure 2 (blue line with shaded confidence intervals). We find a positive slope of ${0.52}_{-0.23}^{+0.38}$, which is consistent with the z ∼ 0 and z ∼ 0.7 relations. However, there is a clear offset to lower [Fe/H]. We do not fit a relation to the z ∼ 2.1 population because of the small sample size.

The observed evolution in [Fe/H] between z = 1.4–2.1 and z < 1 is corroborated by studies of other quiescent galaxies at similar redshifts. At z ∼ 1.4, galaxies in the VANDELS+KMOS survey (Carnall et al. 2022) and two lensed galaxies in the AGEL survey (Zhuang et al. 2023) have [Fe/H] ∼ −0.25, in agreement with the z ∼ 1.4 Heavy Metal galaxies (see Figure 2). At z = 2.1, Kriek et al. (2016) measure [Fe/H] = −0.25 for a massive quiescent galaxy (COSMOS-11494), consistent with the values of the z ∼ 2.1 Heavy Metal galaxies (refer to Figure 2). Interestingly, the only other study at z ∼ 2.1 finds supersolar [Fe/H] for a lensed galaxy (MRG-M0138) at z = 1.98 (Jafariyazani et al. 2020). It was previously unclear whether this large variation in [Fe/H] is typical of the z ∼ 2.1 quiescent galaxy population. Now, with four Heavy Metal galaxies at z ∼ 2.1 with subsolar [Fe/H], there is additional evidence suggesting that typical massive galaxies at z ∼ 2.1 are indeed Fe-deficient and that the Fe abundance of MRG-M0138 may be anomalous. Nevertheless, the sample size is small, and the large uncertainties on the Fe abundances prevent any definitive statement about the z ∼ 2.1 population.

In the middle panel of Figure 2, we see that the Heavy Metal galaxies at z ∼ 1.4 have similar [Mg/Fe] to what is found in the SDSS (z ∼ 0) and LEGA-C (z ∼ 0.7) samples. Meanwhile, the z ∼ 2.1 Heavy Metal galaxies have [Mg/Fe] that is ≈0.2 dex higher than what is found at lower redshifts, though the typical uncertainties are large (0.2 dex) and the sample size is small. We derive a σ–[Mg/Fe] relation for the z ∼ 1.4 Heavy Metal galaxies using the same method as with [Fe/H], resulting in a negative slope of $-{0.15}_{-0.29}^{+0.17}$ (blue line with shaded confidence intervals). This slope is consistent with 0 to within 1σ, and there is no significant offset compared to the z < 1 relations. These results are supported by the other studies of quiescent galaxies at z > 1. The two studies at z ∼ 1.4 find [Mg/Fe] consistent with the LEGA-C and SDSS, while the two at z ∼ 2 find enhanced [Mg/Fe] consistent with the z ∼ 2.1 Heavy Metal galaxies.

Finally, in the right panel of Figure 2, we show galaxy σ_v versus galaxy formation time. This age is determined using the best-fit stellar population age. We find that the formation times of the Heavy Metal galaxies range from 1 Gyr (z_form = 9) to 4 Gyr (z_form = 2) after the Big Bang, with an average formation time of 2.5 Gyr (z_form = 3) and 1.5 Gyr (z_form = 4) for the z ∼ 1.4 and z ∼ 2.1 samples, respectively. We note that the galaxy formation times are calculated using an SSP-equivalent age, and more complex star formation histories would have likely resulted in slightly earlier formation times. The assumption of an SSP, however, does not bias the stellar metallicities or elemental abundance measurements (Serra & Trager 2007). In the next section, we consider the results in Figure 2 in the context of galaxy evolution.

In this section, we investigate the changes of average stellar metallicities, abundance ratios, and ages of massive quiescent galaxies as a function of redshift and discuss implications for their evolution over the past 11 billion yr. In Figure 3, we show the average [Fe/H], [Mg/H], t_form, and [Mg/Fe] for quiescent galaxies in the SDSS, LEGA-C survey, and Heavy Metal Survey as a function of redshift. To allow for a fair comparison between surveys, we only consider the SDSS and LEGA-C galaxies that fall in bins with similar σ_v as the Heavy Metal galaxies. Specifically, we remove three SDSS bins and one LEGA-C bin with σ < 170 km s⁻¹. Thus, the average velocity dispersions of the SDSS and LEGA-C galaxies used in Figure 3 are σ_v = 240 km s⁻¹, with the number of galaxies in each respective sample now being 2990 (SDSS) and 101 (LEGA-C).

Zoom In Zoom Out Reset image size

Examining the top left panel of Figure 3, we observe an increase in the average [Fe/H] from z = 2.1 to z = 0, with distinct offsets of −0.2 and −0.3 dex at z ∼ 1.4 and z ∼ 2.1, respectively, between the Heavy Metal galaxies and the lower-redshift samples. Moving to the upper middle panel, we find less evolution in [Mg/H] compared to [Fe/H]; the z = 1.4 sample is 0.1 dex lower compared to the z = 0.7 and z = 0 galaxies. The z = 2.1 galaxies have no [Mg/H] offset compared to lower-redshift galaxies, although the uncertainties remain substantial. Primarily due to the increase in [Fe/H], the average [Mg/Fe] decreases toward lower redshift, as revealed in the upper right panel. The z = 2.1 sample exhibits [Mg/Fe] elevated by 0.25 dex, while the z = 1.4 sample is only elevated by 0.08 dex. In the bottom left panel, we show that the average t_form changes from t_form = 1.5 to t_form = 4 Gyr from z = 2.1 to z = 0.

In the bottom middle panel of Figure 3, we compare the average [Mg/Fe] as a function of [Fe/H] with a simple closed-box chemical evolution model from Kriek et al. (2016; represented by a black dashed line). The model provides mass-weighted abundances for various different star formation timescales (indicated in Gyr next to the black boxes). The behavior of the model is driven by the differential enrichment timescales of Fe and Mg. While Mg is predominantly synthesized in massive stars and promptly returned to the interstellar medium through core-collapse supernovae (SNe), Fe is produced through the explosions of intermediate-mass stars (Type Ia SNe; SNe Ia) and is released on delayed timescales. Therefore, galaxies with shorter star formation timescales experience less Fe enrichment from SNe Ia prior to quenching, resulting in stellar populations with lower [Fe/H] and higher [Mg/Fe].

The simple chemical evolution model used in this paper is illustrative and is by no means meant to capture the complex feedback processes and star formation histories likely present in the galaxies at hand. More detailed chemical evolution models, including inflows and outflows, along with varied star formation histories, will help decipher the various factors shaping the elemental abundances and stellar metallicities of integrated stellar populations, especially as we move to higher redshifts (e.g., see N. Gountanis 2024, in preparation).

A visual comparison between the observations and the simple chemical evolution model highlights that the inferred star formation timescales of the z ∼ 0 and z ∼ 0.7 samples are longer than the extreme timescales observed for the Heavy Metal galaxies, which are approximately 0.15 and 0.7 Gyr for the z ∼ 2.1 and z ∼ 1.4 samples, respectively. These inferred star formation timescales derived from [Mg/Fe], along with the results for t_form, paint a picture of galaxy formation in which typical massive quiescent galaxies at higher redshift formed earlier and underwent more rapid star formation compared to their counterparts at z ∼ 0. To visually depict this scenario, we present the implied average star formation histories in the lower right panel of Figure 3. We adopt a Gaussian parameterization for the star formation histories, with the mean centered at t_form and the FWHM corresponding to the star formation timescale, which is determined by comparing the abundances to the chemical evolution model. It is important to note that this panel serves as an illustrative tool, and the choice of a Gaussian star formation history is driven by its simplicity rather than a specific physical motivation. Nevertheless, it is evident that the star formation histories of the Heavy Metal galaxies, particularly at z ∼ 2.1, are extremely short compared to the lower-redshift samples. We also note that this panel does not take into account the evolution of the stellar mass function or the spread in formation redshifts and star formation timescales within the sample. Instead, it shows the typical star formation history at each redshift interval.

Before considering how one should interpret the observed evolution, it is important to understand how selection biases influence the results. Unlike the Heavy Metal Survey, the z ∼ 0 and z ∼ 0.7 samples are mass-complete; therefore, the average stellar population properties should reflect the general galaxy population at each redshift. The Heavy Metal sample, however, has a complex selection function (see Kriek et al. 2024). The selection favors quiescent galaxies with lower mass-to-light ratios, indicative of younger ages. As a result, the ages of the Heavy Metal galaxies are likely biased young compared to the true population of quiescent galaxies at z > 1. Accounting for this bias, however, would strengthen the observed evolution in Figure 3; galaxies at z > 1 could have even earlier formation redshifts, and therefore more extreme chemical properties, than what is found in this sample.

One possible explanation for the observed evolution is the growth of the quiescent galaxy population due to the continuous quenching of star-forming galaxies between z = 2.1 and z = 0 (i.e., progenitor bias; van Dokkum & Franx 2001; Carollo et al. 2013; Lilly et al. 2013). Galaxies that quench at later epochs presumably form their stars over longer timescales, resulting in later t_form and more Fe enrichment from SNe Ia (lower [Mg/Fe]) before quenching. Hence, the addition of these galaxies to the quiescent population steadily increases the average t_form and [Fe/H], as observed. In fact, the number of quiescent galaxies with M_* > 10^10.5 M_⊙ has increased tenfold since z = 2, as evident from the evolution of the stellar mass function (e.g., McLeod et al. 2021; Conselice et al. 2022). Furthermore, most of this evolution occurs at z > 0.5 (McLeod et al. 2021), which could explain the pronounced jump in [Fe/H] and t_form from z ∼ 1.4 to z ∼ 0.7.

An alternate explanation is mergers and late-time star formation, which impact the elemental abundances of individual galaxies; when galaxies accrete satellites or form new stars, their integrated elemental abundances are altered by the newly added stars. Indeed, the importance of (dry) minor mergers in the evolution of quiescent galaxies has been posited by observations that they were more compact and had flatter color gradients at earlier times (e.g., van Dokkum et al. 2008; Naab et al. 2009; Trujillo et al. 2009; Taylor et al. 2010; van Dokkum et al. 2010; Suess et al. 2019b, 2020; Miller et al. 2023), as well as by direct pair fraction measurements (Newman et al. 2012). Furthermore, Suess et al. (2023) show that tiny companions with stellar mass ratios of 1:900 compromise some 30% of their total stellar mass budget.

The extent to which mergers, late-time star formation, and progenitor bias contribute to the evolution of the massive quiescent galaxy population remains a topic of ongoing debate (e.g., Bezanson et al. 2009; Hopkins et al. 2009; Carollo et al. 2013; Poggianti et al. 2013; van de Sande et al. 2013; Williams et al. 2017; Damjanov et al. 2019). In the following section, we examine the elemental abundances and stellar ages of individual galaxies to help differentiate between the various evolutionary scenarios; if mergers and/or late-time star formation play only a small role in the evolution of massive quiescent galaxies, then the extreme galaxies observed at z ∼ 2 should persist as rare remnants in the z ∼ 0 population.

In the previous section, we observed significant evolution in the average elemental abundances of massive quiescent galaxies from z ∼ 2.1 to z ∼ 0, suggesting that galaxies at higher redshift are typically metal-poor, have shorter star formation timescales, and formed earlier. In this section, we examine the abundance ratios of individual galaxies over cosmic time to help differentiate between various evolutionary channels.

Figure 4 shows [Fe/H], [Mg/H], and [Mg/Fe] as a function of σ_v (top row) and M_* (bottom row). The stellar masses of the Heavy Metal galaxies are derived in Kriek et al. (2024). We compare the Heavy Metal galaxies to the abundances of individual LEGA-C galaxies (Beverage et al. 2023; gray circles), and to galaxies in SDSS (black diamonds) presented by Zhuang et al. (2023). Additionally, we compare to alf-derived abundances of 41 galaxies from the MASSIVE survey, a volume-limited survey of the most massive galaxies in the local Universe (open circles; Ma et al. 2014; Gu et al. 2022). We also include abundance results from other studies of quiescent galaxies at z ∼ 1.4 and z ∼ 2.1, with the same markers as Figure 2. We note that the SDSS fiber and the LEGA-C and Heavy Metal slits cover similar physical extents within the galaxies (3–4 kpc). However, as galaxies were smaller at earlier times, the Heavy Metal and LEGA-C slits cover 1–1.2 R_e , while SDSS only covers 0.4–0.8 R_e . The MASSIVE spectra were extracted from the inner R_e /8 and therefore represent only the cores of nearby massive early-type galaxies.

Zoom In Zoom Out Reset image size

In the previous section, we proposed two theories to explain the observed abundance differences between the Heavy Metal Survey and z ∼ 0. In the first scenario, the quiescent galaxy population is continuously growing by the addition of newly quenched galaxies, driving the evolution of the observed average abundances ("progenitor bias," e.g., van Dokkum & Franx 2001; Carollo et al. 2013; Poggianti et al. 2013). In this scenario, we expect to find galaxies with properties similar to the Heavy Metal galaxies in the lower-redshift Universe, and thus the populations at higher and lower redshift should overlap. In other words, the Heavy Metal galaxies should constitute the lower-metallicity and high-[Mg/Fe] tail of the z ∼ 0 samples. Indeed, the top panels of Figure 4 show that this is mostly the case for the Heavy Metal galaxies, though there are a few exceptions; most of the Heavy Metal galaxies with high σ_v have significantly lower [Fe/H] and [Mg/H] than any galaxy found at z ∼ 0 and z ∼ 0.7. Furthermore, in the rightmost panel, we find that galaxies at z ∼ 2.1 are all [Mg/Fe]-enhanced compared to the galaxies at lower redshifts.

When considered as a function of stellar mass, the metallicity differences between the lower- and higher-redshift populations become more apparent. At all stellar masses, and particularly at the highest stellar masses, the Heavy Metal galaxies have [Fe/H] and [Mg/H] ≈0.2 dex lower than any individual galaxy in SDSS or LEGA-C. This result poses a problem for the progenitor bias scenario; if these Heavy Metal galaxies passively evolve to z ∼ 0, we would expect to find Fe- and Mg-deficient remnants in the SDSS and MASSIVE populations. Specifically, in this scenario, 10% of the SDSS and MASSIVE galaxies should be direct descendants of z > 2 quiescent galaxies given the rapid evolution of the stellar mass function (McLeod et al. 2021). Thus, the z ∼ 0 galaxies are likely not the direct descendants of the z ∼ 2 Heavy Metal galaxies. It is important to acknowledge, however, that uncertainties are substantial and sample sizes are limited. Thus, the possibility remains that we may have simply missed the rare metal-poor remnants in the z ∼ 0 samples. We also note that the Heavy Metal galaxies were selected to maximize the number of bright quiescent galaxies in one field of view. As such, the environments of these galaxies may not be representative of the broader massive quiescent galaxy population at higher redshift.

As previously mentioned, minor mergers provide an alternate way to explain the observed evolution. In this picture, the z ∼ 1.4–2.1 quiescent galaxies may grow into the cores of today's massive early-type galaxies by accreting low-mass galaxies (e.g., Bezanson et al. 2009; Hopkins et al. 2009; van de Sande et al. 2013). This scenario can successfully explain the observed increase in stellar metallicity and decrease in [Mg/Fe], but only if the accreted galaxies form over long timescales (i.e., lower [Mg/Fe]) and are more metal-rich. There are, however, two problems with this scenario. First, our observations show there is evolution within the inner 1 R_e , while minor mergers are typically accreted onto galaxy outskirts (Naab et al. 2009; Oser et al. 2010). Second, given the mass–metallicity relation, the accreted low-mass galaxies will likely be metal-poor and therefore cannot explain the increase in metallicity (Kirby et al. 2013). Major mergers with galaxies that have longer star formation timescales (and thus lower [Mg/Fe]) solve both of these problems, as more massive galaxies tend to be more metal-rich and can more easily disrupt galaxy cores. However, the observed major merger rate is likely too low to fully account for the observed evolution (Newman et al. 2012; Nipoti et al. 2012; Oser et al. 2012), and major mergers alone cannot explain the size evolution (e.g., Hopkins et al. 2009; van Dokkum et al. 2010; Hilz et al. 2013; Belli et al. 2014). Therefore, a combination of major and minor mergers may be necessary to explain the observed evolution.

Finally, late-time star formation, likely in combination with or triggered by mergers, can increase the metallicities of individual galaxies while decreasing [Mg/Fe]; if the star formation material is enriched in SN Ia products by previous epochs of star formation, then the newly formed stars will be younger, more metal-rich, and more α-enhanced. It should be noted that a burst of central star formation is expected in a wet merger scenario (e.g., Hopkins et al. 2006); thus, the two scenarios together are effective at explaining the absence of metal-poor and [Mg/Fe]-enriched galaxies by z ∼ 0 as observed in Figure 4.

Crucially, all of the above scenarios rely on the assumption that at the same stellar masses, younger galaxies have higher metallicity and lower [Mg/Fe]. For example, the progenitor bias scenario only works if the galaxies that quench at later times are more metal-rich and form over longer timescales. In Figure 5, we search for such trends in [Fe/H], [Mg/H], and [Mg/Fe] as a function of t_form for the SDSS, LEGA-C, and Heavy Metal quiescent galaxies. In each panel, we remove the first-order mass dependence by subtracting the SDSS M_*–abundance relations from all galaxies and then scaling them to the value of the M_*–abundance relations at M_* = 10¹¹ M_⊙ (e.g., Leethochawalit et al. 2019; Zhuang et al. 2023). To assess the level of correlation between quenching time (as traced by t_form) and stellar metallicities, we use the Pearson correlation coefficient. A correlation is deemed "marginal" if the Pearson correlation coefficient is 0.2 ≤ r < 0.4, "moderate" if 0.4 ≤ r < 0.6, and "strong" if r ≥ 0.6. We list the correlation coefficients along with the associated p-value in Figure 5. At fixed stellar mass, we do in fact find that younger SDSS and LEGA-C galaxies exhibit higher [Fe/H] (r = 0.53), marginally higher [Mg/H] (r = 0.26), and moderately lower [Mg/Fe] (r = −0.40), as previously found by Beverage et al. (2023) and Zhuang et al. (2023), respectively (see also Leethochawalit et al. 2018).

Zoom In Zoom Out Reset image size

The trends with [Fe/H] and [Mg/Fe] can be naturally explained by galaxies joining the quiescent galaxy population at later times (higher t_form) having more time to form stars from gas that is rich in SN Ia products (namely, Fe). However, the marginal trend with [Mg/H] is more challenging to explain, as Mg enrichment is independent of the star formation timescale and instead primarily depends on the size of the gas reservoir at the time of star formation quenching. Beverage et al. (2021) propose a model to explain the marginal trend with t_form by invoking rapid gas removal during the quenching process. In summary, if galaxies at earlier epochs, which form earlier and over shorter timescales, quench more abruptly and have larger gas fractions at the time of quenching (e.g., effective AGN outflows), then [Mg/H] is truncated to lower values. Thus, the marginal trend with [Mg/H] may suggest that galaxies at earlier times experience more effective outflows (see also the outflow model proposed by Leethochawalit et al. 2019).

Furthermore, Figure 5 shows that high-z galaxies have a slight systematic offset to lower metallicities compared to local galaxies even when the dependence on M_* and stellar age is removed. This same result was also found by Zhuang et al. (2023) but for a smaller galaxy sample. Thus, this result confirms the finding by Zhuang et al. (2023) and reinforces the need for mergers or late-time star formation to alter the chemical compositions and stellar ages of these extreme individual galaxies by z ∼ 0. Again, we note that comparisons at "fixed formation times" between galaxies at different redshift regimes should be approached with caution; the SSP-equivalent ages can be biased young. As a result, a more complicated star formation history may shift the Heavy Metal galaxies in Figure 5 toward even earlier formation times. As a result, the ages of the Heavy Metal galaxies are likely biased young compared to the true population of quiescent galaxies at z > 1. Nonetheless, correcting for such a bias would not change our interpretation.

Figures 4 and 5 together demonstrate that progenitor bias, minor and major mergers, and late-time star formation are likely all responsible, to some extent, for the observed evolution in the elemental abundances of massive quiescent galaxies over cosmic time. This analysis further shows that the existing story of massive quiescent galaxy assembly—namely, that the z > 1 quiescent galaxies are the progenitor cores of today's massive early-type galaxies—is more complicated than initially anticipated. We arrive at a similar conclusion based on the dynamical masses of the Heavy Metal galaxies in Kriek et al. (2024).

In this paper, we present the stellar metallicities, elemental abundance ratios, and stellar ages of 19 massive quiescent galaxies from the Keck Heavy Metal Survey. The galaxies are divided over two redshift intervals around z ∼ 1.4 and z ∼ 2.1. We fit the ultradeep LRIS+MOSFIRE spectra using a full-spectrum modeling code with variable elemental abundance ratios. Combined, these measurements represent the largest sample of individual quiescent galaxies at z ≳ 1.4 with stellar metallicity and elemental abundance measurements.

We find that massive quiescent galaxies at 1.4 ≲ z ≲ 2.2 have subsolar Fe abundances, ranging from [Fe/H] = −0.5 to −0.1 dex, with typical values of −0.2 and −0.3 dex at z ∼ 1.4 and z ∼ 2.1, respectively. We also find a trend between σ_v and [Fe/H], with a best-fit slope consistent with what is found at lower redshifts. We find a high [Mg/Fe] of ≈0.5 for the z ∼ 2.1 galaxies, while those at z ∼ 1.4 are lower with [Mg/Fe] ∼ 0.3. The ages range from 0.9 to 3.7 Gyr, implying formation redshifts z_form = 2–8.

We examine the evolution of the massive quiescent galaxy population over cosmic time by comparing the average Heavy Metal elemental abundances and ages to those in the LEGA-C survey (z ∼ 0.7) and SDSS (z ∼ 0). The z ∼ 1.4 and z ∼ 2.1 galaxies have [Fe/H] that is 0.15 and 0.23 dex lower than what is observed at z ≲ 0.7. Interestingly, we find no evolution in [Mg/Fe] at z ≲ 1.4, though the z ∼ 2.1 galaxies have 0.2 dex higher [Mg/Fe]. By comparing these results with a simple chemical evolution model, we find that galaxies at higher redshifts formed at progressively earlier epochs and shorter timescales, with the z ∼ 2.1 galaxies forming in only ≈150 Myr at z = 4.

The observed evolution of stellar population properties provides a unique insight into the assembly histories of massive quiescent galaxies. One possible explanation is the growth of the quiescent galaxy population. If the galaxies that join the population at later times are younger, more metal-rich, and less α-enhanced, then the average properties will evolve. However, this scenario cannot fully explain the observations because we do not find the direct descendants of the chemically extreme Heavy Metal galaxies at z ≲ 0.7. Minor mergers also fail to account for the observed evolution; due to the mass–metallicity relation, they cannot increase [Fe/H], in particular in the galaxy centers. Thus, major mergers and/or late-time star formation are required to explain the evolution of chemical properties across cosmic time. These results imply a more complicated assembly history than the minor-merger-driven inside-out growth scenario previously inferred from the evolution in the sizes, color gradients, and central densities of quiescent galaxies. In our accompanying paper, Kriek et al. (2024) arrive at a similar conclusion based on the dynamical masses of the Heavy Metal galaxies.

This study was enabled by ultradeep spectroscopic observations conducted with the most efficient ground-based spectrographs. With up to 16 hr of integration per galaxy per band, these spectra represent some of the deepest spectroscopic observations of individual galaxies at z ≳ 1.4 with coverage of key absorption features. As a result, the Heavy Metal Survey has tripled the number of individual stellar metallicity measurements at z ≳ 1.4. Nonetheless, sample sizes at these redshifts remain small, and the uncertainties on each measurement remain substantial.

In the future, JWST will provide even deeper spectroscopy of massive quiescent galaxies at z ≳ 1.4, extending the existing sample sizes, reaching to lower stellar masses, and enabling the measurements of more elements. Thus, JWST will revolutionize our understanding of the assembly of massive quiescent galaxies over comsic time.

We thank the referee for a very thoughtful and constructive report. We also thank Zhuyun Zhuang and Glenn Van de Ven for useful conversations. We acknowledge support from NSF AAG grants AST-1908748 and 1909942. A.G.B. is supported by the National Science Foundation Graduate Research Fellowship Program under grant Nos. DGE 1752814 and DGE 2146752 and by the H2H8 Association and LKBF. We gratefully acknowledge the efforts and dedication of the Keck Observatory staff for support, especially amid the challenges posed by a global pandemic. We recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has within the indigenous Hawaiian community. We are deeply fortunate to have the opportunity to conduct observations from this mountain.

Software: alf (Conroy & van Dokkum 2012; Conroy et al. 2018), FAST (Kriek et al. 2009).

In Figure 6, we present the LRIS+MOSFIRE spectra for all 19 galaxies from the Heavy Metal Survey (black) along with their uncertainties (gray). The spectra were binned such that each pixel represents 5 Å in the rest frame. We show the corresponding best-fit alf model for each galaxy. The fitting was accomplished prior to binning the spectra. Galaxies with (*) in front of the listed ID numbers are too young for reliable abundance measurements, galaxies with (^∧) had inaccurate photometric redshifts and therefore we miss key absorption features, and galaxies with (#) have unconstrained parameters. All (five) galaxies with symbols in front of their IDs were removed from the main analysis of this paper.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 7, we present corner plots from the full-spectrum alf fitting. The four corner plots correspond to the same four spectra showcased in Figure 1. Next to each corner plot, we provide information about the object: ID, redshift, and S/N (Å⁻¹). We present the four essential parameters: [Fe/H], [Mg/H], [Mg/Fe], and the age of the stellar population. The dashed lines in the plots represent the 1σ uncertainties. Though the uncertainties are large, the results follow a normal distribution. We find there to be a degeneracy between stellar age and metallicity. Target 213931 shows the strongest degeneracy, likely because the observations are a blend of two nearby systems, resulting in the superposition of stellar populations with different ages and metallicities. For the other galaxies, the age–Z degeneracy is considerably less pronounced.

Zoom In Zoom Out Reset image size

In this Appendix, we compute absorption line indices for the Heavy Metal galaxies. The Lick/IDS system, first introduced by Faber et al. (1985) and updated by Worthey et al. (1994), is the most common method for measuring elemental abundances and stellar ages of nearby quiescent galaxies. However, as we continue to push observations to higher redshifts, absorption features become fainter and more contaminated by NIR skylines. As a result, Lick indices and their derived stellar population parameters become highly uncertain. Full-spectrum modeling was, in part, developed to account for these issues. By fitting the entire stellar continuum, noise peaks can be easily downweighted, and all absorption features for a single element can be evaluated simultaneously.

In this Appendix, we check that the absorption line indices of the Heavy Metal galaxies agree with our main results from full-spectrum modeling. Specifically, we compute the <Fe> and Mgb indices, where <Fe > = 0.5 × (Fe5270 + Fe5335). We then compare these indices to simple spectrum models from Thomas et al. (2003).

Given the prevalence of skylines in the Heavy Metal spectra, we must account for the prominent noise peaks. However, by default, the Lick/IDS system does not use the variance spectrum because it biases the equivalent widths. Thus, to overcome skyline contamination, we follow the approach of van der Wel et al. (2021). For each spectrum, we take the best-fitting alf model and compute its χ². We then scale the noise to χ² = 1. Any pixel that is more than 2σ or 20% offset from the model flux is masked. If more than 30% of the pixels in the index wavelength interval (including the pseudocontinuum intervals) are masked, the index is deemed unreliable. Four Heavy Metal galaxies are removed at this step. For the other 10 galaxies, the masked pixels are replaced by a linear interpolation using the good pixels. We note that the alf models are only used to help identify which pixels should be masked.

Next, we correct the absorption line indices to a constant velocity dispersion. This step is important if we want to compare indices between objects and to stellar population models. In the Lick/IDS system, indices are corrected to have zero intrinsic velocity dispersion and to have a standard wavelength-dependent resolution ranging from 8.4 Å to 11.5 Å (Worthey & Ottaviani 1997). We compute correction factors using alf models at a variety of velocity dispersions (e.g., see Davies et al. 1993; Kuntschner 2004). This correction factor is applied after measuring the equivalent widths. For the Heavy Metal galaxies, the correction ranges from C(σ_v ) = 1.0 to 1.5. With the corrected Lick indices in hand, we finally estimate uncertainties by bootstrapping the spectrum and recomputing equivalent widths 10,000 times, taking the 16th and 84th percentiles as the index uncertainty.

In Figure 8, we present <Fe> as a function of Mgb for 10 Heavy Metal galaxies. We compare these indices with a set of 2 Gyr models from Thomas et al. (2003). Drawing conclusions from this figure is challenging, as numerous Heavy Metal galaxies display index measurements with substantial uncertainties that span the entire model grid. This issue is particularly pronounced for galaxies at a redshift z ∼ 2.1. The increased uncertainties at higher redshifts can be attributed to the younger ages of these galaxies, resulting in less pronounced metal lines. Additionally, significant contamination from skylines in the NIR underscores the necessity for full-spectrum modeling, as opposed to index fitting, at higher redshifts.

Zoom In Zoom Out Reset image size

In this Appendix, we test the robustness of our results as a function of S/N. The primary goal of this test is to check for any systematic uncertainties in full-spectrum modeling, especially as we push the limits of this fitting to lower S/Ns. We begin by generating mock Heavy Metal spectra at various S/Ns, fitting them using alf, and finally comparing the best-fit values and their uncertainties to the true values. It is worth noting that Choi et al. (2014) conducted a similar test in their Appendix A. However, there are a few notable differences in our approach. First, alf has undergone significant updates since the work of Choi et al. (2014). Second, while Choi et al. (2014) used a flat noise spectrum for their simulated spectra, we use actual LRIS+MOSFIRE noise spectra. Finally, the Heavy Metal galaxies have much younger stellar ages, and therefore less pronounced absorption features, than what was tested in Choi et al. (2014).

All mock spectra are generated from a single alf model with stellar population parameters reflecting average values from Table 1. We degrade the model spectrum to the desired S/N at rest-frame 5000 Å using the Heavy Metal noise spectrum. In this way, we can replicate the observations, accounting for wavelength-dependent sensitivities. The model spectrum is also resampled to the LRIS+MOSFIRE resolution using spectres (Carnall 2017). We generate 45 realizations at S/Ns = 8, 12, 20, and 40 and fit them using the same method described in Section 2.2. For all realizations, we compute the difference between the best-fit parameter and the known model parameters.

We present the results in Figure 9. The black rectangles show the average parameter offset for the 45 realizations, with the blue error bars representing the 1σ scatter of these offsets. We also show the typical uncertainties on a single realization in red. The gray tick marks at the bottom of each panel represent the S/Ns of the individual Heavy Metal spectra. All but one Heavy Metal spectrum has S/N ≳ 10 Å⁻¹. σ_v is recovered extremely well with no systematic uncertainties at S/N ≳ 12. The stellar ages are also recovered well, with only a slight (∼0.05 dex) systematic offset toward older ages at S/N ≲ 20. The elemental abundance recovery also shows very little systematic uncertainty. There is a −0.1 dex offset in [Fe/H] at S/N = 8, but at the typical S/N of the Heavy Metal galaxies (S/N > 12), there is no appreciable systematic offset to within 0.05 dex. Similarly, [Mg/Fe] shows a slight systematic offset of 0.05 dex at S/N = 8 and S/N = 12, though at larger S/N, there is no appreciable systematic effect.

Zoom In Zoom Out Reset image size

This parameter recovery test confirms that the abundance results corresponding to the Heavy Metal spectra (S/N ≳ 12 Å⁻¹) have systematic uncertainties of ≲0.05 dex and can be constrained to within at least 0.15 dex. Furthermore, the results are mostly consistent with Choi et al. (2014), despite major updates to alf and mock spectra generated using real LRIS+MOSFIRE noise spectra.