# Mass Metallicity Relationship of SDSS Star Forming Galaxies: Population
Synthesis Analysis and Effects of Star Burst Length, Extinction Law, Initial
Mass Function and Star Formation Rate

Eva Sextl  Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians
Universität München, Scheinerstr. 1, 81679 München, Germany  Rolf-Peter
Kudritzki  Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians
Universität München, Scheinerstr. 1, 81679 München, Germany  Institute for
Astronomy, University of Hawaii at Manoa, 2680 Woodlawn Drive, Honolulu, HI
96822, USA  H. Jabran Zahid  Microsoft Research, 14820 NE 36th St, Redmond, WA
98052, USA  I-Ting Ho  Max-Planck-Institute for Astronomy, Königstuhl 17,
D-69117 Heidelberg, Germany

###### Abstract

We investigate the mass-metallicity relationship of star forming galaxies by
analysing the absorption line spectra of $\sim$200,000 galaxies in the Sloan
Digital Sky Survey. The galaxy spectra are stacked in bins of stellar mass and
a population synthesis technique is applied yielding metallicities, ages and
star formation history of the young and old stellar population together with
interstellar reddening and extinction. We adopt different lengths of the
initial starbursts and different initial mass functions for the calculation of
model spectra of the single stellar populations contributing to the total
integrated spectrum. We also allow for deviations of the ratio of extinction
to reddening RV from 3.1 and determine the value from the spectral fit. We
find that burst length and RV have a significant influence on the
determination of metallicities whereas the effect of the initial mass function
is small. RV values are larger than 3.1. The metallicities of the young
stellar population agree with extragalactic spectroscopic studies of
individual massive supergiant stars and are significantly higher than those of
the older stellar population. This confirms galaxy evolution models where
metallicity depends on the ratio of gas to stellar mass and where this ratio
decreases with time. Star formation history is found to depend on galaxy
stellar mass. Massive galaxies are dominated by stars formed at early times.

techniques: spectroscopic — galaxies: abundances — galaxies: evolution —
galaxies: stellar content

††journal: ApJ

##  1 Introduction

Star forming galaxies show a tight relationship between total stellar mass and
average metallicity, the mass-metallicity relationship (“MZR”, see for
instance Lequeux et al. 1979; Tremonti et al. 2004; Kudritzki et al. 2016;
Bresolin et al. 2022). This relationship while evolving and becoming steeper
holds out to redshifts z larger than 3.3 (Vale Asari et al. 2009; Zahid et al.
2014; Genzel et al. 2015; Sanders et al. 2021). In view of the complicated
interplay between complex processes such as gas accretion, star formation,
nucleosynthesis and chemical evolution, stellar and galactic winds, dynamical
evolution and merging the existence of a simple relationship such as the MZR
seems to be a Rosetta stone to understand galaxy formation and evolution. As a
consequence, a large variety of approaches has been developed to use the MZR
as a constraint of galaxy formation and evolution, Davé et al. (2011), Yates
et al. (2012), Dayal et al. (2013), Schaye et al. (2015), Peng et al. (2015),
Spitoni et al. (2017a, b), Pantoni et al. (2019), Lapi et al. (2020), Spitoni
et al. (2020), Trussler et al. (2020), Kudritzki et al. (2021a, b) to name
only a few. However, there is an important problem. Many of the observational
MZR results obtained are based on measurements of the strongest ionized
interstellar medium (ISM) emission lines, which are then converted into oxygen
abundances. These “strong line methods” depend heavily on the calibration
method used. Kewley & Ellison (2008) in an emission line analysis of 50,000
star forming Sloan Digital Sky Survey (SDSS) galaxies demonstrated
convincingly that the systematic uncertainties in oxygen abundances can be as
large as 0.8 dex. As a result, the slope of the MZR can vary significantly for
the same set of observed emission line data depending on the calibration
applied.

An obvious alternative approach is the analysis of low resolution absorption
line spectra of individual stars. Blue and red supergiants (BSG and RSG) as
the brightest stars at visual or infra-red light, respectively, are the ideal
objects for this purpose. Galaxy metallicities obtained by this method have
only small uncertainties of at most 0.1 dex (systematic and random, see Gazak
et al. 2015) and have been used to construct a very reliable MZR, see
Kudritzki et al. 2016; Davies et al. 2017; Bresolin et al. 2022, and
references therein. Unfortunately, the method is limited with respect to
galactic distance. The most distant objects studied so far are the Antennae
galaxies and their super star clusters at 20 Mpc (Lardo et al. 2015).

In view of this limitation a third technique which has been well advanced over
the last two decades has become truly important, the analysis of stellar
absorption line spectra of the integrated stellar population through
population synthesis techniques. The advantage of this technique is that it
reaches out to much larger distances because of the relatively high surface
brightnesses of galaxies. In addition, it can provide information about the
young and old stellar population and the star formation history.

The general idea was pioneered by Tinsley (1968) and Spinrad & Taylor (1972).
With significant further improvements of stellar evolution theory and stellar
spectral modelling the concept was then advanced by many groups, see for
instance Arimoto & Yoshii (1987), Bressan et al. (1996), Reichardt et al.
(2001), Bruzual & Charlot (2003), Cardiel et al. (2003) and references
therein. More recent work used the advantages of increased computer power and
improved the method for applications on massive spectroscopic galaxy surves,
first focussing on spectra of the central regions of galaxies (Gallazzi et al.
2005; Cid Fernandes et al. 2005; Panter et al. 2008; Conroy et al. 2009;
Koleva et al. 2009; Sánchez-Blázquez et al. 2011; Peng et al. 2015; Zahid et
al. 2017) and later on spectra of spatially resolved regions of galaxies
obtained with Integral Field Units (Yoachim et al. 2012; González Delgado et
al. 2015; Schaefer et al. 2017; Parikh et al. 2021).

Given the enormous potential of this technique it is important to carefully
assess the systematic effects arising from the large variety of assumptions
which enter this approach. While this has been done already to a large extend
in previous work, our work presented here is investigating a variety of
additional important aspects, which have not been addressed so far in much
detail. We study the effects of finite star formation burst lengths compared
with instantaneous bursts, the effects of pre-main sequence evolution and
initial mass function and we assess the influence of an interstellar reddening
law characterized by a ratio RV of extinction to reddening deviating from the
standard value of 3.1. We will also investigate the influence of star
formation rate (SFR) on MZR, extinction and ages. For this purpose we start
from our previous work (Zahid et al. 2017) and apply an improved technique to
analyse a sample of 200,000 SDSS star forming galaxies with spectra stacked
into bins of stellar mass. We compare the MZRs obtained in this way for the
different assumptions with the MZR from supergiant spectroscopy, which we use
as benchmark. We will also discuss stellar ages and star formation history.

Figure 1:  Fluxes (arbitrary units) of the stacked SDSS spectra in 25 mass
bins from log M∗ = 8.55 (top) to 10.95 (bottom). For clarity the spectra are
shifted by 0.2 units from one mass bin to the next.  Figure 2:  The grid of
SSB metallicities [Z] and ages for the calculation of model spectra using eq.
(1). Filled red circles: high age resolution grid, green circles: low
resolution grid.

##  2 Observed Spectra and Analysis Technique

The same dataset as described in Zahid et al. (2017) is used. The observations
consist of SDSS spectra of about 200,000 star forming galaxies in the redshift
range 0.027 $\leq$ z $\leq$ 0.25. The galaxies were selected according to the
strengths of their interstellar medium (ISM) emission lines and the redshift
range was chosen to ensure rest-frame wavelength coverage of the [OII]
$\lambda$3727 and $\lambda$7330 emission lines. Since we are interested only
in the average galaxy properties as a function of their stellar mass, we have
stacked the spectra in 25 bins of stellar mass log M∗ (measured in units of
$M_{\odot}\;$) from 8.55 to 10.95 in steps of 0.1 dex. In this way the signal-
to-noise is strongly enhanced for the analysis and ranges from 160 at the low
mass end to 1300 at higher masses (see Table in Zahid et al. (2017)). To avoid
systematic effects as a function of stellar mass caused by different line
widths we have convolved all spectra to the line width FWHM = 330 km s-1 found
for the bin with the highest stellar mass. In addition, the spectra are
normalized by setting the mean flux between 4400 and 4450 Å equal to unity.
Details of galaxy selection, stellar mass determination, stacking and
convolution procedure are given in Zahid et al. (2017).

The spectra of the 25 mass bins are shown in Figure 1. Strikingly, the SED
slope changes from low to high stellar mass and the stellar absorption lines
become stronger. We will show that this is caused by systematic differences in
interstellar extinction, average stellar age and metallicity.

The technique to analyse these spectra in order to extract information about
stellar population metallicity and ages is a modification of the approach
introduced by Zahid et al. (2017). We model the spectrum of the integrated
stellar population Mλ as a linear combination of the spectra of single stellar
bursts (SSB) fλ,i(ti, [Z]i) with age ti and logarithmic metallicity [Z]i = log
Z/Z☉

$M_{\lambda}=D_{\lambda}(R_{V},E(B-V))\sum_{i=1}^{n_{SSB}}b_{i}f_{\lambda,i}(t_{i},[Z]_{i}).$
(1)

The coefficients bi describe the contribution of burst i with age ti and
metallicity [Z]i. D${}_{\lambda}(R_{V},E(B-V))$ accounts for the absorption by
interstellar dust, which depends on interstellar reddening E(B-V) and RV =
AV/E(B-V) the ratio of interstellar extinction to reddening in the V-band. The
value of RV characterizes the reddening law, at least to first order. As is
well known from studies of stellar SEDs in star forming regions, RV can cover
a wide range from 2 to 6 (see, for instance, Urbaneja et al. 2017). We use the
interstellar reddening law by O’Donnell (1994), which is a modification of
Cardelli et al. (1989).

Figure 3:  Luminosities of SSBs as a function of age (yrs) for different
spectral passbands (pink: U, blue: B, green: V, red: R, orange: I). The solid
curves correspond to a burst duration tb = 0.1 Myr, the dashed curves to tb =
10 Myr. The bottom figure shows the logarithm of the ratio of luminosities of
10 Myr to 0.1 Myr bursts at larger ages. A metallicity of [Z]=0 has been used
for this example.  Figure 4:  Same as Figure 3 (top) but the dashed curves
correspond to tb = 1 Myr.

The SSB spectra are calculated using the Flexible Stellar Population Synthesis
code (FSPS, version 3.2) (Conroy et al. 2009; Conroy & Gunn 2010) and adopting
a Chabrier (2003) initial mass function, the MILES stellar spectral library
(Sánchez-Blázquez et al. 2006) and the MESA stellar evolution isochrones
(MIST, Choi et al. 2016; Dotter 2016). We also account for finite time lengths
of the stellar bursts (0.1, 1.0, 10.0 Myr). This will be discussed in the next
section. Finally, our model spectra are normalized in the same way as the
observed spectra, convolved to the equal resolution and interpolated to the
same integer values of wavelength points as the observed spectra.

The modelling of the spectra by eq. (1) does not account for the contribution
of a potential HII-region nebular continuum, which can become important in
cases of extremely high star formation rates with very strong nebular emission
as shown by Cardoso et al. (2022). However, the strength of the hydrogen
nebular emission lines relative to the stellar continuum allows us to estimate
the effects of nebular continuum emission for our analysis. The Hβ emission
equivalent widths in our stacked spectra is $\lesssim$ 20 Å and, thus, the
nebular continuum contribution to the total spectrum is very small. We have
verified this conclusion by estimating the nebular contribution from the
measurement of Hβ emission equivalent width and by subtracting a theoretical
reddened nebular continuum. We have then repeated the analysis with modified
spectrum and found only very small changes with respect to metallicity ($\leq$
0.01 dex) and age ($\leq$ 1%). In our analysis of the SFR dependence of the
stellar MZRs in section 9 below we have encountered a very few cases in the
highest SFR quintile with larger Hβ equivalent widths. But even here the
effects on metallicity were smaller than 0.05 dex.

Figure 5:  Spectra of SSBs at different ages for different burst lengths
(pink: tb = 0.1 Myr, cyan: 1 Myr, orange: 10 Myr). Top: age = 10 Myr, middle:
19 Myrs, bottom: 45 Myrs)

For the selection of metallicities and ages of the SSB spectra we use two
rectangular grids consisting of 9 [Z] values and 52 (high age resolution) and
20 (low age resolution) ages, respectively. The metallicities range from [Z] =
-1.5 to 0.5 increasing in steps of 0.25 dex. The low resolution ages start at
1Myr and increase in logarithmic steps of $\Delta$ log t[Gyr] = 0.25 until
7.08 Gyr and then continue with 8.91, 10.0 and 12.59 Gyr. The high resolution
age grid starts at 0.1 Myr with a step to 1.0 Myr and then continues with
logarithmic steps $\Delta$ log t[Gyr] alternating between 0.05 and 0.1 until
12.59 Gyr are reached. This leads to a total number of SSB nSSB = 468 in the
high age resolution and 180 in the low resolution case. Metallicities and ages
are illustrated in Figure 2.

In the fit of the observed spectra with the model spectra we then adopt a grid
of RV and E(B-V) values. For each pair of RV and E(B-V) we calculate the
wavelength dependent dust attenuation factor D${}_{\lambda}(R_{V},E(B-V))$ and
use the bonded-variables least squares (BVLS) algorithm by Stark & Parker
(1995), translated to IDL by Michele Cappellari111https://www-
astro.physics.ox.ac.uk/~cappellari/idl/bvls.pro to determine the coefficients
bi. With the bi the comparison of the model spectrum with the observations
leads to a $\chi^{2}$ value for each RV, E(B-V) pair and the minimum of
$\chi^{2}$ defines the best fit yielding final values for RV and E(B-V) and
average metallicity [Z] and age tav defined as

$[Z]=\sum_{i=1}^{n_{SSB}}b_{i}[Z]_{i}$ (2)

and

$t_{av}=\sum_{i=1}^{n_{SSB}}b_{i}t_{i}.$ (3)

Introducing a young and old stellar population through t${}_{i}\leq$ 1.6 Gyr
and t${}_{i}\geq$ 1.6 Gyr, respectively, we can calculate corresponding
metallicities [Z]young, [Z]old and ages tyoung, told via

$b_{young}=\sum_{i_{young}}b_{i},b_{old}=\sum_{i_{old}}b_{i}$ (4)

and

$[Z]_{young}={1\over b_{young}}\sum_{i_{young}}b_{i}[Z]_{i}$ (5)
$[Z]_{old}={1\over b_{old}}\sum_{i_{old}}b_{i}[Z]_{i}$ (6) $t_{young}={1\over
b_{young}}\sum_{i_{young}}b_{i}t_{i}$ (7) $t_{old}={1\over
b_{old}}\sum_{i_{old}}b_{i}t_{i}.$ (8)

The choice of 1.6 Gyr for the definition of the young population is motivated
by galaxy chemical evolution models which show only small changes ($\leq$ 0.1
dex) of metallicity during the final 1 to 2 Gyrs of evolution (Spitoni et al.
2017a, b; Kudritzki et al. 2021a, b). Moreover, the molecular gas depletion
time which defines star formation efficiency, is of the same order (1.9 Gyr,
see Leroy et al. 2008). They young and old population fitting is insensitive
to the selected value, for instance, chosing 2.5 instead of 1.6 Gyr does not
change the results.

We note that due to the normalization of our spectra between 4400 and 4450 Å
the [Z] and t-values obtained in this way are effectively B-band luminosity
weighted metallicities and ages. This will become more evident further below.

The results for metallicities, ages, reddening and extinction obtained in this
way will be discussed in the next sections.

Figure 6:  Same as Figure 5 but with a different wavelength range for ages of
45 Myr (top) and 501 Myr (bottom).

##  3 The Effects of Star Formation Burst Length on SSB spectra

Most commonly, the spectra of star bursts are calculated by adopting the
approximation of so-called “instantaneous bursts”, for which the stellar
evolution isochrone at one singular time step is used. However, as is well
known, star formation in star forming regions is not instantaneous but happens
over finite time scales. For instance, as shown by Efremov & Elmegreen (1998),
the star formation time tb of a burst is correlated with the size of the star
forming region and can cover a wide range from 0.1 to 10 Myr in regions where
associations and star clusters form. More recent work on star formation in
molecular clouds confirmed burst lengths of a few million years (Corbelli et
al. 2017; Kim et al. 2021). It is, thus, important to investigate what
influence tb can have on the burst model spectra and the resulting analysis.

If Liso($\lambda,\tau$) is the total luminosity per Angstroem at wavelength
$\lambda$ of an isochrone of age $\tau$, then the luminosity
L($\lambda$,t,t0,tb) of a burst at time t, which started at t0 is

$L(\lambda,t,t_{0},t_{b})={1\over(t-t_{0})}\int_{t_{0}}^{t}L_{iso}(\lambda,\tau)d\tau,$
(9)

for t${}_{0}\leq$ t $\leq$ tb+t0 and

$L(\lambda,t,t_{0},t_{b})={1\over
t_{b}}\int_{t-t_{0}}^{t}L_{iso}(\lambda,\tau)d\tau$ (10)

for t $\geq$ tb+t0.

We adopted constant star formation rate during the burst and
L($\lambda$,t,t0,tb) is the luminosity averaged over the duration of the
burst. Figures 3 and 4 show the effects of burst duration on the luminosities
of SSBs as a function of age in different spectral bands. They are quite
dramatic around the peak at ages of 5-10 Myrs (up to 0.3 dex) and then become
weaker but continue to be noticeable ($\geq$ 0.005 dex) compared with the S/N
of the spectra until values of tb/t $\approx$ 0.02 are reached (see Figures 3
bottom). This is a consequence of the power law behaviour Liso(X,t)$\propto$
t-n for larger values of t which leads to L(X,t,t0,tb) = Liso(X,t)(1+ntb/(2t))
for small values of tb/t. X stands for the spectral passband (U, V, B, R, I).
In the same way the spectra of SSBs show the effects of burst length. Examples
are given in figures 5 and 6. The differences between the spectra of different
burst length become small at ages much larger than the burst length as
discussed above. We note that all luminosities are given in units of the solar
luminosity L☉.

The influence of the adopted burst lengths on the analysis of the observed
spectra will be investigated in the next section.

Figure 7: Analysis of SDSS star forming galaxies. Metallicities of the young
(blue) and old (red) stellar population and average metallicities (green). The
grid of SSB with high time resolution was used. Top: 0.1 Myr burst length,
middle: 1 Myr Myr, bottom: 10 Myr. The MZR of blue and red supergiant stars
(pink and orange stars) obtained from the analysis of individual stars in 17
galaxies is also shown as a benchmark for the young stellar population.

Figure 8:  Same as Figure 7 but for the SSB grid with low time resolution.

Figure 9: Top: SDSS fit metallicities of Figure 7 (high time resolution) with
all burst lengths overplotted (0.1 Myr: solid, 1 Myr: dashed, 10 Myr: dashed
dotted). Middle: same as top but with metallicities of Figure 8 (low time
resolution). Bottom: Assessement of the quality of the fits obtained with the
different sets of SSB. Difference $\Delta\chi^{2}$ =
$\chi^{2}_{min}-\chi^{2}_{min,0}$ for the different SSB spectral model fits as
a function of galaxy stellar mass. $\chi^{2}_{min,0}$ is the value obtained
with SSBs calculated with high time resolution and tb = 0.1 Myr. Circles
correspond to SSB with high time resolution, squares to low resolution, red
color to 0.1 Myr bursts, blue color to 1.0 Myr and green to 10 Myr.  Figure
10:  SSB spectral fit (red dots) of the observed spectrum (blue) for log M∗ =
9.75. Note that for the fit we mask out regions contaminated by ISM emission
or absorption lines. The fit spectrum is calculated with a burst length tb =
0.1 Myr and high time resolution.

##  4 Analysis with Different Burst Lengths and Time Resolution

In this section we use SSBs calculated for different burst lengths and time
resolution to analyse our stacked SDSS spectra. We note that we allow for
variations of RV in our spectral fits. As will be discussed in section 5, the
RV obtained are substantially higher than 3.1 (see Figure 13).

Figures 7 and 8 summarize the results with respect to the MZR. We show
metallicities of the young and old population and average metallicities as a
function of stellar mass obtained when using the different sets of SSBs and
compare with with the MZR of massive supergiant stars, which has been obtained
by quantitative spectroscopy of individual objects in 17 nearby galaxies (see
Bresolin et al. 2022 and references therein). These metallicities are accurate
to 0.1 dex and serve as a benchmark for the young stellar population. The
metallicity uncertainties of the SDSS spectral fits as shown in Figures 7 and
8 are up to 0.07 dex for the young and old population, respectively, and 0.03
dex for the average over all ages. They are obtained from Monte Carlo
simulations with fits of the observed stacked spectra modified by adding
Gaussian noise.

Figures 7 and 8 indicate that burst length and SSB time resolution have a
significant effect of up to 0.4 dex on the determination [Z]young, the
metallicity of the young stellar population. The influence on the average
metallicity [Z] is somewhat smaller (up to 0.2 dex) and the differences for
the old population are below or about 0.1 dex and, thus, of the order of the
fit uncertainties. A comparison of all metallicities obtained with different
burst lengths is given in Figure 9.

The goodness of the fits is assessed in Figure 9 where we compare the
$\chi^{2}$-values obtained with the different sets of SSB for the fit of the
observed SDSS spectra. (We note that number of wavelength points used for the
calculation of $\chi^{2}$ is npix = 2523). Figure 9 shows that the SSB models
with the shortest burst length of 0.1 Myr and high time resolution provide the
best fit, although the high time resolution fits with 1 Myr burst length come
close. It is reassuring that for these two cases the SDSS young stellar
population metallicities are in very good agreement with the metallicities of
our supergiant benchmark sample (see Figure 7) . We take this as an indication
that the population synthesis method is reliable.

Figure 10 shows the fit of the SDSS spectrum at log M∗ = 9.75 as a typical
example.

Figure 11:  Ages tyoung (top), tav (middle) and told (bottom) as defined by
equations (3), (7) and (8) as a function of galactic stellar mass obtaned with
spectral fits suing high time reolution and a burst length tb = 0.1 Myr. Red
values obtained with RV fitted. Blue values with RV = 3.1.  Figure 12:  The
ratio of star formation of young to the old stellar population as a function
of galactic stellar mass. The fit of the observed spectra was carried out with
tb = 0.1 Myr and high time resolution.

The metallicities obtained for the old stellar population [Z]old are
significantly lower than [Z]young. This reflects galactic chemical evolution
and will be discussed below.

We note that Zahid et al. (2017) in their work have used SSB with 10 Myr burst
length and compared their average metallicities [Z] with the supergiant values
finding good agreement. We confirm their result (see Figure 8), but we note
that this comparison is misleading. Since supergiants have ages around 50 Myr,
a comparison with [Z]young would have been more appropriate.

Figure 13:  Top: The ratio of V-band extinction to reddening RV as a function
of stellar mass log M∗. Middle: Reddening E(B-V) as a function of stellar
mass. Red: values obtained with RV as a free parameter and determined as in
the top figure. Blue: values obtained with RV = 3.1. Bottom: Interstellar
V-band extinction AV. Red and blue as in the middle figure. tb = 0.1Myr and
high age resolution were used for the spectra fit.  Figure 14:  Same MZR as in
Figure 7 (tb = 0.1Myr, high age resolution) but with the results of spectral
fits with RV = 3.1 additionally shown as squares.

The average stellar ages tav and the ages of the young and old population
obtained with SSB of tb = 0.1 Myr and the high age resolution grid are
displayed in Figure 11 together with their error bars. The ages between 50 and
300 Myr of the young population confirm that a comparison with supergiant
stars as a benchmark is appropriate. We notice a clear correlation of tyoung
with galactic stellar mass. No correlation is obtained for tav and the age
told of the old population.

Figure 15:  Top: Metallicities obtained with different normalization, where
violet, cyan and blue belong to the young population, pink, orange and red to
the old population and green colors show the average metallicities. Circles
correspond to B-band normalization, squares to V-band, and triangles to
R-band. Middle: Ages of the young population. Blue color: B-band
normalization, green color: V-Band, red color: R-band. Bottom: Same as the
middle plot, but for the old population. Again tb = 0.1Myr and high age
resolution were used for the spectra fit.

The average redshift of our galaxy sample is z=0.08 corresponding to a
cosmological lookback time of 1.05 Gyr (assuming a Hubble constant H0=70.4 km
s-1 Mpc-1 and a flat universe with $\Omega_{\Lambda}$=0.728 and
$\Omega_{m}$=1-$\Omega_{\Lambda}$). With t${}_{old}\approx$ 10 Gyr we then
know that the stars of the old population have been formed on average at a
lookback time of 11 Gyr equivalent to a redshift of z $\approx$ 2.5. Thus, the
metallicities obtained for the old stellar population in Figure 7 and 8
correspond to z $\approx$ 2.5.

The fit of the observed spectrum yields the coefficients bi in eq. (1) for the
individual bursts, which can be related to the relative number contribution Ni
of stars of isochrone i, for which the SSB model spectrum is calculated, via

$b_{i}=N_{i}L_{i}(X).$ (11)

Figure 16:  MZRs obtained by analysis with model spectra assuming different
IMFs. Curves the in upper part of the figure correspond to the young
population, curves in the lower part to the old population. Blue: Chabrier
(2003), pink: Kroupa (2001), green: van Dokkum (2008), orange: Salpeter
(1955).

Li(X) is the luminosity of the isochrone in passband X with an effective
wavelength which corresponds to the wavelength at which the observed spectra
and the SSB spectra are normalized. Since our normalization interval is 4400
to 4450 Å the B-band is the appropriate choice for X (examples for Li(X) as a
function of isochrone age are given in Figures 3 and 4). Eq. (11) explains why
the metallicities and ages of our spectral fits are B-band luminosity weighted
quantities.

Since N${}_{i}\propto\psi_{i}\Delta$ti, where $\Delta$ti is the time interval
between isochrone i and the next older isochrone in our SSB grid, Eq. (11) can
be used to estimate star formation rates $\psi$ for the young and old
population

$\psi_{young}\propto\sum_{i_{young}}{b_{i}\over\Delta
t_{i}L_{i}};\psi_{old}\propto\sum_{i_{old}}{b_{i}\over\Delta t_{i}L_{i}}.$
(12)

Figure 12 shows the ratio of star formation of the young to the old population
as a function of galaxy stellar mass. We see the obvious trend that low mass
galaxies are presently more active than high mass galaxies.

Figure 17:  The chemical evolution of star forming SDSS galaxies. Top:
Metallicities of the young and old stellar population (as in Figure 7 top) and
of supergiant stars compared with B-band luminosity weighted metallicities
predicted by galaxy evolution models (Kudritzki et al. 2021a, b) for very
young stars (cyan) and stars 10 Gyr older (violet). Bottom: Mass of the galaxy
evolution models at lookback time $\Delta$t = 10 Gyr as a function of galaxy
final stellar mass.

##  5 Reddening and the Effects of the Extinction Law

Since a long time there have been clear indications of deviations from the
standard value RV = 3.1 in star forming regions of the Milky Way (see, for
example, Feinstein et al. 1973; Herbst 1976; The et al. 1980; Cardelli et al.
1989). Similar results have been found for the Large Magellanic Cloud, see
Maíz Apellániz et al. (2017), Urbaneja et al. (2017) and Holwerda & Keel
(2013) in their differential IFU SED study of two galaxy pairs conclude that
the canonical extinction law with RV = 3.1 does not fit their data. In
consequence, it is important to keep RV as free parameter for the SSB fits of
the observed spectra as described in section 2.

Figure 13 (top) shows the RV values obtained in our analysis for the case of
tb = 0.1Myr and high age resolution. While the values at lower stellar masses
have large uncertainties, we generally find values larger than 3.1 at all
stellar masses. We believe the reason is that the SDSS observations sample the
central regions of the star forming galaxies where still a lot of star
formation in dense clouds is going on. These natal birthclouds act as an
additional source of attenuation. In consequence, the contribution of the
diffuse ISM with RV = 3.1 is small. We speculate that the reason for the
obvious correlation with stellar mass is that the young stellar population is
younger at lower mass (see Figure 11) and that the density in the star forming
regions is higher. The influence of varying density has been described by Witt
& Gordon in a series of seminal papers (Witt & Gordon 1996; Gordon et al.
1997; Witt & Gordon 2000). As a rule of thumb, an increasing inhomogeneity
with alternating high- and low-density regions in the ISM structure leads to
flattened (’grayer’) attenuation curves. In our case, using a CCM/O’Donnell
extinction law, a shallower slope is equivalent to high values of $R_{V}$. See
also Berlind et al. (1997) for an in depth discussion of this issue in case of
the galaxy NGC2207.

We note that the average value found for the LMC (log M∗ = 9.2) by Urbaneja et
al. (2017), RV = 4.6, is in good agreement with Figure 13. It is important to
note that the difference from the standard reddening law in star forming
galaxies at all masses has important consequences for other aspects of
astrophysics. For instance, the estimate of galaxy luminosities and the
corresponding mass to light ratios can be affected. Distance determinations
based on stellar distance indicators and the resulting cosmological distance
ladder might be influenced as well (see, for instance, Urbaneja et al. 2017;
Falco et al. 1999).

Reddening E(B-V) and V-band extinction AV are also shown in Figure 13 (middle
and bottom, respectively). Again, there is a clear correlation with stellar
mass. Such correlations have already been found in previous work (Brinchmann
et al. 2004; Asari et al. 2007; Garn & Best 2010; Zahid et al. 2013). We note,
however, that the deviation of RV from the standard value 3.1 leads to
significantly higher extinction of about 0.2 magnitudes.

Ignoring a potential variation of RV and restricting to values of 3.1 also
affects the determination of metallicities and the MZR. This is demonstrated
in Figure 14, where results obtained with RV = 3.1 are over-plotted. With
differences up to 0.4 dex the effects are significant for the young and old
population, most importantly in the middle mass range.

The determination of ages tyoung, told, and tav is also affected, if a fixed
value RV = 3.1 is adopted (see Figure 11). The differences can be as large 50%
and 10% for the young and population, respectively, and 15% for the average
stellar age.

Figure 18:  Metallicities of star forming SDSS galaxies stacked in quintiles
of star formation rates at each stellar mass. Yellow, orange, blue, cyan,
violet colors correspond to the different quintiles from lowest to highest
SFR. Top: MZR of the young population, bottom: old population.

##  6 Spectral Normalization

The metallicities and ages obtained in the previous sections are B-band
luminosity averages because we have normalized our spectra by setting the mean
flux between 4400 and 4450 Å to unity. This is illustrated by eq. (11). It is,
therefore, important to investigate how metallicities and ages change
depending on the wavelength interval selected for normalization. For that
purpose we have repeated the analysis with normalizations between 5500 and
5550 Å (V-band) and 6950 and 7000 Å (R-band). The comparison of the
metallicities and ages obtained with the different normalizations are shown in
Figure 15.

The differences between the different normalizations are small. The
metallicities of the young and old population increase gradually with the
normalization shifted towards the red, but the maximum difference is $\sim$
0.1 dex only slightly larger than the metallicity uncertainties discussed in
section 4. Ages also increase by a small amount but the differences are
comparable to the errors shown in Figure 16.

We note that RV, E(B-V), AV and the star formation histories (as displayed in
Figure 12) remain unchanged and do not depend on the spectral normalization
interval.

Figure 19:  Top: average age tav (top figure) of the stellar population as a
function of stellar mass in different bins of star formation rate. Middle and
bottom: weight contributions of the young and old population byoung and bold.
The color code is the same as in Figure 18.

##  7 Initial Mass Function and Pre-Main Sequence Evolution

For our calculation of SSB model spectra we have adopted a Chabrier (2003)
initial mass function (IMF). As discussed by Conroy et al. (2009), there are
viable alternatives such as the IMFs by Kroupa (2001) and van Dokkum (2008).
In order to investigate the influence of the choice of the IMF on our analysis
we have additionally calculated model spectra with these two alternatives and
then repeated the analysis. The resulting MZRs are shown in Figure 16 for high
time resolution and tb = 0.1Myr. While the differences in metallicity arising
from the choices of the different IMFs are noticeable, they are at most 0.05
dex and, thus, not a major source of systematic uncertainties. We have also
included calculations with the original IMF suggested by Salpeter (1955),
which is a simple power law without a modification at the low mass end. Here
the differences are larger, 0.1 dex for the young and up to 0.2 dex for the
old stellar population, respectively. Figure 16 confirms the results found in
previous work (Cid Fernandes et al. 2005; Wilkinson et al. 2017).

The MESA stellar isochrones used in the calculation of our spectra include
pre-main sequence evolution towards the zero-age mains sequence (ZAMS). As a
consequence, young isochrones for which the more massive stars are leaving the
main sequence have a contribution from stars with lower mass, which are still
moving towards the main sequence and have luminosities significantly higher
than the ZAMS at the same stellar mass (see Choi et al. 2016 and Dotter 2016).
Because of the IMF these lower mass pre-main sequence objects could, thus,
influence the integrated isochrone spectra and lead to differences when
compared with isochrones just starting at the ZAMS. Since such isochrones have
also been used in population synthesis diagnostic work, we have carried out an
additional test with modified MESA isochrones, which start at the ZAMS. As it
turns out, the MZRs obtained with model spectra calculated in this way are
very similar to the ones where pre-main sequence evolution is included. The
reason is the high luminosity of the more massive star on the main sequence or
moving away, which have a dominating influence on the total spectrum. The
inclusion of the pre-main sequence phase is obviously of minor importance for
the results of our spectral diagnostics.

Figure 20:  The ratio of star formation rate of young to the old stellar
population as a function of stellar mass for the galaxies in different bins of
star formation rate. The color code is the same as in Figure 18 but only three
bins are shown for the sake of clarity.

##  8 A Comparison with Lookback Galaxy Evolution Models

Figure 7 and 8 indicate a strong evolution of stellar metallicities as a
function of age. The metallicities of the old population which is about 10 Gyr
older are significantly lower and the difference is anti-correlated with
galaxy stellar mass. This can be understood with the help of galaxy evolution
models using the standard framework of chemical evolution. Kudritzki et al.
(2021a, b) have introduced ’lookback’ models where the competing processes of
gas infall from the cosmic web and the circumgalactic medium, galactic winds
and star formation are assumed to lead to a redshift dependent power law
relationship between galactic gas mass and stellar mass, Mg =
A(z)M${}_{*}^{\beta}$ and stars are formed along the observed star formation
main sequences with star formation rates $\psi$ = $\psi_{0}$(z)Mδ. With these
assumptions and the usual chemical evolution equation metallicity [Z] becomes
a simple analytical function of the ratio of stellar to gas mass M∗/Mg and
decreases when this ratio becomes smaller. This is very similar to the
classical closed box model of chemical evolution (Searle & Sargent 1972) but
quantitatively the result is substantially different because of the power law
dependence between gas and stellar mass.

Figure 17 (top) shows the metallicities of the lookback evolution model for
each galaxy with different stellar mass for the young stars and stars 10 Gyr
older overplotted to the result of our SDSS spectral analysis and the
individual supergiants stars in galaxies. We note that these metallicities are
luminosity weighted in the B-band. Details of the model are described in
Kudritzki et al. (2021a, b). We note that $\beta$ = 0.5 is adopted and
$\psi_{0}$(z) has been increased by a factor of three to be in better
agreement with the star formation rates of the SDSS galaxies at redshift
z=0.08. The relative redshift dependence $\psi$(z) remains the same.

The agreement of the lookback models with the observations is good. The model
reproduces the shape of the young stellar population MZR well and also the
difference of the metallicities between the young stars and the old stars. We
note that the galaxy stellar masses at lookback time $\Delta$t = 10 Gyr are
significantly smaller than the final stellar masses. This is shown in the
bottom part of Figure 17. The difference in stellar mass together with
dependence of metallicity on M∗/Mg is the reason why the decline of the
metallicity of the old population with decreasing stellar mass is steeper than
for the younger population. Because of this mass difference the red curve in
Figure 17 (top) does not represent a mass metallicity relationship at high
redshift (z $\approx$ 2.5), but instead, describes the chemical evolution of
these galaxies.

##  9 Star Formation Rate Dependence of Metallicity and Ages

From the study of HII region emission lines it has been found that MZRs
contain a dependence on star formation rates as a third parameter (Mannucci et
al. 2010; Yates et al. 2012; Sanders et al. 2021). Galaxies with higher star
formation rates tend to have lower metallicities. Our analysis of stellar
spectra provides the opportunity to investigate this effect with respect to
stellar populations. For that purpose we follow the approach by Zahid et al.
(2017) and divide the sample of our 200,000 SDSS spectra in quintiles of star
formation rate. We then stack the spectra of these quintiles at each stellar
mass and run our analysis on the stacked spectra of the quintiles. Since the
spectra in the highest SFR quintile show strong HII emission lines for the
lowest galaxy stellar masses, we include the effects of nebular continuum
emission in our analysis as described in section 2. We use SSB with high age
resolution and a short burst length tb = 0.1Myr for the analysis.

The resulting MZRs for the young and old population are given in Figure 18. We
see a clear correlation with SFR for both the young population at high masses
and the old population at intermediate masses. The metallicity errors are of
the same order as shown in Figure 7 and described in section 4.

We have also determined population ages tyoung, told, and $t_{av}$. While
tyoung and told seem independent of star formation rate, the average age of
the total population tav as shown in Figure 19 reveals a very clear
dependence. The average population in galaxies with a lower star formation
rate is significantly older. Obviously, at the lower star formation rates it
takes the galaxies much longer to build up their stellar population. In
consequence, the contribution by the old stellar population bold to the
spectral fit is much higher, as is also shown in Figure 19. (Note that tav =
byoungtyoung \+ boldtold and that the second term dominates).

The stronger contribution of the older population is also reflected in the
star formation history expressed by the ratio $\psi_{young}$/ $\psi_{old}$ of
star formation rates of the young to the old population. As Figure 20
demonstrates, this ratio is significantly smaller for galaxies with presently
lower star formation rates.

##  10 Summary and Conclusions

In this section we summarize and discuss the major results of this work. In a
first step we have investigated the influence of star formation burst lengths
tb on the results of the population synthesis analysis. As observations show
(see section 3), tb depends on the size of the star forming region and can
range from a fraction of a Myr to many Myr. We find that the assumption of
burst length influences the spectra of stellar isochrones significantly
depending, of course, on isochrone age t and tb. For ages t with tb/t $\leq$
0.02 the signatures of burst length are clearly visible in the spectra. This
affects the results of the population synthesis analysis, in particular, the
metallicities of the young stellar population. Judging the quality of the
spectral fits by means of their minimum $\chi^{2}$ value we find that burst
spectra calculated with tb = 0.1 and 1.0 Myr give the best fit. Fits with 10
Myr bursts are worse.

Since we model the spectrum of the integrated stellar population as a sequence
of bursts with different ages, the time resolution of this sequence is
important. We have, therefore, tested the influence of the time resolution
using two grids with different resolution of isochrone age steps. We find that
the higher time resolution leads to lower $\chi^{2}$ values of the spectral
fits.

The analysis of the galaxy spectra with respect to metallicity and age
requires a determination of interstellar reddening and extinction. Since it is
well know that in the star forming regions of galaxies reddening laws deviate
from the standard law of the diffuse ISM characterized by the RV = 3.1, we
determine RV as an independent additional parameter. We find values clearly
larger than 3.1 indicating a strong contribution of dense star forming regions
to the reddening. The RV values are largest for low mass galaxies. For these
galaxies the average age of the young population is lowest. In consequence,
the contribution by stars originating from compact regions of star formation
might be strongest. Allowing for deviations of RV from 3.1 has a significant
influence on the determination of metallicity (up to 0.4 dex), population ages
(10 to 50% and, of course, interstellar extinction (0.2 magnitudes). We note
that this striking general difference from the standard reddening law in star
forming galaxies has important consequences for other aspects of astrophysics
such as, for instance, distance determinations and the cosmological distance
ladder or the estimate of galaxy luminosities.

The metallicity of the young population of our star forming SDSS galaxies is
in good agreement with massive supergiant stars where the metallicities are
obtained through accurate spectral analysis of a large sample of galaxies in
the local universe. We take this as a crucial confirmation that our population
synthesis analysis technique is reliable. The average age we derive for the
young population is between 50 and 300 Myrs increasing with galactic stellar
mass. It is, thus, in the same range as the supergiant ages and, consequently,
we expect the metallicities to be similar.

We also determine the average metallicities and ages of the old population.
The metallicities are significantly smaller than for the young population. The
age difference to the young population is about 10 Gyr. Thus, we are seeing
the results of chemical evolution over 10 Gyr. We compare these results with
galaxy evolution look back models and find good agreement. Both, the mass-
metallicity relationship of the young population and the metallicities of the
old population are well reproduced. This supports a galaxy evolution model
where metallicity depends on the ratio of gas mass to stellar mass and where
this ratio decreases with time.

For this comparison we calculated B-band luminosity weighted metallicities of
our models because our observed and model spectra are normalized to unity
between 4400 to 4450 Å and our analysis effectively provides B-band luminosity
weighted metallicities and ages. However, we have checked the effects of
normalization by introducing normalizations between 5500 and 5550 and 6950 and
7000 Å and found that they are small. Changes due to the different
normalizations are comparable to the fit uncertainties.

From the weight contributions to the total spectrum of isochrones at different
ages we can also estimate star formation rates. We calculate the ratio of star
formation rate of the young and old population and find an anti-correlation
with galaxy stellar mass. Low mass galaxies show a much higher star formation
contribution of the younger population. On the other hand, more massive
galaxies are dominated by stars that formed at early times. This is in
agreement with the observations of “galaxy downsizing” as originally detected
by Cowie et al. (1996) and followed up in comprehensive galaxy survey work
(see, for instance, Pilyugin & Thuan 2011; Thomas et al. 2019 and references
therein).

For the calculation of our model spectra we have assumed a Chabrier (2003)
initial mass function. In order to test the influence of this assumption we
repeat the analysis using the IMFs by Kroupa (2001), van Dokkum (2008) and
Salpeter (1955). The differences with respect to metallicity are very small
($\leq$ 0.05 dex) for the first three IMFs, but are larger for the Salpeter
IMF (0.1 dex and 0.2 dex for the young and old population, respectively).

Finally, we investigate the dependence of metallicity and ages on star
formation rate. We stack the spectra of the 200,000 SDSS galaxies in quintiles
of SFR at each stellar mass and repeat the analysis for each quintile. We find
a correlation of metallicity with SFR for higher and intermediate stellar
masses for the young and old population. In addition, the average age of the
total stellar population shows a clear dependence on SFR. The population in
galaxies with lower SFR is significantly older and the contribution of the old
population to the integrated spectrum is much larger. In the same way, the
ratio of SFR of the young to the old population is smaller for galaxies with
smaller SFR.

In summary, the results provide important implications for galaxy evolution
and underline the power of spectroscopic population synthesis analysis
techniques.

This work was initiated and supported by the Munich Excellence Cluster Origins
funded by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) under Germany’s Excellence Strategy EXC-2094 390783311. It is a
pleasure to thank our colleague Andi Burkert for inspiring and stimulating
discussion along this project.

## References

  * Arimoto & Yoshii (1987) Arimoto, N., & Yoshii, Y. 1987, A&A, 173, 23
  * Asari et al. (2007) Asari, N. V., Cid Fernandes, R., Stasińska, G., et al. 2007, MNRAS, 381, 263, doi: 10.1111/j.1365-2966.2007.12255.x
  * Berlind et al. (1997) Berlind, A. A., Quillen, A. C., Pogge, R. W., & Sellgren, K. 1997, AJ, 114, 107, doi: 10.1086/118457
  * Bresolin et al. (2022) Bresolin, F., Kudritzki, R.-P., & Urbaneja, M. A. 2022, ApJ, 940, 32, doi: 10.3847/1538-4357/ac9584
  * Bressan et al. (1996) Bressan, A., Chiosi, C., & Tantalo, R. 1996, A&A, 311, 425, doi: 10.48550/arXiv.astro-ph/9602032
  * Brinchmann et al. (2004) Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151, doi: 10.1111/j.1365-2966.2004.07881.x
  * Bruzual & Charlot (2003) Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000, doi: 10.1046/j.1365-8711.2003.06897.x
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245, doi: 10.1086/167900
  * Cardiel et al. (2003) Cardiel, N., Gorgas, J., Sánchez-Blázquez, P., et al. 2003, A&A, 409, 511, doi: 10.1051/0004-6361:20031096
  * Cardoso et al. (2022) Cardoso, L. S. M., Gomes, J. M., Papaderos, P., et al. 2022, A&A, 667, A11, doi: 10.1051/0004-6361/202243856
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763, doi: 10.1086/376392
  * Choi et al. (2016) Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102, doi: 10.3847/0004-637X/823/2/102
  * Cid Fernandes et al. (2005) Cid Fernandes, R., Mateus, A., Sodré, L., Stasińska, G., & Gomes, J. M. 2005, MNRAS, 358, 363, doi: 10.1111/j.1365-2966.2005.08752.x
  * Conroy & Gunn (2010) Conroy, C., & Gunn, J. E. 2010, ApJ, 712, 833, doi: 10.1088/0004-637X/712/2/833
  * Conroy et al. (2009) Conroy, C., Gunn, J. E., & White, M. 2009, ApJ, 699, 486, doi: 10.1088/0004-637X/699/1/486
  * Corbelli et al. (2017) Corbelli, E., Braine, J., Bandiera, R., et al. 2017, A&A, 601, A146, doi: 10.1051/0004-6361/201630034
  * Cowie et al. (1996) Cowie, L. L., Songaila, A., Hu, E. M., & Cohen, J. G. 1996, AJ, 112, 839, doi: 10.1086/118058
  * Davé et al. (2011) Davé, R., Finlator, K., & Oppenheimer, B. D. 2011, MNRAS, 416, 1354, doi: 10.1111/j.1365-2966.2011.19132.x
  * Davies et al. (2017) Davies, B., Kudritzki, R.-P., Lardo, C., et al. 2017, ApJ, 847, 112, doi: 10.3847/1538-4357/aa89ed
  * Dayal et al. (2013) Dayal, P., Ferrara, A., & Dunlop, J. S. 2013, MNRAS, 430, 2891, doi: 10.1093/mnras/stt083
  * Dotter (2016) Dotter, A. 2016, ApJS, 222, 8, doi: 10.3847/0067-0049/222/1/8
  * Efremov & Elmegreen (1998) Efremov, Y. N., & Elmegreen, B. G. 1998, MNRAS, 299, 588, doi: 10.1046/j.1365-8711.1998.01819.x
  * Falco et al. (1999) Falco, E., Impey, C., Kochanek, C., et al. 1999, ApJ, 523, 617, doi: 10.1086/307758
  * Feinstein et al. (1973) Feinstein, A., Marraco, H. G., & Muzzio, J. C. 1973, A&AS, 12, 331
  * Gallazzi et al. (2005) Gallazzi, A., Charlot, S., Brinchmann, J., White, S. D. M., & Tremonti, C. A. 2005, MNRAS, 362, 41, doi: 10.1111/j.1365-2966.2005.09321.x
  * Garn & Best (2010) Garn, T., & Best, P. N. 2010, MNRAS, 409, 421, doi: 10.1111/j.1365-2966.2010.17321.x
  * Gazak et al. (2015) Gazak, J. Z., Kudritzki, R., Evans, C., et al. 2015, ApJ, 805, 182, doi: 10.1088/0004-637X/805/2/182
  * Genzel et al. (2015) Genzel, R., Tacconi, L. J., Lutz, D., et al. 2015, ApJ, 800, 20, doi: 10.1088/0004-637X/800/1/20
  * González Delgado et al. (2015) González Delgado, R. M., García-Benito, R., Pérez, E., et al. 2015, A&A, 581, A103, doi: 10.1051/0004-6361/201525938
  * Gordon et al. (1997) Gordon, K. D., Calzetti, D., & Witt, A. N. 1997, ApJ, 487, 625, doi: 10.1086/304654
  * Herbst (1976) Herbst, W. 1976, ApJ, 208, 923, doi: 10.1086/154681
  * Holwerda & Keel (2013) Holwerda, B. W., & Keel, W. C. 2013, A&A, 556, A42, doi: 10.1051/0004-6361/201321455
  * Kewley & Ellison (2008) Kewley, L. J., & Ellison, S. L. 2008, ApJ, 681, 1183, doi: 10.1086/587500
  * Kim et al. (2021) Kim, J., Chevance, M., Kruijssen, J. M. D., et al. 2021, MNRAS, 504, 487, doi: 10.1093/mnras/stab878
  * Koleva et al. (2009) Koleva, M., Prugniel, P., Bouchard, A., & Wu, Y. 2009, A&A, 501, 1269, doi: 10.1051/0004-6361/200811467
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231, doi: 10.1046/j.1365-8711.2001.04022.x
  * Kudritzki et al. (2016) Kudritzki, R. P., Castro, N., Urbaneja, M. A., et al. 2016, ApJ, 829, 70, doi: 10.3847/0004-637X/829/2/70
  * Kudritzki et al. (2021a) Kudritzki, R.-P., Teklu, A. F., Schulze, F., et al. 2021a, ApJ, 922, 274, doi: 10.3847/1538-4357/ac32cf
  * Kudritzki et al. (2021b) —. 2021b, ApJ, 922, 274, doi: 10.3847/1538-4357/ac32cf
  * Lapi et al. (2020) Lapi, A., Pantoni, L., Boco, L., & Danese, L. 2020, ApJ, 897, 81, doi: 10.3847/1538-4357/ab9812
  * Lardo et al. (2015) Lardo, C., Davies, B., Kudritzki, R. P., et al. 2015, ApJ, 812, 160, doi: 10.1088/0004-637X/812/2/160
  * Lequeux et al. (1979) Lequeux, J., Peimbert, M., Rayo, J. F., Serrano, A., & Torres-Peimbert, S. 1979, A&A, 500, 145
  * Leroy et al. (2008) Leroy, A. K., Walter, F., Brinks, E., et al. 2008, AJ, 136, 2782, doi: 10.1088/0004-6256/136/6/2782
  * Maíz Apellániz et al. (2017) Maíz Apellániz, J., Trigueros Páez, E., Bostroem, A. K., Barbá, R. H., & Evans, C. J. 2017, in Highlights on Spanish Astrophysics IX, 510–510
  * Mannucci et al. (2010) Mannucci, F., Cresci, G., Maiolino, R., Marconi, A., & Gnerucci, A. 2010, MNRAS, 408, 2115, doi: 10.1111/j.1365-2966.2010.17291.x
  * O’Donnell (1994) O’Donnell, J. E. 1994, ApJ, 422, 158, doi: 10.1086/173713
  * Panter et al. (2008) Panter, B., Jimenez, R., Heavens, A. F., & Charlot, S. 2008, MNRAS, 391, 1117, doi: 10.1111/j.1365-2966.2008.13981.x
  * Pantoni et al. (2019) Pantoni, L., Lapi, A., Massardi, M., Goswami, S., & Danese, L. 2019, ApJ, 880, 129, doi: 10.3847/1538-4357/ab2adc
  * Parikh et al. (2021) Parikh, T., Thomas, D., Maraston, C., et al. 2021, MNRAS, 502, 5508, doi: 10.1093/mnras/stab449
  * Peng et al. (2015) Peng, Y., Maiolino, R., & Cochrane, R. 2015, Nature, 521, 192, doi: 10.1038/nature14439
  * Pilyugin & Thuan (2011) Pilyugin, L. S., & Thuan, T. X. 2011, ApJ, 726, L23, doi: 10.1088/2041-8205/726/2/L23
  * Reichardt et al. (2001) Reichardt, C., Jimenez, R., & Heavens, A. F. 2001, MNRAS, 327, 849, doi: 10.1046/j.1365-8711.2001.04768.x
  * Salpeter (1955) Salpeter, E. E. 1955, ApJ, 121, 161, doi: 10.1086/145971
  * Sánchez-Blázquez et al. (2011) Sánchez-Blázquez, P., Ocvirk, P., Gibson, B. K., Pére, I., & Peletier, R. F. 2011, MNRAS, 415, 709, doi: 10.1111/j.1365-2966.2011.18749.x
  * Sánchez-Blázquez et al. (2006) Sánchez-Blázquez, P., Peletier, R. F., Jiménez-Vicente, J., et al. 2006, MNRAS, 371, 703, doi: 10.1111/j.1365-2966.2006.10699.x
  * Sanders et al. (2021) Sanders, R. L., Shapley, A. E., Jones, T., et al. 2021, ApJ, 914, 19, doi: 10.3847/1538-4357/abf4c1
  * Schaefer et al. (2017) Schaefer, A. L., Croom, S. M., Allen, J. T., et al. 2017, MNRAS, 464, 121, doi: 10.1093/mnras/stw2289
  * Schaye et al. (2015) Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521, doi: 10.1093/mnras/stu2058
  * Searle & Sargent (1972) Searle, L., & Sargent, W. L. W. 1972, ApJ, 173, 25, doi: 10.1086/151398
  * Spinrad & Taylor (1972) Spinrad, H., & Taylor, B. J. 1972, ApJ, 171, 397, doi: 10.1086/151291
  * Spitoni et al. (2020) Spitoni, E., Calura, F., Mignoli, M., et al. 2020, A&A, 642, A113, doi: 10.1051/0004-6361/202037879
  * Spitoni et al. (2017a) Spitoni, E., Gioannini, L., & Matteucci, F. 2017a, A&A, 605, A38, doi: 10.1051/0004-6361/201730545
  * Spitoni et al. (2017b) Spitoni, E., Vincenzo, F., & Matteucci, F. 2017b, A&A, 599, A6, doi: 10.1051/0004-6361/201629745
  * Stark & Parker (1995) Stark, P., & Parker, R. 1995, Computational Statistics, 10, 129–141
  * The et al. (1980) The, P. S., Bakker, R., & Tjin A Djie, H. R. E. 1980, A&A, 89, 209
  * Thomas et al. (2019) Thomas, R., Le Fèvre, O., Zamorani, G., et al. 2019, A&A, 630, A145, doi: 10.1051/0004-6361/201935813
  * Tinsley (1968) Tinsley, B. M. 1968, ApJ, 151, 547, doi: 10.1086/149455
  * Tremonti et al. (2004) Tremonti, C. A., Heckman, T. M., Kauffmann, G., et al. 2004, ApJ, 613, 898, doi: 10.1086/423264
  * Trussler et al. (2020) Trussler, J., Maiolino, R., Maraston, C., et al. 2020, MNRAS, 491, 5406, doi: 10.1093/mnras/stz3286
  * Urbaneja et al. (2017) Urbaneja, M. A., Kudritzki, R.-P., Gieren, W., et al. 2017, AJ, 154, 102, doi: 10.3847/1538-3881/aa79a8
  * Vale Asari et al. (2009) Vale Asari, N., Stasińska, G., Cid Fernandes, R., et al. 2009, MNRAS, 396, L71, doi: 10.1111/j.1745-3933.2009.00664.x
  * van Dokkum (2008) van Dokkum, P. G. 2008, ApJ, 674, 29, doi: 10.1086/525014
  * Wilkinson et al. (2017) Wilkinson, D. M., Maraston, C., Goddard, D., Thomas, D., & Parikh, T. 2017, MNRAS, 472, 4297, doi: 10.1093/mnras/stx2215
  * Witt & Gordon (1996) Witt, A. N., & Gordon, K. D. 1996, ApJ, 463, 681, doi: 10.1086/177282
  * Witt & Gordon (2000) —. 2000, ApJ, 528, 799, doi: 10.1086/308197
  * Yates et al. (2012) Yates, R. M., Kauffmann, G., & Guo, Q. 2012, MNRAS, 422, 215, doi: 10.1111/j.1365-2966.2012.20595.x
  * Yoachim et al. (2012) Yoachim, P., Roškar, R., & Debattista, V. P. 2012, ApJ, 752, 97, doi: 10.1088/0004-637X/752/2/97
  * Zahid et al. (2014) Zahid, H. J., Dima, G. I., Kudritzki, R.-P., et al. 2014, ApJ, 791, 130, doi: 10.1088/0004-637X/791/2/130
  * Zahid et al. (2017) Zahid, H. J., Kudritzki, R.-P., Conroy, C., Andrews, B., & Ho, I. T. 2017, ApJ, 847, 18, doi: 10.3847/1538-4357/aa88ae
  * Zahid et al. (2013) Zahid, H. J., Yates, R. M., Kewley, L. J., & Kudritzki, R. P. 2013, ApJ, 763, 92, doi: 10.1088/0004-637X/763/2/92