# Modeling the Ages and Chemical Abundances of Elliptical Galaxies

Nicole Marcelina Gountanis David H. Weinberg Department of Astronomy and
Center for Cosmology and AstroParticle Physics, The Ohio State University,
Columbus, OH 43210, USA  Aliza G. Beverage Astronomy Department, University of
California, Berkeley, CA 94720, USA  Nathan R. Sandford Department of
Astronomy and Astrophysics, University of Toronto, 50 St. George Street,
Toronto, ON M5P 0A2, Canada  Charlie Conroy Center for Astrophysics—Harvard &
Smithsonian, Cambridge, MA 02138, USA  Mariska Kriek Leiden Observatory,
Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands

###### Abstract

Spectroscopic studies of elliptical galaxies show that their stellar
population ages, mean metallicity, and $\alpha$-enhancement traced by $[{\rm
Mg/Fe}]$ all increase with galaxy stellar mass or velocity dispersion. We use
one-zone galactic chemical evolution (GCE) models with a flexible star
formation history (SFH) to model the age, $[{\rm Mg}/{\rm H}]$, and $[{\rm
Mg/Fe}]$ inferred from simple stellar population (SSP) fits to observed
ellipticals at $z\sim 0$ and $z\sim 0.7$. We show that an SSP fit to the
spectrum computed from a full GCE model gives ages and abundances close to the
light-weighted, logarithmically averaged values of the composite stellar
population, ${\langle{\rm age}\rangle}$, ${\langle[{\rm Mg}/{\rm H}]\rangle}$,
and ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$. With supernova Mg and Fe yields
fixed to values motivated by Milky Way stellar populations, we find that
predicted ${\langle[{\rm Mg}/{\rm H}]\rangle}$-${\langle{\rm age}\rangle}$ and
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$-${\langle{\rm age}\rangle}$ relations
are surprisingly insensitive to SFH parameters: older galaxies have higher
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$, but the detailed form of the SFH has
limited impact. The star formation efficiency and outflow efficiency affect
the early and late evolution of ${\langle[{\rm Mg}/{\rm H}]\rangle}$,
respectively; explaining observed trends requires higher star formation
efficiency and lower outflows in more massive galaxies. With core collapse
supernova yields calibrated to the plateau $[{\rm Mg}/{\rm Fe}]_{\rm
cc}\approx 0.45$ observed in many Milky Way studies, our models underpredict
the observed ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ ratios of ellipticals by
0.05-0.1 dex. Increasing the core collapse yield ratio to $[{\rm Mg}/{\rm
Fe}]_{\rm cc}=0.55$ improves the agreement, though the models still lie below
the data. We discuss potential resolutions of this discrepancy, including the
possibility that many ellipticals terminate their star formation with a self-
enriching, terminating burst that reduces the light-weighted age and boosts
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$.

##  1 Introduction

Population synthesis modeling of the spectra of elliptical galaxies provides
insights into the ages and chemical abundances of their stellar populations.
At $z\sim 0$, ellipticals of higher stellar mass or higher velocity dispersion
exhibit systematically older ages, higher metallicities, and larger
enhancements of $\alpha$-elements relative to iron (Thomas et al., 2005;
Schiavon, 2007; Graves et al., 2009; Smith et al., 2009; Johansson et al.,
2012; Conroy et al., 2014; Greene et al., 2019). Recent studies of early-type
galaxies at $z\sim 0.7$ show near-solar metallicity ($[{\rm Fe/H}]\sim 0$) and
strong $\alpha$-enhancement ($[{\rm Mg/Fe}]\sim 0.2-0.3$), but the trends with
velocity dispersion are weaker than those at $z\sim 0$ (Beverage et al., 2021,
2023). The elevated $[{\rm Mg/Fe}]$ ratios of high-redshift ellipticals and
the most massive low-redshift ellipticals are interpreted as a sign of rapid
star formation, with many stars forming before Type Ia supernovae (SNIa) have
had time to make large contributions to the abundance of iron-peak elements.

Traditional observational analyses fit the integrated spectrum of an early-
type galaxy as a simple stellar population (SSP) of specified age and chemical
composition, but of course, the stars form over time from an interstellar
medium (ISM) with steadily evolving elemental abundances. In this paper, we
use parameterized galactic chemical evolution (GCE) models to predict the
light-weighted age, $[{\rm Mg}/{\rm H}]$, and $[{\rm Mg/Fe}]$ of observed
elliptical populations. We also use these models to test how well the SSP fit
recovers the true light-weighted quantities of the full population. We focus
on Mg as a well-measured element that is expected to come entirely from core-
collapse supernovae (CCSN), and Fe as a well-measured element with
contributions from both CCSN and SNIa. We use $[{\rm Mg}/{\rm H}]$ rather than
$[{\rm Fe/H}]$ as our primary measure of overall metallicity because our
predictions for it do not depend on our assumptions about SNIa, and the two
quantities are measured comparably well from observed spectra. We adopt one-
zone GCE models that treat the star-forming gas reservoir as chemically
homogeneous at a given time, an approximation that we consider adequate for
predicting integrated properties, though it does not allow us to model radial
gradients in age or abundance.

The population-averaged Mg and Fe yields of CCSN and the population-averaged
Fe yield and delay time distribution (DTD) of SNIa are crucial inputs to our
models. We adopt values that are empirically motivated by observations of
Milky Way stellar abundances and cosmic SNIa rates, implicitly assuming that
the stellar scale astrophysics that determines these quantities is similar in
elliptical galaxies despite their distinctive star formation histories. We aim
to address three broad questions:

  1. 1.

Given these yields motivated by independent observations, can GCE models with
reasonable parameter values reproduce the ages, metallicities, and
$\alpha$-enhancements of observed ellipticals at $z\sim 0$ and $z\sim 0.7$?

  2. 2.

If so, what do the observed properties imply about the star formation
histories, star formation efficiencies, and outflows of these galaxies?

  3. 3.

If not, what modifications might be required to the stellar yields or other
model assumptions, and what do they imply about the distinctive features of
elliptical galaxies?

We rely principally on analytic GCE models (Weinberg et al. 2017, hereafter
WAF) because they are sufficient for our purposes here. We provide a Python
implementation of these analytic models, fanCE111fanCE (pronounced fan-C-E)
for Flexible Analytic Chemical Evolution, so that others can easily adapt them
to interpret other galaxies or stellar populations or explore alternative
assumptions (see Appendix A). It would be straightforward to extend our
approach using numerical models that allow more flexible assumptions such as
bursty star formation histories and metallicity-dependent yields.

While the high metallicities and $\alpha$-enhancements of ellipticals are
usually attributed to “rapid star formation,” this phrase is imprecise. In our
models, there are four distinct timescales in the star formation history: the
time at which star formation begins, the duration of an initial linear rise,
the $e$-folding timescale of a subsequent decline, and the time at which star
formation ceases entirely. The star formation efficiency, ${\rm
SFE}\equiv\dot{M}_{*}/M_{\rm gas}\equiv\tau_{*}^{-1}$, defines another
important timescale. One goal of our modeling is to tease apart the role of
these distinct timescales in governing ages and chemical enrichment. While our
flexible parameterization of star formation history (SFH) allows a variety of
behaviors, all our models assume a sharp final transition from star formation
to quiescence.

Stellar population trends in the massive quiescent galaxy population were
first identified in the low-redshift universe using Lick indices (Worthey,
1994; Thomas, 1999; Trager et al., 2000; Thomas et al., 2005). These trends
have since been confirmed using a variety of methods, including full-spectrum
fitting with variable abundance patterns (Walcher et al., 2009; Conroy et al.,
2014). By stacking thousands of SDSS spectra, Conroy et al. (2014) measure the
average elemental abundances of 11 different elements as a function of
velocity dispersion. Their analysis reproduces the well-established trends at
$z\sim 0$ with high precision and adds many individual elements.

One way to separate the star-formation histories from the assembly histories
of massive quiescent galaxies is by linking the $z\sim 0$ galaxies to their
progenitors at higher redshift. However, measuring stellar population
properties and chemical compositions becomes increasingly challenging at
higher redshift, relying on already faint absorption lines that are shifted to
near-infrared wavelengths. Studies of large samples of quiescent galaxies at
$z\sim 0.7$ find that their abundance patterns look remarkably similar to
today’s massive elliptical galaxies (Beverage et al., 2023). However, the
lower velocity dispersion quiescent galaxy population at $z\sim 0.7$ appears
to have both higher [$\alpha$/Fe] and earlier formation redshifts compared to
the $z\sim 0$ population (Leethochawalit et al., 2018; Beverage et al., 2023).
This finding shows that the $z\sim 0$ ellipticals are not simply the passively
evolved descendants of the $z\sim 0.7$ quiescent population.

Galaxies in observational samples such as Conroy et al. (2014) and Beverage et
al. (2023) are selected by a combination of early-type (bulge-dominated)
morphology and photometric or spectroscopic diagnostics that imply little or
no ongoing star formation. Throughout this paper we use the term “ellipticals”
as shorthand for this combination of characteristics, even though many
galaxies in these samples would not have a strictly elliptical morphological
classification. From the point of view of chemical evolution, the critical
feature of these galaxies is their quiescent star formation history, with
minimal or zero star formation for at least the last Gyr prior to the epoch at
which they are observed. Because the galaxies in our observational samples
have a minimum velocity dispersion $\sigma\geq 80\,{\rm km}\,{\rm s}^{-1}$, we
sometimes refer to them as “massive ellipticals” or “massive quiescent
galaxies” to distinguish them from dwarf elliptical or dwarf spheroidal
populations, which may have different evolutionary histories.

Full GCE modeling allows us to go beyond the empirical trends of SSP fits to
the physical processes that drive them — nucleosynthesis, SFH, SFE, feedback,
and quenching. Section 2 describes our observational samples and SSP fitting
tools and summarizes tests of our basic methodology, predicting SSP quantities
from the light-weighted average quantities of composite stellar populations.
Section 3 describes our GCE models, with particular attention to our flexible
SFH parameterization, our choice of empirically motivated Mg and Fe yields,
and the impact of model parameters on the evolution of the mean age, $[{\rm
Mg}/{\rm H}]$, and $[{\rm Mg/Fe}]$ of stellar populations. In Section 4 we
compare our model predictions to the observed trends. Our models struggle to
produce $[{\rm Mg/Fe}]$ values as high as the observed values, so we also
consider variant models with reduced SNIa Fe production or a terminating burst
of star formation. Section 5 provides further discussion of nucleosynthetic
yields and identifies directions of future investigation. We summarize our
conclusions in Section 6. Appendix A gives instructions for accessing fanCE;
Appendix B gives details of our treatment of terminating bursts of star
formation; and Appendix C tests the relation between SSP fits and the light-
weighted average quantities of composite stellar populations.

##  2 Observational data

Table 1: Stellar age and abundance results from stacked SDSS ($z\sim 0$) and LEGA-C ($z\sim 0.7$) spectra in bins of velocity dispersion. All values are from Beverage et al. (2023).  $\sigma$ ($\,{\rm km}\,{\rm s}^{-1}$) | Age (Gyr) | [Fe/H] | [Mg/H] | [Mg/Fe] |
---|---|---|---|---|---
|  | SDSS |  |  |
87 | $5.5^{+0.1}_{-0.1}$ | $-0.138^{+0.005}_{-0.005}$ | $0.007^{+0.005}_{-0.005}$ | $0.145^{+0.004}_{-0.004}$ |
112 | $6.0^{+0.1}_{-0.1}$ | $-0.109^{+0.005}_{-0.005}$ | $0.065^{+0.005}_{-0.005}$ | $0.174^{+0.004}_{-0.004}$ |
138 | $6.7^{+0.6}_{-0.1}$ | $-0.085^{+0.008}_{-0.020}$ | $0.106^{+0.007}_{-0.019}$ | $0.191^{+0.005}_{-0.004}$ |
170 | $7.9^{+0.1}_{-0.1}$ | $-0.072^{+0.005}_{-0.006}$ | $0.155^{+0.006}_{-0.005}$ | $0.227^{+0.004}_{-0.004}$ |
200 | $9.4^{+0.1}_{-0.1}$ | $-0.064^{+0.004}_{-0.005}$ | $0.188^{+0.005}_{-0.005}$ | $0.252^{+0.004}_{-0.004}$ |
246 | $10.3^{+0.1}_{-0.1}$ | $-0.059^{+0.005}_{-0.005}$ | $0.240^{+0.005}_{-0.005}$ | $0.299^{+0.004}_{-0.004}$ |
299 | $10.7^{+0.2}_{-0.2}$ | $-0.012^{+0.007}_{-0.007}$ | $0.292^{+0.007}_{-0.007}$ | $0.305^{+0.005}_{-0.005}$ |
|  | LEGA-C |  |  |
147 | $3.0^{+0.2}_{-0.2}$ | $-0.10^{+0.03}_{-0.03}$ | $0.18^{+0.03}_{-0.03}$ | $0.28^{+0.03}_{-0.03}$ |
176 | $3.8^{+0.3}_{-0.3}$ | $-0.09^{+0.03}_{-0.03}$ | $0.21^{+0.02}_{-0.03}$ | $0.30^{+0.03}_{-0.03}$ |
207 | $3.6^{+0.4}_{-0.3}$ | $-0.05^{+0.03}_{-0.03}$ | $0.22^{+0.03}_{-0.03}$ | $0.27^{+0.03}_{-0.02}$ |
247 | $4.4^{+0.3}_{-0.3}$ | $-0.03^{+0.02}_{-0.02}$ | $0.28^{+0.03}_{-0.02}$ | $0.31^{+0.03}_{-0.02}$ |

For our observational comparisons in Section 4 we use the stellar ages,
[Mg/H], and [Mg/Fe] reported by Beverage et al. (2023), who measure these
quantities for massive quiescent galaxies in bins of velocity dispersion using
full-spectrum absorption line fitting. Beverage et al. (2023) examine two
galaxy samples, one at $z\sim 0$ and one at $z\sim 0.7$.

The $z\sim 0$ sample was originally selected by Conroy et al. (2014) from SDSS
Data Release 7 (Abazajian et al., 2009). Galaxies were chosen from a narrow
redshift interval $0.025<z<0.06$ and identified as quiescent using emission
line fluxes. Structural properties were also used to identify galaxies with
early-type morphologies. The sample was split into seven evenly spaced bins of
velocity dispersion centered at $\sigma=87-299\,{\rm km}\,{\rm s}^{-1}$.

The $z\sim 0.7$ sample was selected by Beverage et al. (2023) from the LEGA-C
survey. LEGA-C is a deep spectroscopic survey of massive galaxies at
intermediate redshifts ($0.5\lesssim z\lesssim 1.0$) conducted with VIMOS on
the VLT (van der Wel et al., 2016; Straatman et al., 2018; van der Wel et al.,
2021). Galaxies were identified as quiescent using their rest-frame $UVJ$
colors (Wuyts, 2007; Muzzin et al., 2013). The sample was limited to galaxies
with a rest-frame wavelength coverage of key absorption features (H$\beta$,
MgI, FeI), translating to a redshift selection of $z\leq 0.75$. The $z\sim
0.7$ galaxies were split equally among four bins of velocity dispersion
centered at $\sigma=147-247\,{\rm km}\,{\rm s}^{-1}$.

The average stellar population properties of the $z\sim 0$ and $z\sim 0.7$
galaxies in bins of velocity dispersion were derived using different methods.
The spectra of the $z\sim 0$ galaxies were combined (“stacked”) to produce a
high signal-to-noise (S/N) spectrum in each velocity dispersion bin before
fitting. For the $z\sim 0.7$ sample, each galaxy spectrum was fit
individually, and the individual parameter estimates were combined to obtain
high-precision values in each velocity dispersion bin. This hierarchical
Bayesian method, unlike the stacking method, avoids the need for smoothing and
interpolating the spectra and polynomial fitting, all of which can introduce
correlated noise and diminish the strength of absorption features. While the
hierarchical Bayesian modeling method is distinctly different from fitting
stacked spectra, Beverage et al. (2023) show that the two methods produce
similar average properties. More details and comparisons between methods can
be found in Section 3 of Beverage et al. (2023).

###  2.1 Stellar Population Fitting

The stellar population fitting is accomplished using the full-spectrum
absorption line fitting code, alf (Absorption-Line Fitter) (Conroy & van
Dokkum, 2012a; Conroy et al., 2018). alf generates simple stellar population
(SSP) models with variable elemental abundances for the optical-NIR. The code
employs an MCMC wrapper (Foreman-Mackey et al., 2013) to estimate the
posteriors of the stellar population properties and elemental abundances.

The SSP model ingredients include metallicity-dependent MIST isochrones (Choi
et al., 2016) and the MILES and Extended-IRTF empirical spectral libraries
(Sánchez-Blázquez et al., 2006; Villaume et al., 2017). To alter the
abundances of the 19 individual elements in these models, alf uses
metallicity- and age-dependent synthetic response functions to determine the
fractional change in spectra due to the variation of each element. alf models
the continuum-normalized spectrum, so it does not incorporate constraints from
the overall shape of the spectral energy distribution (SED).

Conroy et al. (2014) also used alf to measure the stellar population
properties of the SDSS stacked spectra used here. However, since the
publication of Conroy et al. (2014), the SSP models have been improved, with
metallicity- and age-dependent response functions, updated isochrones, and
improved interpolation of the empirical spectral libraries (Villaume et al.,
2017; Conroy et al., 2018). Therefore in this analysis, we adopt the updated
SDSS measurements from Beverage et al. (2023).

Figure 1: (Top) Stacked SDSS spectrum ($z\sim 0$) of a velocity dispersion bin
(black) and the best-fit SSP spectrum (red) computed with alf. The blue curve
shows the alf model with $[{\rm Fe/H}]$ raised from its best-fit value of
$\sim 0.05$ to 0.0, a change of $3\times$ the fit uncertainty in $[{\rm
Fe/H}]$. Other parameters are held fixed. The lower panel plots the fractional
residuals between the two SSP models and the observed spectrum. (Bottom)
similar to the top, but for an individual LEGA-C galaxy spectrum at $z\sim
0.7$. The blue curve again shows a model with $[{\rm Fe/H}]$ perturbed by
$3\sigma$ from the best-fit value. Gray bands show the statistical noise per
spectral pixel. In the top panels, spectra are binned such that one pixel is
$\simeq 1$ Å in the rest-frame for both SDSS and LEGA-C.

The best-fit stellar population properties for the $z\sim 0$ and $z\sim 0.7$
populations are tabulated in Table 1. Figure 1 shows the observed spectrum and
best-fit model spectrum for one velocity dispersion stack from the SDSS (top)
and an individual LEGA-C galaxy (bottom). We also show the spectrum for a
model in which we increase the total metallicity, [Z/H] by 3$\sigma$ while
keeping other parameters fixed. The $\chi^{2}/{\rm dof}$ for each model is
shown in the upper left. The residuals and $\chi^{2}/{\rm dof}$ for the
perturbed model are larger, as expected, but the differences are subtle. The
response of the predicted spectrum to changes in $[{\rm Fe/H}]$ or $[{\rm
Mg}/{\rm H}]$ is spread over much of the spectrum, as shown by Conroy et al.
(2014, 2018), so a high S/N allows precise constraints even though the changes
are small. Each velocity dispersion bin of the $z\sim 0.7$ sample contains
30-40 galaxies, so the effective S/N for determining the mean stellar
properties is much higher than that of the individual LEGA-C spectrum shown
here.

###  2.2 SSP Fits vs. Light-Weighted Averages

The age and abundance measurements in Table 1 are derived by modeling the
observed spectrum as produced by a single stellar population, even though the
light is a composite of populations with different ages and abundances. This
approach is a long-standing practice in the study of elliptical galaxies
(e.g., Worthey, 1994; Thomas et al., 2003; Vazdekis et al., 2010; Conroy et
al., 2014). With the flexible GCE models used in this paper, it would be
possible to compute composite spectra for each model and compare them to data
directly in spectral space. However, in this paper we adopt the simpler
procedure of computing the light-weighted average quantities from models and
comparing them to the empirical values from SSP fits.

Younger stellar populations have higher luminosity per unit stellar mass, so
to compute model averages that can be compared to SSP fits we must include the
impact of varying $L/M$ ratio. In Figure 2, orange curves show the V-band
$L_{\mathrm{V}}/M$ predicted by alf as a function of population age for three
different $[{\rm Fe/H}]$ values. The blue curve shows an approximation that we
use to weight stellar populations in our models,

$w(t)_{\mathrm{light}}=L_{\mathrm{V}}/M=1.5\left(\frac{t_{\mathrm{obs}}-t}{\mathrm{Gyr}}+0.1\right)^{-\beta}\frac{L_{\odot}}{M_{\odot}},$
(1)

with $\beta=0.90$. Here, $t_{\mathrm{obs}}$ is the age of the universe at the
galaxy’s observed redshift ($t_{\mathrm{obs}}=7$ and 14 Gyr for $z=0$ and
$z=0.7$, respectively), and $t_{\mathrm{obs}}-t$ is the age of the stellar
population formed at time $t$. We add 0.1 to avoid divergence at
$t=t_{\mathrm{obs}}$, but our model predictions are insensitive to this
choice, and the youngest stellar populations in the models we compare to data
in Section 2 have ages of 1-2 Gyr. The metallicity dependence of
$L_{\mathrm{V}}/M$ is small compared to the age dependence, so we ignore it in
our light-weighting. We have checked that our predictions are insensitive to
the precise value of $\beta$, though ignoring light-weighting entirely
($\beta=0$) gives biased results.

We investigate the accuracy of our light-weighted average method by performing
a series of simple tests similar to those of Serra & Trager (2007). We
describe these tests and their results in Appendix C. In brief, we use alf to
compute the spectrum of a composite stellar population (CSP) comprised of two
SSPs that have different $[{\rm Mg}/{\rm H}]$, $[{\rm Mg/Fe}]$, or age. We
then fit this composite spectrum with alf assuming a single population, and we
compare the SSP parameters to the light-weighted average of the two input
populations. We find excellent agreement between the averaged and SSP-fit
values of $[{\rm Mg}/{\rm H}]$ and $[{\rm Mg/Fe}]$, and our tests show that it
is more accurate to average these logarithmic quantities than to average
(Mg/H) and (Mg/Fe) before taking logarithms. Similarly, we find that the mean
log(age) agrees with the SSP value, and that linear age averaging is less
accurate. For large age differences between the two populations, the log(age)
average slightly overestimates the SSP value, e.g., by 0.1 dex for populations
that differ by 0.7 dex. Our findings are consistent with those of Serra &
Trager (2007).

In addition to these two-population tests, we conduct an end-to-end test using
continuous GCE models in Section 4.2 (see Figure 13). These tests confirm the
viability of our approach, with the caveat that SSP ages are slightly younger
($\sim 1$ Gyr) than the light-weighted mean ages for the oldest galaxies
(${\langle{\rm age}\rangle}>10$ Gyr) at $z\sim 0$.

Figure 2: $L_{\mathrm{V}}/M$ as a function of stellar population age. The
V-band luminosity is used to compute light-weighted quantities from the
SFH+GCE models. The solid orange line shows the true $L_{\mathrm{V}}/M$, as
measured from the Conroy et al. (2014) SSP models, while the blue line shows
the adopted approximation (Eq. 1). The dashed/dotted transparent orange lines
correspond to SSP models with varying [Fe/H] = (-1, -0.5, -0.2). The adopted
light-weighing function appropriately captures the true $L_{\mathrm{V}}/M$ and
is robust to variations in [Fe/H].

##  3 Galactic chemical evolution models

###  3.1 Star Formation History

To interpret the observational data described in Section 2, we use one-zone
models with accretion and outflows treated using the analytic solutions of
WAF. We adopt a “rise-fall” star formation history (Johnson et al., 2021b),

$\rm{SFH}\equiv\dot{M}_{\star}(t)\propto(1-e^{-t/\tau_{1}})\cdot
e^{-t/\tau_{2}},$ (2)

with parameters $\tau_{1}$ and $\tau_{2}$ describing the rising and falling
SFH timescales, respectively. For modeling elliptical galaxies, we consider
this form superior to the one-parameter exponential SFH, in which
$\dot{M}_{\star}(t)$ begins immediately at its maximum value, or the linear-
exponential (a.k.a. “delayed tau”) model, in which the timescale of the rising
phase is locked to that of the declining phase. The application of WAF
solutions to this SFH is described in the Appendix of Weinberg et al. (2024),
hereafter WGJT24. Although the $\tau_{1}$ parameter provides flexibility in
the early SFH, we do not want to require that the initially linear form
already applies at $t=0$ (i.e., at extremely high redshift). We therefore
allow $t_{\mathrm{start}}$ to be another adjustable parameter, with
$t_{\mathrm{start}}=0.5$ Gyr as our reference choice and $t_{\mathrm{start}}=$
0 and 1 Gyr as alternatives. We take $t_{\mathrm{obs}}=14$ Gyr and
$t_{\mathrm{obs}}=7$ Gyr as the age of the Universe at $z=0$ and $z=0.7$,
respectively.

Figure 3: Examples of star formation histories for different $\tau_{1}$,
$\tau_{2}$, $t_{\mathrm{start}}$ (see Equation 2). The solid, orange curve
represents our reference SFH with $\tau_{1}=2$ Gyr, $\tau_{2}=8$ Gyr, and
$t_{\mathrm{start}}=0.5$ Gyr. The solid blue and dashed orange curves show the
impact of longer rising timescales ($\tau_{1}=4$ Gyr) and shorter declining
timescales ($\tau_{2}=4$ Gyr), respectively. The dotted orange curve shows the
reference SFH with an earlier start of star formation ($t_{\mathrm{start}}=0$
Gyr). The SFR normalization gives the fraction of the galaxy’s stellar mass
formed per Gyr. Our elliptical galaxy SFH models are parameterized by
$\tau_{1}$, $\tau_{2}$, $t_{\rm start}$, and the time $t_{\rm cut}$ at which
star formation cuts off abruptly.

Figure 3 illustrates SFHs for a variety of parameter choices. For most
calculations in this section, we will consider variations about a reference
model with $t_{\mathrm{start}}=0.5$ Gyr, $\tau_{1}=2$ Gyr, and $\tau_{2}=8$
Gyr (solid orange curve). For $\tau_{1}\ll\tau_{2}$, the SFH reaches a maximum
at $t\approx t_{\mathrm{start}}+\tau_{1}$, so 2.5 Gyr ($z\approx$ 2.6) in our
reference case, before settling into exponential decline at $t\gg\tau_{1}$.
Changing $\tau_{2}$ to 4 Gyr leads to a sharper maximum and a more rapid
decline, while changing $\tau_{1}$ to 4 Gyr leads to a slower rise and a
later, broader maximum. Changing $t_{\mathrm{start}}$ to 0 simply translates
the reference SFH in time. For our elliptical models, we also assume that the
SFH truncates sharply at a time $t_{\mathrm{cut}}$. Our model SFHs have four
adjustable timescale parameters — $t_{\mathrm{start}}$, $\tau_{1}$,
$\tau_{2}$, and $t_{\mathrm{cut}}$ — and $t_{\mathrm{cut}}$ is the one with
the largest impact on the predicted mean stellar age and abundances. The form
of Equation 2 allows a much wider range of behavior than shown in Figure 3,
including a linearly rising SFH, a pure exponential, and a rise followed by a
constant SFR plateau. The parameter choices here represent values we consider
plausible for massive ellipticals.

###  3.2 Yields

Nucleosynthetic yields play a critical role in any GCE calculation.
Theoretically predicted yields have uncertainties associated with supernova
models, black hole formation, and the stellar IMF. We instead base our choice
of yields on the empirical analysis of WGJT24, who anchor the overall scale of
yields to the Rodríguez et al. (2023) determination of the mean Fe yield of
CCSN. This determination implies (with some uncertainties described in WGJT24)
an IMF-averaged CCSN yield

$\frac{y^{\rm{cc}}_{\rm{Fe}}}{Z_{\rm{Fe,\odot}}}=0.345,$ (3)

where $Z_{\rm{Fe,\odot}}=0.00137$ is the solar Fe mass fraction and
$y^{\rm{cc}}_{\rm{Fe}}$ is the mass of iron produced per unit mass of star
formation (and is therefore dimensionless). For the IMF-averaged Mg yield we
take

$\frac{y^{\rm{cc}}_{\rm{Mg}}}{Z_{\rm{Mg,\odot}}}=0.973\cdot 10^{[{\rm Mg}/{\rm
Fe}]_{\rm cc}-0.45}.$ (4)

In combination with Equation 3, Equation 4 implies that pure CCSN enrichment
produces a Mg/Fe ratio that corresponds to $[{\rm Mg/Fe}]=[{\rm Mg}/{\rm
Fe}]_{\rm cc}$. Following WGJT24, we adopt $[{\rm Mg}/{\rm Fe}]_{\rm cc}=0.45$
as our reference value based on the abundance plateau observed for low
metallicity Milky Way disk and halo stars. We also consider a higher CCSN
yield ratio, corresponding to $[{\rm Mg}/{\rm Fe}]_{\rm cc}=0.55$, which
implicitly assumes that stars with $[{\rm Mg/Fe}]=0.45$ already have a
significant fractional Fe contribution from SNIa (Maoz & Graur, 2017; Conroy
et al., 2022; Mason et al., 2023). We assume that the IMF-averaged Mg and Fe
yields are metallicity-independent.

For the delay time distribution (DTD) of SNIa we adopt a sum of two
exponentials that accurately approximates the $t^{-1.1}$ power-law form
advocated by Maoz & Graur (2017). For the minimum delay time we take
$t_{\mathrm{d,min}}=0.15$ Gyr as our reference value and consider
$t_{\mathrm{d,min}}=0.05$ Gyr as an alternative that increases the amount of
early SNIa enrichment. The choice of the exponential parameters is described
by WAF. The most important SNIa parameter is the time-integrated Fe yield,
which we set with the condition

$y_{\rm Fe}^{\rm Ia}=y_{\rm Fe}^{\rm cc}\left(10^{\Delta[\alpha/{\rm
Fe}]}-1\right)\mu^{-1}~{},$ (5)

where $\Delta[\alpha/{\rm Fe}]$ is the drop between the CCSN plateau and the
late-time equilibrium $[\alpha/{\rm Fe}]$, and $\mu\approx 1$ is a parameter
that accounts for the impact of SFH on this equilibrium. Following WGJT24 (see
their Section 3.2), we set $\Delta[\alpha/{\rm Fe}]=[{\rm Mg}/{\rm Fe}]_{\rm
cc}$ and $\mu=1.1$, so that a GCE model with a slowly declining SFH evolves to
$[{\rm Mg/Fe}]\approx 0$ at late times. The implied $y_{\rm Fe}^{\rm Ia}$ is
also consistent with empirical estimates of the SNIa rate and a mean Fe yield
of $\approx 0.7$ $M_{\odot}$ per supernova.

Combining Equations 3-5 with our adopted solar abundances gives

$\\{y_{\rm Mg}^{\rm cc},y_{\rm Fe}^{\rm cc},y_{\rm Fe}^{\rm
Ia}\\}=\\{6.52,4.73,7.82\\}\times 10^{-4}$ (6)

for models with $[{\rm Mg}/{\rm Fe}]_{\rm cc}=0.45$ and

$\\{y_{\rm Mg}^{\rm cc},y_{\rm Fe}^{\rm cc},y_{\rm Fe}^{\rm
Ia}\\}=\\{8.21,4.73,10.95\\}\times 10^{-4}$ (7)

for models with $[{\rm Mg}/{\rm Fe}]_{\rm cc}=0.55$.

Details aside, the core assumption of our models is that we can take yields
motivated by empirical studies of the Milky Way and general observations of
supernova rates and galaxy properties, and apply them to the specific case of
elliptical galaxy evolution. This assumption is reasonable if the stellar
scale astrophysics that governs properties such as the stellar initial mass
function (IMF) and supernova yields is decoupled from the galactic scale
astrophysics that governs star formation history and morphology. However, it
is not guaranteed to hold. We return to this issue in Section 5.

###  3.3 Star formation efficiency and outflows

The other parameters of our GCE models are the SFE timescale, $\tau_{*}\equiv
M_{g}/\dot{M}_{\star}$, and the outflow mass loading factor,
$\eta\equiv\dot{M}_{\mathrm{out}}/\dot{M}_{\star}$. The SFE timescale has an
important impact on the rate of metal enrichment at early times, and $\eta$
governs the metallicity at late times by regulating the degree of metal loss
from the ISM. In the molecular gas of the present-day star-forming galaxies, a
typical value of $\tau_{*}$ is 2 Gyr, with substantial variations (Leroy et
al., 2008; Sun et al., 2023). Studies of higher redshift galaxies suggest that
this timescale decreases as $(1+z)^{-0.5}$ (Tacconi et al., 2018), and the SFE
may be higher than the average in galaxies that evolve to become ellipticals.
We therefore adopt $\tau_{*}=1$ Gyr as our reference value and consider values
as short as 0.5 Gyr.

It has long been conjectured that the galaxy mass-metallicity relation is
driven primarily by a mass dependence of $\eta$ (Finlator & Davé, 2008;
Peeples & Shankar, 2011; Davé et al., 2012; Zahid et al., 2012; Chartab et
al., 2023), with lower values of $\eta$ at higher stellar and halo mass. There
is a degeneracy between the overall scale of yields and the values of $\eta$
required to reproduce the observed mass-metallicity relation (see discussions
by WGJT24). For our adopted yields, the implied values of $\eta$ are
relatively low. We take $\eta=0.3$ as a reference value, but in line with
previous work, we will assume that $\eta$ decreases with increasing
$M_{\star}$ to reproduce the observed mass-metallicity trend.

The WAF analytic solutions assume metallicity-independent yields, time-
independent values of $\tau_{*}$ and $\eta$, instantaneous enrichment of CCSN
products, and an SNIa DTD that is exponential or (in our case) a sum of
exponentials. These solutions also assume that a fraction $r$ of a stellar
population’s mass is immediately returned to the ISM with its birth
metallicity, where the recycling fraction $r\approx 0.4$ for a Kroupa (2001)
IMF. We have used the numerical GCE code VICE (Johnson & Weinberg, 2020) to
compare results for instantaneous recycling and full, time-dependent recycling
for our typical parameter choices, and we find that changes to the abundance
evolution relevant to this paper are negligible, consistent with the small
impact found by WAF.

For a single exponential SFH, with $\tau_{1}=0$ Gyr and $\dot{M}_{*}\propto
e^{-t/\tau_{2}}$, the time evolution of the ISM mass fraction of a pure CCSN
element such as Mg follows

$Z_{{\rm Mg}}={y_{\rm Mg}^{\rm cc}\over
1+\eta-r-\tau_{*}/\tau_{2}}\left(1-e^{-t/\bar{\tau}}\right)$ (8)

with

$\bar{\tau}={\tau_{*}\over 1+\eta-r-\tau_{*}/\tau_{2}}~{}.$ (9)

For our more general SFH, the abundance still approaches the same late-time
equilibrium

$Z_{\rm Mg,eq}={y_{\rm Mg}^{\rm cc}\over 1+\eta-r-\tau_{*}/\tau_{2}}~{},$ (10)

but the time dependence is more complex. From these equations, it is clear
that $\eta$ plays a key role in governing late-time metallicity.

At early times, $t\ll\bar{\tau}$, Equation 8 gives $Z_{\mathrm{Mg}}=y_{\rm
Mg}^{\rm cc}(t/\tau_{*})$, independent of $\eta$. A longer $\tau_{*}$ leads to
slower metallicity evolution. A shorter $\tau_{2}$ at fixed $\tau_{*}$ leads
to a higher equilibrium abundance but a slower approach to that equilibrium
because of the $-\tau_{*}/\tau_{2}$ term in the denominators of Equations
8-10. Because our observable quantities are calibrated to solar abundance,
e.g. $[{\rm Mg}/{\rm H}]=\log(Z_{{\rm Mg}}/Z_{{\rm Mg},\odot})$, it is the
scaled yields $y_{\mathrm{Mg}}/Z_{{\rm Mg},\odot}$ and
$y_{\mathrm{Fe}}/Z_{{\rm Fe},\odot}$ that matter, which is why we specify
these scaled values in Equations 3 and 4.

Table 2 summarizes our model parameters and reference choices of these
parameters.

Table 2: GCE Model Parameters Parameter | Description | Reference Value
---|---|---
Smooth SFH |  |
$\tau_{*}$ | =$M_{g}$/$\dot{M}_{\star}$, (SFE) timescale | 1 Gyr
$\eta$ | =$\dot{M}_{out}/\dot{M}_{\star}$, mass loading factor | 0.3
$t_{\rm start}$ | star formation start time | 0.5 Gyr
$\tau_{1}$ | SFH rise timescale | 2 Gyr
$\tau_{2}$ | SFH decline timescale | 8 Gyr
$t_{\mathrm{cut}}$ | star formation cutoff time | …
$y_{\rm Mg}^{\rm cc}/Z_{\mathrm{Mg},\odot}$ | IMF-integrated CCSN Mg yield | 0.973
$y_{\rm Fe}^{\rm cc}/Z_{\mathrm{Fe},\odot}$ | IMF-integrated CCSN Fe yield | 0.345
$y_{\rm Fe}^{\rm Ia}/Z_{\mathrm{Fe},\odot}$ | time-integrated SNIa Fe yield | 0.571
$t_{{\rm d,min}}$ | minimum delay time for SNIa | 0.15 Gyr
$r$ | mass recycling parameter | 0.4
Terminating Burst Models |  |
$F_{D}$ | accreted gas/existing gas | 0, 3
$f_{g}$ | fraction of gas consumed in burst | 0.05, 0.3
$\eta^{\prime}$ | burst mass loading factor | 0
$r^{\prime}$ | burst mass recycling parameter | 0.3
Our GCE models utilize a smooth SFH defined by four parameters: $t_{\rm
start}$, $\tau_{1}$, $\tau_{2}$, and $t_{\rm cut}$. Other parameters of our
GCE models describe the star formation efficiency timescale $\tau_{*}$,
outflows $\eta$, and supernova yields. In addition, our burst model adds a
terminating burst of star formation to a smooth SFH at $t=t_{\rm cut}$ and
contains parameters that describe gas accretion $F_{D}$ and the fraction of
gas consumed in the burst $f_{g}$.

###  3.4 Light vs. Mass Weighting

As discussed in Section 2, our tests (Figure 13, Appendix C) show that the
inferred stellar population properties and elemental abundances from SSP-
equivalent measurements are similar to the light-weighted, log-averaged
abundances and ages in our GCE models. The log-averaged abundance ratio
$\langle[\mathrm{X/H}]\rangle$ of element $X$ is defined as

$\langle[\mathrm{X/H}]\rangle=\frac{\int_{0}^{t_{\mathrm{cut}}}\dot{M}_{\star}(t)\cdot[\mathrm{X/H}](t)\cdot
w(t)\,dt}{\int_{0}^{t_{\mathrm{cut}}}\dot{M}_{\star}(t)\cdot w(t)\,dt},$ (11)

where $\dot{M}_{\star}(t)$ is the SFR, $[{\rm X}/{\rm H}](t)$ is the gas-phase
abundance, and $w(t)$ is the light-weighting function described by Equation 1.
This function assigns more weight to recently formed stars to account for the
fact that younger stars are more luminous. For the mass-weighted quantities
shown in some plots, we set $w(t)=1$.

From ${\langle[{\rm Mg}/{\rm H}]\rangle}$ and ${\langle[{\rm Fe}/{\rm
H}]\rangle}$ we compute ${\langle[{\rm Mg}/{\rm Fe}]\rangle}={\langle[{\rm
Mg}/{\rm H}]\rangle}-{\langle[{\rm Fe}/{\rm H}]\rangle}$. We determine the
light-weighted, log-averaged age from

$\langle\log_{10}(\mathrm{age})\rangle=\frac{\int_{0}^{t_{\mathrm{cut}}}\dot{M}_{\star}(t)\cdot\log_{10}(t_{\mathrm{obs}}-t)\cdot
w(t)\,dt}{\int_{0}^{t_{\mathrm{cut}}}\dot{M}_{\star}(t)\cdot w(t)\,dt}.$ (12)

For brevity, we use the notation ${\langle{\rm age}\rangle}\equiv
10^{\langle\log_{10}(\mathrm{age})\rangle}$ throughout this work.

Figure 4: Influence of light-weighting on averaged stellar quantities for two
example SFHs, shown by the blue and orange curves in panel (a). The legends
apply to all panels: all models have $\tau_{1}=1$ Gyr, $\tau_{*}=1$ Gyr, and
$\eta=0.3$, blue and orange curves show models with $\tau_{2}=$ 4 and 12 Gyr,
respectively, and solid and dashed curves (sometimes indistinguishably) show
the light- and mass-weighted quantities, respectively. Panel (b) shows the
light- and mass-weighted ${\langle{\rm age}\rangle}$ vs. time. Panels (c) and
(e) show the evolution of gas phase $[{\rm Mg}/{\rm H}]$ and $[{\rm Mg/Fe}]$.
Panels (d) and (f) show the corresponding mean stellar quantities, with light-
weighting producing slightly lower ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ at
late times. Panels (g) and (h) plot ${\langle[{\rm Mg}/{\rm H}]\rangle}$ and
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ against ${\langle{\rm age}\rangle}$,
which nearly cancels out the impact of light-weighting on ${\langle[{\rm
Mg}/{\rm Fe}]\rangle}$.

Figure 4 shows the influence of light-weighting for two example SFH models. We
consider models with a quickly declining SFH ($\tau_{2}=4$ Gyr, shown in blue)
and a slowly declining SFH ($\tau_{2}=12$ Gyr, shown in orange) for fixed
$\tau_{1}=1$ Gyr, $\tau_{*}=1$ Gyr, and $\eta=0.3$. The normalized SFR in
Figure 4a gives the fraction of the galaxy’s stellar mass formed per Gyr. A
galaxy with a quickly declining SFH forms most of its stars early on, while a
galaxy with a slowly declining SFH forms stars more steadily over time. As a
result, the stellar ${\langle{\rm age}\rangle}$ at the time of observation is
younger for galaxies with a slowly declining SFH versus a quickly declining
SFH, as shown by the difference between the blue and orange dashed curves in
Figure 4b. Light-weighting, represented by the solid curves, gives more weight
to younger stars, lowering the ${\langle{\rm age}\rangle}$ for both models by
a couple of Gyr at late times, with deviations from the mass-weighted
${\langle{\rm age}\rangle}$ rapidly increasing at times close to
$t_{\mathrm{obs}}$.

Figures 4c and 4e show the evolution of gas-phase $[{\rm Mg}/{\rm H}]$ and
$[{\rm Mg/Fe}]$ for each model. $[{\rm Mg}/{\rm H}]$ initially increases over
time as Mg is produced by CCSN, then plateaus as the production of Mg reaches
an equilibrium with the removal of Mg by stellar lockup and outflows. $[{\rm
Mg/Fe}]$ starts out at the CCSN yield ratio, here $[{\rm Mg}/{\rm Fe}]_{\rm
cc}=0.45$, then decreases as SNIa begin contributing Fe. $[{\rm Mg/Fe}]$ is
lower for a galaxy with a rapidly declining SFR because the ratio of CCSN to
delayed SNIa is lower at late times. Figures 4d and 4f show the corresponding
mean stellar quantities, which have smaller differences between models. Light-
weighting has minimal impact on ${\langle[{\rm Mg}/{\rm H}]\rangle}$ because
the gas-phase $[{\rm Mg}/{\rm H}]$ is flat at late times. However, light-
weighting results in lower ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ at late times
because it gives higher weight to the populations formed with the lowest
$[{\rm Mg/Fe}]$. Figures 4g and 4h plot ${\langle[{\rm Mg}/{\rm H}]\rangle}$
and ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ versus ${\langle{\rm age}\rangle}$.
Somewhat coincidentally, the impact of light-weighting on ${\langle{\rm
age}\rangle}$ and on ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ nearly cancels, so
mass-weighted and light-weighted tracks are nearly identical.

In Figure 5 we explore the influence of start time on evolution for our
reference model with $\tau_{1}=2$ Gyr, $\tau_{2}=8$ Gyr, $\tau_{*}=1$ Gyr, and
$\eta=0.3$. Incorporating a nonzero $t_{\mathrm{start}}$ is nearly identical
to simply shifting the $t_{\mathrm{start}}=0$ curve to later $t$ or younger
${\langle{\rm age}\rangle}$ by $t_{\mathrm{start}}$, despite the effects of
light-weighting. This near-perfect equivalence is convenient as it allows us
to treat different $t_{\mathrm{start}}$ values as simple horizontal
translations along the time or ${\langle{\rm age}\rangle}$ axis, without
recomputing the models.

Figure 5: Influence of start time on evolution for our reference model with
$\tau_{1}=2$ Gyr, $\tau_{2}=8$ Gyr, $\tau_{*}=1$ Gyr, and $\eta=0.3$. The
solid orange and blue curves in each panel show models with
$t_{\mathrm{start}}$ = 0 and 1 Gyr, respectively. Orange dashed curves, nearly
superposed on the blue curves, show the $t_{\mathrm{start}}=0$ curves shifted
by 1 Gyr to later $t$ or lower ${\langle{\rm age}\rangle}$. The shifted SFR is
renormalized to sum to 1. This 1 Gyr translation captures the impact of
$t_{\mathrm{start}}=1$ almost exactly, despite the effects of light-weighting.

###  3.5 Impact of Parameters on Observables

Figure 6: Influence of star formation history on the evolution of gas phase
([Mg/H], [Mg/Fe]) and stellar (${\langle[{\rm Mg}/{\rm H}]\rangle}$,
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$) abundances and the light-weighted mean
stellar age, ${\langle{\rm age}\rangle}$. The panel layout is the same as in
Figure 4. Solid blue, orange, and green curves show models with $\tau_{2}=8$
Gyr and $\tau_{1}=$ 0, 2, and 4 Gyr, respectively. Dashed and dotted orange
curves show models with $\tau_{1}=2$ Gyr and $\tau_{2}=$ 4 and 12 Gyr,
respectively. An early start of star formation (blue) or a rapid exponential
decline (dashed orange) leads to higher ${\langle{\rm age}\rangle}$ at a given
$t$ and slightly lower ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ at a given
${\langle{\rm age}\rangle}$. All models have $\tau_{*}=1$ Gyr, $\eta=0.3$, and
$t_{\mathrm{start}}=0.5$ Gyr.

We now investigate the impact of model parameters on the observables
${\langle[{\rm Mg}/{\rm H}]\rangle}$, ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$,
and ${\langle{\rm age}\rangle}$ at different epochs. We begin our analysis
with the influence of star formation history on age and abundance evolution
for several SFH models. Our results are presented in Figure 6 with the same
panel layout as Figure 4. However, we have changed the $x$-axis label in the
upper panels from time to $t_{\rm cut}$ to emphasize the role of this
parameter in our predictions: we simply halt model evolution at $t=t_{\rm
cut}$ and read the values of the observables there. Evolving a model to
$t=14\,{\rm Gyr}$ allow us to make predictions for all possible choices of
$t_{\rm cut}$.

Figure 6a shows the selected SFH models, with the reference model represented
by the solid orange curve. These SFHs are identical to those in Figure 3 but
with the addition of a slowly declining SFR ($\tau_{2}=12$ Gyr) and a strictly
declining SFR with high initial star formation ($\tau_{1}=0$ Gyr). As
expected, we find that models with higher SFR at later times, corresponding to
longer $\tau_{1}$ and $\tau_{2}$, have younger, light-weighted ${\langle{\rm
age}\rangle}$ at $z=0$.

Figure 6c shows the time evolution of the gas $[{\rm Mg}/{\rm H}]$. The
different SFH models attain almost the same equilibrium $[{\rm Mg}/{\rm H}]$;
however, the pure exponential model reaches equilibrium sooner because of its
higher early star formation. The time evolution of the gas $[{\rm Mg/Fe}]$,
shown in Figure 6e, has a stronger dependence on the SFH. As already seen in
Figure 4, a rapidly declining SFH leads to lower $[{\rm Mg/Fe}]$ at late
times. Curves for the other models are similar, but the pure exponential again
approaches equilibrium more quickly. Figures 6d and 6f show the time evolution
of the stellar ${\langle[{\rm Mg}/{\rm H}]\rangle}$ and ${\langle[{\rm
Mg}/{\rm Fe}]\rangle}$, respectively. Averaging over the stellar population
greatly reduces the model differences; in particular, while a sharply
declining SFH produces lower $[{\rm Mg/Fe}]$ at late times, there are few
stars formed at these abundances. Overall, the impact of SFH on the predicted
time evolution of ${\langle[{\rm Mg}/{\rm H}]\rangle}$ and ${\langle[{\rm
Mg}/{\rm Fe}]\rangle}$ is remarkably small, at least over the range of
histories considered here.

Figures 6g and 6h plot the mean stellar abundance ratios against ${\langle{\rm
age}\rangle}$, yielding model predictions that can be compared to the SSP
observables discussed in Section 2. The predicted ${\langle[{\rm Mg}/{\rm
H}]\rangle}$-${\langle{\rm age}\rangle}$ tracks are insensitive to the SFH,
except that the exponential SFH, which has its maximum SFR at
$t=t_{\mathrm{start}}$, achieves higher ${\langle[{\rm Mg}/{\rm H}]\rangle}$
at the oldest ages. For ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ the model curves
are more separated as a function of ${\langle{\rm age}\rangle}$ than they were
as a function of $t_{\rm cut}$ because both the $\tau_{1}=0$ and $\tau_{2}=4$
Gyr models have older ages at a given time (Figure 6b). Nonetheless, the
influence of SFH on the ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$-${\langle{\rm
age}\rangle}$ trajectory remains small. In other words, while $\tau_{1}$,
$\tau_{2}$, and $t_{\rm cut}$ all influence ${\langle{\rm age}\rangle}$,
combinations of these parameters that give the same ${\langle{\rm
age}\rangle}$ also give nearly the same ${\langle[{\rm Mg}/{\rm H}]\rangle}$
and ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$.

Figure 7 shows the influence of SFE and outflows on the mean abundance
evolution for models with the reference SFH (solid orange curves in Figure 6a
and 6b). Reducing (increasing) $\tau_{*}$ raises (lowers) ${\langle[{\rm
Mg}/{\rm H}]\rangle}$ at early times and thus in galaxies with high
${\langle{\rm age}\rangle}$ caused by early truncation of star formation.
Conversely, reducing (increasing) $\eta$ raises (lowers) ${\langle[{\rm
Mg}/{\rm H}]\rangle}$ at late times and thus in galaxies with younger
${\langle{\rm age}\rangle}$, while having little impact at early times. This
behavior is consistent with our previous discussion: high SFE (short
$\tau_{*}$) leads to more rapid growth of metallicity, but the late time
equilibrium is governed by the degree of metal loss in outflows. In contrast
to ${\langle[{\rm Mg}/{\rm H}]\rangle}$, the influence of $\tau_{*}$ and
$\eta$ on the ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$-$t_{\rm cut}$ and
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$-${\langle{\rm age}\rangle}$ trajectories
is barely discernible. Although SFE and outflows have a well-known impact on
evolutionary tracks in [$\alpha$/Fe] - [Fe/H] (see, e.g., Figure 3 of Andrews
et al. (2017)), the evolution of [$\alpha$/Fe] with time is determined mainly
by yields and by the SNIa DTD. With these quantities specified, the prediction
of ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ vs. ${\langle{\rm age}\rangle}$ is
remarkably robust within our family of models, with little sensitivity to the
SFH parameters (Figure 6) and even less to $\tau_{*}$ and $\eta$.

Figure 7: Influence of SFE and outflows on mean abundance evolution for models
with the same SFH ($\tau_{1}=2$ Gyr and $\tau_{2}=8$ Gyr, solid orange curve
in Figure 6a). The solid orange curve shows our reference choices of
$\tau_{*}=1$ Gyr and $\eta=0.3$. Solid blue and green curves show models with
high and low SFE, respectively, while dashed and dotted curves show models
with no outflows and high outflows, respectively. Corresponding ($\tau_{*}$,
$\eta$) values are given in the legend. Changing $\tau_{*}$ affects the early
evolution of ${\langle[{\rm Mg}/{\rm H}]\rangle}$, while changing $\eta$
alters the late evolution. Neither parameter has much impact on ${\langle[{\rm
Mg}/{\rm Fe}]\rangle}$ evolution. Figure 8: Influence of the $[{\rm Mg}/{\rm
Fe}]_{\rm cc}$, $t_{\rm start}$, and $\tau_{\rm Ia}$ on models with the
reference values of SFH, SFE, and outflow parameters ($\tau_{1}=2$ Gyr,
$\tau_{2}=8$ Gyr, $\tau_{*}=1$ Gyr and $\eta=0.3$). Solid orange curves,
identical to those in Figure 7, adopt our reference choices of
$t_{\mathrm{d,min}}=0.15$ Gyr, $t_{\mathrm{start}}=0.5$ Gyr, and $y_{\rm
Mg}^{\rm cc}/Z_{{\rm Mg},\odot}=0.973$, corresponding to $[{\rm Mg}/{\rm
Fe}]_{\rm cc}=0.45$ in Equation 4. The solid blue curve shows a model with
boosted $y_{\rm Mg}^{\rm cc}$, corresponding to $[{\rm Mg}/{\rm Fe}]_{\rm
cc}=0.55$. The dotted green and dashed orange curves show models with short
$t_{\mathrm{d,min}}$ and early $t_{\mathrm{start}}$, respectively, with
parameter values given in the legend. Boosting the Mg yield (high $[{\rm
Mg}/{\rm Fe}]_{\rm cc}$) raises ${\langle[{\rm Mg}/{\rm H}]\rangle}$ and
${\langle[{\rm Fe}/{\rm H}]\rangle}$ by 0.1 dex at all times and boosts
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ by 0.1 dex at early times. Changing
$t_{\mathrm{d,min}}$ only affects the Fe abundance, slightly raising
${\langle[{\rm Fe}/{\rm H}]\rangle}$ above the reference value at all times
and altering ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ significantly at early and
intermediate times, but having little impact on the late evolution. Changing
$t_{\mathrm{start}}$ affects the early evolution of ${\langle[{\rm Mg}/{\rm
H}]\rangle}$, ${\langle[{\rm Fe}/{\rm H}]\rangle}$, and ${\langle[{\rm
Mg}/{\rm Fe}]\rangle}$ by shifting the curves in $t$ and ${\langle{\rm
age}\rangle}$ (see Figure 5), but it has little impact on the late evolution.

Figure 8 explores the impact of the Mg yield, the SNIa minimum time delay, and
the star formation start time, with other parameters fixed to their reference
values. As discussed in Section 3.2, we adopt an empirically estimated $y_{\rm
Fe}^{\rm cc}$ and compute $y_{\rm Mg}^{\rm cc}$ based on the plateau value of
$[{\rm Mg}/{\rm Fe}]_{\rm cc}$ taken to represent pure CCSN enrichment. Our
reference choice of $y_{\rm Mg}^{\rm cc}=0.973$ $Z_{{\rm Mg},\odot}$
corresponds to $[{\rm Mg}/{\rm Fe}]_{\rm cc}=0.45$. Boosting $y_{\rm Mg}^{\rm
cc}$ by 0.1 dex, with $[{\rm Mg}/{\rm Fe}]_{\rm cc}=0.55$, raises
${\langle[{\rm Mg}/{\rm H}]\rangle}$ by 0.1 dex at all times. It also boosts
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ by 0.1 dex at early times, as expected.
However, our empirical normalization of $y_{\rm Fe}^{\rm Ia}$ is designed to
produce $[{\rm Mg/Fe}]\approx 0$ at late times for a non-truncated SFH, and
this boosted $y_{\rm Fe}^{\rm Ia}$ largely erases the impact of boosted
$y_{\rm Mg}^{\rm cc}$ on ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ for $t>6$ Gyr
and ${\langle{\rm age}\rangle}<10$ Gyr.

By definition, $t_{\mathrm{d,min}}$ has no impact on ${\langle[{\rm Mg}/{\rm
H}]\rangle}$. However, reducing $t_{\mathrm{d,min}}$ from 0.15 Gyr to 0.05 Gyr
depresses ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ noticeably ($\sim 0.05$ dex)
over a wide range of time and ${\langle{\rm age}\rangle}$. Although we hold
the time-integrated $y_{\rm Fe}^{\rm Ia}$ fixed, a shorter
$t_{\mathrm{d,min}}$ with a $t^{-1.1}$ DTD leads to substantially more rapid
Fe enrichment, and the signature of this faster enrichment persists for
several Gyr. As previously shown (Figure 5), changing $t_{\mathrm{start}}$ is
effectively equivalent to shifting ${\langle[{\rm Mg}/{\rm H}]\rangle}$ and
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ horizontally along the time and
${\langle{\rm age}\rangle}$ axes. Because these trajectories are steep at
early times and flat at late times, the impact on the predicted abundances is
significant at the oldest ages but small at ${\langle{\rm age}\rangle}<10$
Gyr. The evolution of ${\langle[{\rm Fe}/{\rm H}]\rangle}$ follows directly
from ${\langle[{\rm Mg}/{\rm H}]\rangle}$ and ${\langle[{\rm Mg}/{\rm
Fe}]\rangle}$. In this figure only, we show ${\langle[{\rm Fe}/{\rm
H}]\rangle}$ evolution in separate panels. While the tracks resemble those of
${\langle[{\rm Mg}/{\rm H}]\rangle}$, the ${\langle[{\rm Fe}/{\rm H}]\rangle}$
tracks are less flat at late times and young ages because of the Fe
contribution from delayed SNIa.

Figure 9: Abundances as a function of $t_{\rm form}\equiv t_{\rm
obs}-{\langle{\rm age}\rangle}$ for three example models at $z=0$ and $z=0.7$.
In SSP fits, $t_{\rm obs}-\mathrm{age}$ is referred to as $t_{\rm form}$
because it represents the time (relative to the big bang) at which model stars
are assumed to form. As a function of $t_{\rm form}$, the predicted light-
weighted abundances are nearly identical for models evolved to $z=0$ (solid
curves) or to $z=0.7$ (dotted curves).

A similar analysis can be done for galaxies at $z=0.7$ by setting
$t_{\mathrm{obs}}=7$ Gyr in Equations 1 and 12. Figure 9 compares the
${\langle[{\rm Mg}/{\rm H}]\rangle}$ and ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$
trajectories at $z=0$ (solid curves) and $z=0.7$ (dashed curves) for three
example models. For this figure only, we plot the model abundances as a
function of $t_{\rm form}\equiv t_{\rm obs}-{\langle{\rm age}\rangle}$ rather
than ${\langle{\rm age}\rangle}$ itself. For a stellar population that forms
instantaneously at a single epoch, as assumed in an SSP fit, $t_{\rm form}$ is
simply the time at which the stars formed, relative to the big bang at $t=0$.
With linear, mass-weighted age averaging, predictions for $z=0$ and $z=0.7$
galaxies with the same model parameters would be identical as a function of
$t_{\rm form}$ because the evolutionary tracks are the same as a function of
time. Figure 9 shows that predictions are nearly identical as a function of
$t_{\rm form}$ even with logarithmic, light-weighted averaging. Equivalently,
the ${\langle[{\rm Mg}/{\rm H}]\rangle}-{\langle{\rm age}\rangle}$ and
${\langle[{\rm Mg}/{\rm Fe}]\rangle}-{\langle{\rm age}\rangle}$ predictions
for galaxies at $z=0$ and $z=0.7$ are simply shifted in age by $t_{\rm
obs}(z=0)-t_{\rm obs}(z=0.7)=7\,{\rm Gyr}$, if we adopt the same model
parameters.

Figure 10: Summary of parameter impact on reference model abundances at
different age cuts. The red, black, and green, dotted curves show the
reference model with t${}_{\mathrm{start}}=$ 0, 0.5, and 1 Gyr, respectively.
The vertical lines show the possible abundance range from changing one
parameter in the t${}_{\mathrm{start}}=0.5$ Gyr reference model at a given
age. Parameters are represented by different colors and are provided in the
legend. Circles mark the reference values, while squares and triangles
represent lower and higher parameter values, respectively. (One can read the
points as an “arrow” that follows the order of parameter values in the
legend.) While some parameters have a significant impact at some epochs, none
of these parameters boost ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ substantially
above the reference model curve.

Figure 10 summarizes the impact of changing each parameter in the reference
model for galaxies at $z=0$ and $z=0.7$. The red, black, and green dotted
curves show the reference model with $t_{\mathrm{start}}=$ 0, 0.5, and 1 Gyr,
respectively. The vertical lines show the abundance shifts from changing one
parameter in the $t_{\mathrm{start}}=0.5$ Gyr reference model at a given
${\langle{\rm age}\rangle}$. Parameters are represented by different colors
and are provided in the legend. Circles mark the reference values, while the
squares and triangles represent the lowest and highest parameter values,
respectively.

The physics underlying the changes shown in Figure 10 has already been
discussed above, but the presentation here helps compare their magnitude. For
${\langle[{\rm Mg}/{\rm H}]\rangle}$, the values of $\eta$ and the CCSN Mg
yield have the largest impact. The SFE timescale $\tau_{*}$ has an important
impact at the oldest ages but a smaller impact at young ages. In the context
of our comparison to observed abundance trends in Section 4, the most
important implication of Figure 10 is that none of the changes we consider
substantially boost ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ relative to the
reference model. Increasing $[{\rm Mg}/{\rm Fe}]_{\rm cc}$ boosts
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ for the oldest galaxies, but the effect
diminishes at younger ages. Several parameter changes can decrease
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ relative to the reference model: a
shorter $t_{{\rm d,min}}$, a shorter rise time $\tau_{1}$, or a longer decline
time $\tau_{2}$. Shifting $t_{\rm start}$ from 0.5 Gyr to 1 Gyr boosts
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ at a given ${\langle{\rm age}\rangle}$
by roughly the same amount as a 0.1-dex bump in $[{\rm Mg}/{\rm Fe}]_{\rm
cc}$. Increasing $t_{\rm start}$ reduces the ${\langle[{\rm Mg}/{\rm
H}]\rangle}$ of the oldest galaxies, as there is less time for enrichment. For
the $z=0$ galaxies, changing $t_{\rm start}$ from 0.5 Gyr to 1 Gyr is
qualitatively similar to extending the rise time $\tau_{1}$ from 2 Gyr to 4
Gyr. However, for the $z=0.7$ galaxies the longer rise time produces a minimal
change in ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ at a given ${\langle{\rm
age}\rangle}$ while shifting $t_{\rm start}$ has a substantial impact.

##  4 Modeling observed elliptical populations

###  4.1 Representative Models

We use our results from Figure 10 to select representative models that
reproduce the observed ${\langle[{\rm Mg}/{\rm H}]\rangle}$ and ${\langle[{\rm
Mg}/{\rm Fe}]\rangle}$ of low, medium, and high velocity dispersion galaxies
at $z=0$ and galaxies at $z=0.7$. These models are described by star formation
histories with rising timescales of $\tau_{1}=2$ and 4 Gyr, declining
timescales of $\tau_{2}=8$ and 12 Gyr, and a star formation start time of
$t_{\rm start}=0.5$ Gyr. The parameter that has the strongest impact on
${\langle{\rm age}\rangle}$ is $t_{\rm cut}$, the time at which star formation
rapidly truncates. Figure 11 illustrates the relationship between $t_{\rm
cut}$ and the predicted ${\langle{\rm age}\rangle}$ for the representative
models (see Figure 12) at $z=0$ and $z=0.7$. The $z=0.7$ models display nearly
identical tracks in ${\langle{\rm age}\rangle}$-$t_{\rm cut}$. However, there
are differences as large as $\sim 1$ Gyr in ${\langle{\rm age}\rangle}$
between the $z=0$ models for a given $t_{\rm cut}$ when $t_{\rm cut}\gtrsim
4.5$ Gyr. Here, models having a higher star formation rate at late times
(longer $\tau_{1}$ and/or longer $\tau_{2}$) have younger ${\langle{\rm
age}\rangle}$ at fixed $t_{\rm cut}$. The horizontal lines mark the observed
${\langle{\rm age}\rangle}$ for the lowest (youngest) and highest (oldest)
velocity dispersion bins at each redshift (see Figure 12). The intersections
of the horizontal lines and model curves indicate the truncation time required
for each model to obtain the observed ${\langle{\rm age}\rangle}$. More
massive, high-$\sigma$ ellipticals are older than less massive, low-$\sigma$
ellipticals and therefore cease star formation at an earlier $t_{\rm cut}$.

Figure 12 compares our representative models to the observed ages and
abundances of galaxies stacked in bins of velocity dispersion, $\sigma$.
Figure 10 shows that once the yields are fixed, there is not much leverage to
boost ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$. With empirical yields described
in Section 3.2, there are no combinations of $\tau_{1}$, $\tau_{2}$,
$\tau_{*}$, and $\eta$ that can reproduce the observations. Because our models
generally underproduce the observed ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$, we
adopt $[{\rm Mg}/{\rm Fe}]_{\rm cc}=0.55$ as our reference value in this
section rather than 0.45 as in Section 3. While Milky Way stellar observations
generally favor a lower $[{\rm Mg/Fe}]$ plateau (see WGJT24), the value of
$[{\rm Mg}/{\rm Fe}]_{\rm cc}=0.55$ remains possible, especially if one allows
for a significant SNIa Fe contribution to stars with $-1.5<[{\rm Fe/H}]<-0.5$
as argued by Maoz & Graur (2017), Conroy et al. (2022), and Mason et al.
(2023). The dotted gray curves in Figure 12 show the reference model that is
the central curve from Figure 10; the most important difference from the
models in this section is the lower $[{\rm Mg}/{\rm Fe}]_{\rm cc}=0.45$.

Figure 11: Relation between ${\langle{\rm age}\rangle}$ and $t_{\mathrm{cut}}$
for representative models (see Figure 12) at redshift $z=0$ and $z=0.7$. The
$z=0.7$ models display overlapping tracks in ${\langle{\rm
age}\rangle}$-$t_{\rm cut}$ while the $z=0$ models have distinguishable tracks
for $t_{\rm cut}\gtrsim 4$ Gyr. However, the pink and orange models at $z=0$
are nearly identical. The horizontal lines mark the observed ${\langle{\rm
age}\rangle}$ for the lowest (youngest) and highest (oldest) velocity
dispersion bins at each redshift. The intersection of the horizontal lines and
model curves indicate the value of $t_{\rm cut}$ required to obtain the
observed ${\langle{\rm age}\rangle}$. Figure 12: Comparison between
representative models and the observed ages and abundances of galaxies stacked
in bins of velocity dispersion, $\sigma$. Observed properties are represented
by blue points and contain SDSS (Conroy et al., 2014) and LEGA-C (Beverage et
al., 2023) data for galaxies at $z=0$ and $z=0.7$, respectively. The points
increase in $\sigma$ moving from light to dark blue. The solid curves show
representative models with parameters chosen from values presented in Figure
10 to reproduce the observed ${\langle[{\rm Mg}/{\rm H}]\rangle}$ and
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ of low, medium, and high velocity
dispersion galaxies at $z=0$ and galaxies at $z=0.7$. All of these models have
$t_{\rm start}=0.5$ Gyr. Reproducing the ${\langle{\rm age}\rangle}$ of
high-$\sigma$ ellipticals requires $\eta\approx 0$ and high SFE, while
low-$\sigma$ ellipticals can accommodate moderate outflows and lower SFE. At
$z=0.7$, the high ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ requires $\eta\approx
0$ and high SFE for all $\sigma$-bins. All representative models have a high
Mg yield, with $[{\rm Mg}/{\rm Fe}]_{\rm cc}=0.55$, but even so, they under-
predict the observed ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ at $z=0$ for all
but the lowest-$\sigma$ galaxies. The dotted gray curve shows the reference
model displayed by the solid orange curves in Figures 6-8 and by the central
dotted curves in Figure 10. The reference model has $[{\rm Mg}/{\rm Fe}]_{\rm
cc}=0.45$, which reduces the predicted ${\langle[{\rm Mg}/{\rm H}]\rangle}$
and ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ and worsens the agreement with
observations.

The blue points in Figure 12a show the ${\langle[{\rm Mg}/{\rm H}]\rangle}$
and ${\langle{\rm age}\rangle}$ of the observed SDSS galaxy bins at
$z=0$.222More precisely, values shown by blue points are the SSP fits, which
we refer to with the same notation ${\langle{\rm age}\rangle}$, ${\langle[{\rm
Mg}/{\rm H}]\rangle}$, and ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ that we use
for light-weighted model averages. Velocity dispersion increases when moving
from the light to dark blue points, corresponding to ${\langle[{\rm Mg}/{\rm
H}]\rangle}$ and ${\langle{\rm age}\rangle}$ increasing with increasing
$\sigma$ (Conroy et al., 2014). The ${\langle[{\rm Mg}/{\rm H}]\rangle}$ and
${\langle{\rm age}\rangle}$ of the low-$\sigma$ bins are best represented by a
slowly rising and declining SFH ($\tau_{1}=4$ and $\tau_{2}=12$ Gyr) with mild
outflows ($\eta=0.3$) and SFE typical of local molecular gas disks
($\tau_{*}=2$ Gyr) (green curve). Longer $\tau_{1}$, $\tau_{2}$, and
$\tau_{*}$ allow for higher late-time star formation that results in young
${\langle{\rm age}\rangle}$. Mild outflows balance the late Mg production and
push ${\langle[{\rm Mg}/{\rm H}]\rangle}$ closer to Solar.

The representative model for the mid-$\sigma$ bins (pink curve) is similar to
the low-$\sigma$ model but with a more rapid SFR rise ($\tau_{1}=2$ Gyr) and
higher SFE ($\tau_{*}=1$ Gyr) that produces more stars at earlier times. This
increases the ${\langle{\rm age}\rangle}$ and raises the ${\langle[{\rm
Mg}/{\rm H}]\rangle}$ by producing more Mg earlier and quenching before all
the SNIa explode. The high-$\sigma$ bins (blue curve) are best represented by
a rapidly rising and declining SFH ($\tau_{1}=2\,{\rm Gyr}$ and
$\tau_{2}=8\,{\rm Gyr}$) with high SFE ($\tau_{*}=0.5$ Gyr) and no outflows
($\eta=0$). Shorter $\tau_{1}$, $\tau_{2}$, and $\tau_{*}$ result in the
majority of stars being formed at early times, increasing the ${\langle{\rm
age}\rangle}$ and producing significant amounts of Mg before early quenching.
No outflows results in all the Mg being retained, driving ${\langle[{\rm
Mg}/{\rm H}]\rangle}$ to $\sim 0.25$ dex above Solar. The high SFE is needed
to achieve the high observed ${\langle[{\rm Mg}/{\rm H}]\rangle}$ at such old
ages.

The parameter choices for these three representative models are not unique,
but they reflect our expectation that the oldest and highest-$\sigma$ galaxies
have high star formation efficiency and rapid early star formation, and they
ascribe the trend of ${\langle[{\rm Mg}/{\rm H}]\rangle}$ with $\sigma$ to the
lower outflow efficiency in deeper potential wells. The mid-$\sigma$ model is
similar to the reference model from Section 3 (gray dotted curve), but it
predicts higher ${\langle[{\rm Mg}/{\rm H}]\rangle}$ because $[{\rm Mg}/{\rm
Fe}]_{\rm cc}=0.55$ implies a higher Mg yield $y_{\rm Mg}^{\rm cc}$.

The blue points in Figure 12c show the ${\langle[{\rm Mg}/{\rm H}]\rangle}$
and ${\langle{\rm age}\rangle}$ of the observed LEGA-C galaxy bins at $z=0.7$,
with velocity dispersion increasing from the light blue to the dark blue
points. At this redshift, ${\langle[{\rm Mg}/{\rm H}]\rangle}$ is nearly
constant across all bins of $\sigma$. The blue and pink representative models
at $z=0.7$ have the same parameters as the blue and pink representative models
at $z=0$. The orange model is similar to the blue model but with a longer rise
timescale $\tau_{1}=4$ Gyr, chosen because it slightly boosts ${\langle[{\rm
Mg}/{\rm Fe}]\rangle}$ for the youngest galaxies at this redshift. As shown in
Figure 10, ${\langle[{\rm Mg}/{\rm H}]\rangle}$ at $z=0.7$ is insensitive to
$\tau_{1}$ for all ${\langle{\rm age}\rangle}$, so the blue and orange models
have nearly identical ${\langle[{\rm Mg}/{\rm H}]\rangle}$-${\langle{\rm
age}\rangle}$ trajectories. The blue and orange models retain all of the Mg
they produce ($\eta=0$) and can match the high observed ${\langle[{\rm
Mg}/{\rm H}]\rangle}$ at $z=0.7$. Outflows with $\eta=0.3$ cause the pink
model to undershoot the observed ${\langle[{\rm Mg}/{\rm H}]\rangle}$ by $\sim
0.2$ dex.

Figures Figure 12b and Figure 12d compare the representative models to the
observed ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ and ${\langle{\rm age}\rangle}$
of the SDSS and LEGA-C bins, respectively. At $z=0$, the green model predicts
the observed ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ of the low-$\sigma$
galaxies and nearly reproduces the observed values at older ages. The pink and
blue models, which reproduce the ${\langle[{\rm Mg}/{\rm H}]\rangle}$ of mid-
and high-$\sigma$ galaxies, respectively, underpredict the observed
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ for all ages. At $z=0.7$, the blue,
orange, and pink models make nearly identical predictions and are reasonably
consistent with the observed ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ of the
LEGA-C galaxies within the estimated uncertainties.

Figure 12 demonstrates the difficulty of reproducing observed abundance trends
in ellipticals with smooth star formation histories and yields that are
empirically motivated by abundance trends in the Milky Way. Our reference
model with $[{\rm Mg}/{\rm Fe}]_{\rm cc}=0.45$ underpredicts the observed
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ values by $\sim 0.08$ dex at $z=0$ and
$\sim 0.05$ dex at $z=0.7$. For $[{\rm Mg}/{\rm Fe}]_{\rm cc}=0.55$, we find
models that come close to the observed ${\langle[{\rm Mg}/{\rm
H}]\rangle}$-${\langle{\rm age}\rangle}$ and ${\langle[{\rm Mg}/{\rm
Fe}]\rangle}$-${\langle{\rm age}\rangle}$ trends at $z=0.7$ (blue and orange
curves in panels c and d) and underpredict the observed ${\langle[{\rm
Mg}/{\rm Fe}]\rangle}$ by $\sim 0.02-0.05$ dex for low-$\sigma$ ellipticals at
$z=0$ (green curves in panels a and b). For the higher-$\sigma$ galaxies,
achieving super-solar ${\langle[{\rm Mg}/{\rm H}]\rangle}$ at older ages
requires higher SFE, which leads to a more rapid drop in ${\langle[{\rm
Mg}/{\rm Fe}]\rangle}$. For the mid-$\sigma$ galaxies, the gap between our
representative model (pink curve) and the observed ${\langle[{\rm Mg}/{\rm
Fe}]\rangle}$ is $\sim 0.06-0.08$ dex, and for the oldest, high-$\sigma$
ellipticals (blue curve) it is $\sim 0.07$ dex. We have made parameter choices
that tend to maximize ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ at a given
${\langle{\rm age}\rangle}$ subject to the constraints imposed by
${\langle[{\rm Mg}/{\rm H}]\rangle}$, and one can see from Figure 10 that many
alternative choices would worsen the agreement with the observed abundance
trends.

Gaps of $\sim 0.05$ dex are small enough that they might be explained by
systematic errors in the SSP fits themselves, e.g., because of imperfections
in the stellar libraries used for the spectral synthesis. However, since we
have already pushed our parameter choices to reduce discrepancies to this
level, we will examine possible physical resolutions, arising from lower SNIa
yields (Section 4.3) or from terminating bursts of star formation (Section
4.4). Even if systematic errors can adequately explain the gaps between our
model predictions and the binned-average data, other mechanisms may be needed
to explain individual galaxies that lie above the mean. Before considering
these mechanisms, we re-examine the relationship between SSP fits and full
evolutionary models, now guided by the models that come closest to matching
the data. We return to a broader discussion of our results, including the
choice of CCSN yields and the effects of the stellar IMF, in Section 5.

###  4.2 Test of SSP Fits

Figure 13: Mock-recovery test results comparing the light-weighted model
predictions (lines) to the best-fit SSP-equivalent quantities derived from
mock CSP spectra (points). CSPs were generated using the pink and blue GCE
models from Figure 12 for $z=0$ (top row) and $z=0.7$ (bottom row). At various
cutoff times, mock spectra with 20 noise realizations were fit using the
method described in Section 2. Error bars represent the 16th and 84th
percentiles of the best-fit SSP-equivalent quantities for the 20 realizations
of each CSP. This test shows that SSP-equivalent quantities can be well-
approximated by light-weighted values, to within 0.02 dex for ${\langle[{\rm
Mg}/{\rm H}]\rangle}$ and ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$, with a small
bias in ${\langle{\rm age}\rangle}$ for the oldest models.

Our two-population tests in Appendix C, summarized in Section 2, show that an
SSP fit recovers the light-weighted, logarithmically averaged age, $[{\rm
Mg}/{\rm H}]$, and $[{\rm Mg/Fe}]$ of a composite stellar population. Now that
we have guidance on the model parameters relevant to observed ellipticals, we
test the consistency between SSP fits and light-weighted averages for models
with extended SFHs and continuously evolving metal content.

We construct composite stellar population spectra corresponding to the pink
and blue GCE models presented in Figure 12. To do so we compute SSP spectra at
each time step based on the evolving gas-phase [Mg/H] and [Mg/Fe] and the
look-back time $t_{\mathrm{obs}}-t$. For each choice of truncation time
$t_{\rm cut}$, we combine all SSPs up to $t_{\rm cut}$, weighting each by the
mass fraction formed in each time step. The appropriate light weighting at
each wavelength follows directly from the spectral synthesis.

We simulate observational data from these CSP spectra by adding Gaussian noise
corresponding to S/N = 40 Å-1, which is typical of individual LEGA-C and SDSS
spectra. These mock observations are then fit using the same full-spectrum SSP
fitting method applied to the observations as described in Section 2.

In Figure 13, smooth curves show the models’ light-weighted ${\langle{\rm
age}\rangle}$, ${\langle[{\rm Mg}/{\rm H}]\rangle}$, and ${\langle[{\rm
Mg}/{\rm Fe}]\rangle}$, as a function of $t_{\rm cut}$, with
$t_{\mathrm{obs}}=14$ Gyr ($z=0$, top row) and $t_{\mathrm{obs}}=7$ Gyr
($z=0.7$, bottom row). Points show the mean SSP-fitting results from 20 mock
spectra noise realizations at selected values of $t_{\rm cut}$. Error bars
show 16th-84th percentile range of the 20 realizations. The uncertainty in the
mean prediction is smaller by $\approx\sqrt{20}$, and the statistical errors
on the observational data points in Figure 12 are often smaller than this
because more than 20 galaxies contribute to an individual $\sigma$-bin.

Figure 13 demonstrates good agreement between light-weighted model quantities
and SSP fits for realistic SFH and enrichment histories, as expected from the
two-population tests. The model averages typically overpredict the SSP
${\langle[{\rm Mg}/{\rm H}]\rangle}$, but by $\lesssim 0.02$ dex. Agreement in
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ is similarly at the 0.01-0.02 dex level
in most cases. The logarithmically averaged ${\langle{\rm age}\rangle}$
typically agrees with the SSP values to $\lesssim 0.5$ Gyr. However, for the
oldest models (earliest $t_{\rm cut}$) at $z=0$, the difference rises to 1-1.5
Gyr. Accounting for this effect would noticeably improve the agreement between
models and data for the most massive $z=0$ galaxies, as one can see by
mentally translating the blue and pink curves in Figure 12 leftward by 1-1.5
Gyr. The differences for the oldest galaxies may be a consequence of
incorporating CSPs at $[{\rm Fe/H}]<-0.7$ that rely on empirical stellar
libraries in a poorly constrained metallicity regime (Villaume et al., 2017).
More accurate libraries in this regime might reduce the difference between
light-weighted evolutionary histories and SSP fits.

###  4.3 SNIa Yields

Figure 14: Influence of SNIa Fe yields. (Top panels) The solid blue and orange
curves compare models with the reference SNIa Fe yield
($y_{\mathrm{Fe}}^{\mathrm{Ia}}/y_{\mathrm{Fe}}^{\mathrm{cc}}=1.65$) and a
reduced SNIa Fe yield
($y_{\mathrm{Fe}}^{\mathrm{Ia}}/y_{\mathrm{Fe}}^{\mathrm{cc}}=1.00$),
respectively. Shortening $t_{\mathrm{d,min}}$ for the reduced yield model
(dashed orange curve) results in lower ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$
at early times but similar ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ at the ages
of the observed galaxies. (Bottom panels) The blue and green curves compare
models with a $t^{-1.1}$ power-law DTD and a single exponential DTD,
respectively. The solid green curve has $\tau_{\rm Ia}=1.5$ Gyr while the
dotted green curve has $\tau_{\rm Ia}=3$ Gyr. All models are based on the
representative model shown by the blue curve in Figure 12 ($\tau_{1}=2$ Gyr,
$\tau_{2}=8$ Gyr, $\tau_{*}=0.5$ Gyr and $\eta=0$), but with the reference
[Mg/Fe]${}_{\mathrm{cc}}=0.45$. A lower SNIa Fe yield or preferential loss of
SNIa Fe in galactic outflows would enable our elliptical galaxy models to
reproduce observed ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ ratios.

Although we can obtain marginal agreement with the elliptical galaxy
observables by adopting a high $[{\rm Mg}/{\rm Fe}]_{\rm cc}=0.55$, most Milky
Way observations favor $[{\rm Mg}/{\rm Fe}]_{\rm cc}\approx 0.45$ or (in the
case of APOGEE) even lower. Another parameter that we have not yet explored is
the SNIa Fe yield. Our normalization condition (Equation 5) is chosen to
produce near-solar $[{\rm Mg/Fe}]$ in the ISM at late times, but galaxies that
become massive ellipticals could have a different $y_{\rm Fe}^{\rm Ia}$ from
disk galaxies like the Milky Way.

For $[{\rm Mg}/{\rm Fe}]_{\rm cc}=0.45$, Equation 5 gives $y_{\rm Fe}^{\rm
Ia}/y_{\rm Fe}^{\rm cc}=1.65$. The blue solid curves in Figure 14 show
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ vs. ${\langle{\rm age}\rangle}$ in a
model with these yields and other parameters ($\tau_{1}=2$ Gyr, $\tau_{2}=8$
Gyr, $\tau_{*}=0.5$ Gyr and $\eta=0$) matched to the blue-curve model of
Figure 12. The solid orange curves show the impact of reducing the SNIa yield
to $y_{\rm Fe}^{\rm Ia}=y_{\rm Fe}^{\rm cc}$, with other yields unchanged.
This reduced $y_{\rm Fe}^{\rm Ia}$ model is remarkably successful at
reproducing the observed ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$-${\langle{\rm
age}\rangle}$ relations at both $z=0$ and $z=0.7$. If we also reduce $t_{{\rm
d,min}}$ from 0.15 Gyr to 0.05 Gyr (dashed orange curves), the agreement with
the data is somewhat worse, but still better than that of the original
reference model.

Reducing $y_{\rm Fe}^{\rm Ia}$ by a factor of 1.65 is equivalent to reducing
the normalization of the SNIa DTD in ellipticals by the same factor. However,
observations suggest that the SNIa rate in galaxy clusters, which have a high
fraction of elliptical and lenticular galaxies, is higher than the rate in
star-forming field galaxies, not lower (Maoz & Graur, 2017). Because the
stellar populations of ellipticals are old, observed rates probe the late
phases of the delay time distribution, so this inference of a higher time-
integrated rate depends on assuming that the DTD form is the same in cluster
and field galaxies. It is also possible that the
$y_{\mathrm{Fe}}^{\mathrm{Ia}}/y_{\mathrm{Fe}}^{\mathrm{cc}}$ in ellipticals
is the same as that in our reference model, or higher, but that ellipticals
preferentially lose SNIa products in galactic winds relative to CCSN products,
even at epochs when they are actively forming stars and enriching their ISM.
Such preferential ejection could arise if CCSN ejecta are trapped in dense
star-forming gas while SNIa ejecta are released in lower-density gas that is
more easily vented from the galaxy. Nonetheless, while these scenarios are
possible, there seems to be no good empirical or theoretical argument for
reduced $y_{\mathrm{Fe}}^{\mathrm{Ia}}/y_{\mathrm{Fe}}^{\mathrm{cc}}$ in
elliptical galaxies other than the high ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$
ratios.

Dubay et al. (2024) find that GCE models are more successful at reproducing
Milky Way observations if the SNIa DTD has fewer short-delay supernovae than
predicted by the $t^{-1.1}$ power-law advocated by, e.g., Maoz & Graur (2017)
and Wiseman et al. (2021). Strolger et al. (2020) argue that an exponential-
like DTD fits observed SNIa rates as well as the power-law model (see also
Palicio et al. 2024). Dubay et al. show (in their Appendix B) that an
exponential DTD with a timescale $\tau_{\rm Ia}\approx 1.5\,{\rm Gyr}$ is
similar to the theoretically motivated single-degenerate DTD of Greggio
(2005), and a number of GCE studies have adopted this exponential form for the
DTD (e.g., Schönrich & Binney 2009; Andrews et al. 2017; WAF). Relative to the
$t^{-1.1}$ power-law with $t_{{\rm d,min}}=0.15\,{\rm Gyr}$, the $\tau_{\rm
Ia}=1.5\,{\rm Gyr}$ exponential DTD has a similar median delay time ($\approx
1\,{\rm Gyr}$) but fewer SNIa at $t\ll 1\,{\rm Gyr}$. The solid green curve in
Figure 14 adopts the same parameters and $y_{\rm Fe}^{\rm Ia}$ as the blue
curve but employs a single exponential DTD with $\tau_{\rm Ia}=1.5\,{\rm
Gyr}$. For the oldest ages, the predicted ${\langle[{\rm Mg}/{\rm
Fe}]\rangle}$ is slightly higher than that of the reference model, which uses
a sum of two exponentials to approximate a $t^{-1.1}$ power-law. However, at
younger ages the predicted ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ is slightly
lower because the exponential DTD has more intermediate-time SNIa than the
long-tailed, power-law DTD. Thus, simply changing the form of the DTD to a
single exponential does not improve agreement with the elliptical galaxy
observations.

Dubay et al. (2024) also consider an exponential DTD with $\tau_{\rm
Ia}=3\,{\rm Gyr}$. We have repeated our calculations for this DTD (dotted
green curve), and in this case, we get much better agreement with the data at
$z=0.7$ and $z=0$. Like the findings of Dubay et al. (2024), Figure 14
highlights the importance of pinning down both the form and the normalization
of the SNIa DTD, and testing its universality across different types of
galaxies.

###  4.4 Terminating Bursts of Star Formation

Our models assume flexible but smooth star formation histories. Bursts of star
formation produce temporary boosts in [Mg/Fe] because the rate of
$\alpha$-element production from massive stars increases during the burst
(WAF; Johnson & Weinberg 2020). If a burst occurs at the beginning or middle
of the SFH, then the Fe from SNIa progenitors formed in the burst will
eventually be incorporated into stars, so we expect the impact on
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ to be small. A burst that terminates the
galaxy’s SFH is different because the newly produced Mg can be incorporated
into stars during the burst (self-enrichment), but the SNIa Fe production
occurs after star formation has ceased. It has long been known that some
galaxies, variously referred to as “E+A”, “K+A”, or “post-starburst”, undergo
bursts of star formation before abruptly transitioning to quiescence (Dressler
& Gunn, 1983; Zabludoff et al., 1996; Balogh et al., 1999; Leung et al., 2023;
Park et al., 2023). However, it is unclear whether most early type galaxies
end their SFH with a final burst.

To estimate the impact of a terminating burst, we compute the mean metallicity
of stars formed in the burst using a “leaky box” chemical evolution model
(Hartwick, 1976), which extends “closed box” evolution (Talbot & Arnett, 1971)
by allowing outflow. We implicitly assume that the burst duration is long
compared to the $\sim 30\,{\rm Myr}$ lifetime of massive stars, so that newly
produced CCSN elements are incorporated into later generations within the
burst, but short compared to the $\sim 1\,{\rm Gyr}$ timescale over which
substantial SNIa Fe production would occur. A terminating burst could be
triggered by a galaxy merger that brings a substantial amount of fresh gas,
which dilutes the metallicity of the ISM but increases the fuel available to
make stars. Our use of the leaky box model implies that any further gas
accretion during the burst is negligible compared to the initial supply. While
our calculation is idealized, it captures the distinctive ways that a
terminating burst affects the mean stellar abundances of a galaxy and gives an
approximate value of the quantitative impact on observables.

We give details of this calculation in Appendix B. The parameters of a burst
in our model are $f_{g}$, the ratio of gas remaining at the end of the burst
to gas at the start of the burst, and $F_{D}$, the ratio of accreted gas to
gas present in the galaxy at the end of smooth evolution. For representative
values of $f_{g}$ we consider $f_{g}=0.3$ and $f_{g}=0.05$, corresponding to
bursts that consume 70% and 95% of the initial gas supply. We assume that
outflows during the burst are negligible, regardless of the value of $\eta$
assumed during smooth evolution. This choice maximizes the impact of the burst
because all metals formed during the burst are retained.

For the dilution factor we consider $F_{D}=0$, corresponding to a burst that
only consumes the pre-existing gas in the ISM, and $F_{D}=3$, corresponding to
an accretion event that increases the gas supply by a factor $1+F_{D}=4$. For
example, $F_{D}=3$ could describe a galaxy with $10^{11}M_{\odot}$ of stars
and $10^{10}M_{\odot}$ of remaining gas that accretes a gas rich companion
with $3\times 10^{10}M_{\odot}$ of gas, triggering a burst. We ignore the
impact of accreted stars. We assume that accretion lowers the ISM abundance of
both Mg and Fe by a factor $(1+F_{D})$, which is a reasonable approximation
for a minor merger in which the accreted system has much lower metallicity
than the main galaxy. The mass and abundances of stars formed during the burst
are computed from Equation (B8) and Equation (B9), respectively.

The impact of a burst on observables is enhanced by light-weighting, as the
stellar population formed during the burst is younger and brighter than the
stellar populations formed during smooth evolution. We assume that the burst
duration is long enough for self-enrichment but short compared to other GCE
timescales. We therefore assign all stellar mass from the burst to the
timestep ($\Delta t=0.02\,{\rm Gyr}$) that immediately follows $t_{\rm
burst}=t_{\rm cut}$. Our procedure should be accurate for burst durations up
to $\sim 0.3\,{\rm Gyr}$, but for longer burst durations one could not ignore
SNIa enrichment.

Figure 15: Impact of a burst of star formation at $t_{\mathrm{burst}}=t_{\rm
cut}=8$ Gyr on age and abundance evolution. The reference model ($\tau_{1}=2$
Gyr, $\tau_{2}=8$ Gyr, $\tau_{*}=1$ Gyr, $\eta=0.3$, solid orange curve in
Figure 4) before the burst is shown in gray. The colored points in panel (a)
represent the fraction of total stars formed in the burst, ranging from $\sim
10\%$ ($F_{D}=0$, $f_{g}=0.3$) to $\sim 40\%$ ($F_{D}=3$, $f_{g}=0.05$); the
colored curves differ only in normalization, such that the area under the
curve plus the final point sum to one. Points in the other panels show the
post-burst age or abundance. The colors indicate different starburst models,
with corresponding ($F_{D}$, $f_{g}$) values given in the legend. A burst of
star formation decreases ${\langle{\rm age}\rangle}$ while increasing
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$.

Figure 15 shows the effects of a final burst of star formation at
$t_{\mathrm{burst}}=8$ Gyr on the reference model for different burst
parameters. The curves in Figure 15a show the smooth SFH as before, while
points at 8 Gyr show the fraction of the galaxy’s stars that form in the final
burst. Burst models that have accretion of pristine gas ($F_{D}=3$) and use a
large fraction of gas in the burst ($f_{g}=0.05$) form a higher fraction of
stars ($\sim$40%) during the burst. In comparison, burst models with no gas
accretion ($F_{D}=0$) and a lower fraction of gas used in the burst
($f_{g}=0.3$) form $\sim$10% of their total stars during the burst. Light-
weighting makes the final burst have a significant impact on ${\langle{\rm
age}\rangle}$, especially for late $t_{\mathrm{burst}}$. The gray curve in
Figure 15b shows age evolution for the reference model before the burst, with
the colored points indicating the post-burst values. A final burst of star
formation lowers ${\langle{\rm age}\rangle}$ by as much as 1 Gyr from the pre-
burst value, with models that form a higher fraction of stars during the burst
having younger, final ${\langle{\rm age}\rangle}$.

Points in Figures 15c and 15e show the mean $[{\rm Mg}/{\rm H}]$ and $[{\rm
Mg/Fe}]$ of stars formed during the burst, computed as described in Appendix
B. If there is no accretion ($F_{D}=0$), $[{\rm Mg}/{\rm H}]$ increases due to
self-enrichment from CCSN producing Mg during the burst. Models that convert a
higher fraction of gas into stars during the burst ($f_{g}=0.05$) produce more
CCSN and therefore show a larger increase in $[{\rm Mg}/{\rm H}]$. If the ISM
is diluted by pristine gas immediately before the burst ($F_{d}=3$), then the
$[{\rm Mg}/{\rm H}]$ of burst stars is lower, and it even lies below the
$[{\rm Mg}/{\rm H}]$ at the end of the smooth evolution if only a small
fraction of gas is converted into stars ($f_{g}=0.3$). show the mean $[{\rm
Mg}/{\rm H}]$ and $[{\rm Mg/Fe}]$ of stars formed during the burst, computed
as described in Appendix B. If there is no accretion ($F_{D}=0$), $[{\rm
Mg}/{\rm H}]$ increases due to self-enrichment from CCSN producing Mg during
the burst. Models that convert a higher fraction of gas into stars during the
burst ($f_{g}=0.05$) produce more CCSN and therefore show a larger increase in
$[{\rm Mg}/{\rm H}]$. If the ISM is diluted by pristine gas immediately before
the burst ($F_{d}=3$), then the $[{\rm Mg}/{\rm H}]$ of burst stars is lower,
and it even lies below the $[{\rm Mg}/{\rm H}]$ at the end of the smooth
evolution if only a small fraction of gas is converted into stars
($f_{g}=0.3$).

We assume that the final bursts are long enough to incorporate CCSN products
produced by stars formed in the burst, but short enough that there is no Fe
enrichment from SNIa. As a result, all burst models increase $[{\rm Mg/Fe}]$.
Models with larger burst fractions show larger $[{\rm Mg/Fe}]$ increases, so
in contrast to $[{\rm Mg}/{\rm H}]$ the $F_{D}=3$ models have the largest
boost of $[{\rm Mg/Fe}]$. Figures 15d and 15f show the final light-weighted
${\langle[{\rm Mg}/{\rm H}]\rangle}$ and ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$
of the entire galaxy. Relative to Figures 15c and 15e, abundance changes are
reduced because the burst comprises only a fraction of the galaxy’s stars.
Figures 15g and 15h plot ${\langle[{\rm Mg}/{\rm H}]\rangle}$ and
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ against ${\langle{\rm age}\rangle}$.
Overall, adding a final burst of star formation increases ${\langle[{\rm
Mg}/{\rm Fe}]\rangle}$ and decreases ${\langle{\rm age}\rangle}$, with both
changes moving the models in the right direction to reproduce the observed
data.

Figure 16: Comparison between representative starburst models and observed
ages and abundances for galaxies stacked in bins of velocity dispersion. The
gray curve is similar to the pink, representative model shown in Figure 12
($\tau_{1}=2$ Gyr, $\tau_{2}=12$ Gyr, $\tau_{*}=1$ Gyr), but with a lower Mg
yield ($[{\rm Mg}/{\rm Fe}]_{\rm cc}=0.45$) and with no outflows ($\eta=0$).
The colored lines show the effects of a burst of star formation occurring at
$t_{\mathrm{burst}}=5.30,$ 7.70, and 10.52 Gyr at $z=0$ and
$t_{\mathrm{burst}}=2.02,$ 3.50, and 5.16 Gyr at $z=0.7$, with endpoints
marking the post-burst value. The colors indicate different starburst models,
with corresponding ($F_{D}$, $f_{g}$) values given in the legend.
Incorporating a substantial burst of star formation right before quenching can
reproduce the observed, elevated ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ ratios.
The bursts typically boost ${\langle[{\rm Mg}/{\rm H}]\rangle}$ as well, but
the effect is reduced if the burst is fueled by freshly accreted, low-
metallicity gas (high $F_{D}$).

Figure 16 provides a comparison between representative starburst models and
the observed ages and abundances for galaxies stacked in bins of velocity
dispersion, $\sigma$. We introduce bursts of star formation at
$t_{\mathrm{burst}}=5.30,$ 7.70, and 10.52 Gyr for the $z=0$ comparison and
$t_{\mathrm{burst}}=2.04,$ 3.50, and 5.16 Gyr for the $z=0.7$ comparison. As a
baseline smooth model, we adopt parameters similar to those of the central
model (pink curve) shown in Figure 12 ($\tau_{1}=2$ Gyr, $\tau_{2}=12$ Gyr,
$\tau_{*}=1$ Gyr), but we adopt the lower Mg yield favored by Milky Way disk
data ($[{\rm Mg}/{\rm Fe}]_{\rm cc}=0.45$), and we set outflows to $\eta=0$ to
compensate for this lower yield. The colored lines show the effects of the
burst models at each $t_{\mathrm{burst}}$, with the endpoints marking the
post-burst value.

Figure 16 shows that adding a substantial final burst of star formation to a
reasonable smooth star formation history can bridge the gap between predicted
and observed ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ ratios, bringing models
with $[{\rm Mg}/{\rm Fe}]_{\rm cc}=0.45$ within the range of the data. The
effects of the burst on ${\langle{\rm age}\rangle}$ and on ${\langle[{\rm
Mg}/{\rm Fe}]\rangle}$ are about equally important. While a burst that
consumes the existing gas (e.g., $F_{D}=0$, $f_{g}=0.05$) can resolve the
tension with data at $z=0.7$, fully resolving the tension at $z=0$ requires a
larger burst with accreted gas. In these $F_{D}=3$ models, the ${\langle[{\rm
Mg}/{\rm H}]\rangle}$ at a given age is similar to or even below that of the
smooth model, as dilution outweighs the effect of self-enrichment. The models
shown in Figure 16 fall short of the observed ${\langle[{\rm Mg}/{\rm
H}]\rangle}$ of the high-$\sigma$ galaxies. However, this discrepancy could be
resolved by adopting a shorter $\tau_{*}$ and $\tau_{2}$ for this population
and, if necessary, boosting the Mg yield above the level assumed here (see
Figure 12).

If self-enrichment during terminating bursts is indeed the explanation for the
high observed ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ relative to our model
predictions, it implies that these large final bursts of star formation are
the norm in elliptical galaxy evolution, not the exception. Population
synthesis models that isolate late star formation can test this conjecture.
Furthermore, since not all galaxies will experience bursts of the same
magnitude, we expect correlations between the fraction of stars formed
immediately before quiescence and the mean abundances of the stellar
population.

To illustrate this point, we construct models with varying values of $F_{D}$
and $f_{g}$, and $t_{\mathrm{burst}}$ chosen to produce a specified
${\langle{\rm age}\rangle}$, so we can predict correlations that might be seen
within ellipticals of a given age or velocity dispersion. The smooth models
are the blue, pink, and green representative models presented in Figure 12 but
with lower Mg yields, $[{\rm Mg}/{\rm Fe}]_{\rm cc}=0.45$, and with outflows
set to $\eta=0$ for the pink model to boost ${\langle[{\rm Mg}/{\rm
H}]\rangle}$. In all cases, the burst is assumed to terminate the galaxy’s SFH
at $t_{\mathrm{burst}}=t_{\rm cut}$. We choose $F_{D}=0$, 1, 2, or 3, and for
each value, we generate models with randomly chosen values of $f_{g}$ in the
range 0.05-0.5 and randomly chosen values of $t_{\mathrm{burst}}$. We then
select $(f_{g},t_{\mathrm{burst}})$ combinations that give ${\langle{\rm
age}\rangle}$ in a narrow range, corresponding to low-$\sigma$, mid-$\sigma$,
and high-$\sigma$ ellipticals at $z=0$ and to mid-$\sigma$ ellipticals at
$z=0.7$.

Figure 17: Correlation of burst stellar fraction and abundance differences for
models with terminating bursts of star formation at $z=0$ (top) and $z=0.7$
(bottom). The quantities $\delta{\langle[{\rm Mg}/{\rm Fe}]\rangle}$ and
$\delta{\langle[{\rm Mg}/{\rm H}]\rangle}$ are the changes in average
abundances relative to the model with no burst. We consider models with
dilution factor $F_{D}=0$, 1, 2, and 3, and for each $F_{D}$ the points trace
out a range of final gas fraction from $f_{g}=0.5$ (smallest
$\delta{\langle[{\rm Mg}/{\rm Fe}]\rangle}$) to $f_{g}=0.05$ (largest
$\delta{\langle[{\rm Mg}/{\rm Fe}]\rangle}$). The smooth model parameters of
the blue, pink, and green sequences correspond to those of the blue, pink, and
green models in Figure 12, except here we adopt $[{\rm Mg}/{\rm Fe}]_{\rm
cc}=0.45$ and set $\eta=0$ for the pink model. In each sequence we choose
$t_{\mathrm{burst}}$ to give ${\langle{\rm age}\rangle}\approx 10.24$ Gyr
(blue), 8 Gyr (pink), and 6 Gyr (green) at $z=0$ or ${\langle{\rm
age}\rangle}\approx 3.76$ Gyr at $z=0.7$.

Figures 17a and 17c correlate the fraction of stars formed during the burst
with $\delta{\langle[{\rm Mg}/{\rm Fe}]\rangle}$, the deviation of the final
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ from the smooth model ${\langle[{\rm
Mg}/{\rm Fe}]\rangle}$ at the corresponding ${\langle{\rm age}\rangle}$, for
the model elliptical populations at $z=0$ and $z=0.7$, respectively. Each
curve shows the results for the different-age populations at a fixed $F_{D}$,
with a lower final gas fraction $f_{g}$ corresponding to a higher burst
stellar fraction and larger $\delta{\langle[{\rm Mg}/{\rm Fe}]\rangle}$. For
$F_{D}=0$, even the $f_{g}=0.05$ models have only small increments of
$\delta{\langle[{\rm Mg}/{\rm Fe}]\rangle}<0.035$ dex, at either redshift.
However, burst models with substantial accreted gas fractions (large $F_{D}$)
can have larger $\delta{\langle[{\rm Mg}/{\rm Fe}]\rangle}$, with large
fractions of stars formed during the burst.

Figures 17b and 17d compare $\delta{\langle[{\rm Mg}/{\rm H}]\rangle}$ to
$\delta{\langle[{\rm Mg}/{\rm Fe}]\rangle}$ for the model elliptical
populations at $z=0$ and $z=0.7$, respectively. For $F_{D}=0$, models that
consume more gas (lower $f_{g}$) have both higher ${\langle[{\rm Mg}/{\rm
Fe}]\rangle}$ and higher ${\langle[{\rm Mg}/{\rm H}]\rangle}$, with a similar
correlation for the three population ages. However, this relationship is more
complicated for $F_{D}>0$ because of the dilution effect on ${\langle[{\rm
Mg}/{\rm H}]\rangle}$. At fixed $\delta{\langle[{\rm Mg}/{\rm Fe}]\rangle}$,
ellipticals with larger $F_{D}$ (darker blue curves) have more negative
$\delta{\langle[{\rm Mg}/{\rm H}]\rangle}$ than other ellipticals of similar
${\langle{\rm age}\rangle}$.

Figure 17 shows that raising ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ as much as
0.05 dex above the smooth model prediction requires a substantial fraction of
stars to be formed in a terminating burst of star formation. This can be
achieved by adding a significant amount of fresh gas during the burst through
a major accretion event. Smaller changes in ${\langle[{\rm Mg}/{\rm
Fe}]\rangle}$ can be obtained by forming 10-20% of stars in a terminating
burst with no gas accretion. A test of this terminating burst hypothesis is
that galaxies with higher ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ at a given
${\langle{\rm age}\rangle}$ or velocity dispersion should show signatures of
substantial late star formation. However, there is no clear predicted
correlation between deviations in ${\langle[{\rm Mg}/{\rm H}]\rangle}$ and
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ because of the impact of dilution on
${\langle[{\rm Mg}/{\rm H}]\rangle}$.

##  5 Discussion

###  5.1 Yields and IMF

Our calculations show that the nucleosynthetic yields and SNIa DTD are the key
factors governing the predicted ${\langle[{\rm Mg}/{\rm
H}]\rangle}$-${\langle{\rm age}\rangle}$ and ${\langle[{\rm Mg}/{\rm
Fe}]\rangle}$-${\langle{\rm age}\rangle}$ relations of ellipticals. Once we
choose the yields and DTD, the SFH has only a small impact on the predicted
relations (see Figures 6gh and 10), within the fairly wide range of behavior
allowed by our smooth SFH parameters. The choices of SFE timescale $\tau_{*}$
and outflow efficiency $\eta$ have an important impact on the early and late
evolution, respectively, of ${\langle[{\rm Mg}/{\rm H}]\rangle}$, but they
have limited impact on ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ (Figure 7). The
high SFE needed to get early enrichment in ${\langle[{\rm Mg}/{\rm
H}]\rangle}$ tends to drive down ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$, making
it difficult to explain the high observed values of both quantities in the
oldest, high-$\sigma$ ellipticals (Figure 10, Figure 12).

For our fiducial choice of yields, Equation 6, we are unable to reproduce the
high observed ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ ratios in our $z\sim 0$
and $z\sim 0.7$ samples, falling short by 0.05-0.1 dex (dotted gray curve in
Figure 12). As discussed in Section 3, the overall normalization of these
yields is motivated by the Rodríguez et al. (2023) determination of the mean
Fe yield of CCSN; the ratio $y_{\rm Mg}^{\rm cc}/y_{\rm Fe}^{\rm cc}$
corresponding to $[{\rm Mg}/{\rm Fe}]_{\rm cc}=0.45$ is motivated by the
observed plateau in $[{\rm Mg/Fe}]$ for low metallicity stars in the Milky Way
disk and halo; and the ratio $y_{\rm Fe}^{\rm Ia}/y_{\rm Fe}^{\rm cc}$ is
motivated by the expectation that models evolve to $[{\rm Mg/Fe}]\approx 0$ at
late times to explain solar abundances in MW disk stars.

Nearly all stars in APOGEE with $[{\rm Fe/H}]>-1.5$ have $[{\rm Mg/Fe}]<0.45$
(see Figure 1 of WGJT24), and a plateau at $[{\rm Mg}/{\rm Fe}]_{\rm cc}\leq
0.45$ is also supported by most of the data compiled by Kobayashi et al.
(2020) (see their Figure 12) and by the Sgr dwarf and Milky Way halo stars
shown by Sestito et al. (2024) (see their Figure 6). However, based on a
sample of halo stars from the H3 survey extending to $[{\rm Fe/H}]=-3$, Conroy
et al. (2022) argue that the true core collapse plateau lies at $[{\rm
Mg}/{\rm Fe}]_{\rm cc}\approx 0.6$, and that the roughly flat trend for stars
with $-1.5\lesssim[{\rm Fe/H}]\lesssim-0.5$ already includes a substantial
SNIa contribution (see also Maoz & Graur 2017). With our higher choice of
$[{\rm Mg}/{\rm Fe}]_{\rm cc}=0.55$, giving the yields in Equation 7, our
smooth SFH models still underproduce the observed ${\langle[{\rm Mg}/{\rm
Fe}]\rangle}$, but at a level that is possibly within the systematic
uncertainties of the SSP abundance measurements (Figure 12). Unfortunately,
massive star and CCSN models are not robust enough to settle the choice of
$[{\rm Mg}/{\rm Fe}]_{\rm cc}$ theoretically, as illustrated by the fact that
state-of-the-art models (e.g., Sukhbold et al. 2016) miss the solar [O/Mg]
ratio by 0.4 dex even though O and Mg are produced in stars of similar mass
(Griffith et al., 2021b).

The empirical evidence for $[{\rm Mg}/{\rm Fe}]_{\rm cc}\approx 0.45$ comes
from stars with $[{\rm Fe/H}]<-0.5$. Our models assume $[{\rm Mg}/{\rm
Fe}]_{\rm cc}$ is independent of metallicity, which is consistent with
theoretical predictions (e.g., Nomoto et al. (2013); Andrews et al. (2017)).
However, these theoretical predictions could be incorrect, and the yield
ratios at the $[{\rm Mg}/{\rm H}]\approx 0$ metallicities that matter most for
our models could be higher than the plateau observed for lower metallicity
stars.

Another possible explanation of the difference between our models and the data
is that the IMF-averaged Mg and Fe yields in ellipticals are simply different
than those in the Milky Way. For CCSN, the IMF-averaged yields could be
different because of a difference in the high-mass slope of the IMF, perhaps
connected to star formation in environments of high ambient pressure. A top-
heavy IMF would increase the Mg and Fe yields $y_{\rm Mg}^{\rm cc}$ and
$y_{\rm Fe}^{\rm cc}$, making it easier to reproduce the observed
${\langle[{\rm Mg}/{\rm H}]\rangle}$ and ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$
in the oldest galaxies and reducing the need for a high SFE to explain this
rapid enrichment. Using the CCSN yields of Chieffi & Limongi (2013), Griffith
et al. (2021a) find that changing the high-mass IMF slope from -2.3 (Kroupa,
2001) to -2.0 slightly decreases the predicted $[{\rm Mg}/{\rm Fe}]_{\rm cc}$,
by 0.03 dex. We used VICE (Johnson & Weinberg, 2020) to perform a similar
calculation for yields from Sukhbold et al. (2016) with different choices of
explosion landscape. Here we find that changing the high-mass IMF slope from
-2.3 to -2.0 increases $[{\rm Mg}/{\rm Fe}]_{\rm cc}$, but only by 0.02-0.03
dex. Both of these calculations suggest that the IMF shape has little leverage
on the CCSN Mg/Fe ratio. However, this dependence merits further investigation
with other yield models.

Raising the CCSN yields relative to SNIa yields (higher $y_{\rm Fe}^{\rm
cc}/y_{\rm Fe}^{\rm Ia}$) with a top-heavy IMF would help to boost
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ even if $[{\rm Mg}/{\rm Fe}]_{\rm cc}$
is unchanged, similar to the effect shown in Figure 14. Black hole formation
can also have a large impact on IMF-averaged yields (Griffith et al. 2021b;
WGJT24), even if the IMF itself is unchanged. However, while black hole
formation may well depend on metallicity, there is no obvious reason it should
be different in ellipticals vs. the Milky Way at a given metallicity.

Much of the discussion about IMF differences in ellipticals has focused on
low-mass behavior, in particular on observations suggesting that high-$\sigma$
ellipticals have more low mass ($M<0.5$ $M_{\odot}$) stars than predicted by a
Kroupa (2001) IMF (van Dokkum & Conroy, 2010, 2012; Conroy & van Dokkum,
2012b; Conroy et al., 2017). Because these star have lifetimes longer than the
age of the Universe, the main impact of increasing their relative numbers is
to lower all yields — $y_{\rm Mg}^{\rm cc}$, $y_{\rm Fe}^{\rm cc}$, and
$y_{\rm Fe}^{\rm Ia}$ — by the same factor. A bottom-heavy IMF thus makes it
more difficult to explain rapid enrichment, but it probably has little impact
on ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ predictions. If the increased numbers
extend above the $\sim 1$ $M_{\odot}$ turnoff, then the SNIa rate could be
boosted by increased white dwarf formation, possibly explaining the high SNIa
rate inferred in galaxy clusters (Maoz & Graur, 2017). An increased $y_{\rm
Fe}^{\rm Ia}/y_{\rm Fe}^{\rm cc}$ would worsen the discrepancy with observed
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$, opposite to the reduced $y_{\rm
Fe}^{\rm Ia}/y_{\rm Fe}^{\rm cc}$ models considered in Section 4.3.

###  5.2 Future Directions

Our results suggest several directions for future investigations. One is to
fit model parameters directly to spectra instead of comparing the light-
weighted ${\langle{\rm age}\rangle}$, ${\langle[{\rm Mg}/{\rm H}]\rangle}$,
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ to SSP fit values. The tests presented
in Figure 13 indicate that predicting SSP quantities is an effective shortcut,
but the residual inaccuracies in this process may still be enough to affect
the assessment of models. Fitting directly to observed spectra requires
computing model spectra for continuous sequences of stellar populations as we
have done in Figure 13, but this process could be made reasonably efficient.
For a smooth SFH and $[{\rm Mg}/{\rm Fe}]_{\rm cc}=0.45$, we would not expect
to obtain statistically acceptable fits to the data, since we cannot reproduce
the observed ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$-${\langle{\rm age}\rangle}$
relation, but inspecting the discrepancies of best-fit models in spectral
space might better reveal whether these discrepancies could arise from
systematic uncertainties in the spectral modeling.

Another promising direction is to use broad-band SED measurements (photometry
or spectrophotometry) in addition to the continuum-normalized absorption
spectra modeled by alf. Codes such as Prospector (Johnson et al., 2021a) and
BagPipes (Carnall et al., 2018) use the broad-band SED to derive SFH
constraints that go beyond effective population age defined by SSP fits. We
might therefore hope to constrain our SFH parameters with the broad-band SED
and allow absorption spectra to better constrain the yields and GCE
parameters. In the context of ellipticals, it would be especially valuable to
use the SED and spectra combination to probe the cessation of star formation
and thus get insights into the physics of quenching. Our models assume a sharp
cutoff at $t_{\rm cut}$, but they allow an exponential decline before that,
and the same ${\langle{\rm age}\rangle}$ can be obtained with a rapid
exponential decline and late cutoff or a slow decline and early cutoff. The
predicted ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$-${\langle{\rm age}\rangle}$
tracks are surprisingly insensitive to the difference (Figure 6), but the
combination with SED-fitting might distinguish such models. Most importantly,
our modeling in Section 4.4 suggests that terminating bursts of star formation
may play an important role in boosting ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$
relative to the smooth-SFH predictions. Broad-band SED modeling may be able to
confirm or refute the hypothesis that typical ellipticals terminate their star
formation with a final burst.

The combined spectra-SED constraints become still more promising if they can
be applied on a galaxy-by-galaxy basis in addition to the high-precision
constraints for binned samples. Our preliminary investigations imply that
galaxies at fixed $M_{\star}$ exhibit intrinsic scatter in SSP age, $[{\rm
Mg}/{\rm H}]$, and $[{\rm Mg/Fe}]$ comparable to the statistical uncertainties
for individual SDSS or LEGA-C spectra (A. Beverage et al., in preparation). If
terminating bursts are a common phenomenon, we expect correlations between
deviations from the mean $[{\rm Mg/Fe}]$ and the fraction of stars formed at
the end of the galaxy’s SFH, as illustrated in Figure 17. The correlation of
scatter in age, $[{\rm Mg}/{\rm H}]$, and $[{\rm Mg/Fe}]$ at fixed $M_{\star}$
or $\sigma$ can provide further insights into the termination of star
formation and the factors that govern abundance ratios.

HST and JWST studies show that there are already substantial populations of
quiescent galaxies present at $z=1-5$. The combination of ground-based and
JWST spectroscopy allows the first assessments of stellar metallicities and
$[{\rm Mg/Fe}]$ ratios in such systems (e.g., Beverage et al., 2024a, b; Kriek
et al., 2024; Slob et al., 2024; Carnall et al., 2024). For massive star-
forming galaxies, one can combine stellar $[{\rm Fe/H}]$ measurements with gas
phase [O/H] measurements to infer the $\alpha$-enhancement (e.g., Steidel et
al., 2016; Topping et al., 2020; Kashino et al., 2022). In a recent study,
Stanton et al. (2024) applied this technique to massive star-forming galaxies
at $z\approx 3.5$, finding [O/Fe]$=0.425\pm 0.026$ for the average value of
this mixed gas-stellar abundance ratio. These high redshift observations
provide another direction to extend our modeling, and probing to early epochs
may help break remaining degeneracies between yields and other GCE parameters.
Applying the WAF one-zone models to their $z\approx 3.5$ data set, Stanton et
al. (2024) find that explaining their high observed [O/Fe] ratios requires a
high oxygen yield corresponding to [O/Fe]${}_{\mathrm{cc}}\approx 0.7$, a
conclusion similar to but even stronger than the one inferred here from the
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ ratios of $z\sim 0$ and $z\sim 0.7$
ellipticals.

The other natural direction is to extend this work to model other element
abundances. With SDSS and LEGA-C samples, we can obtain robust measurements
for the mean abundances of up to a dozen elements at $z\sim 0$ and six
elements at $z\sim 0.7$ (Conroy et al., 2014; Beverage et al., 2023). These
include elements with different ratios of CCSN/SNIa enrichment (e.g., Si, Ca,
Ni, Mn) and elements in which the delayed contribution comes from asymptotic
giant branch (AGB) stars rather than SNIa (e.g., C, N, Sr, Ba). For all these
elements, we would like to know whether their abundances in ellipticals can be
explained by adopting stellar yields that reproduce Milky Way data combined
with the GCE parameters inferred from Mg and Fe. The differing CCSN/SNIa
ratios and the distinct DTD of AGB enrichment can further test the GCE models
themselves.

##  6 Conclusions

Massive quiescent galaxy populations exhibit well-defined trends of stellar
age, $[{\rm Mg}/{\rm H}]$, and $[{\rm Mg/Fe}]$ with galaxy velocity dispersion
or stellar mass. We have constructed one-zone GCE models with flexible star
formation histories to interpret the trends measured by Beverage et al.
(2023), who fit SSP models to the spectra of galaxies from the SDSS (at
$z\approx 0$) and from LEGA-C (at $z\approx 0.7$). We predict SSP quantities
from the light-weighted, logarithmically averaged quantities in the GCE model,
denoted ${\langle{\rm age}\rangle}$, ${\langle[{\rm Mg}/{\rm H}]\rangle}$, and
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$. Our tests with synthetic spectra from
composite stellar populations show that this prediction is generally accurate,
though the light-weighted ${\langle{\rm age}\rangle}$ slightly exceeds the SSP
age (by $\sim 0.5$-2 Gyr) for the oldest galaxies (Figures 13, 18).

The components of our GCE models are the star formation history (SFH), the
star formation efficiency timescale $\tau_{*}=M_{\rm gas}/\dot{M}_{*}$, the
outflow efficiency $\eta=\dot{M}_{\rm out}/\dot{M}_{*}$, and the IMF-averaged
Mg and Fe yields from CCSN and SNIa. Following WGJT24, we normalize $y_{\rm
Fe}^{\rm cc}$ using the Rodríguez et al. (2023) determination of the mean
observed Fe yield of CCSN, choose $y_{\rm Mg}^{\rm cc}/y_{\rm Fe}^{\rm cc}$
based on the plateau of $[{\rm Mg/Fe}]$ seen in metal-poor Milky Way stars,
and choose $y_{\rm Fe}^{\rm Ia}/y_{\rm Fe}^{\rm cc}$ so that GCE models with a
continuous, slowly declining SFH evolve to $[{\rm Mg/Fe}]\approx 0$ at late
times. We adopt $[{\rm Mg}/{\rm Fe}]_{\rm cc}=0.45$ as our default value of
the plateau, but we consider models with a higher $y_{\rm Mg}^{\rm cc}/y_{\rm
Fe}^{\rm cc}$ corresponding to $[{\rm Mg}/{\rm Fe}]_{\rm cc}=0.55$.

With our adopted SFH parameterization (Equation 2), star formation begins at
$t_{\rm start}$, rises approximately linearly to a maximum at
$t\approx\tau_{1}$, then declines exponentially on a timescale $\tau_{2}$ and
truncates sharply at a cutoff time $t_{\rm cut}$, where the galaxy transitions
to quiescence. Because our parameterization already allows a variety of rise
times $\tau_{1}$, we fix $t_{\rm start}=0.5\pm 0.5$ Gyr, assuming that star
formation in massive ellipticals commences at $z\gtrsim 6$. While all four SFH
parameters affect the mean age and abundances, the trajectories of
${\langle[{\rm Mg}/{\rm H}]\rangle}$ vs. ${\langle{\rm age}\rangle}$ and
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ vs. ${\langle{\rm age}\rangle}$ are not
sensitive to the individual values of $\tau_{1}$, $\tau_{2}$, and $t_{\rm
cut}$, and the value of $t_{\rm start}$ simply translates these trajectories
in time (Figure 5). Thus, once the yields, $\tau_{*}$, and $\eta$ are
specified, the choice of SFH has surprisingly little leverage on the predicted
${\langle[{\rm Mg}/{\rm H}]\rangle}$ and ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$
ratios for galaxies of a specified ${\langle{\rm age}\rangle}$.

Figure 10 summarizes the impact of GCE model parameters on ${\langle[{\rm
Mg}/{\rm H}]\rangle}$-${\langle{\rm age}\rangle}$ and ${\langle[{\rm Mg}/{\rm
Fe}]\rangle}$-${\langle{\rm age}\rangle}$ trajectories. Beyond the linear
scaling of Mg abundance with $y_{\rm Mg}^{\rm cc}$, the predicted
${\langle[{\rm Mg}/{\rm H}]\rangle}$ depends mainly on $\tau_{*}$ in old
galaxies and on $\eta$ in younger galaxies: higher SFE (lower $\tau_{*}$)
leads to more rapid early enrichment, while higher $\eta$ leads to greater
metal loss and lower metallicity at later times. Several parameters affect the
predicted ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ at fixed ${\langle{\rm
age}\rangle}$, but only at the $\sim 0.02$-dex level, and it is generally
easier to decrease the predicted ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ than to
increase it. Our fiducial models adopt a double-exponential DTD for SNIa that
approximates a $t^{-1.1}$ power-law with a minimum delay time $t_{{\rm
d,min}}=150$ Myr. If we reduce $t_{{\rm d,min}}$ to 50 Myr, predicted
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ ratios drop by $\sim 0.05$ dex because
of the more rapid onset of SNIa Fe enrichment.

Figure 12 compares a representative set of GCE models to the SDSS and LEGA-C
data. With reasonable parameter choices, the models can reproduce the observed
${\langle[{\rm Mg}/{\rm H}]\rangle}$ values, with the SDSS trend of higher
${\langle[{\rm Mg}/{\rm H}]\rangle}$ at higher $\sigma$ explained mainly by
decreasing outflow efficiency (lower $\eta$) in deeper potential wells. Even
with $\eta=0$, reproducing the super-solar ${\langle[{\rm Mg}/{\rm
H}]\rangle}$ observed in higher-$\sigma$ SDSS galaxies and all LEGA-C galaxies
requires $y_{\rm Mg}^{\rm cc}$ to be $\sim 0.1$ dex higher than our default
value, and a high SFE ($\tau_{*}\sim 0.5$ Gyr) is needed for the oldest SDSS
galaxies. The higher $y_{\rm Mg}^{\rm cc}$ is within the uncertainty estimated
by WGJT24, so we do not consider the high observed ${\langle[{\rm Mg}/{\rm
H}]\rangle}$ to be problematic for our models. However, a bottom-heavy IMF in
the most massive ellipticals (van Dokkum & Conroy, 2010; Conroy et al., 2017)
would lower the $y_{\rm Mg}^{\rm cc}$ in these systems, making their high
${\langle[{\rm Mg}/{\rm H}]\rangle}$ more difficult to explain; shifting from
a Kroupa (2001) IMF to a Salpeter-like IMF lowers $y_{\rm Mg}^{\rm cc}$ by
0.14 dex.

With our fiducial CCSN yield ratio of $[{\rm Mg}/{\rm Fe}]_{\rm cc}=0.45$, our
models underpredict the observed ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ values
at a given ${\langle{\rm age}\rangle}$ by $\sim 0.08$ dex at $z=0$ and $\sim
0.05$ dex at $z=0.7$. With the higher ratio $[{\rm Mg}/{\rm Fe}]_{\rm
cc}=0.55$, the discrepancy with the data is smaller, though once we choose
parameters to match ${\langle[{\rm Mg}/{\rm H}]\rangle}-{\langle{\rm
age}\rangle}$ the models still consistently underpredict the observed
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ (Figure 12bd). At $z=0.7$ the difference
is comparable to the $1\sigma$ statistical uncertainties in the observed
values, but at $z=0$ it is a statistically significant 0.05-0.07 dex in most
velocity dispersion bins. Spectral fits imply that the observed stellar
populations of ellipticals are old, but they are still young enough to allow
ample time for SNIa Fe enrichment. Within the fairly flexible framework of our
models, that SNIa enrichment should have depressed ${\langle[{\rm Mg}/{\rm
Fe}]\rangle}$ below the observationally inferred values.

The tension between models and the data is mild enough that it might be
explained by systematic uncertainties in the SSP spectral fits, e.g., from
imperfect stellar libraries. Reducing the observational estimates of
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ by $\sim 0.05$ dex or increasing the
estimated ages by 1-2 Gyr would align the data with our $[{\rm Mg}/{\rm
Fe}]_{\rm cc}=0.55$ models. Reducing ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ by
0.1 dex would remove the conflict with our $[{\rm Mg}/{\rm Fe}]_{\rm cc}=0.45$
models.

Alternatively, this tension could be a sign that the IMF-averaged yields in
elliptical galaxies are systematically different from the empirical estimates
of WGJT24, which are motivated by global statistics of CCSN and SNIa rates and
by abundance trends in the Milky Way. Agreement with the observed $[{\rm
Mg/Fe}]$ would be improved by increasing $y_{\rm Mg}^{\rm cc}/y_{\rm Fe}^{\rm
cc}$ (above $[{\rm Mg}/{\rm Fe}]_{\rm cc}=0.55$), or by decreasing $y_{\rm
Fe}^{\rm Ia}/y_{\rm Fe}^{\rm cc}$. For example, Figure 14 shows that a model
with $[{\rm Mg}/{\rm Fe}]_{\rm cc}=0.45$ and $y_{\rm Fe}^{\rm Ia}/y_{\rm
Fe}^{\rm cc}$ reduced from 1.65 to 1.0 can reproduce the ${\langle[{\rm
Mg}/{\rm Fe}]\rangle}$-${\langle{\rm age}\rangle}$ trends from SDSS and
LEGA-C. However, observations suggest a higher SNIa DTD normalization in
galaxy clusters relative to the field (Maoz & Graur, 2017), which would imply
a higher $y_{\rm Fe}^{\rm Ia}/y_{\rm Fe}^{\rm cc}$ in ellipticals, not lower.

Agreement between models and data can also be improved by adding a self-
enriching burst of star formation to the end of a galaxy’s smooth SFH,
immediately before the transition to quiescence. This terminating burst
reduces the light-weighted age, and it increases both ${\langle[{\rm Mg}/{\rm
H}]\rangle}$ and ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ because self-enrichment
in a short burst is dominated by CCSN. We have developed a simple, approximate
model for the impact of these terminating bursts (Figures 15 and 16, Appendix
B). Bursts that consume just the gas that remains after the smooth star
formation can boost ${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ at fixed
${\langle{\rm age}\rangle}$ by $\sim 0.03$ dex, while also boosting
${\langle[{\rm Mg}/{\rm H}]\rangle}$ by $\sim 0.02-0.08$ dex. If the burst is
fueled by accretion of a metal-poor, gas-rich companion, then the increase of
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$ can be larger, while ${\langle[{\rm
Mg}/{\rm H}]\rangle}$ may increase, stay flat, or even drop because of the
dilution of the pre-burst ISM. The importance of terminating bursts could be
probed observationally by looking for correlations between residuals from the
${\langle[{\rm Mg}/{\rm Fe}]\rangle}$-${\langle{\rm age}\rangle}$ relation and
spectral signatures of late star formation, going beyond the SSP-fitting
techniques used here.

Our results suggest several directions for future work. One is to compare our
model star formation histories to SED-based observational estimates, including
investigating the frequency and magnitude of terminating bursts. A second is
to apply our models to emerging observational results at higher redshift,
which should include progenitors of some of the oldest low-$z$ quiescent
populations. A third is to extend our analysis to other elements whose
abundance can be estimated by full spectral fitting, including elements that
have substantial contributions from intermediate-mass ABG stars. Our results
also highlight the importance of sharpening empirical constraints on IMF-
averaged yields, in the Milky Way and in other local galaxies.

Elliptical galaxies differ from disk galaxies in their star formation
histories and potentially in their outflows and other galactic scale
phenomena. Stellar yields are governed by stellar-scale astrophysics,
including the IMF, the physics of supernova explosions and black hole
populations, and binary populations that may affect massive star evolution and
the SNIa DTD. Building on the framework developed here, we hope to learn
whether galactic scale effects fully explain the distinctive abundance
patterns of ellipticals, or whether these patterns also imply different
stellar scale astrophysics in the progenitors of massive quiescent galaxies.

## Acknowledgments

We thank Adam Wheeler, Todd Thompson, Liam Dubay, and Jennifer Johnson for
helpful discussions. This work is supported by NSF grant AST-2307621. NMG
acknowledges the hospitality of the Kavli Institute for Theoretical Physics
(KITP), supported by NSF PHY-2309135, during the completion of this work. This
work is dedicated to the loving memory of Marcelina Kroenke. Niech spoczywa w
pokoju.

##  Appendix A fanCE

fanCE (Flexible Analytic Chemical Evolution) is an analytic one-zone chemical
evolution code. It implements the Weinberg, Andrews, & Freudenburg (2017)
analytic solution for a constant, exponential declining, or linear-exponential
(delayed tau) star formation history. It also implements a flexible, two-
parameter star formation history $\propto(1-e^{-t/\tau_{1}})\cdot
e^{-t/\tau_{2}}$, which rises linearly on a timescale $\tau_{1}$ and then
falls exponentially on a timescale $\tau_{2}$ for $\tau_{1}>\tau_{2}>0$.

We provide Python implementations of fanCE via
https://github.com/nmgountanis/fanCE so that others can easily adapt them to
interpret other galaxies or stellar populations or explore alternative
assumptions. The source code for fanCE can be downloaded from GitHub and
installed by running:

cd <path_to_installation>

git clone https://github.com/nmgountanis/fanCE

cd fanCE

pip install -e .

If you make use of this code, please cite Weinberg et al. (2017) for the
analytic solutions and Gountanis et al. (2024, this paper) for the code
implementation. The reference for the default yield choices is Weinberg et al.
(2024).

##  Appendix B Mean metallicity of burst stars

Our goal is to compute the mean metallicity of stars formed during a burst
described by a “leaky box” chemical evolution model (Hartwick 1976; Binney &
Merrifield 1998, section 5.3; WAF, section 3.5), which generalizes “closed
box” evolution (Talbot & Arnett, 1971) to allow outflow. These models assume
instantaneous enrichment and instantaneous recycling of stellar mass loss. For
notational simplicity, we will consider a single element of mass fraction $Z$
and yield $y$. The metallicity evolution of the leaky box model is

$Z(t)={y\over 1+\eta^{\prime}-r^{\prime}}\left[-\ln f_{g}(t)\right]~{},$ (B1)

where

$f_{g}(t)=M_{g}(t)/M_{g,i}$ (B2)

is the ratio of the gas mass to the initial gas mass. We have denoted the
outflow mass-loading by $\eta^{\prime}$ to emphasize that it may differ from
the value $\eta$ that applies during pre-burst evolution, and we have denoted
the recycling factor by $r^{\prime}$ because a burst may not last long enough
for the full recycling factor $r$ to be applicable. If the initial gas
metallicity is non-zero (as it is in our models), then the enrichment
described by Equation (B1) is added to the initial abundance $Z_{i}$ at the
start of the burst.

Because there is, by assumption, no ongoing gas accretion during the burst,
mass conservation implies

$M_{g}(t)+(1-r^{\prime})M_{*,{\rm form}}(t)+\eta^{\prime}M_{*,{\rm
form}}(t)=M_{g,i}~{},$ (B3)

where

$M_{*,{\rm form}}(t)=\int_{0}^{t}\dot{M}_{*}(t^{\prime})dt^{\prime}$ (B4)

is the mass of stars formed and multiplying by $(1-r^{\prime})$ gives the
current stellar mass. Dividing by $M_{g,i}$ and rearranging yields

$f_{*}(t)\equiv{M_{*,{\rm form}}(t)\over M_{g,i}}={1-f_{g}(t)\over
1+\eta^{\prime}-r^{\prime}}~{}.$ (B5)

For small $\eta^{\prime}$ and $f_{g}$, the value of $f_{*}$ can exceed one
because recycled gas can be formed into new stars. Stars formed at time
$t^{\prime}$ have metallicity $Z(t^{\prime})$, so the mean metallicity of all
stars formed up to time $t$ is

$\langle
Z_{*}\rangle={\int_{0}^{t}Z(t^{\prime})\dot{M}_{*}(t^{\prime})dt^{\prime}\over\int_{0}^{t}\dot{M}_{*}(t^{\prime})dt^{\prime}}={\int_{0}^{f_{*}}Z(f_{*}^{\prime})df_{*}^{\prime}\over\int_{0}^{f_{*}}df_{*}^{\prime}}={\int_{1}^{f_{g}}Z(f_{g}^{\prime})df_{g}^{\prime}\over\int_{1}^{f_{g}}df_{g}^{\prime}}~{},$
(B6)

where the last equality uses the fact that
$df_{*}^{\prime}/df_{g}^{\prime}=-(1+\eta^{\prime}-r^{\prime})^{-1}$ is
constant and therefore cancels between the numerator and denominator. Applying
Equation (B1) and evaluating the integral yields the final result

$\langle Z_{*}\rangle={y\over 1+\eta^{\prime}-r^{\prime}}\left[1+{f_{g}\ln
f_{g}\over 1-f_{g}}\right]~{}.$ (B7)

The factor in $[...]$ grows from $\frac{1}{2}(1-f_{g})$ at early times, when
$(1-f_{g})\ll 1$, to unity at late times, when $f_{g}\ll 1$.

Our terminating burst models are parameterized by a dilution factor $F_{D}$
and a final gas fraction $f_{g}$. (Recall that $f_{g}$ is the ratio of the gas
mass at the end of the burst to the gas mass at the start of the burst, not to
the stellar mass.) At the end of smooth evolution the gas mass is
$\tau_{*}\dot{M}_{*}$, so the initial gas mass for the burst is
$M_{g,i}=(1+F_{D})\tau_{*}\dot{M}_{*}$. We assume that the accreted gas is
much lower metallicity than that in the parent elliptical, so that the gas
phase abundance of element X is reduced by a factor $(1+F_{D})$. We therefore
compute the impact of a burst on our elliptical galaxy models by adding, in
the final timestep, a mass of stars

$M_{\rm burst}=\tau_{*}\dot{M}_{*}(1+F_{D}){1-f_{g}\over
1+\eta^{\prime}-r^{\prime}}$ (B8)