# Toward a Better Understanding of Cosmic Chronometers:
Stellar Population Properties of Passive Galaxies at Intermediate Redshift

Nicola Borghi Dipartimento di Fisica e Astronomia “Augusto Righi”–Università
di Bologna, via Piero Gobetti 93/2, I-40129 Bologna, Italy  INAF—Osservatorio
di Astrofisica e Scienza dello Spazio di Bologna, via Piero Gobetti 93/3,
I-40129 Bologna, Italy  Michele Moresco Dipartimento di Fisica e Astronomia
“Augusto Righi”–Università di Bologna, via Piero Gobetti 93/2, I-40129
Bologna, Italy  INAF—Osservatorio di Astrofisica e Scienza dello Spazio di
Bologna, via Piero Gobetti 93/3, I-40129 Bologna, Italy  Andrea Cimatti
Dipartimento di Fisica e Astronomia “Augusto Righi”–Università di Bologna, via
Piero Gobetti 93/2, I-40129 Bologna, Italy  INAF—Osservatorio Astrofisico di
Arcetri, Largo E. Fermi 5, I-50125, Firenze, Italy  Alexandre Huchet Faculté
des sciences d’Orsay—Université Paris-Saclay, 15 rue Georges Clemenceau,
F-91405 Orsay cedex I, France  Salvatore Quai Department of Physics and
Astronomy, University of Victoria, 3800 Finnerty Road, Victoria, BC V8P 5C2,
Canada  Lucia Pozzetti INAF—Osservatorio di Astrofisica e Scienza dello Spazio
di Bologna, via Piero Gobetti 93/3, I-40129 Bologna, Italy Nicola Borghi
<EMAIL_ADDRESS>
(Accepted October 20, 2021)

###### Abstract

We take advantage of the publicly available LEGA-C spectroscopic survey to
measure the stellar population properties of 140 individual massive and
passive galaxies at $z\sim 0.7$. We develop and publicly release PyLick, a
flexible Python code to measure UV to near-IR spectral indices. With PyLick we
study the $\mathrm{H/K}$ ratio as a new diagnostic based on the pseudo-Lick Ca
II H and K indices and find that a cut in $\mathrm{H/K}<1.1$ can be used
jointly with other criteria to select (or verify the purity of) samples of
passive galaxies. By combining photometric and spectroscopic criteria, we
select a reliable sample of passively evolving galaxies. We constrain single-
burst stellar ages, metallicities $\mathrm{[Z/H]}$, and $\mathrm{[\alpha/Fe]}$
with an optimized set of Lick indices, exploring in detail the robustness of
our measurement against different combinations. Even without imposing
cosmological priors, the derived ages follow a clear trend compatible with the
expected cosmological aging of the Universe. We observe no significant
redshift evolution for the metal abundance with respect to the values derived
at $z=0$, with median $\mathrm{[Z/H]}=0.08\pm 0.18$ and
$\mathrm{[\alpha/Fe]}=0.13\pm 0.11$. Finally, we analyze the relations between
$\log\mathrm{age}$, $\mathrm{[Z/H]}$, $\mathrm{[\alpha/Fe]}$ and the stellar
velocity dispersion, finding slopes of ($0.5\pm 0.1$), ($0.3\pm 0.2$), and
($0.2\pm 0.1$), respectively; the small scatter of $<0.2$ dex points to rather
homogeneous and short star formation histories. Overall, these results confirm
and extend low-redshift findings of a mass-downsizing evolution. This work
further strengthens the possibility of selecting pure samples of passive
galaxies to be exploited reliably as cosmic chronometers to place independent
cosmological constraints.

Galaxy evolution (594), Galaxy abundances (574), Galaxy ages (576),
Observational cosmology (1146)

††software: Astropy (Astropy Collaboration et al., 2018); ChainConsumer
(Hinton, 2016); emcee (Foreman-Mackey et al., 2013); LOESS (Cappellari et al.,
2013a); LtsFit (Cappellari et al., 2013b); Matplotlib (Hunter, 2007); Numpy
(Harris et al., 2020); Scipy (Virtanen et al., 2020); and Topcat (Taylor,
2005).

monthyearday []

##  1 Introduction

The advent of deep spectroscopic surveys led to considerable progress in
understanding galaxy formation and evolution over different cosmic epochs. The
reason is that they enable us to combine two main observational techniques:
look-back statistical studies of galaxy properties over different cosmic
times, and the archaeological reconstruction of galaxy properties and star
formation histories (SFHs) from spectroscopic data.

Galaxies with no or negligible levels of star formation (hereafter “passive”
or “quiescent”) are ideal laboratories to perform these studies because their
stellar population is relatively simple to model (for a detailed review, see
Renzini, 2006). They are mostly found after the peak epoch of galaxy assembly
(the so-called ‘cosmic noon’), especially below $z\sim 1.5$, when they
dominate the stellar-mass density (Muzzin et al., 2013; Ilbert et al., 2013),
and the number density of the most massive ones ($M_{\star}\gtrsim
10^{11}M_{\odot}$) remains almost constant between $z=1$ and $z=0$ (Pozzetti
et al., 2010; Moresco et al., 2013). On the contrary, at higher redshift, they
constitute only a minor population (Cimatti et al., 2004; Daddi et al., 2005).
Recent spectroscopic observations confirmed the existence of a few of these
systems up to $z\sim 4$ (Tanaka et al., 2019; Valentino et al., 2020; Santini
et al., 2021) when the Universe was only $\sim 2$ Gyr old, requiring an early
intense star formation followed by a complete decline (known as quenching).

Many different methods are used to investigate their physical properties. An
example is the study of the photometric spectral energy distribution (SED;
e.g., Pacifici et al. 2016). Photometric observations can span a very large
range of wavelengths (usually from UV to FIR), allowing a broad view of
various physical quantities, such as stellar mass, star formation rate,
stellar age, and metallicity. However, it is not possible to find precise
constraints due to strong degeneracies, e.g., age-metallicity (Worthey et al.,
1994). To break them, high signal-to-noise spectroscopy is needed. The well-
established method of Lick indices (Faber, 1973; Burstein et al., 1984;
Worthey et al., 1994) involves the use of a set of absorption features, each
one having a unique sensitivity to stellar ages and element abundances.
Another main approach is based on the use of the entire spectral information
(full spectral fitting; e.g., Cid Fernandes et al., 2005; Conroy, 2013). All
of these techniques independently contributed to establishing the downsizing
scenario (first introduced by Cowie et al., 1996), according to which galaxy
mass plays a major role in galaxy formation and evolution: massive galaxies
are found to have evolved earlier and over shorter time scales than less
massive ones. In the local Universe, this is supported by positive scaling
relations between stellar ages, metallicities ($\rm[Z/H]$), and
$\alpha$-element abundances ($\rm[\alpha/Fe]$) with galaxy mass (dynamical and
stellar) found in large, high-quality spectroscopic samples of early-type
galaxies (Kauffmann et al., 2003; Gallazzi et al., 2005, 2006; Thomas et al.,
2005, 2010; Treu et al., 2005; Conroy et al., 2014; McDermid et al., 2015;
Scott et al., 2017; Siudek et al., 2017; Gallazzi et al., 2021).

The stellar metallicity, i.e. the amount of metals locked in stars, can shed
light on the evolutionary stage of the stellar population, as well as on the
external mechanisms that modify the chemical content of the interstellar
medium. For instance, Peng et al. (2015) showed that local quiescent and star-
forming galaxies form two distinct relations in the stellar-mass–metallicity
plane. The difference between their shape can be interpreted either as an
evidence for ‘strangulation’ (i.e. the lack of new gas supply Peng et al.,
2015; Trussler et al., 2021), or for shorter formation timescales coupled with
galactic winds (Spitoni et al., 2017) as the main mechanism driving galaxy
quenching. The relative $\alpha$-element abundance with respect to iron is
another key parameter. Core-collapsing massive stars are the main producers of
$\alpha$-elements (O, Mg, Si, Ca, Ti), polluting the interstellar medium over
relatively short timescales ($<100$ Myr). On the other hand, iron-peak
elements (Fe, Cr) are primarily produced in Type Ia supernovae, which pollute
the interstellar medium over longer timescales (100 Myr–2 Gyr). For this
reason, the mean stellar $\rm[\alpha/Fe]$ has been traditionally used as a SFH
timescale diagnostic (Thomas et al., 2005; de La Rosa et al., 2011). In
addition to stellar population properties, the environment can in principle
have a role in shaping galaxy evolution. However, its effects seem weak once
the correlation with mass is removed (Thomas et al., 2010; La Barbera et al.,
2014; McDermid et al., 2015; Trussler et al., 2021). To explore them in
greater detail, very large samples of galaxies are needed (e.g., Bluck et al.,
2020).

At higher redshift, studies of stellar population properties are more
challenging and require deep near-infrared spectroscopy. For this reason, they
are mostly limited to samples of a few up to dozens galaxies (e.g., Jørgensen
& Chiboucas, 2013; Gallazzi et al., 2014; Lonoce et al., 2015, 2020; Belli et
al., 2019; Carnall et al., 2019; Kriek et al., 2019; Tacchella et al., 2021;
Beverage et al., 2021), and/or require the stacking of different galaxy
spectra (e.g., Choi et al., 2014; Onodera et al., 2015). In particular,
detailed studies of galaxy ages and chemical abundances for individual
galaxies at intermediate redshift have been presented in very few works and
usually by assuming the age of a $\Lambda$CDM universe as the maximum age
allowed in the stellar population analysis.

Beyond galaxy evolution studies, massive and passive galaxies encode valuable
information about the underlying cosmological framework. In fact, as first
proposed by Jimenez & Loeb (2002), it has been demonstrated that these objects
can be used as cosmic chronometers to trace the differential age evolution of
the Universe $\mathrm{d}t/\mathrm{d}z$, and to provide in this way a
cosmology-independent estimate of the expansion rate of the Universe
$H(z)=-1/(1+z)\,\mathrm{d}z/\mathrm{d}t$. While we refer to the literature for
a detailed discussion of the method and of the systematics involved (Moresco
et al., 2012b; Moresco, 2015; Moresco et al., 2016, 2020), we underline here
that there are two key ingredients to be met. First, the best cosmic
chronometers that have been studied are very massive and passively evolving
galaxies; therefore, it is fundamental to accurately select these objects by
carefully excluding any possible hints of star-forming or young outliers.
Second, it is crucial that in the differential age estimate $\mathrm{d}t$ no
cosmological assumption is made. With the advantage of providing a cosmology-
independent estimate of $H(z)$, the cosmic chronometers method has also been
considered in several cosmological studies to place constraints on various
cosmological models and parameters (Moresco et al., 2012a; Seikel et al.,
2012; Capozziello et al., 2014; Valkenburg et al., 2014; Sapone et al., 2014;
Nunes et al., 2016; Solà et al., 2017; Moresco & Marulli, 2017; L’Huillier &
Shafieloo, 2017; Yang et al., 2018; Haridasu et al., 2018; Gómez-Valent &
Amendola, 2018; Lin et al., 2020, 2021), with particular benefits over more
standard cosmological probes (see e.g., Vagnozzi et al., 2021).

In this work, we take advantage of the deep spectroscopic observations of the
second data release of the Large Early Galaxy Astrophysics Census (LEGA-C DR2;
van der Wel et al. 2016; Straatman et al. 2018) at $0.6\lesssim z\lesssim 1$
to infer physical properties of individual passive galaxies without relying on
any cosmological model or assumption. Our analysis is based on Pylick, a
flexible Python tool to measure absorption features that we publicly release.
It includes a wide set of indices already defined in the literature and also
new diagnostics introduced and explored in this work for identifying passively
evolving systems. Stellar population properties of each individual galaxy are
derived with a Bayesian approach, adopting the simple stellar population (SSP)
models of Thomas et al. (2011). The analysis of their trends with redshift and
stellar velocity dispersion will allow us to understand the individual and
median properties of passive galaxies over different cosmic epochs and explore
the underlying cosmology.

The work is organized as follows. Section 2 gives an overview of the dataset,
selection process, spectral index measurements, and main observational
properties. Section 3 presents background information on the models and the
stellar population analysis. The main science results are presented and
discussed in Section 4. A summary is presented in Section 5.

In some cases, as reference values and mostly for illustrative purposes, we
will use some theoretical relations based on cosmological models; for these
cases, we adopt a ‘737 cosmology’ (with $H_{0}=70$ km s-1 Mpc-1,
$\Omega_{m}=0.3$, and $\Omega_{\Lambda}=0.7$).

##  2 Passive galaxy sample

###  2.1 Data

The data used in this study are drawn from the second data release of LEGA-C,
a recently completed ESO Public Spectroscopic Survey targeting $\sim$3000
$K_{s}$-band-selected galaxies at $0.6\lesssim z\lesssim 1$ in the COSMOS
field. The observations were carried out with the Visible Multi-Object
Spectrograph (VIMOS) on the VLT at Paranal Observatory. The flux-calibrated
spectra span a wavelength range of $6300<\lambda<8800$ Å with a spectral
resolution of $R\sim 3500$, and a median signal-to-noise ratio (S/N) of $\sim
15$ per pixel (0.6 Å). Spectra are obtained with $1\arcsec$ wide slits,
corresponding to $\sim 7$ kpc at these redshifts.

A catalog of spectroscopic measurements for DR2 has been publicly released,
comprising spectroscopic redshift, flux measurements for the main emission
lines (Balmer lines, [O ii]$\lambda$3727, [O iii]$\lambda$4959, [O
iii]$\lambda$5007), velocity dispersion determinations, as well as
measurements for a sample of 14 Lick/IDS indices (Straatman et al., 2018). In
our analysis, we consider the galaxy spectra, as well as the measurements of
redshift ($z$), observed stellar velocity dispersion ($\sigma_{\star}$), and
[O ii]$\lambda$3727 emission-line flux. For the main analysis, we do not use
the spectral indices measurements provided in the catalog, but instead
determine our own line strengths from the spectra after matching the stellar
population models resolution. This also allows us to extend the measurements
to redder indices up to $\sim 5000$ Å (see Section 2.3). The LEGA-C set of
Lick indices will be used to validate our measurement code PyLick in A.

We cross-match the LEGA-C sample with the COSMOS2015 catalog (Laigle et al.,
2016) using a search radius of $1\arcsec$, to complement the spectroscopic
information with photometric data, including bands NUV, $r$, and $J$, as well
as stellar masses $M_{\star}$ and (specific) star formation rates (${\rm
sSFR=SFR}/M_{\star}$) derived through SED fitting. Finally, we use
morphological information from the Zurich Estimator of Structural Types
catalog (Scarlata et al., 2007), based on principal component analysis of the
surface brightness profiles. Our parent sample is selected requiring good
quality spectra (see Straatman et al. 2018 for spectra quality flags) and
available NUV, $r$, and $J$ absolute magnitudes. In this way, we end up with
1622 sources.

Figure 1: Distribution of the main properties of subsequently refined samples
of passive galaxies in the LEGA-C survey. Upper panels: NUVrJ and EW[O
ii]$-$$\mathrm{D_{n}4000}$ diagnostic diagrams with selection criteria adopted
in Ilbert et al. (2013) and Mignoli et al. (2009), respectively. We use black
colors for the criteria also adopted in this work, while gray lines are for
illustrative purposes only. The hatched region ($\pm 0.05$ mag with respect to
the NUVrJ cut) indicates the location of green-valley galaxies (Davidzon et
al., 2017). Note that only 65/350 passive galaxies have detected [O
ii]$\lambda$3727 line. Lower panels: distribution in redshift ($z$), observed
stellar velocity dispersion ($\sigma_{\star}$), stellar mass ($M_{\star}$),
and star formation rate (SFR). Arrows represent the median values of the
various subsamples. Table 1: Median Properties of the Various Subsamples of
Passive Galaxies Defined

Sample | $N$ | $\left\langle\right.\\!\\!{z}\\!\\!\left.\right\rangle$ | $\left\langle\right.\\!\\!{\sigma_{\star}}\\!\\!\left.\right\rangle$ | $\left\langle\right.\\!\\!{\log M_{\star}}\\!\\!\left.\right\rangle$ | $\left\langle\right.\\!\\!{\log{\rm SFR}}\\!\\!\left.\right\rangle$ | $\left\langle\right.\\!\\!{\rm H/K}\\!\\!\left.\right\rangle$ | $\left\langle\right.\\!\\!{\mathrm{D_{n}4000}}\\!\\!\left.\right\rangle$ | $\left\langle\right.\\!\\!{\mathrm{H\delta_{A}}}\\!\\!\left.\right\rangle$ | $\left\langle\right.\\!\\!{\mathrm{G4300}}\\!\\!\left.\right\rangle$ | $\left\langle\right.\\!\\!{\mathrm{Fe4383}}\\!\\!\left.\right\rangle$ | $\left\langle\right.\\!\\!{\mathrm{C_{2}4668}}\\!\\!\left.\right\rangle$
---|---|---|---|---|---|---|---|---|---|---|---
|  |  | $\rm km\,s^{-1}$ | $M_{\odot}$ | $M_{\odot}\,\mathrm{yr^{-1}}$ |  |  | Å | Å | Å | Å
parent | 1622 | 0.752 | 165 | 10.75 | 0.50 | 1.166 | 1.426 | 2.941 | 2.357 | 2.583 | 3.256
photometric | 658 | 0.758 | 201 | 10.98 | -0.80 | 0.969 | 1.654 | 0.692 | 4.286 | 3.594 | 4.935
spectrophotometric | 485 | 0.732 | 202 | 10.95 | -1.30 | 0.957 | 1.681 | 0.462 | 4.387 | 3.741 | 4.996
bona fide passive | 350 | 0.735 | 206 | 10.95 | -1.46 | 0.957 | 1.690 | 0.347 | 4.443 | 3.832 | 5.097

Note. — Data from: LEGA-C DR2 (3–4, Straatman et al. 2018); COSMOS2015 (5–6,
Laigle et al. 2016); this work (7–12, Section 2.3).

###  2.2 Selection criteria

In the literature, multiple methods have been proposed to separate ‘passive’
from ‘star-forming’ galaxies, including a morphological selection of
spheroidal systems (following the original separation by Hubble, 1936), cuts
on color–color diagrams (e.g., UVJ, Williams et al. 2009; NUVrJ, Ilbert et al.
2013) or on a color–mass diagram (e.g., Peng et al., 2010), SED fitting (e.g.,
Ilbert et al., 2009), and sSFR criteria (e.g., Pozzetti et al., 2010).

However, these different criteria do not perfectly overlap (Renzini, 2006). In
fact, selections based on a single criterion are not stringent enough to
reduce the contamination from star-forming outliers, with a percentage of
contamination up to 10%–30% depending on the considered criterion. On the
other hand, a combination of different criteria, maximizing the overlap of
complementary information (photometric and spectroscopic), is significantly
more effective in selecting a pure sample (Franzetti et al., 2007; Moresco et
al., 2013). In this analysis, with the aim of studying the physical properties
of cosmic chronometers, we are interested in having the purest possible
sample, minimizing as much as possible the eventual residual contamination
from star-forming outliers. We decided to combine different complementary
cuts, as outlined below.

  * •

NUVrJ selection. Photometric passive galaxies are selected using the rest-
frame $\rm{NUV}-r)$ and $(r-J)$ colors following the criterion proposed by
Ilbert et al. (2013): $(\mathrm{NUV}-r)>3~{}(r-J)+1$ and
$(\mathrm{NUV}-r)>3.1$. These colors have been demonstrated to be extremely
sensitive to reveal objects with recent (1–100 Myr) star formation episodes
(blue $\mathrm{NUV}-r$ colors) even if they are dust-obscured (redder $r-J$
colors), and are therefore optimized to safely separate quiescent and star-
forming galaxies (Arnouts et al., 2007; Ilbert et al., 2015; Davidzon et al.,
2017). With this cut, 658 sources are selected. We refer to these as the
photometric passive sample.

  * •

Emission-line cut. We further restrict our sample by excluding galaxies with a
significant [O ii]$\lambda$3727 emission line, a tracer of photoionized gas,
which is typically considered an indicator of ongoing star formation111Low-
ionization nuclear emission-line regions (LINERs) and ionization from old
stars can also be responsible for [O ii] and [O iii] emission-lines (e.g., Yan
et al., 2006; Singh et al., 2013; Cimatti et al., 2019), but because we aim to
sample purity we still exclude these sources.. In particular, we exclude those
galaxies that have an $\text{EW[O\>{ii}]}>5$ Å. The threshold is found to
separate well star-forming and passive galaxies in previous spectroscopic
surveys at similar redshift (e.g., Mignoli et al., 2009). Combining this and
the previous cut, we obtain 485 galaxies that we refer to as the
spectrophotometric passive sample.

  * •

Visual inspection. The sample is further refined by visually inspecting all
the remaining spectra. We remove galaxies with clearly strong [O
ii]$\lambda$3727 and/or [O iii]$\lambda$5007 lines, obtaining typical
$\mathrm{S/N}<3$ in their EWs. The latter is crucial to spectroscopically
characterize galaxies at $z\lesssim 0.65$ for which [O ii] is not available in
LEGA-C spectra. Based upon these criteria, we obtain a final sample of 350
bona fide passive galaxies.

In Figure 1 we show the distribution of LEGA-C galaxies in two diagnostic
diagrams (NUVrJ, Ilbert et al. 2013; EW[O ii]$-$$\mathrm{D_{n}4000}$ Mignoli
et al. 2009), redshift, $\sigma_{\star}$, $M_{\star}$, and SFR, across the
parent, photometric, spectrophotometric and bona fide passive subsamples. The
median values for the main quantities are quoted in Table 1. Note that objects
without [OII] detection are not displayed in the upper-right diagram.

Overall, the LEGA-C galaxies present two distinctive peaks in the redshift
distribution at $z\sim 0.7$ and $z\sim 0.9$, with very few galaxies at $z>1$.
They form two separate populations in the NUVrJ plot, with a blue sequence
reaching low $(\mathrm{NUV}-r)\sim 1.1$ colors and a red cloud that
constitutes the photometric passive sample. Although a NUVrJ-only criterion
drastically reduces the presence of star-forming systems in the sample
($\left\langle\right.\\!\\!{\rm\log sSFR/yr}\\!\\!\left.\right\rangle=-11.8$),
about one-third of objects still have a significant [O ii] emission.
Therefore, a spectroscopic selection, here performed by combining a cut on the
EW followed by a careful inspection of all the remaining [O ii] and/or [O iii]
contributions, is fundamental to ensure the purity of the sample. It is
interesting to note that this does not only remove the tail of bluer
$(\mathrm{NUV}-r)$ galaxies associated with the green-valley region, but also
systems with redder colors. The final bona fide passive sample has a median
redshift of $\left\langle\right.\\!\\!{z}\\!\\!\left.\right\rangle=0.735$.
Passive galaxies are located toward the high-$\sigma_{\star}$ and $M_{\star}$
tails of the parent distribution. In particular, the median $\sigma_{\star}$
($\log M_{\star}/M_{\odot}$) increases from 164.5 $\rm km\,s^{-1}$ (10.75) to
205.7 $\rm km\,s^{-1}$ (10.95) moving from the parent to the bona fide passive
sample, and most of the bona fide passive galaxies ($85\%$) have $\log
M_{\star}/M_{\odot}>10.6$. With respect to the spectrophotometric, this sample
has a SFR lower by $0.16$ dex (with a median uncertainty of 0.18 dex).
Finally, we note that the passive sample has a median specific star formation
rate of $\left\langle\right.\\!\\!{\rm\log
sSFR/yr}\\!\\!\left.\right\rangle=-12.1$, with only 15 galaxies ($\sim
4~{}\%$) reaching $>-11$, a value commonly adopted to classify “passive”
galaxies (see Pozzetti et al., 2010).

Figure 2: Median composite spectrum of the 350 selected passive galaxies.
Upper panel: number of stacked galaxies at each wavelength. Lower panel:
Median composite spectrum and associated 1$\sigma$ error. Gray boxes show
central regions of the 22 Lick and pseudo-Lick indices measured in this work;
narrower indices in overlapping regions ($\mathrm{H\delta_{F}}$ ,
$\mathrm{H\gamma_{F}}$ , and $\mathrm{Mg}\,b$ ) are shown with darker shading.
Horizontal segments identify blue and red bandpasses where
$\mathrm{D_{n}4000}$ is computed.

We measure the composite spectrum of the bona fide passive sample by
normalizing each rest-frame spectrum to the median flux between 4200 and 4400
Å, and interpolating it onto a common grid (3000–5500 Å,
$\Delta\lambda=0.35~{}\text{\AA}~{}\rm{pix^{-1}}$), considering only pixels
with no problematic spectral flags. For each pixel, a $\sigma$-clipping is
applied to reject fluxes $F_{i}$ deviating more than $4\sigma$ from the mean,
minimizing the impact of potential night sky emission-line residuals. The
composite spectrum is obtained by computing the median flux and the associated
uncertainty, given by the normalized median absolute deviation
(NMAD222$\mathrm{NMAD}=1.4826\cdot\mathrm{median}(\left|F_{i}-\mathrm{median}(F_{i})\right|)$,
see Hoaglin et al. (1983)) divided by the square root of the number of objects
at each pixel. This technique provides a composite spectrum robust against
imperfections and not biased toward galaxies with higher S/N. The resulting
S/N per pixel measured around $4000$ Å rest frame is $\sim 230$.

In Figure 2 we show the median composite spectrum of the 350 passive galaxies
selected in this work. The spectrum shows features typical of a passively
evolving population. Qualitatively, we note a red continuum, indicative of an
old stellar population; a significant 4000 Å discontinuity; the absence of a
Balmer break (3650 Å); a $\mathrm{Ca\>\textsc{ii}\>K}$ absorption-line deeper
than the associated $\mathrm{Ca\>\textsc{ii}\>H}$ (Section 2.3.1); and the
presence of several metallic indices (e.g., $\mathrm{G4300}$ ,
$\mathrm{Fe4383}$ , $\mathrm{Fe4531}$ , and $\mathrm{Mg_{2}}$ ). It is very
important to notice that even at the very high S/N of the stacked spectrum, no
emission line is detectable, confirming the robustness of the selection. On
the contrary, the composite spectrum of galaxies excluded by visual inspection
has significant [O ii] and [O iii] emission-lines.

In conclusion, the distribution of SED-fitting-derived stellar masses and SFR,
as well as the analysis of the composite spectrum, confirms the reliability of
our selection of passive galaxies. These systems at $0.6<z<1$ show no
detectable evidence of recent star formation. The presence of a possible
underlying young component will be further assessed by studying spectral
absorption features (Section 2.4).

###  2.3 Measuring spectral indices with PyLick

In its complete version, the Lick system consists of 25 indices in the rest-
frame wavelength interval $4000-6000$ Å (Worthey & Ottaviani, 1997; Trager et
al., 1998). Each index requires the definition of a central region
($\lambda_{c1}$, $\lambda_{c2}$) and two other regions located toward the red
($\lambda_{r1}$, $\lambda_{r2}$) and blue ($\lambda_{b1}$, $\lambda_{b2}$) of
the central one to estimate a reference pseudo-continuum level. Following the
approach of the Lick group, the strenghts of atomic $I_{a}$ and molecular
$I_{m}$ indices are calculated using the following equations:

$\displaystyle
I_{a}(\text{[\AA]})=\int_{\lambda_{c_{1}}}^{\lambda_{c_{2}}}\left(1-\frac{F(\lambda)}{F_{c}(\lambda)}\right)\;\mathrm{d}{\lambda},$
(1) $\displaystyle
I_{m}([\mathrm{mag}])=-2.5\log\left(\frac{1}{\lambda_{c_{2}}-\lambda_{c_{1}}}\int_{\lambda_{c_{1}}}^{\lambda_{c_{2}}}\frac{F(\lambda)}{F_{c}(\lambda)}\;\mathrm{d}{\lambda}\right),$
(2)

where $F(\lambda)$ and $F_{c}(\lambda)$ are the spectrum flux and the local
pseudo-continuum, respectively. The latter is commonly derived through linear
interpolation:

$F_{c}(\lambda)=F_{b}\frac{\lambda_{r}-\lambda}{\lambda_{r}-\lambda_{b}}+F_{r}\frac{\lambda-\lambda_{b}}{\lambda_{r}-\lambda_{b}},$
(3)

where $\lambda_{b,r}$ and $F_{b,r}$ are the central wavelengths and mean
fluxes of the lateral blue and red regions.

Table 2: Spectral Indices Properties for the 350 Passive Galaxies

Index | $z^{(a)}$ | Range(b) | $\left\langle\right.\\!\\!{I}\\!\\!\left.\right\rangle^{(c)}$ |  $\left\langle\right.\\!\\!{S/N}\\!\\!\left.\right\rangle$(c)
---|---|---|---|---
$\mathrm{Ca\>\textsc{ii}\>K}$ | $>0.65$ | 4.928–9.157 | 7.167 | 22.40
$\mathrm{Ca\>\textsc{ii}\>H}$ | $>0.65$ | 5.472–8.340 | 6.891 | 40.93
$\mathrm{D4000}$ | $>0.65$ | 1.732–2.028 | 1.881 | 236.61
$\mathrm{D_{n}4000}$ | $>0.65$ | 1.517–1.846 | 1.690 | 172.13
$\mathrm{H\delta_{A}}$ | All | –1.438–3.826 | 0.347 | 3.44
$\mathrm{H\delta_{F}}$ | All | 0.630–3.516 | 1.709 | 10.37
$\mathrm{CN_{1}}$ | All | –0.039–0.098 | 0.035 | 5.34
$\mathrm{CN_{2}}$ | All | –0.005–0.145 | 0.074 | 9.21
$\mathrm{Ca4227}$ | All | 0.290–1.701 | 0.922 | 7.69
$\mathrm{G4300}$ | All | 2.605–5.818 | 4.443 | 18.29
$\mathrm{H\gamma_{A}}$ | All | –4.225–0.431 | –2.729 | 9.86
$\mathrm{H\gamma_{F}}$ | All | –1.121–2.348 | –0.105 | 3.64
$\mathrm{Fe4383}$ | All | 1.675–5.778 | 3.832 | 11.87
$\mathrm{Ca4455}$ | $<0.96$ | 0.219–2.366 | 1.395 | 8.54
$\mathrm{Fe4531}$ | $<0.92$ | 1.438–4.599 | 2.958 | 11.16
$\mathrm{C_{2}4668}$ | $<0.88$ | 1.423–7.925 | 5.097 | 13.94
$\mathrm{H\beta}$ | $<0.77$ | 1.175–3.490 | 1.989 | 13.72
$\mathrm{Fe5015}$ | $<0.74$ | 1.641–7.171 | 4.737 | 13.98
$\mathrm{Mg_{1}}$ | $<0.68$ | 0.036–0.145 | 0.077 | 23.94
$\mathrm{Mg_{2}}$ | $<0.68$ | 0.079–0.280 | 0.205 | 46.88
$\mathrm{Mg}\,b$ | $<0.69$ | 1.789–4.827 | 3.274 | 17.08
$\mathrm{Fe5270}$ | $<0.64$ | 0.628–4.054 | 2.549 | 15.96
$\mathrm{Fe5335}$ | $<0.63$ | 0.425–3.553 | 2.448 | 13.02
$\mathrm{Fe5406}$ | $<0.62$ | 0.144–2.462 | 1.573 | 8.67

Note. — $(a)$ Expected redshift coverage within $0.6<z<1$ of the VIMOS HR red
spectrograph; $(b)$ Computed between $5^{th}-95^{th}$ percentiles; $(c)$
Median index value and signal-to-noise ratio. Index units are angström for all
indices except: $\mathrm{CN_{1}}$ , $\mathrm{CN_{2}}$ , $\mathrm{Mg_{1}}$ ,
$\mathrm{Mg_{2}}$ (mag); $\mathrm{D4000}$ , $\mathrm{D_{n}4000}$ (dex).

One of the strongest features in passive galaxy spectra is the 4000 Å
discontinuity (D4000). At bluer wavelengths, the flux suddenly declines due to
the accumulation of a large number of spectral lines that are present in
stellar types cooler than G0. This property makes it a good age tracer (e.g.,
Kauffmann et al., 2003; Moresco et al., 2011). A discontinuity index can be
quantified as the ratio of the average flux density in two regions redwards
and bluewards of 4000 Å:

$\mathrm{D4000}=\frac{\lambda_{b2}-\lambda_{b1}}{\lambda_{r2}-\lambda_{r1}}\frac{\int_{\lambda_{r1}}^{\lambda_{r2}}\lambda^{2}F(\lambda)\;\mathrm{d}{\lambda}}{\int_{\lambda_{b1}}^{\lambda_{b2}}\lambda^{2}F(\lambda)\;\mathrm{d}{\lambda}}$
(4)

where $\lambda_{b1},\lambda_{b2}=3750,3950$ Å and
$\lambda_{r1},\lambda_{r2}=4050,4250$ Å in the original definition by Bruzual
(1983). Alternatively, Balogh et al. (1999) introduced a narrower definition
of the $\mathrm{D4000}$ ($\mathrm{D_{n}4000}$ , estimated over the ranges
3850–3950 and 4000–4100 Å), which is optimized to be less sensitive to
reddening effects.

To perform all the measurements of absorption features, we developed a
dedicated Python code named PyLick. It is a public333The code is available at
https://gitlab.com/mmoresco/pylick/ flexible Python library developed for the
fast and accurate estimation of the main spectral indices introduced in the
literature from the UV to the near-IR. A more detailed description of the main
elements of the code, as well as validation tests against available LEGA-C DR2
data, can be found in A. In this work, we will focus on atomic and molecular
(Lick) indices and $\mathrm{D4000}$ . Because spectral indices are defined on
fixed wavelength intervals, the measurement can be affected by systematics if
spectral broadening effects (e.g., instrumental resolution, stellar velocity
dispersion) are not properly taken into account. In general, this results in
lower index values compared to intrinsic ones. To address this issue, we
convolve each rest-frame spectrum ($\mathrm{FWHM}_{\mathit{spec}}\approx 1.3$
Å at $z=0.7$) to match the resolution of the models that will be implemented
($\mathrm{FWHM}_{\mathit{mod}}=2.5$ Å), by using a Gaussian kernel with a
standard deviation of
$\sigma=(\mathrm{FWHM}_{\mathit{model}}^{2}-\mathrm{FWHM}_{\mathit{spec}}^{2})^{1/2}/2.355$.
After the measurement, we calibrate indices to zero-velocity dispersion
following the approach described by Carson & Nichol (2010). The procedure to
derive the correction coefficients is described in detail in Huchet et al.
(2022, in preparation). Briefly, we measure indices on SSP spectra broadened
at different $\sigma_{\star}$ generated from the 2016 version of Bruzual &
Charlot (2003) models using the MILES stellar spectral library (Falcón-Barroso
et al., 2011) and Chabrier (2003) initial mass function444The models are
available at http://www.bruzual.org/bc03/Updated_version_2016/ and fit each
index–$\sigma_{\star}$ relation with a fourth-order polynomial.

At the end of this process, we obtain a set of 24 spectral indices, namely:

  1. -

Balmer indices: $\mathrm{H\delta_{A}}$ , $\mathrm{H\delta_{F}}$ ,
$\mathrm{H\gamma_{A}}$ , $\mathrm{H\gamma_{F}}$ , $\mathrm{H\beta}$ ;

  2. -

Iron-dominated indices: $\mathrm{Fe4383}$ , $\mathrm{Fe4531}$ ,
$\mathrm{Fe5015}$ , $\mathrm{Fe5270}$ , $\mathrm{Fe5335}$ ;

  3. -

Molecular indices: $\mathrm{CN_{1}}$ , $\mathrm{CN_{2}}$ , $\mathrm{Mg_{1}}$ ,
$\mathrm{Mg_{2}}$ ;

  4. -

Other Lick indices available in the LEGA-C spectral range: $\mathrm{Ca4227}$ ,
$\mathrm{G4300}$ , $\mathrm{Fe4531}$ , $\mathrm{C_{2}4668}$ , $\mathrm{Mg}\,b$
;

  5. -

4000 Å discontinuity indices: $\mathrm{D4000}$ , $\mathrm{D_{n}4000}$ ;

  6. -

Two recently defined pseudo-Lick indices: $\mathrm{Ca\>\textsc{ii}\>K}$ ,
$\mathrm{Ca\>\textsc{ii}\>H}$ (Section 2.3.1).

Table 3: Definitions of Pseudo-Lick $\mathrm{Ca\>\textsc{ii}\>K}$ and $\mathrm{Ca\>\textsc{ii}\>H}$ Indices Introduced in Fanfani (2019) Index | Central | Blue | Red
---|---|---|---
| Å | Å | Å
$\mathrm{Ca\>\textsc{ii}\>K}$ | 3925.65 – 3945 | 3845 – 3880 | 3950 – 3954
$\mathrm{Ca\>\textsc{ii}\>H}$ | 3959.40 – 3978 | 3950 – 3954 | 3983 – 3993

This dataset extends LEGA-C DR2 public catalog from $\mathrm{H\beta}$ to
$\mathrm{Fe5406}$ indices and, in particular, Mg ones, commonly used as
proxies to study the $\alpha$ enhancement. For each index, we compute redshift
coverage, 5th–95th percentile range, median value, and median S/N. These
values are presented in Table 2.

Trivially, bluer (redder) indices are available only at higher (lower)
redshift. This information, combined with the upper panel of Figure 2, gives
an idea of the spectral and redshift coverage of this study. Most of the
indices are available only for a narrow rest-frame wavelength range between
3900–4500 Å. However, a significant number of galaxies ($\sim 200$) still
share a common wavelength range between 3700–4900 Å. In this case, the
redshift coverage is reduced to $z\lesssim 0.9$. The median S/N of the
measured indices is $>10$, except for those with a signal $\sim 0$ (Balmer and
CN indices), or those defined on a narrow central region ($\mathrm{Ca4227}$ ,
$\mathrm{Ca4455}$ )

Figure 3: Measurement of pseudo-Lick indices $\mathrm{Ca\>\textsc{ii}\>K}$ and
$\mathrm{Ca\>\textsc{ii}\>H}$ of a bona fide passive galaxy (violet) and a
NUVrJ star-forming galaxy (blue, with a vertical offset of $+0.35$). Vertical
bands show the regions where these indices are defined (red pseudo-continuum
of $\mathrm{Ca\>\textsc{ii}\>H}$ and blue pseudo-continuum of
$\mathrm{Ca\>\textsc{ii}\>K}$ are overlapped). Vertical dotted lines are the
centroids of $\mathrm{Ca\>\textsc{ii}\>K}$ , $\mathrm{Ca\>\textsc{ii}\>H}$ ,
and H$\epsilon$ lines. Dashed segments show the estimated pseudo-continuum
flux used as a reference to compute the indices values (filled regions).

####  2.3.1 The H/K ratio

Rose (1984) firstly proposed to study stellar populations using the ratio of
the minimum fluxes of two nearby lines. Specifically,
$\mathrm{Ca\>\textsc{ii}\>H}$ over $\mathrm{Ca\>\textsc{ii}\>K}$ ratio was
found to correlate with the starburst ages in post-starburst galaxies. As a
matter of fact, H$\epsilon$ absorption line ($\lambda=3970.1$ Å), which is
deeper in the presence of young A- and B-type stars, overlaps to
$\mathrm{Ca\>\textsc{ii}\>H}$ ($\lambda=3968.5$ Å), while
$\mathrm{Ca\>\textsc{ii}\>K}$ ($\lambda=3933.7$ Å) remains relatively
uncontaminated. The nomenclature H/K is thus a compact notation for
$(\mathrm{Ca\>\textsc{ii}\>H}+{\rm H}\epsilon)/\mathrm{Ca\>\textsc{ii}\>K}$
and can be used as a stellar population diagnostic. Even the presence of a
fraction of a young ($\lesssim 1$ Gyr) stellar component can significantly
alter this spectral feature (Longhetti et al., 1999; Lonoce et al., 2014;
Moresco et al., 2018). In the literature, this value has been usually measured
as the ratio of minimum fluxes in H and K lines (e.g., Rose, 1985; Leonardi &
Rose, 1996; Longhetti et al., 1999; Lonoce et al., 2014; Moresco et al.,
2018):

$\left|\mathrm{H/K}\right|_{\mathrm{min}}=\frac{F_{\mathrm{min}}({\rm
H})}{F_{\mathrm{min}}({\rm K})}.$ (5)

This method does not rely on measuring integrated quantities (differently from
Lick indices), hence it is relatively independent of changes in spectral
broadening, like those due to the instrumental resolution and stellar velocity
dispersion. However, measurements of minima can be strongly biased by the
presence of noise peaks, especially in low-S/N regimes. For this reason, we
adopt here a hybrid approach recently introduced by Fanfani (2019), where H/K
is computed as a ratio of two pseudo-Lick indices,

$\mathrm{H/K}=\frac{I_{\rm H}}{I_{\rm K}},$ (6)

using the $\mathrm{Ca\>\textsc{ii}\>K}$ and $\mathrm{Ca\>\textsc{ii}\>H}$
index passbands listed in Table 3 that have also been included in PyLick. We
measure and $\sigma_{\star}$-correct these indices to derive H/K values and
associated uncertainties for the entire parent sample. Two representative
examples are shown in Figure 3. For a typical passive population, the H line
is less deep than the K line, therefore
$\left|\mathrm{H/K}\right|_{\mathrm{min}}>1$. The so-called H/K inversion can
already arise for a contribution of $\sim 5$ % in mass of a young stellar
population (with ages $<200$ Myr; Moresco et al. 2018). The two quantities
defined in Eqs. 5 and 6 have an inverse relationship, but there is not a
strict one-to-one correspondence. From the parent sample, we find an
equivalent inversion value for the $\mathrm{H/K}$ between 1.2 and 1.5;
therefore, the passive regime corresponding to
$\left|\mathrm{H/K}\right|_{\mathrm{min}}>1$ can be safely defined at
$\mathrm{H/K}<1.2$.

Figure 4: Distribution of LEGA-C galaxies as a function of the main spectral
indices measured in this work. Different colors are used for the parent
(gray), photometric (yellow), spectrophotometric (pink) and bona fide passive
(violet) subsamples. Arrows represent the median values.

###  2.4 Observed spectral features

The distributions of the main measured absorption features for the various
subsamples of passive galaxies are presented in Figure 4, and the
corresponding median values are quoted in Table 1.

Overall, the bona fide passive galaxies tend to segregate with respect to the
parent sample median. We observe typical $\rm{H/K}$ ratios of $0.96\pm 0.08$
and very few systems with $\rm{H/K}>1.1$, while parent sample galaxies reach
values well above $1.2$. It is important to stress here that even without
applying any cut on the $\rm{H/K}$ ratio, our selection criteria allows us to
discard the majority of the tail at high H/K, leaving us a sample with values
characteristic of a passive population. This is important and independent
evidence supporting the purity of our selection and the negligible
contamination from underlying young stellar populations in the passive sample.
Another interesting result that we report for the first time is that the
$\mathrm{Ca\>\textsc{ii}\>K}$ line itself shows a clear bimodality, with
passive galaxies having $\mathrm{Ca\>\textsc{ii}\>K}>5$ Å. This is not the
case for $\mathrm{Ca\>\textsc{ii}\>H}$ , because as mentioned above, the
H$\epsilon$ line strengthens the index in younger populations. We note here
that $\mathrm{Ca\>\textsc{ii}\>K}$ is likely to be one of the main drivers of
the D4000 bimodality as it is included in the blue passband of the D4000
index, while the $\mathrm{Ca\>\textsc{ii}\>H}$ shows a much flatter trend as a
function of redshift. The well-known $\mathrm{D_{n}4000}$ and H$\delta$
bimodalities are already studied in detail in the local universe (Kauffmann et
al., 2003; Siudek et al., 2017) and in LEGA-C data (Wu et al., 2018).
Photometrically selected massive and passive galaxies already populate
high-$\mathrm{D_{n}4000}$ and low-$\mathrm{H\delta_{A}}$ tails, indicating
relatively old stellar populations. The addition of a spectroscopic criterion
removes a significant number of remaining low-$\mathrm{D_{n}4000}$ and
high-$\mathrm{H\delta_{A}}$ galaxies. Among bona fide passive galaxies, only
11 ($3\%$) have $\mathrm{D_{n}4000}<1.5$, and only 28 ($8\%$) have
$\mathrm{H\delta_{A}}>2.5$ Å. However, we underscore here that no additional
cuts have been applied to minimize selection biases. Interestingly, we observe
a slight bimodality also for $\mathrm{G4300}$ , which measures the optical CH
band (also known as G-band) and is very sensitive to the carbon abundance
(Tripicco & Bell, 1995; Korn et al., 2005). The bona fide passive galaxies
have a relatively high $\mathrm{G4300}\gtrsim 2.5$ Å. Finally, they are
characterized by high $\mathrm{Fe4383}$ , a primary indicator of stellar
metallicity, with typical values of
$\left\langle\right.\\!\\!{\mathrm{Fe4383}}\\!\\!\left.\right\rangle\sim 3.8$
Å.

Figure 5: Binned relationships between the main spectral indices and redshift
for the bona fide passive sample divided into two velocity dispersion regimes
($\sigma_{\star}=215$ $\rm km\,s^{-1}$ as threshold). Each bin contains on
average $N\sim 30$ galaxies. Error bars on the $x$-axis represent the bin
width, while those on the $y$-axis are errors associated to the mean values.

We further explore the bona fide passive sample by looking at
index–$z$–$\sigma_{\star}$ trends. For this purpose, we firstly divide
galaxies into two stellar velocity dispersion bins using $\sigma_{\star}=215$
$\rm km\,s^{-1}$ as the threshold value. Then, we divide each subsample into
four or three redshift bins depending on the redshift coverage, considering
intervals $\Delta z\sim 0.08-0.1$. All the bins have on average $N\sim 30$
objects and in all cases $N>10$. Finally, we compute the mean of the redshift
and corresponding index in each bin, along with its associated errors. The
resulting trends for the main spectral features are presented in Figure 5. We
also quantify these trends within individual galaxies by using the
nonparametric Spearman rank correlation test.

In general, we observe that Balmer indices correlate moderately with redshift
with a Spearman $\rho\sim 0.4$ ($p$-value $10^{-8}$). At fixed $z$,
high-$\sigma_{\star}$ galaxies have weaker Balmer indices. Similar behavior is
observed for the $\mathrm{Ca\>\textsc{ii}\>K}$ index and - with an opposite
trend - for the H/K ratio. $\mathrm{D_{n}4000}$ anticorrelates moderately with
redshift with a Spearman $\rho\sim-0.3$ ($p$-value $10^{-5}$). At fixed $z$,
there is a clear separation between the two $\sigma_{\star}$ regimes. Note
that these relations could be (and have been) used in the cosmic chronometers
framework to constrain the expansion history of the Universe, once the
D4000–age relation is carefully calibrated (Moresco et al., 2011). Temporarily
ignoring metallicity effects that will be discussed in the following, these
observed trends are fully consistent with the mass-downsizing scenario, i.e.
more massive galaxies formed earlier and in short timescales, and then
experienced a passive evolution. Further support to this idea is given by iron
indices, which are stronger for high-$\sigma_{\star}$ galaxies but do not show
significant evolution over $\sim 2$ Gyr of cosmic time. This finding supports
the scenario in which this population has already exhausted its gas reservoir,
hence not being able to significantly evolve their metallic content, and is
evolving passively as a function of cosmic time.

The absence of a correlation in redshift for $\mathrm{C_{2}4668}$ and
$\mathrm{G4300}$ , instead, shows that these indices cannot be used as age
indicators throughout cosmic time. However, the segregation in
$\sigma_{\star}$ is still consistent with the idea that they are sensitive to
stellar population age in the very first Gyr after the formation. Another
explanation is that the segregation is due to different stellar metallicities
and abundances. Even in that case, they could be crucial to determine galaxy
ages indirectly by breaking the age–metallicity degeneracy.

It is also interesting to notice that $\mathrm{Fe4383}$ shows a very flat
behavior as a function of redshift, for both $\sigma_{\star}$ bins. This
trend, confirmed also by the analysis of the stellar metallicity presented
below, supports the scenario for which these systems formed the majority of
their mass (and hence metallicity) at high redshifts and in very concentrated
episodes of star formation, exhausting their gas reservoirs and therefore not
being able to further change their metal content.

Qualitatively similar trends are obtained also maximizing the number of
galaxies per bin instead of using bins of constant $\Delta z$. Given the
inhomogeneous redshift coverage of LEGA-C, this latter method helps to
strengthen bin statistics but has the cost of losing leverage in redshift.

###  2.5 Morphology

We analyze the morphological content of the bona fide passive sample using the
Zurich Estimator of Structural Types (ZEST) classification (Scarlata et al.,
2007), based on principal component analysis of the surface brightness
profiles. We find that the great majority (71%) of galaxies are classified as
E/S0, 27% as intermediate (where the contribution of the bulge is similar to
that of the disk), 2% (six galaxies) as irregular, and none show a late-type
morphology. Among the six galaxies classified as irregular, only two appear in
the final sample of constrained galaxies. Their morphology is not
significantly disturbed and their spectra are typical of old stellar
populations, therefore we still include them in the bona fide passive sample
and verify that their exclusion does not affect our results. Within the bona
fide passive sample, we observe no significant correlation between different
morphological classes and $\sigma_{\star}$. Interestingly, a similar
percentage of E/S0 types (72%) has also been found by Moresco et al. (2013) in
their sample of $\sim 17000$ zCOSMOS galaxies. The two works share 127
galaxies (37% of the bona fide passive sample) and adopt the same
morphological catalog. However, selection criteria are different because the
authors make use of the best-fit SED matching a local E-S0 template and
different color cuts, along with spectroscopic cuts based on both [O ii] and
H$\alpha$ emission.

The presence of passive systems with non-purely early-type morphologies has
been already discussed in the literature (e.g., Dressler et al., 1999;
Pozzetti et al., 2010). This result can be explained by considering the
existence of a class of objects for which the morphological transformation
lasts longer than changes in stellar population, i.e. galaxy colors redden
before the galaxy reaches an early-type morphology. It is therefore improper
to treat galaxies without an early-type morphology as contaminants.

In conclusion, physical properties derived from SED fitting (Figure 1) and
observed spectral properties of the three different subsamples considered
(Figures 2, 4, and 5), confirm the existence of a population of passive
systems characterized by lower/absent star formation, older stellar
populations, and higher metallicities with respect to the parent galaxy
population, and consistent with the mass-downsizing scenario. These
observational data also support that the bona fide selection adopted is able
to maximize the purity of the sample, providing a sample of massive and
passive galaxies, with a negligible (if any) contamination by star-forming
outliers. The analysis of each individual galaxy will add more granularity to
this picture.

##  3 Analysis

With this work, we are exploring the capability of studying cosmic
chronometers with Lick indices, which allow us to compare with several
literature studies (e.g., Thomas et al., 2010; Gallazzi et al., 2014; Onodera
et al., 2015). The stellar ages derived with this method are light weighted
and can be biased toward lower values because young stars outshine the older
ones (Conroy, 2013), while metallicity and abundance are found to be less
affected (Serra & Trager, 2007). When a young stellar population is present,
mass weighted age - which can be derived with full spectral fitting codes -
can provide a better estimate of the integrated galaxy SFH. From the analysis
of the H/K ratio (Sect. 2.3.1) we can exclude a contamination from young (200
Myr–1 Gyr) stars in our bona fide passive galaxies, therefore also the
difference between light- and mass-weighted is expected to be negligible.

To more quantitatively assess the robustness for the CC method, in Kang et al.
(2022, in preparation) we perform a full spectral fitting analysis using the
publicly available BAGPIPES code (Carnall et al., 2018) to determine the ages,
metallicities, and SFHs of the bona fide galaxies selected in this paper.
While we refer to that paper for the detailed description of the methods
adopted and results obtained, we anticipate here that the slope of the
$\mathrm{age}-z$ relation is very robust and in perfect agreement between our
two studies, with a percentage difference $\lesssim 5\%$ ($\lesssim
0.2\sigma$).

###  3.1 Stellar population model

In this work, we adopt the Thomas, Maraston, & Johansson (2011) models
(hereafter TMJ11), which make predictions of the Lick indices by varying SSP-
like stellar ages ($\mathrm{age}$), metallicities ($\rm[Z/H]$), and
$\alpha$-element enhancements ($\rm[\alpha/Fe]$) of the stellar population.
These models have been widely used in the literature for similar Lick indices
studies (e.g., Johansson et al., 2012; Jørgensen & Chiboucas, 2013; Onodera et
al., 2015; Scott et al., 2017; Lonoce et al., 2020). They are based on the
evolutionary synthesis code of Maraston (2005) and element response functions
from Korn et al. (2005), and are carefully calibrated with galactic globular
cluster data. The main ingredients are single-burst SFH, Salpeter (1955)
initial mass function, MILES empirical stellar libraries (Sánchez-Blázquez et
al., 2006), and Cassisi et al. (1997) stellar evolutionary tracks. In this
work, we use the models provided at a MILES resolution of 2.5 Å (Beifiori et
al., 2011).

The original grid spans the following parameter space:
$0.1<\mathrm{age}/\mathrm{Gyr}<15$, $-2.25<\hbox{$\rm[Z/H]$}<0.67$, and
$-0.3<\hbox{$\rm[\alpha/Fe]$}<0.5$ with a total of 480 grid points, each one
corresponding to the prediction for a single SSP. We perform a three-
dimensional linear interpolation of the grid to reach a resolution of
$\Delta\mathrm{age}=0.1$ Gyr, and
$\Delta\hbox{$\rm[Z/H]$}=\Delta\hbox{$\rm[\alpha/Fe]$}=0.01$ dex. This process
allows us to achieve higher precision in parameter estimation. We verified
that our approach does not introduce any systematic bias, since by considering
different grid choices the final results are always fully compatible with each
other within $1\sigma$.

These models assume an instantaneous burst of star formation. While this
assumption might in principle lead to a significant underestimation of the
global stellar age in mixed populations, we note here that our selection
criteria were chosen to obtain a sample with a very concentrated SFH and to
minimize the contamination from a significant residual star formation, as
confirmed by the analysis of various indicators discussed in Section 2. In
particular, for this population we expect the SFHs to be extremely coeval,
with very small durations ($\tau\lesssim 0.3$ Gyr, if modeled with a delayed
exponential SFH), as confirmed from a parallel analysis (Borghi et al., 2022).

###  3.2 MCMC Analysis

To compare the measured absorption features to TMJ11 models, we develop a
fully Bayesian analysis pipeline. We assume that the uncertainties on indices
are well-determined, Gaussianly distributed, and independent. A set of modeled
indices, which are a function of parameters
$\boldsymbol{\theta}=(\mathrm{age},\hbox{$\rm[Z/H]$},\hbox{$\rm[\alpha/Fe]$})$,
can therefore be fitted to the observed ones using the log-likelihood function

$\ln\left(\mathcal{L}\right)=k-\frac{1}{2}\sum_{i}^{\rm
N_{ind}}\left(\frac{I_{i}-I^{\mathrm{mod}}_{i}(\boldsymbol{\theta})}{\sigma_{i}}\right)^{2},$
(7)

where $k$ is a constant, $I^{\mathrm{mod}}_{i}(\boldsymbol{\theta})$ the model
prediction for the $i$th observed index $I_{i}$, and $\sigma_{i}$ its
uncertainty. The posterior probability distributions of $\boldsymbol{\theta}$
are explored using the affine-invariant ensemble sampler emcee (Foreman-Mackey
et al., 2019). Chains are initialized with 200 walkers randomly scattered
around the center of the parameter space. Each walker performs at least 2000
model realizations. Parameters and uncertainties are defined as the median,
and 16th and 84th percentiles of the marginalized posterior distributions.

We use flat priors that span the entire parameter space allowed from the
models. An important point that we stress here is that we do not assume any
cosmological prior for galaxy ages. This is a crucial point, and a difference
with respect to other similar works, to avoid introducing cosmological biases
in the age determination and keep the results cosmological-independent.

After the analysis, we carefully assess the convergence of each chain. A
detailed description is given in B.

###  3.3 Index combination

The choice of the index set to analyze must be carefully addressed because,
given the number of the measured indices, there are more than one million
possible combinations. The relative sensitivity of different indices to
different stellar population parameters and abundances is not identical
(Tripicco & Bell, 1995; Korn et al., 2005; Lee et al., 2009). Balmer and
$\mathrm{D_{n}4000}$ indices are better suited to constrain ages, Fe-dominated
indices to measure the Fe abundance and total stellar metallicity, while Mg
indices to estimate the $\alpha$ elements abundance. It is therefore intuitive
that different index sets provide different constraints.

Because the ultimate goal of this work is to avoid any effect that could
potentially introduce biases in the $\mathrm{age}-z$ relation, we choose to
fit galaxies considering always the same set of indices. This means, on one
hand, that using indices covering a wide redshift range will significantly
reduce the number of galaxies analyzed; on the other hand, to maximize the
number of analyzable galaxies the set of indices has to be chosen to span a
small wavelength range (from 3600 to 4900 Å; see the upper panel of Figure 2
and Table 2). The tradeoff between these two aspects combined has the cost of
reducing by about one-third the number of galaxies measured but guarantees the
crucial advantage of providing a homogeneous analysis.

Table 4: Extract from the catalog of measured stellar ages, metallicities
$\rm[Z/H]$, and $\rm[\alpha/Fe]$ and associated $16^{\mathit{th}}-84^{th}$
percentile uncertainties

ID |  |  | $\mathrm{age}$ | $\rm[Z/H]$ | $\rm[\alpha/Fe]$
---|---|---|---|---|---
[MMS2013] |  |  | Gyr | dex | dex
133240 |  |  | $3.08_{-0.64}^{+0.64}$ | $+0.26_{-0.09}^{+0.12}$ | $+0.16_{-0.06}^{+0.06}$
133783 |  |  | $1.90_{-0.08}^{+0.07}$ | $+0.33_{-0.03}^{+0.05}$ | $+0.14_{-0.03}^{+0.03}$
134169 |  |  | $2.79_{-0.36}^{+0.36}$ | $+0.16_{-0.07}^{+0.08}$ | $+0.24_{-0.05}^{+0.05}$
139772 |  |  | $3.05_{-0.12}^{+0.14}$ | $+0.12_{-0.02}^{+0.02}$ | $+0.23_{-0.03}^{+0.03}$
205742 |  |  | $3.07_{-0.27}^{+0.39}$ | $-0.01_{-0.04}^{+0.06}$ | $+0.04_{-0.03}^{+0.03}$
206573 |  |  | $3.30_{-0.19}^{+0.29}$ | $-0.02_{-0.03}^{+0.02}$ | $+0.28_{-0.03}^{+0.03}$
207825 |  |  | $1.94_{-0.07}^{+0.12}$ | $+0.16_{-0.05}^{+0.01}$ | $+0.16_{-0.02}^{+0.02}$
…

Note. — Full table is available online.

For these reasons, we decide to use the following set of Lick indices:
$\mathrm{H\delta_{A}}$ , $\mathrm{CN_{1}}$ , $\mathrm{CN_{2}}$ ,
$\mathrm{Ca4227}$ , $\mathrm{G4300}$ , $\mathrm{H\gamma_{A}}$ ,
$\mathrm{H\gamma_{F}}$ , $\mathrm{Fe4383}$ , $\mathrm{Fe4531}$ , and
$\mathrm{C_{2}4668}$ . They are chosen among those that are calibrated against
globular cluster data in TMJ11, but excluding those redder than
$\mathrm{H\beta}$ , because, as mentioned above, they would not allow us to
obtain a statistically meaningful sample. We also exclude $\mathrm{H\beta}$
because it can be biased by the presence of a residual emission component
(e.g., Concas et al., 2017). The final set spans a narrow wavelength range in
the optical regime, from 4000 to 4800 Å. We verified on a few available
galaxies that the inclusion of $\mathrm{Mg}\,b$ does not significantly change
the results. We also performed extensive tests with very different sets of
indices. A detailed discussion is presented in C. Given the available data, we
find that this is the optimal set, because it maximizes both the spectral
coverage and the number of galaxies for which we obtain constraints, and it
also provides a good balance between age-, metallicity-, and
$\alpha$-sensitive indices.

In the end, we obtain robust constraints for 140 galaxies over 199 that have
been analyzed (Table 4), after removing the galaxies with nonconverging fits
(see B). These galaxies are at
$\left\langle\right.\\!\\!{z}\\!\\!\left.\right\rangle=0.70$, and have
spectral $\left\langle\right.\\!\\!{\rm S/N}\\!\\!\left.\right\rangle=26.4$
per pixel. Typical uncertainties are $\pm 0.33$ Gyr in $\mathrm{age}$ and $\pm
0.05$ dex in $\rm[Z/H]$ and $\rm[\alpha/Fe]$.

##  4 Results and discussion

In this section, we present the results obtained from the analysis of the
previously defined H/K ratio and for the physical properties of selected
passive galaxies in the LEGA-C DR2 survey. We start by discussing the H/K
ratio and the correlation with commonly used diagnostics in Section 4.1. In
Section 4.2 we explore scaling relations of stellar population parameters
versus $\sigma_{\star}$ and $M_{\star}$. In Section 4.3 we discuss trends with
redshift, with a focus on the age–redshift relation. Finally, in Section 4.4
we present the median binned parameters–redshift relations.

We underline here below the main points of our analysis to be kept in mind
when comparing our results with those of previous analyses in the same field.

  1. (i)

Our galaxies are passive. Similar studies are mostly targeting morphological
early-type galaxies. These samples may contain galaxies with a nonnegligible
level of star-formation ($\sim 20\%$ at the current stellar masses; Moresco et
al. 2013);

  2. (ii)

We do not use cosmological priors. We also verified that, especially for
low-S/N galaxies, the use of a cosmological upper limit for the galaxy ages
may produce an apparent convergence of Markov Chain Monte Carlo (MCMC) chains
at higher ages toward the prior, introducing a bias in the sample (B);

  3. (iii)

We adopt single-burst SFHs. In this case, the time of the onset of star-
formation $t_{\mathrm{form}}$ coincides with the time of the main (and only)
star-formation event $t_{\rm peak}$, and with the time of quenching $t_{\rm
quench}$. A proper comparison with other data sets where an extended SFH is
assumed should carefully take into account how ages are defined and intrinsic
degeneracies between SFH parameters. For instance, by assuming an
$\mathrm{SFR}(t)\propto\exp((t-\mathrm{age})/\tau)$ (tau model), a positive
correlation between the timescale $\tau$ and the galaxy age (time since
$t_{\mathrm{form}}$) can be present, with values of $\Delta\tau\sim
0.3~{}\Delta\mathrm{age}$ (Borghi et al., 2022).

Before applying the cosmic chronometer method, points (i) and (ii) should be
carefully considered, but it is important to stress that this method relies on
differential - not absolute - galaxy ages. Therefore, by assuming a more
extended SFH for the entire population of these galaxies (i.e. a vertical
offset in the $\mathrm{age}-z$ plane), the final $H(z)$ measurement is not
affected (see Section 3.1).

Figure 6: Four well-known diagnostic diagrams (NUVrJ, UVJ, EW[O
ii]–$\mathrm{D_{n}4000}$ , SFR–$M_{\star}$) color-coded by
$\mathrm{H/K\;ratio}$ using the LOESS method (see text). Each panel shows the
parent sample galaxies extracted from LEGA-C DR2 with a S/N of the
$\mathrm{H/K\;ratio}$ higher than 3. Black borders identify bona fide passive
galaxies. Dashed lines are the criteria to separate star-forming from passive
galaxies taken from the literature (Ilbert et al., 2013; Williams et al.,
2009; Mignoli et al., 2009; Spilker et al., 2018, respectively). We use black
colors for the criteria also adopted in this work, while gray lines are for
illustrative purposes only.

###  4.1 The H/K ratio as a stellar population diagnostic

In Figure 6 we show the distribution of the entire parent sample in four
widely used diagnostic diagrams, namely NUVrJ, UVJ, EW[O
ii]–$\mathrm{D_{n}4000}$ , and SFR–$M_{\star}$, color-coded by the H/K ratio.
To better capture the mean trends of H/K, we perform locally weighted
regression (LOESS). We use the LOESS package of Cappellari et al. (2013a)
based on the two-dimensional algorithm of Cleveland & Devlin (1988), with a
linear local approximation and a regularization factor $f=0.5$. Selected bona
fide passive galaxies are highlighted with black borders. To quantify these
trends, we use Spearman correlation coefficients.

The diagrams show a clear and significant correlation between the H/K ratio
and other diagnostic tools despite its small dynamic range ($0.8\lesssim{\rm
H/K}\lesssim 2$, and a median error of $\sim 0.14$). In more detail, we
observe a strong correlation with (NUV–$r$) ($\rho=-0.72$;
$\text{$p$-value}=10^{-164}$) and with ($U-V$) ($\rho=-0.72$;
$\text{$p$-value}=10^{-178}$). The correlation with ($r-J$) and ($V-J$) colors
is weaker but still significant ($\rho=-0.25$; $\text{$p$-value}=10^{-19}$ and
$10^{-15}$, respectively). Interestingly, we find that a selection based on
the threshold value $\rm{H/K}<1.1$ can reproduce the NUVrJ selection with
$17\%$ incompleteness and $19\%$ contamination.

As shown in the third panel, an H/K cut does not exclude galaxies with a
significant [O ii] emission. For instance, $19\%$ of the plotted galaxies have
${\rm H/K}<1$ and $\text{EW[O\>{ii}]}>5$ Å. Unfortunately, the lack of other
spectral features as H$\alpha$, [N ii], and [S ii], due to the limited
wavelength coverage of the current data set, does not allow us to investigate
the nature of these sources. For these systems, the combination of multiple
indicators is still needed to obtain a pure sample of passive galaxies.

A strong correlation is also observed for the sSFR ($\rho=0.70$;
$\text{$p$-value}=10^{-161}$). We find that a threshold on H/K values of ${\rm
H/K}<1.1$ can reproduce a ${\rm\log sSFR/yr}<-11$ cut with $15\%$
incompleteness and $16\%$ contamination.

These are remarkable results if we compare the data needed for the two
different selections, and the range of wavelengths spanned. On one hand, a
wide photometric coverage is needed for a reliable estimate of a NUVrJ
diagram, SFR, or $M_{\star}$ (typically from the UV to the near-IR), with an
accuracy increasing with the number of available photometric points; it is
therefore not always available in many surveys. On the other hand, we have a
feature defined over a window of only about $150$ Å for which deep rest-frame
optical spectroscopy is needed. The H/K can therefore play a key role in the
selection of pure samples of passive galaxies in future wide-field grism
surveys such as Euclid (Laureijs et al., 2011) and the Roman Space Telescope
(Spergel et al., 2015). Another advantage of this diagnostic is the mild
dependence on spectral resolution. Differences between H/K values obtained on
8 and $2.5$ Å FWHM spectra are $\lesssim 4\%$. Performing the same analysis on
individual $\mathrm{Ca\>\textsc{ii}\>K}$ and $\mathrm{Ca\>\textsc{ii}\>H}$
indices yields $\sim 10\%$ systematic differences.

Figure 7: Distribution of mean stellar ages, metallicities, and $\alpha$
enhancements as a function of stellar velocity dispersion for individual
massive and passive galaxies selected in LEGA-C DR2 (violet points). Gray
contours represent early-type galaxies in the local Universe from SDSS/MOSES
(Thomas et al., 2010), and enclose $1\sigma$, $2\sigma$, and $5\sigma$
regions. Dashed lines and shaded regions are robust linear fits and associated
$2\sigma$ scatter regions, for $2.1<\log\sigma_{\star}<2.5$. The resulting
slopes and their uncertainties are annotated in the bottom left. Black arrows
represent the age evolution expected for a passive SSP.

Concerning the 140 bona fide passive galaxies, we note here that no
correlation between stellar population parameters (especially age) and H/K is
present. This result is expected given the fact that no significant
contribution from a young stellar component is present, i.e. the H$\epsilon$
line deepening effect becomes negligible in the passive regime.

In summary, young ($\sim 200$ Myr) stellar populations whose light is
predominantly due to A- and B-type stars (${\rm H/K}>1.1$) are characterized
by higher UV fluxes, lower D4000, a higher EW[O ii], and dominate the star
formation main sequence of LEGA-C galaxies. The H/K lowers when approaching
the quiescence criteria, but the study of timescales and the interplay with
stellar population parameters for the whole population of galaxies will
require a further assessment.

###  4.2 Physical parameters versus $\sigma_{\star}$ and $M_{\star}$

Figure 7 shows stellar population parameters as a function of the observed
stellar velocity dispersion for passive galaxies at $z\sim 0.7$. As a local
control sample, we consider the SDSS/MOSES data set of Thomas et al. (2010) of
morphologically selected early-type galaxies (ETG) at $z\sim 0.05$, since
similar models and analysis techniques are adopted. Their sample is $\sim 20$
times larger than ours and spans a wider range in $\log\sigma_{\star}$
(1.7–2.5). For a fair comparison, the local sample is limited to comparable
stellar velocity dispersion values, i.e. $2.1<\log\sigma_{\star}~{}[\rm
km\,s^{-1}]<2.5$. To this catalog we add galaxy stellar-mass estimates from
MPA-JHU DR 8 (Kauffmann et al., 2003). To study $\log\mathrm{age}$,
$\rm[Z/H]$, and $\rm[\alpha/Fe]$ versus $\log\sigma_{\star}$ and $\log
M_{\ast}$ we perform robust linear regression with the least trimmed squares
(LTS) algorithm (Rousseeuw, 1984) and measure their Spearman coefficients.
Results are quoted in Table 5.

Table 5: Coefficients for the Scaling Relations in the bona fide Passive
Sample and Associated Spearman Coefficients and $p$-values

y | $a\pm{\rm err}(a)$ | $b\pm{\rm err}(b)$ | rms |  $\rho$ ($p$-value)
---|---|---|---|---
$y=a\;\log\sigma_{\star}+b$
$\log\mathrm{age}/\mathrm{Gyr}$ | $0.5\pm 0.1$ | $-0.6\pm 0.3$ | $0.1$ |  $0.2$ $(0.02)$
$\rm[Z/H]$ | $0.3\pm 0.2$ | $-0.5\pm 0.4$ | $0.2$ |  $0.2$ $(0.15)$
$\rm[\alpha/Fe]$ | $0.2\pm 0.1$ | $-0.4\pm 0.3$ | $0.1$ |  $0.2$ $(0.03)$
$y=a\;\log(M_{\star}/10^{11}M_{\odot})+b$
$\log\mathrm{age}/\mathrm{Gyr}$ | $0.19\pm 0.04$ | $0.47\pm 0.01$ | $0.1$ |  $0.3$ $(<0.01)$
$\rm[Z/H]$ | $0.03\pm 0.05$ | $0.07\pm 0.02$ | $0.2$ |  $0.0$ $(0.63)$
$\rm[\alpha/Fe]$ | $0.02\pm 0.03$ | $0.14\pm 0.01$ | $0.1$ |  $0.1$ $(0.50)$

Note. — Linear fits are obtained with the LtsFit routine (Cappellari et al.,
2013b).

Figure 8: Formation redshifts as a function of galaxy stellar mass for 140
bona fide passive galaxies selected in LEGA-C DR2 (violet points). We compare
our measurements with literature data of massive quiescent galaxies in the
local (Thomas et al., 2010) and intermediate-redshift (Jørgensen & Chiboucas,
2013; Choi et al., 2014; Gallazzi et al., 2014; Onodera et al., 2015)
Universe. Lines are obtained by performing robust linear fits (see text).

The general trends in $\sigma_{\star}$ are consistent with a stellar
population that experienced a passive evolution from $z\sim 0.7$ to $z\sim
0.05$. In our $z\sim 0.7$ sample we find shallower relations between the
stellar population parameters and $\sigma_{\star}$. This is expected from the
lower statistics and the different selection criteria adopted. Moreover, the
dynamic range available in $\mathrm{age}$ decreases with increasing redshift
because the Universe gets younger.

We first consider SSP-equivalent stellar ages. Remarkably, the difference of
5.5 Gyr between the two samples is perfectly in agreement with the age
evolution of the Universe, confirming that the stellar populations of this
population of galaxies experienced a pure passive evolution within the assumed
reference cosmology. Contamination from young stellar populations in the
high-$z$ (low-$z$) sample would produce a larger (smaller) offset. We find
mild correlations between $\log\mathrm{age}$ vs. $\log\sigma_{\star}$ (with a
slope of $0.5\pm 0.1$) and $\log\mathrm{age}$ vs. $\log M_{\star}$ (with a
slope of $0.19\pm 0.04$). To facilitate the comparison between different
works, we convert galaxy ages to formation redshifts $z_{\mathrm{form}}$, and
study the relation with $M_{\star}$, obtaining

$z_{\mathrm{form}}=\left(0.40\pm
0.05\right)\>\log_{10}\left(\frac{M_{\star}}{10^{11}\,M_{\odot}}\right)+\left(1.46\pm
0.02\right),$ (8)

with an intrinsic scatter of $\sim 0.24$. This means that galaxies with higher
stellar mass ($\log M_{\star}/M_{\odot}=11.3$) formed their stars at
$z_{\mathrm{form}}\sim 1.6$, while less massive ones ($\log
M_{\star}/M_{\odot}=10.7$) formed their stars at $z_{\mathrm{form}}\sim 1.3$.
Interestingly, the formation epoch of this population of passive galaxies is
found to take place right after the peak of the cosmic star formation rate
density ($z\sim 2$; Madau & Dickinson, 2014). We also find few (22, $\sim
16\%$ of the passive sample) very massive ($\log M_{\star}/M_{\odot}>11$)
galaxies with a $z_{\mathrm{form}}>2.5$ up to 5.

In Figure 8 our results are compared with the existing literature where
similar analysis methods are used. We find very good agreement with the work
of Jørgensen & Chiboucas (2013), who studied spectra of $\sim 80$ cluster
galaxies at $z=0.5-0.9$ comparing observed Lick indices with TMJ11 models.
They found formation redshifts of $z_{\mathrm{form}}\approx 1.24$ and $1.95$
for stellar masses of $\log M_{\star}/M_{\odot}\approx 10.6$ and $11.4$,
respectively. We also find an excellent agreement with the work of Choi et al.
(2014), who analyzed stacked spectra of sSFR-selected passive galaxies at
comparable redshifts and masses with the full spectral fitting method. By
assuming single-burst SFHs, the authors found typical formation epochs of
$z_{\mathrm{form}}\sim 1.5$. Finally, Gallazzi et al. (2014) analyzed ages and
stellar metallicities for $\sim 70$ between star-forming and quiescent
galaxies at $z\sim 0.7$ using age-sensitive $\mathrm{D_{n}4000}$ ,
$\mathrm{H\beta}$ , and $\mathrm{H\delta_{A}}$ +$\mathrm{H\gamma_{A}}$ , and
metal-sensitive ${[\rm Mg_{2}Fe]}$ and ${\rm[MgFe]^{\prime}}$ indices.
Rewriting Equation 8 in terms of the formation time $t_{\mathrm{form}}$ versus
$\log M_{\star}$, we obtain a trend of $-1.26\pm 0.27$ Gyr per decade in mass.
Remarkably, this is in excellent agreement with the trend observed by Carnall
et al. (2019) of $-1.48\pm 0.37$ Gyr per decade in mass, using a very
different approach (full spectral fitting) and assuming more extended (double
power-law) SFHs.

We now move to the analysis of mean stellar metallicities. Our results show no
significant evolution in $\rm[Z/H]$ for passive galaxies since $z\sim 0.7$,
with a median offset of 0.05 dex comparable to the median uncertainty. This
confirms and statistically strengthens earlier results from Gallazzi et al.
(2014) at similar redshifts, and it confirms several works that found similar
metallicities up to $z\approx 2$ (Onodera et al., 2012, 2015; Citro et al.,
2016; Estrada-Carpenter et al., 2019). On the other hand, recent studies
report a significant evolution of $\Delta\rm{[Z/H]}>0.1$ dex. An example is
the work by Beverage et al. (2021) based on a sample of 65 LEGA-C quiescent
galaxies analyzed with a full spectral fitting code. The authors found no
evolution in $\rm{[Mg/Fe]}$ values, but a $\rm{\Delta[Fe/H]}$ (hence
$\rm{\Delta[Z/H]}$) of about 0.2 dex with respect to local
$\log(M_{\star}/M_{\odot})=11$ quiescent stacks. We also study $\rm[Z/H]$
versus $\log\sigma_{\star}$ and $\log M_{\star}$, finding no significant
correlations. Again, this can be attributed to our stricter selection criteria
of the most passive systems, which are known to have shallower
metallicity–mass relations (see Peng et al., 2015; Gallazzi et al., 2014,
2021). Indeed, when dividing them into two $\sigma_{\star}$ bins, we find that
the typical $\rm[Z/H]$ of $\sigma_{\star}>215~{}\rm km\,s^{-1}$ systems is 0.1
dex higher than low-$\sigma$ ones, consistently with the downsizing scenario.

Figure 9: Left panels: distribution of single-burst stellar ages,
metallicities, and $\alpha$-enhancements as a function of redshift for 140
bona fide passive galaxies selected in LEGA-C DR2 (violet points). The error
bars are obtained as the 16th and 84th percentiles of the marginalized
posterior distributions. In the age–redshift panel, we shade the parameter
space not allowed for a general cosmology (gray solid), as well as the
formation redshift assuming a pure passive evolution (dotted lines). Right
panels: median binned relations obtained by dividing the sample into two
$\sigma_{\star}$ regimes, with $\sigma_{\star}=215$ $\rm km\,s^{-1}$ as a
threshold. Each bin contains 5–40 objects. Error bars on the $y$-axis are
errors associated to the median values, while those on the $x$-axis indicate
the bin widths.

Finally, we analyze the mean stellar $\alpha$ enhancements. The
$\rm[\alpha/Fe]$ values of our galaxies are $\sim 0.1$ dex lower than the
local sample. However, as discussed in C, this could be completely explained
by the fact that we could not use Mg indices in our baseline analysis. For the
first time, we identify a positive trend between $\rm[\alpha/Fe]$ and
$\log\sigma_{\star}$ with a slope of ($0.2\pm 0.1$) for a large population of
individual passive galaxies at $z\sim 0.7$. This confirms the trends observed
with stacks of quiescent galaxies at $0.1<z<0.7$ by Choi et al. (2014). We do
not find a significant trend in $\log M_{\star}$. On one hand, local early-
type galaxies show stronger correlations with the gravitational potential well
(of which $\sigma_{\star}$ is a tracer) than with $M_{\star}$ (e.g., Thomas et
al., 2005; Barone et al., 2018). On the other hand, uncertainties in the
$M_{\star}$ derivation may play a role in washing out these relations.

In conclusion, the sample of selected passive galaxies at $z\sim 0.7$ shows
trends in $\mathrm{age}$, $\rm[Z/H]$, and $\rm[\alpha/Fe]$ versus
$\sigma_{\star}$ in agreement with those expected from a passively evolving
population. The age offset, as well as the lack of a significant offset
between the typical values of $\rm[Z/H]$ and $\rm[\alpha/Fe]$ is evidence that
these systems should have formed their stars on short timescales, depleting
the great majority of their gas reservoirs, and experienced passive evolution
since then.

###  4.3 Physical parameters versus $z$

Figure 9 shows the derived stellar population parameters as a function of
redshift. Interestingly, despite not having imposed any cosmological prior on
the age of the galaxies, our derived ages are in all cases in agreement with a
generic cosmological model, never exceeding the age of the Universe at any
redshift. Moreover, the upper envelope of the distribution follows the
decrease expected from the aging of the Universe. At all redshifts, only a few
passive galaxies have ages $\lesssim 2$ Gyr owing to the stringent selection
criteria. The median value is
$\left\langle\right.\\!\\!{\mathrm{age}}\\!\\!\left.\right\rangle=3.01$ Gyr
with a $\pm 0.97$ Gyr $1\sigma$ scatter. This implies a formation time of
$t_{\mathrm{form}}\sim 4$ Gyr after the Big Bang, corresponding to a formation
redshift of $z_{\mathrm{form}}\sim 1.5$ as previously discussed (Section 4.2).

Stellar metallicities have solar or slightly supersolar values,
$\left\langle\right.\\!\\!{\hbox{$\rm[Z/H]$}}\\!\\!\left.\right\rangle=0.08$
dex with a $\pm 0.18$ $1\sigma$ scatter, spanning a very narrow range if
compared to the initial parameter space ($-2.25<\hbox{$\rm[Z/H]$}<0.67$).
Differently from stellar ages, they show no sign of evolution even within the
redshift range explored in this work.

We find slightly supersolar $\rm[\alpha/Fe]$ values,
$\left\langle\right.\\!\\!{\hbox{$\rm[\alpha/Fe]$}}\\!\\!\left.\right\rangle=0.13$
dex with a $\pm 0.11$ $1\sigma$ scatter. In particular, 124 (89%) galaxies
have $\hbox{$\rm[\alpha/Fe]$}>0$, pointing to very short formation timescales,
i.e. before Type Ia supernova explosions can pollute the interstellar medium
with a relatively high amount of iron-peak elements. As for the metallicities,
these star formation timescales indicators do not significantly evolve over
the redshift probed.

Overall, these results confirm the passive evolution within the $\sim 1.6$ Gyr
interval of cosmic time explored in this work shedding light on the
granularity of the physical properties and SFHs. The uniformly small scatter
in $\rm[Z/H]$ and $\rm[\alpha/Fe]$, $\lesssim 0.20$ dex at fixed
$\sigma_{\star}$ or $z$, confirms the large homogeneity of the sample and puts
strong constraints on the duration of the chemical assembly of these systems.

###  4.4 Median binned relations

To provide a comprehensive picture with the results discussed in Section 4.2
and 4.3, we bin $\mathrm{age}$, $\rm[Z/H]$, and $\rm[\alpha/Fe]$ in
$\sigma_{\star}$ and in $z$. Binning in $\sigma_{\star}$ instead of
$M_{\star}$ has the benefit of avoiding model-dependent effects introduced by
SED-fit modeling, other than using an observational quantity. The galaxies
were firstly binned into two $\sigma_{\star}$ using
$\left\langle\right.\\!\\!{\sigma_{\star}}\\!\\!\left.\right\rangle=215$ $\rm
km\,s^{-1}$ (approximately equivalent to $\left\langle\right.\\!\\!{\log
M_{\star}/M_{\odot}}\\!\\!\left.\right\rangle\approx 11$) as a threshold, then
into four equally spaced redshift bins, with $\Delta z=0.075$ from $z=0.6$ to
$z=0.9$. Bins have from $\sim 5$, at higher $z$, to $\sim 40$, at lower $z$,
objects. To each bin, we assign a mean $z$ value and a median ($\mathrm{age}$,
$\rm[Z/H]$, $\rm[\alpha/Fe]$) with associated uncertainty. Results are shown
in the right panels of Figure 9.

Clearly, the median properties of the analyzed passive galaxies follow a
downsizing pattern. At each cosmic epoch, stellar populations hosted in
galaxies with higher mass are older, more metal rich and $alpha$ enhanced.
This suggests that their formation occurred at earlier times, with a
difference of $\Delta\mathrm{age}\approx 0.5$ Gyr, and on shorter time scales
with respect to less massive ones. We note that these trends were already
qualitatively confirmed from the analysis of the main absorption indices
(Section 2.4).

Last but not least, it is remarkable that for each $\sigma_{\star}$ regime we
find a clear, almost parallel, age–redshift relation. The study of their
differential evolution will allow us to perform cosmological studies using the
cosmic chronometer approach (Borghi et al., 2022).

##  5 Conclusions

In this work, we take advantage of the public Data Release 2 of the LEGA-C
spectroscopic survey to place constraints on the stellar population properties
of individual massive and passive galaxies at $0.6<z<0.9$. Based on a robust
spectral analysis of Lick indices, our aim is to characterize this population
and to explore the reliability of using these galaxies as cosmic chronometers.
Our main results are summarized below.

  1. 1.

We select a pure sample of 350 passive galaxies at $z\sim 0.7$ combining a
photometric NUVrJ selection, a spectroscopic EW[O ii] cut, and a careful
visual inspection of individual spectra to further remove galaxies with
significant emission-lines (Figure 1). As confirmed by the stacked spectrum
(Figure 2), no underlying emission line components are present in the sample,
confirming its high purity. Selected passive galaxies have a median observed
velocity dispersion of
$\left\langle\right.\\!\\!{\sigma_{\star}}\\!\\!\left.\right\rangle=206$ $\rm
km\,s^{-1}$, stellar mass of $\left\langle\right.\\!\\!{\log
M_{\star}/M_{\odot}}\\!\\!\left.\right\rangle=10.95$, and very low specific
star formation rate $\left\langle\right.\\!\\!{\rm\log
sSFR/yr}\\!\\!\left.\right\rangle=-12.1$. Most of them have an early-type
morphology, but there is also a nonnegligible percentage of systems (about
one-third) with an intermediate morphology.

  2. 2.

We develop, validate, and publicly release PyLick, a flexible Python tool to
measure absorption features, implementing several different index definitions
(see A). This allows us to measure spectral indices over a wide wavelength
range in LEGA-C data, extending the current public catalog of Lick indices by
Straatman et al. (2018) and enabling a more detailed exploration of the
dependence of our results on different index combinations.

  3. 3.

We introduce the H/K ratio, a new diagnostic feature defined as the ratio of
pseudo-Lick indices $\mathrm{Ca\>\textsc{ii}\>K}$ and
$\mathrm{Ca\>\textsc{ii}\>H}$ (Figure 3). We verify that it is an excellent
tracer of potential contamination of the sample due to star-forming or young
populations, confirming that our sample is compatible with no or negligible
contamination, with $\left\langle\right.\\!\\!{\rm
H/K}\\!\\!\left.\right\rangle=0.96\pm 0.08~{}(1\sigma)$. Moreover, a selection
based on ${\rm H/K}<1.1$ is found to reproduce a NUVrJ selection (Ilbert et
al., 2013) or a ${\rm\log sSFR/yr}<-11$ cut (Figure 6) with a small percentage
of incompleteness ($\sim 15\%$) or contamination ($\sim 15\%$) while requiring
a much narrower wavelength range.

  4. 4.

Using an optimized combination of Lick indices (namely, $\mathrm{H\delta_{A}}$
, $\mathrm{CN_{1}}$ , $\mathrm{CN_{2}}$ , $\mathrm{Ca4227}$ , $\mathrm{G4300}$
, $\mathrm{H\gamma_{A}}$ , $\mathrm{H\gamma_{F}}$ , $\mathrm{Fe4383}$ ,
$\mathrm{Fe4531}$ , and $\mathrm{C_{2}4668}$ ), we measure single-burst
stellar $\mathrm{age}$, $\rm[Z/H]$, and $\rm[\alpha/Fe]$ for 140 passive
galaxies, without assuming cosmological priors on the maximum value of
$\mathrm{age}$ as a function of redshift. We also performed an extended
analysis to assess the impact of different choices of indices, verifying that
our findings are robust against the choice of a different combination of
indices (C).

  5. 5.

We find trends between $\log\mathrm{age}$, $\rm[Z/H]$, $\rm[\alpha/Fe]$, and
the stellar velocity dispersion consistent with those expected from a
passively evolving population, with slopes of ($0.5\pm 0.1$), ($0.3\pm 0.2$),
and ($0.2\pm 0.1$), respectively (Figure 7). This analysis shows, for the
first time using individual galaxies, that a relation between $\rm[\alpha/Fe]$
(a star formation timescale proxy) and $\sigma_{\star}$ is already in place at
$z\sim 0.7$. Moreover, the age difference of 5.5 Gyr between our sample and
local ETGs can be entirely accounted for by a pure passive evolution. Assuming
a standard $\Lambda$CDM cosmology, the relation between formation redshifts
and galaxy stellar masses is found to agree with several previous analyses
(Figure 8), confirming that this population of massive galaxies forms at
$z_{\mathrm{form}}\sim 1.3(1.6)$ at masses $\log
M_{\star}/M_{\odot}=10.7(11.3)$, after the peak of the cosmic star formation
rate density.

  6. 6.

Even if we do not impose any cosmological prior to the age of the population,
the obtained $\mathrm{age}$–redshift evolution is consistent with a
$\Lambda$CDM universe (Figure 9). The stellar $\rm[Z/H]$ and $\rm[\alpha/Fe]$
do not evolve significantly over $z$ and their values are slightly supersolar,
$\left\langle\right.\\!\\!{\hbox{$\rm[Z/H]$}}\\!\\!\left.\right\rangle=0.08\pm
0.18\,(1\sigma)$ dex, and supersolar,
$\left\langle\right.\\!\\!{\hbox{$\rm[\alpha/Fe]$}}\\!\\!\left.\right\rangle=0.13\pm
0.11\,(1\sigma)$ dex, compatibly with their local counterparts.

  7. 7.

Finally, the analysis of median binned relations confirms the downsizing
scenario and the passive nature of this population. Remarkably, we obtain two
clear nearly parallel $\mathrm{age}-z$ relations for both the higher
($\sigma_{\star}\approx 230~{}\rm km\,s^{-1}$) and the lower
($\sigma_{\star}\approx 200~{}\rm km\,s^{-1}$) mass regimes. This difference
of $\Delta\mathrm{age}\approx 0.5$ Gyr can be interpreted as a delay in
formation time between the two, with later formation epochs for the population
of less massive galaxies.

Overall, our analysis of individual galaxies confirms the existence of a
population of passively evolving galaxies at intermediate redshift that
follows a downsizing pattern. In a subsequent paper, we will make use of their
median $\mathrm{age}-z$ relation to derive constraints on cosmological
parameters and obtain a direct measure of the Hubble parameter $H(z)$ using
the cosmic chronometers approach. In this context, we also plan to study in
more detail the effect of assuming different SFHs. In this way, we will be
able to study for the first time at these redshifts the detailed stellar
population properties of passive galaxies and their underlying cosmology,
jointly.

In the more distant future, dedicated large spectroscopic surveys would be
crucial to extend the present analysis, providing at the same time unique
opportunities for gravitational wave astronomy and cosmology (Palmese et al.,
2019).

We thank the anonymous referee for the constructive comments and suggestions
that helped in improving this paper. This work is based on data products from
observations made with ESO Telescopes at the La Silla Paranal Observatory
under program ID 194.AF2005(A-N). We thank the LEGA-C team for making their
dataset public, Daniel Thomas for supplying higher-resolution SSP models than
those available publicly, and Nicholas Scott for useful suggestions for the
analysis of spectral indices. N.B. and M.M. acknowledge support from MIUR,
PRIN 2017 (grant 20179ZF5KS). M.M., A.C., and L.P. acknowledge the grants ASI
n.I/023/12/0 and ASI n.2018-23-HH.0. A.C. acknowledges the support from grant
PRIN MIUR 2017 - 20173ML3WW_001.

## References

  * Arnouts et al. (2007) Arnouts, S., Walcher, C. J., Le Fèvre, O., et al. 2007, A&A, 476, 137, doi: 10.1051/0004-6361:20077632
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Balogh et al. (1999) Balogh, M. L., Morris, S. L., Yee, H. K. C., Carlberg, R. G., & Ellingson, E. 1999, ApJ, 527, 54, doi: 10.1086/308056
  * Barone et al. (2018) Barone, T. M., D’Eugenio, F., Colless, M., et al. 2018, ApJ, 856, 64, doi: 10.3847/1538-4357/aaaf6e
  * Beifiori et al. (2011) Beifiori, A., Maraston, C., Thomas, D., & Johansson, J. 2011, A&A, 531, A109, doi: 10.1051/0004-6361/201016323
  * Belli et al. (2019) Belli, S., Newman, A. B., & Ellis, R. S. 2019, ApJ, 874, 17, doi: 10.3847/1538-4357/ab07af
  * Beverage et al. (2021) Beverage, A. G., Kriek, M., Conroy, C., et al. 2021, ApJ, 917, L1, doi: 10.3847/2041-8213/ac12cd
  * Bluck et al. (2020) Bluck, A. F. L., Maiolino, R., Piotrowska, J. M., et al. 2020, MNRAS, 499, 230, doi: 10.1093/mnras/staa2806
  * Borghi et al. (2022) Borghi, N., Moresco, M., Cimatti, A., & et al. 2022, arXiv e-prints, arXiv:2106.14894.  https://arxiv.org/abs/2106.14894
  * Bruzual (1983) Bruzual, G. 1983, ApJ, 273, 105, doi: 10.1086/161352
  * Bruzual & Charlot (2003) Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000, doi: 10.1046/j.1365-8711.2003.06897.x
  * Burstein et al. (1984) Burstein, D., Faber, S. M., Gaskell, C. M., & Krumm, N. 1984, ApJ, 287, 586, doi: 10.1086/162718
  * Capozziello et al. (2014) Capozziello, S., Farooq, O., Luongo, O., & Ratra, B. 2014, Phys. Rev. D, 90, 044016, doi: 10.1103/PhysRevD.90.044016
  * Cappellari et al. (2013a) Cappellari, M., McDermid, R. M., Alatalo, K., et al. 2013a, MNRAS, 432, 1862, doi: 10.1093/mnras/stt644
  * Cappellari et al. (2013b) Cappellari, M., Scott, N., Alatalo, K., et al. 2013b, MNRAS, 432, 1709, doi: 10.1093/mnras/stt562
  * Cardiel et al. (1998) Cardiel, N., Gorgas, J., Cenarro, J., & Gonzalez, J. J. 1998, A&AS, 127, 597, doi: 10.1051/aas:1998123
  * Carnall et al. (2018) Carnall, A. C., McLure, R. J., Dunlop, J. S., & Davé, R. 2018, MNRAS, 480, 4379, doi: 10.1093/mnras/sty2169
  * Carnall et al. (2019) Carnall, A. C., McLure, R. J., Dunlop, J. S., et al. 2019, MNRAS, 490, 417, doi: 10.1093/mnras/stz2544
  * Carson & Nichol (2010) Carson, D. P., & Nichol, R. C. 2010, MNRAS, 408, 213, doi: 10.1111/j.1365-2966.2010.17151.x
  * Cassisi et al. (1997) Cassisi, S., Castellani, M., & Castellani, V. 1997, A&A, 317, 108.  https://ui.adsabs.harvard.edu/abs/1997A&A...317..108C
  * Cenarro et al. (2001) Cenarro, A. J., Cardiel, N., Gorgas, J., et al. 2001, MNRAS, 326, 959, doi: 10.1046/j.1365-8711.2001.04688.x
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763, doi: 10.1086/376392
  * Choi et al. (2014) Choi, J., Conroy, C., Moustakas, J., et al. 2014, ApJ, 792, 95, doi: 10.1088/0004-637X/792/2/95
  * 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
  * Cimatti et al. (2019) Cimatti, A., Fraternali, F., & Nipoti, C. 2019, Introduction to Galaxy Formation and Evolution: From Primordial Gas to Present-Day Galaxies (Cambridge University Press).  https://ui.adsabs.harvard.edu/abs/2019igfe.book.....C
  * Cimatti et al. (2004) Cimatti, A., Daddi, E., Renzini, A., et al. 2004, Nature, 430, 184, doi: 10.1038/nature02668
  * Cirasuolo et al. (2020) Cirasuolo, M., Fairley, A., Rees, P., et al. 2020, The Messenger, 180, 10, doi: 10.18727/0722-6691/5195
  * Citro et al. (2016) Citro, A., Pozzetti, L., Moresco, M., & Cimatti, A. 2016, A&A, 592, A19, doi: 10.1051/0004-6361/201527772
  * Cleveland & Devlin (1988) Cleveland, W. S., & Devlin, S. J. 1988, Journal of the American Statistical Association, 83, 596, doi: 10.1080/01621459.1988.10478639
  * Concas et al. (2017) Concas, A., Pozzetti, L., Moresco, M., & Cimatti, A. 2017, MNRAS, 468, 1747, doi: 10.1093/mnras/stx453
  * Conroy (2013) Conroy, C. 2013, ARA&A, 51, 393, doi: 10.1146/annurev-astro-082812-141017
  * Conroy et al. (2014) Conroy, C., Graves, G. J., & van Dokkum, P. G. 2014, ApJ, 780, 33, doi: 10.1088/0004-637X/780/1/33
  * Cowie et al. (1996) Cowie, L. L., Songaila, A., Hu, E. M., & Cohen, J. G. 1996, AJ, 112, 839, doi: 10.1086/118058
  * Daddi et al. (2005) Daddi, E., Renzini, A., Pirzkal, N., et al. 2005, ApJ, 626, 680, doi: 10.1086/430104
  * Davidzon et al. (2017) Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, A&A, 605, A70, doi: 10.1051/0004-6361/201730419
  * de La Rosa et al. (2011) de La Rosa, I. G., La Barbera, F., Ferreras, I., & de Carvalho, R. R. 2011, MNRAS, 418, L74, doi: 10.1111/j.1745-3933.2011.01146.x
  * Dressler et al. (1999) Dressler, A., Smail, I., Poggianti, B. M., et al. 1999, ApJS, 122, 51, doi: 10.1086/313213
  * Estrada-Carpenter et al. (2019) Estrada-Carpenter, V., Papovich, C., Momcheva, I., et al. 2019, ApJ, 870, 133, doi: 10.3847/1538-4357/aaf22e
  * Faber (1973) Faber, S. M. 1973, ApJ, 179, 731, doi: 10.1086/151912
  * Falcón-Barroso et al. (2011) Falcón-Barroso, J., Sánchez-Blázquez, P., Vazdekis, A., et al. 2011, A&A, 532, A95, doi: 10.1051/0004-6361/201116842
  * Fanelli et al. (1992) Fanelli, M. N., O’Connell, R. W., Burstein, D., & Wu, C.-C. 1992, ApJS, 82, 197, doi: 10.1086/191714
  * Fanfani (2019) Fanfani, V. 2019, PhD thesis, University of Bologna.  https://amslaurea.unibo.it/17900/
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Foreman-Mackey et al. (2019) Foreman-Mackey, D., Farr, W., Sinha, M., et al. 2019, , 4, 1864, doi: 10.21105/joss.01864
  * Franzetti et al. (2007) Franzetti, P., Scodeggio, M., Garilli, B., et al. 2007, A&A, 465, 711, doi: 10.1051/0004-6361:20065942
  * Gallazzi et al. (2014) Gallazzi, A., Bell, E. F., Zibetti, S., Brinchmann, J., & Kelson, D. D. 2014, ApJ, 788, 72, doi: 10.1088/0004-637X/788/1/72
  * Gallazzi et al. (2006) Gallazzi, A., Charlot, S., Brinchmann, J., & White, S. D. M. 2006, MNRAS, 370, 1106, doi: 10.1111/j.1365-2966.2006.10548.x
  * 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
  * Gallazzi et al. (2021) Gallazzi, A. R., Pasquali, A., Zibetti, S., & Barbera, F. L. 2021, MNRAS, 502, 4457, doi: 10.1093/mnras/stab265
  * Gómez-Valent & Amendola (2018) Gómez-Valent, A., & Amendola, L. 2018, J. Cosmology Astropart. Phys, 2018, 051, doi: 10.1088/1475-7516/2018/04/051
  * Haridasu et al. (2018) Haridasu, B. S., Luković, V. V., Moresco, M., & Vittorio, N. 2018, J. Cosmology Astropart. Phys, 2018, 015, doi: 10.1088/1475-7516/2018/10/015
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hinton (2016) Hinton, S. R. 2016, , 1, 00045, doi: 10.21105/joss.00045
  * Hoaglin et al. (1983) Hoaglin, D. C., Mosteller, F., & Tukey, J. W. 1983, Understanding robust and exploratory data anlysis (Wiley).  https://ui.adsabs.harvard.edu/abs/1983ured.book.....H
  * Hogg & Foreman-Mackey (2018) Hogg, D. W., & Foreman-Mackey, D. 2018, ApJS, 236, 11, doi: 10.3847/1538-4365/aab76e
  * Hubble (1936) Hubble, E. 1936, The Realm of the Nebulae (Yale University Press)
  * Hunter (2007) Hunter, J. D. 2007, , 9, 90, doi: 10.1109/MCSE.2007.55
  * Ilbert et al. (2009) Ilbert, O., Capak, P., Salvato, M., et al. 2009, ApJ, 690, 1236, doi: 10.1088/0004-637X/690/2/1236
  * Ilbert et al. (2013) Ilbert, O., McCracken, H. J., Le Fèvre, O., et al. 2013, A&A, 556, A55, doi: 10.1051/0004-6361/201321100
  * Ilbert et al. (2015) Ilbert, O., Arnouts, S., Le Floc’h, E., et al. 2015, A&A, 579, A2, doi: 10.1051/0004-6361/201425176
  * Jimenez & Loeb (2002) Jimenez, R., & Loeb, A. 2002, ApJ, 573, 37, doi: 10.1086/340549
  * Johansson et al. (2012) Johansson, J., Thomas, D., & Maraston, C. 2012, MNRAS, 421, 1908, doi: 10.1111/j.1365-2966.2011.20316.x
  * Jørgensen & Chiboucas (2013) Jørgensen, I., & Chiboucas, K. 2013, AJ, 145, 77, doi: 10.1088/0004-6256/145/3/77
  * Kauffmann et al. (2003) Kauffmann, G., Heckman, T. M., White, S. D. M., et al. 2003, MNRAS, 341, 33, doi: 10.1046/j.1365-8711.2003.06291.x
  * Korn et al. (2005) Korn, A. J., Maraston, C., & Thomas, D. 2005, A&A, 438, 685, doi: 10.1051/0004-6361:20042126
  * Kriek et al. (2019) Kriek, M., Price, S. H., Conroy, C., et al. 2019, ApJ, 880, L31, doi: 10.3847/2041-8213/ab2e75
  * La Barbera et al. (2014) La Barbera, F., Pasquali, A., Ferreras, I., et al. 2014, MNRAS, 445, 1977, doi: 10.1093/mnras/stu1626
  * Laigle et al. (2016) Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24, doi: 10.3847/0067-0049/224/2/24
  * Laureijs et al. (2011) Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv e-prints, arXiv:1110.3193.  https://ui.adsabs.harvard.edu/abs/2011arXiv1110.3193L
  * Lee et al. (2009) Lee, H.-c., Worthey, G., Dotter, A., et al. 2009, ApJ, 694, 902, doi: 10.1088/0004-637X/694/2/902
  * Leonardi & Rose (1996) Leonardi, A. J., & Rose, J. A. 1996, AJ, 111, 182, doi: 10.1086/117772
  * L’Huillier & Shafieloo (2017) L’Huillier, B., & Shafieloo, A. 2017, J. Cosmology Astropart. Phys, 2017, 015, doi: 10.1088/1475-7516/2017/01/015
  * Lin et al. (2021) Lin, W., Chen, X., & Mack, K. J. 2021, ApJ, 920, 159, doi: 10.3847/1538-4357/ac12cf
  * Lin et al. (2020) Lin, W., Mack, K. J., & Hou, L. 2020, ApJ, 904, L22, doi: 10.3847/2041-8213/abc894
  * Longhetti et al. (1999) Longhetti, M., Bressan, A., Chiosi, C., & Rampazzo, R. 1999, A&A, 345, 419.  https://ui.adsabs.harvard.edu/abs/1999A&A...345..419L
  * Lonoce et al. (2014) Lonoce, I., Longhetti, M., Saracco, P., Gargiulo, A., & Tamburri, S. 2014, MNRAS, 444, 2048, doi: 10.1093/mnras/stu1593
  * Lonoce et al. (2020) Lonoce, I., Maraston, C., Thomas, D., et al. 2020, MNRAS, 492, 326, doi: 10.1093/mnras/stz3404
  * Lonoce et al. (2015) Lonoce, I., Longhetti, M., Maraston, C., et al. 2015, MNRAS, 454, 3912, doi: 10.1093/mnras/stv2150
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415, doi: 10.1146/annurev-astro-081811-125615
  * Maraston (2005) Maraston, C. 2005, MNRAS, 362, 799, doi: 10.1111/j.1365-2966.2005.09270.x
  * Maraston et al. (2009) Maraston, C., Nieves Colmenárez, L., Bender, R., & Thomas, D. 2009, A&A, 493, 425, doi: 10.1051/0004-6361:20066907
  * McDermid et al. (2015) McDermid, R. M., Alatalo, K., Blitz, L., et al. 2015, MNRAS, 448, 3484, doi: 10.1093/mnras/stv105
  * Mignoli et al. (2009) Mignoli, M., Zamorani, G., Scodeggio, M., et al. 2009, A&A, 493, 39, doi: 10.1051/0004-6361:200810520
  * Moresco (2015) Moresco, M. 2015, MNRAS, 450, L16, doi: 10.1093/mnrasl/slv037
  * Moresco et al. (2011) Moresco, M., Jimenez, R., Cimatti, A., & Pozzetti, L. 2011, J. Cosmology Astropart. Phys, 2011, 045, doi: 10.1088/1475-7516/2011/03/045
  * Moresco et al. (2020) Moresco, M., Jimenez, R., Verde, L., Cimatti, A., & Pozzetti, L. 2020, ApJ, 898, 82, doi: 10.3847/1538-4357/ab9eb0
  * Moresco et al. (2018) Moresco, M., Jimenez, R., Verde, L., et al. 2018, ApJ, 868, 84, doi: 10.3847/1538-4357/aae829
  * Moresco & Marulli (2017) Moresco, M., & Marulli, F. 2017, MNRAS, 471, L82, doi: 10.1093/mnrasl/slx112
  * Moresco et al. (2012a) Moresco, M., Verde, L., Pozzetti, L., Jimenez, R., & Cimatti, A. 2012a, J. Cosmology Astropart. Phys, 2012, 053, doi: 10.1088/1475-7516/2012/07/053
  * Moresco et al. (2012b) Moresco, M., Cimatti, A., Jimenez, R., et al. 2012b, J. Cosmology Astropart. Phys, 2012, 006, doi: 10.1088/1475-7516/2012/08/006
  * Moresco et al. (2013) Moresco, M., Pozzetti, L., Cimatti, A., et al. 2013, A&A, 558, A61, doi: 10.1051/0004-6361/201321797
  * Moresco et al. (2016) —. 2016, J. Cosmology Astropart. Phys, 2016, 014, doi: 10.1088/1475-7516/2016/05/014
  * Muzzin et al. (2013) Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJ, 777, 18, doi: 10.1088/0004-637X/777/1/18
  * Nunes et al. (2016) Nunes, R. C., Pan, S., & Saridakis, E. N. 2016, Phys. Rev. D, 94, 023508, doi: 10.1103/PhysRevD.94.023508
  * Onodera et al. (2012) Onodera, M., Renzini, A., Carollo, M., et al. 2012, ApJ, 755, 26, doi: 10.1088/0004-637X/755/1/26
  * Onodera et al. (2015) Onodera, M., Carollo, C. M., Renzini, A., et al. 2015, ApJ, 808, 161, doi: 10.1088/0004-637X/808/2/161
  * Pacifici et al. (2016) Pacifici, C., Kassin, S. A., Weiner, B. J., et al. 2016, ApJ, 832, 79, doi: 10.3847/0004-637X/832/1/79
  * Palmese et al. (2019) Palmese, A., Graur, O., Annis, J. T., et al. 2019, BAAS, 51, 310.  https://ui.adsabs.harvard.edu/abs/2019BAAS...51c.310P
  * Peng et al. (2015) Peng, Y., Maiolino, R., & Cochrane, R. 2015, Nature, 521, 192, doi: 10.1038/nature14439
  * Peng et al. (2010) Peng, Y.-j., Lilly, S. J., Kovač, K., et al. 2010, ApJ, 721, 193, doi: 10.1088/0004-637X/721/1/193
  * Pozzetti et al. (2010) Pozzetti, L., Bolzonella, M., Zucca, E., et al. 2010, A&A, 523, A13, doi: 10.1051/0004-6361/200913020
  * Renzini (2006) Renzini, A. 2006, ARA&A, 44, 141, doi: 10.1146/annurev.astro.44.051905.092450
  * Rose (1984) Rose, J. A. 1984, AJ, 89, 1238, doi: 10.1086/113618
  * Rose (1985) —. 1985, AJ, 90, 1927, doi: 10.1086/113898
  * Rousseeuw (1984) Rousseeuw, P. J. 1984, Journal of the American Statistical Association, 79, 871, doi: 10.1080/01621459.1984.10477105
  * Roy (2020) Roy, V. 2020, Annual Review of Statistics and Its Application, 7, 387, doi: 10.1146/annurev-statistics-031219-041300
  * Salpeter (1955) Salpeter, E. E. 1955, ApJ, 121, 161, doi: 10.1086/145971
  * 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
  * Santini et al. (2021) Santini, P., Castellano, M., Merlin, E., et al. 2021, A&A, 652, A30, doi: 10.1051/0004-6361/202039738
  * Sapone et al. (2014) Sapone, D., Majerotto, E., & Nesseris, S. 2014, Phys. Rev. D, 90, 023012, doi: 10.1103/PhysRevD.90.023012
  * Scarlata et al. (2007) Scarlata, C., Carollo, C. M., Lilly, S., et al. 2007, ApJS, 172, 406, doi: 10.1086/516582
  * Scott et al. (2017) Scott, N., Brough, S., Croom, S. M., et al. 2017, MNRAS, 472, 2833, doi: 10.1093/mnras/stx2166
  * Seikel et al. (2012) Seikel, M., Yahya, S., Maartens, R., & Clarkson, C. 2012, Phys. Rev. D, 86, 083001, doi: 10.1103/PhysRevD.86.083001
  * Serra & Trager (2007) Serra, P., & Trager, S. C. 2007, MNRAS, 374, 769, doi: 10.1111/j.1365-2966.2006.11188.x
  * Singh et al. (2013) Singh, R., van de Ven, G., Jahnke, K., et al. 2013, A&A, 558, A43, doi: 10.1051/0004-6361/201322062
  * Siudek et al. (2017) Siudek, M., Małek, K., Scodeggio, M., et al. 2017, A&A, 597, A107, doi: 10.1051/0004-6361/201628951
  * Solà et al. (2017) Solà, J., Gómez-Valent, A., & de Cruz Pérez, J. 2017, , 32, 1750054, doi: 10.1142/S0217732317500547
  * Spergel et al. (2015) Spergel, D., Gehrels, N., Baltay, C., et al. 2015, arXiv e-prints, arXiv:1503.03757.  https://ui.adsabs.harvard.edu/abs/2015arXiv150303757S
  * Spilker et al. (2018) Spilker, J., Bezanson, R., Barišić, I., et al. 2018, ApJ, 860, 103, doi: 10.3847/1538-4357/aac438
  * Spinrad et al. (1997) Spinrad, H., Dey, A., Stern, D., et al. 1997, ApJ, 484, 581, doi: 10.1086/304381
  * Spitoni et al. (2017) Spitoni, E., Vincenzo, F., & Matteucci, F. 2017, A&A, 599, A6, doi: 10.1051/0004-6361/201629745
  * Straatman et al. (2018) Straatman, C. M. S., van der Wel, A., Bezanson, R., et al. 2018, ApJS, 239, 27, doi: 10.3847/1538-4365/aae37a
  * Tacchella et al. (2021) Tacchella, S., Conroy, C., Faber, S. M., et al. 2021, arXiv e-prints, arXiv:2102.12494.  https://ui.adsabs.harvard.edu/abs/2021arXiv210212494T
  * Tamura et al. (2016) Tamura, N., Takato, N., Shimono, A., et al. 2016, in Ground-based and Airborne Instrumentation for Astronomy VI, ed. C. J. Evans, L. Simard, & H. Takami, Vol. 9908, 99081M, doi: 10.1117/12.2232103
  * Tanaka et al. (2019) Tanaka, M., Valentino, F., Toft, S., et al. 2019, ApJ, 885, L34, doi: 10.3847/2041-8213/ab4ff3
  * Taylor (2005) Taylor, M. B. 2005, in , ed. P. Shopbell, M. Britton, & R. Ebert, Vol. 347, 29.  https://ui.adsabs.harvard.edu/abs/2005ASPC..347...29T
  * Thomas et al. (2005) Thomas, D., Maraston, C., Bender, R., & Mendes de Oliveira, C. 2005, ApJ, 621, 673, doi: 10.1086/426932
  * Thomas et al. (2011) Thomas, D., Maraston, C., & Johansson, J. 2011, MNRAS, 412, 2183, doi: 10.1111/j.1365-2966.2010.18049.x
  * Thomas et al. (2010) Thomas, D., Maraston, C., Schawinski, K., Sarzi, M., & Silk, J. 2010, MNRAS, 404, 1775, doi: 10.1111/j.1365-2966.2010.16427.x
  * Trager et al. (1998) Trager, S. C., Worthey, G., Faber, S. M., Burstein, D., & González, J. J. 1998, ApJS, 116, 1, doi: 10.1086/313099
  * Treu et al. (2005) Treu, T., Ellis, R. S., Liao, T. X., et al. 2005, ApJ, 633, 174, doi: 10.1086/444585
  * Tripicco & Bell (1995) Tripicco, M. J., & Bell, R. A. 1995, AJ, 110, 3035, doi: 10.1086/117744
  * Trussler et al. (2021) Trussler, J., Maiolino, R., Maraston, C., et al. 2021, MNRAS, 500, 4469, doi: 10.1093/mnras/staa3545
  * Vagnozzi et al. (2021) Vagnozzi, S., Loeb, A., & Moresco, M. 2021, ApJ, 908, 84, doi: 10.3847/1538-4357/abd4df
  * Valentino et al. (2020) Valentino, F., Tanaka, M., Davidzon, I., et al. 2020, ApJ, 889, 93, doi: 10.3847/1538-4357/ab64dc
  * Valkenburg et al. (2014) Valkenburg, W., Marra, V., & Clarkson, C. 2014, MNRAS, 438, L6, doi: 10.1093/mnrasl/slt140
  * van der Wel et al. (2016) van der Wel, A., Noeske, K., Bezanson, R., et al. 2016, ApJS, 223, 29, doi: 10.3847/0067-0049/223/2/29
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Williams et al. (2009) Williams, R. J., Quadri, R. F., Franx, M., van Dokkum, P., & Labbé, I. 2009, ApJ, 691, 1879, doi: 10.1088/0004-637X/691/2/1879
  * Worthey (1998) Worthey, G. 1998, PASP, 110, 888, doi: 10.1086/316203
  * Worthey et al. (1994) Worthey, G., Faber, S. M., Gonzalez, J. J., & Burstein, D. 1994, ApJS, 94, 687, doi: 10.1086/192087
  * Worthey & Ottaviani (1997) Worthey, G., & Ottaviani, D. L. 1997, ApJS, 111, 377, doi: 10.1086/313021
  * Wu et al. (2018) Wu, P.-F., van der Wel, A., Gallazzi, A., et al. 2018, ApJ, 855, 85, doi: 10.3847/1538-4357/aab0a6
  * Yan et al. (2006) Yan, R., Newman, J. A., Faber, S. M., et al. 2006, ApJ, 648, 281, doi: 10.1086/505629
  * Yang et al. (2018) Yang, W., Mukherjee, A., Di Valentino, E., & Pan, S. 2018, Phys. Rev. D, 98, 123527, doi: 10.1103/PhysRevD.98.123527

##  Appendix A Measuring spectral indices with PyLick

Absorption-line features produced by atomic and molecular absorbers in stellar
photospheres can unlock crucial information about the physical properties of
stellar populations. Here we present PyLick, a flexible tool to measure
spectral indices and associated uncertainties in galaxy spectra that includes
different index definitions and measuring methods. The code is entirely
written in Python language and object oriented. We developed this tool in the
view of upcoming observational facilities to facilitate the detailed analysis
of high-resolution spectroscopy (e.g., VLT/MOONS, Cirasuolo et al. 2020;
Subaru/PFS, Tamura et al. 2016), and a quicker analysis of lower-resolution
spectroscopic surveys (e.g., Euclid, Laureijs et al., 2011).

PyLick is currently based on five modules: “io” implements the built-in
methods for spectra I/O and preliminary analysis; “indices” loads the index
library including passband definitions; “measure” contains the methods to
derive indices values and errors; “plot” deals with data visualization.
Ultimately, the module “analysis” contains two classes: Galaxy, which is
optimized to analyze a single spectrum, and Catalog, optimized to perform the
analysis of a bulk set of spectra.

In the current version, 54 indices are already defined in the index library:

  * •

UV line indices: 11 far-UV and 8 mid-UV (Fanelli et al. 1992, see also
Maraston et al. 2009);

  * •

Line break indices B2640, B2900 (Spinrad et al., 1997);

  * •

The MgUV index (Daddi et al., 2005);

  * •

4000 Å discontinuity indices: D4000 and Dn4000 (Bruzual, 1983; Balogh et al.,
1999, respectively);

  * •

Lick indices (Worthey & Ottaviani, 1997; Trager et al., 1998);

  * •

CaII H and CaII K “pseudo-Lick” indices (Fanfani, 2019);

  * •

Generic indices CaT, PaT, and CaT∗ (Cenarro et al., 2001).

However, new indices can be easily introduced by defining the wavelength
regions of interest and the measuring method in a custom library file.

Different measuring methods are included: lick_atomic (Equation 1),
lick_molecular (Equation 2), break_nu (Equation 4), and break_lb (as Equation
4, but integrating over $F(\lambda)\mathrm{d}\lambda$, instead of
$F(\nu)\mathrm{d}\nu$). More specific methods are directly defined in the
“measure” module. This is the case for the $\rm{Mg_{UV}}$ index that traces
the absorption bump present at 2640–2850 Å (Daddi et al., 2005):

$\rm{Mg_{UV}}=\frac{2\int^{2725}_{2625}F(\lambda)\rm{d}\lambda}{\int^{2625}_{2525}F(\lambda)\rm{d}\lambda+\int^{2825}_{2725}F(\lambda)\rm{d}\lambda},$
(A1)

and for the calcium triplet index

$\rm{CaT^{\ast}}=\rm{CaT}-0.93\,\rm{PaT},$ (A2)

which traces the strength of the CaII lines ($\rm{CaT}$) corrected from the
contamination by Paschen lines ($\rm{PaT}$), as presented in detail in Cenarro
et al. (2001).

Index errors are evaluated following the signal-to-noise method proposed by
Cardiel et al. (1998). The code is also able to handle bad pixels, which
should be passed as a Boolean array using the spec_mask argument. The user can
choose a bad-to-total pixel ratio (we considered BPR=0.15 in this work) above
which the measurement is not performed. Otherwise, a zeroth- or first-order
interpolation is done over the bad pixels prior to the measurement. Finally,
publication-quality figures can be produced with the plotting routines (e.g.,
Figures 2 and 3).

PyLick has been released as an open-source project and is available on GitLab
with extensive documentation and notebook examples 555The code is available at
https://gitlab.com/mmoresco/pylick/.

Figure 10: Differences between the indices values measured with PyLick and
those published in LEGA-C DR2 (Straatman et al., 2018), quantified as
$\eta=(I_{\mathrm{PyLick}}-I_{\mathrm{LEGA-C}})/I_{\mathrm{LEGA-C}}$, except
for the $\mathrm{D_{n}4000}$ break and the molecular indices $\mathrm{CN_{1}}$
and $\mathrm{CN_{2}}$ , for which the absolute differences
$(I_{\mathrm{PyLick}}-I_{\mathrm{LEGA-C}})$ are displayed. Gray bands
correspond to the $\eta=\pm 5\times 10^{-4}$ regions. Note: $y$-axis limits
are set to 5th–95th percentiles for illustrative purposes.

###  A.1 Validation of the code with LEGA-C data

To confirm the reliability and robustness of PyLick, we compare the measured
indices and errors obtained on the unconvolved LEGA-C spectra
$I_{\mathrm{PyLick}}$, with those released in the LEGA-C DR2 catalog
$I_{\mathrm{LEGA-C}}$ (Straatman et al., 2018). The comparison is performed
within the passive sample to minimize differences due to emission-line
subtraction performed in LEGA-C DR2 (see Straatman et al., 2018). Differences
are computed as

$\eta(I)=\frac{I_{\mathrm{PyLick}}-I_{\mathrm{LEGA-C}}}{\left|I_{\mathrm{LEGA-C}}\right|}.$
(A3)

The same analysis is performed for uncertainties of the indices (the notation
$\eta(\sigma)$ will be used). For a fair comparison, we multiplied our
uncertainties for the same coefficients applied in the LEGA-C DR2 pipeline
(see Straatman et al., 2018, Tab. 3). Figure 10 shows the results for values
of the indices, sorting each index available in LEGA-C DR2 by increasing
wavelength.

Figure 11: Examples of contour plots that may generate unreliable (red) and
reliable (blue) results. Axes cover the full parameter space allowed from the
models. Contours enclose 1 and 2$\sigma$ confidence regions.

Overall, we find excellent agreement with existing data, with a typical
$\eta(I)$ of $\sim 5\times 10^{-5}$ ($\sim 10^{-4}$ $1\sigma$ scatter).
Differences are lower for higher S/N indices $\mathrm{G4300}$ ,
$\mathrm{Fe4383}$ , $\mathrm{Fe4531}$ , and $\mathrm{C_{2}4668}$
($|\eta(I)|<10^{-5}$), and outliers are mostly galaxies with lower-quality
spectra. There is no trend between relative differences and the indices
values, and the distribution is qualitatively Gaussian.

We find good agreement also for the uncertainties, with a typical
$\eta(\sigma)$ of $\sim 10^{-2}$ ($\sim 10^{-1}$ $1\sigma$ scatter). Larger
deviations are seen for indices with values $\sim 0$ (Balmer indices,
$\mathrm{CN_{1}}$ ), for which measured uncertainties are $\sim 10\%$ lower.
This difference can be due to different methods to estimate formal errors.

We note that at the current resolution ($R\sim 3500$), the method used to
interpolate the spectra can introduce higher discrepancies than those observed
before. In particular, using $0^{th}$ order interpolation the typical scatter
in $\eta(I)$ increases up to $\sim 10^{-2}$. Overall, these results confirm
the reliability of PyLick to measure indices values and formal errors from
observed spectra.

##  Appendix B Assessing convergence and reliability of MCMC posterior
distributions

In a Bayesian analysis, it is crucial to determine whether MCMC chains are
reproducing with sufficient accuracy the target posterior distribution.
However, there is not an established standard to assess the convergence (Hogg
& Foreman-Mackey, 2018; Roy, 2020). A possibility is to take into account the
autocorrelation time $\tau_{\mathrm{int}}$ analysis, where
$\tau_{\mathrm{int}}$ quantifies how many steps are needed to generate
independent samples. We consider a chain to be formally converged when
$\tau_{\mathrm{int}}$ for each parameter is greater than 1/100th of the chain
size. In this work, analyzed galaxies typically require $\sim 7000$ steps.
However, sample averages derived from formally converged chains could still be
unreliable. This is the case when posterior distributions:

  1. 1.

are skewed toward the priors (in this work, the limits of the parameter space
allowed from the models);

  2. 2.

are not predictive (mainly due to high age-metallicity degeneracy);

  3. 3.

are multimodal.

In Figure 11 we show contour plots representative for the three categories and
a typical contour plot for a good fit. We note that the age-metallicity
degeneracy follows the so-called “3/2 rule” (Worthey, 1998), where an increase
(decrease) of log age by a factor of 2, when accompanied by a decrease
(increase) of stellar metallicity $\rm[Z/H]$ by a factor of 3, can reproduce
the same set of observed indices. Multimodal distributions are obtained for
less than 4% of the analyzed galaxies and follow the direction of the age-
metallicity degeneracy. As they do not recur for the same galaxy when slightly
different sets of indices are considered, they should be considered as a
special case of (2) and to be due to intrinsic degeneracies rather than a real
complexity in the stellar population. We check joint and marginal
distributions of all the 199 analyzed galaxies, and we flag and exclude
galaxies belonging to these three categories.

In Figure 12 we show the distributions of the main parameters before and after
this process. The great majority ($\sim 85\%$) of excluded galaxies have
overall spectral $\left\langle\right.\\!\\!{\rm
S/N}\\!\\!\left.\right\rangle<25$. However, a cut of this kind applied a
priori would have halved the final sample, excluding also good converged fit.
Low-$\sigma_{\star}$ galaxies are preferentially excluded because they have
relatively lower S/Ns. Distributions of stellar population parameters do not
change significantly before and after this process. The posterior distribution
of ages for galaxies with $\mathrm{age}>8$ Gyr tend to be heavily skewed
toward the 15 Gyr prior and are also characterized by low metallicities
($\hbox{$\rm[Z/H]$}<-0.4$) and relatively high $\alpha$/Fe values.

Figure 12: Distribution of overall spectral signal-to-noise ratio
$\left\langle\right.\\!\\!{\rm S/N}\\!\\!\left.\right\rangle_{spec}$, stellar
velocity dispersion $\sigma_{\star}$, stellar population parameters
($\mathrm{age}$, $\rm[Z/H]$, $\rm[\alpha/Fe]$), and signal-to-noise ratios of
individual indices for 199 analyzed passive galaxies in LEGA-C DR2 (gray
histograms). Violet (black hatched) histograms represent galaxies included
(excluded) after inspecting the MCMC posterior distributions.

We observe that S/Ns of individual indices are a good indicator to separate
included and excluded galaxies. For this purpose, we use indices with
relatively higher S/N (see Table 2), namely, $\mathrm{G4300}$ ,
$\mathrm{Fe4383}$ , $\mathrm{Fe4531}$ , and $\mathrm{C_{2}4668}$ . The great
majority ($\sim 90\%$) of excluded galaxies have at least one index with ${\rm
S/N<10}$. A cut of this kind would still exclude $\sim 29\%$ of galaxies with
reliable constraints. Anyway, we note that these cuts are not universal
because they vary for different index combinations, and individual S/Ns of
spectral indices vary for different spectral resolutions.

In conclusion, we find that the inspection of MCMC posterior distributions is
an efficient procedure to detect unreliable constraints–mostly due to low
indices S/N–while maximizing the number of galaxies for which we obtain
reliable constraints.

##  Appendix C Different sets of indices

In this section we illustrate the results obtained with different index sets.
We explore the huge number of viable combinations ($\sim 1$ million) following
three approaches:

  * •

Maximize the information to be fitted. Using a higher number of indices should
provide more stable results. However, indices should be calibrated, and able
to disentangle degeneracies giving equal weight to each model parameter.

  * •

Use already-proposed index combinations. In particular we include redder Mg
and Fe indices as done by previous works.

  * •

Use a small, essential set to break the existing degeneracy between
parameters. Spectral indices are sensitive to variations in $\mathrm{age}$,
$\rm[Z/H]$, $\rm[\alpha/Fe]$, but the relative sensitivity to these parameters
is not identical. One can choose a small combination of four to six indices
based upon their different sensitivity (e.g., for element abundances see
Tripicco & Bell, 1995; Korn et al., 2005; Lee et al., 2009).

As discussed in the text, we find an optimal combination that maximizes the
spectral coverage, the number of constrained galaxies, and the precision of
the constraints: $\mathrm{H\delta_{A}}$ , $\mathrm{CN_{1}}$ ,
$\mathrm{CN_{2}}$ , $\mathrm{Ca4227}$ , $\mathrm{G4300}$ ,
$\mathrm{H\gamma_{A}}$ , $\mathrm{H\gamma_{F}}$ , $\mathrm{Fe4383}$ ,
$\mathrm{Fe4531}$ , and $\mathrm{C_{2}4668}$ (hereafter baseline). In Table 6,
we report the other main sets of indices that we analyzed, and the differences
in the derived parameters with respect to the baseline combination. We also
analyzed many other sets ($\sim 50$), but they do not add significant
information to this study.

Table 6: Examples of Analyzed Index Combinations and Their Definitions.

Combo ID |  $\mathrm{H\delta_{A}}$ |  $\mathrm{H\delta_{F}}$ |  $\mathrm{CN_{1}}$ |  $\mathrm{CN_{2}}$ |  $\mathrm{Ca4227}$ |  $\mathrm{G4300}$ |  $\mathrm{H\gamma_{A}}$ |  $\mathrm{H\gamma_{F}}$ |  $\mathrm{Fe4383}$ |  $\mathrm{Ca4455}$ |  $\mathrm{Fe4531}$ |  $\mathrm{C_{2}4668}$ |  $\mathrm{H\beta}$ |  $\mathrm{Fe5015}$ |  $\mathrm{Mg_{1}}$ |  $\mathrm{Mg_{2}}$ |  $\mathrm{Mg}\,b$ |  $\mathrm{Fe5270}$ |  $\mathrm{Fe5335}$ |  $\mathrm{Fe5406}$ | Combo ID | N (in common) | $\Delta\mathrm{age}~{}(\sigma)$ | $\Delta\hbox{$\rm[Z/H]$}~{}(\sigma)$ | $\Delta\hbox{$\rm[\alpha/Fe]$}~{}(\sigma)$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
baseline | $\blacksquare$ |  | $\blacksquare$ | $\blacksquare$ | $\blacksquare$ | $\blacksquare$ | $\blacksquare$ | $\blacksquare$ | $\blacksquare$ |  | $\blacksquare$ | $\blacksquare$ |  |  |  |  |  |  |  |  |  | 140 (140) | – | – | –
1 | $\blacksquare$ |  |  |  | $\blacksquare$ | $\blacksquare$ | $\blacksquare$ | $\blacksquare$ | $\blacksquare$ |  | $\blacksquare$ | $\blacksquare$ |  |  |  |  |  |  |  |  |  | 105 (95) | 0.04 (0.06) | -0.04 (0.42) | -0.05 (0.64)
2 | $\blacksquare$ |  |  | $\blacksquare$ | $\blacksquare$ | $\blacksquare$ | $\blacksquare$ |  | $\blacksquare$ |  | $\blacksquare$ | $\blacksquare$ |  |  |  |  |  |  |  |  |  | 131 (115) | -0.42 (0.63) | 0.03 (0.33) | -0.06 (0.76)
3 | $\blacksquare$ |  | $\blacksquare$ | $\blacksquare$ | $\blacksquare$ | $\blacksquare$ | $\blacksquare$ | $\blacksquare$ | $\blacksquare$ |  | $\blacksquare$ | $\blacksquare$ |  |  |  |  | $\blacksquare$ |  |  |  |  | 39 (39) | 0.05 (0.08) | 0.00 (0.02) | 0.08 (0.92)
4 | $\blacksquare$ |  |  |  | $\blacksquare$ | $\blacksquare$ | $\blacksquare$ | $\blacksquare$ | $\blacksquare$ |  | $\blacksquare$ | $\blacksquare$ |  |  | $\blacksquare$ | $\blacksquare$ | $\blacksquare$ | $\blacksquare$ | $\blacksquare$ | $\blacksquare$ |  | 11 (8) | 0.25 (0.36) | -0.01 (0.09) | 0.09 (1.15)
5 |  |  |  | $\blacksquare$ | $\blacksquare$ | $\blacksquare$ |  | $\blacksquare$ | $\blacksquare$ |  |  |  |  |  |  |  |  |  |  |  |  | 133 (98) | -0.17 (0.25) | 0.01 (0.08) | -0.13 (1.61)
6 |  |  |  | $\blacksquare$ | $\blacksquare$ | $\blacksquare$ |  |  | $\blacksquare$ | $\blacksquare$ |  | $\blacksquare$ |  |  |  |  |  |  |  |  |  | 119 (102) | -1.11 (1.54) | 0.11 (0.95) | -0.06 (0.69)

Note. — For each combination we report the number of constrained galaxies,
along with simple and $1\sigma$ differences in the derived parameters with
respect to the baseline combination.

We find no significant systematic differences in the derived parameters when
small changes to the baseline combination are applied, i.e. by adding or
removing one to two indices. In particular, we focus here on the removal of CN
indices, as the nitrogen abundance is not a free parameter in TMJ11 models
(Combo 1), and on the removal of those indices that sample twice the same
spectral region (Combo 2). In the first case, we obtain constraints for fewer
galaxies with respect to the baseline combination, but with an overall
excellent agreement. In the second case, we constrain about the same number of
galaxies, obtaining lower ages but still in agreement with the baseline set.
This is likely due to the removal of $\mathrm{H\gamma_{F}}$ that reduces the
weight of age-sensitive features, therefore producing less reliable age
estimates.

Figure 13: Differences between the parameters obtained with and without
$\mathrm{Mg}\,b$ (Combo 3 vs. baseline). Median (dashed line) and $\pm$NMAD
values (shaded region) are annotated in the upper corners.

Index combinations discussed above lack redder indices, such as the Mg
indices, traditionally used as $\alpha$-abundance indicators. Therefore, we
repeat the analysis including $\mathrm{Mg}\,b$ (Combo 3). Amongst the 59 with
such relatively large spectral coverage, we obtain constraints for 39
galaxies. Differences are shown in greater detail in Figure 13. While ages and
$\rm[Z/H]$ are in overall good agreement, the inclusion of another
$\alpha$-sensitive index suggests that $\rm[\alpha/Fe]$ values derived with
Combo 1 may be underestimated by $\sim 0.1$ dex. However, given the small
statistical significance, we do not correct for this offset. A similar
discussion can be done by extending the analysis at all the redder indices
(Combo 4, which is also the same combination used in Onodera et al. 2015).

Finally, we note here that although minimal sets of $N=4$ indices in the
wavelength region between 4000 and 4600 Å allow us to analyze a large number
of galaxies ($\sim 300$), we do not find a relevant set to place constraints
on more than one-third of them. This situation is improved when $N=5-7$ (e.g.,
Combo 5), but results show an overall stronger age–-metallicity degeneracy
with respect to $N\geq 8$. It is also interesting to note that even if we do
not include Balmer indices (Combo 6), we still obtain constraints for more
than 100 galaxies. But we stress again that results are less reliable because
we removed most of the age-sensitive features.

In summary, we find that a blind choice of index combinations can lead to
less-robust results. This happens when a combination is unbalanced toward one
or more parameters of the fit, but also if indices are measured on spectra
where the sky subtraction was imperfect. After an extensive analysis, we
demonstrate that within the limited statistics and wavelength coverage of
current data, results do not show significant systematic differences.