# The stellar metallicities of massive quiescent galaxies at
$\mathbf{1.0<z<1.3}$ from KMOS+VANDELS

A. C. Carnall1, R. J. McLure1, J. S. Dunlop1, M. Hamadouche1, F. Cullen1, D.
J. McLeod1, R. Begley1, R. Amorin2,3, M. Bolzonella4, M. Castellano5, A.
Cimatti6,7, F. Fontanot8,9, A. Gargiulo10, B. Garilli10, F. Mannucci7, L.
Pentericci5, M. Talia4,6, G. Zamorani4, A. Calabro5, G. Cresci7, N. P. Hathi11

1 SUPA, Institute for Astronomy, University of Edinburgh, Royal Observatory,
Edinburgh EH9 3HJ, UK
2 Instituto de Investigación Multidisciplinar en Ciencia y Tecnología,
Universidad de La Serena, Raúl Bitrán 1305, La Serena, Chile
3 Departamento de Física y Astronomía, Universidad de La Serena, Av. Juan
Cisternas 1200 Norte, La Serena, Chile
4 INAF - OAS Bologna, Via P. Gobetti 93/3, I-40129, Bologna, Italy
5 INAF - Osservatorio Astronomico di Roma, Via Frascati 33, I-00078
Monteporzio Catone, Italy
6 University of Bologna, Department of Physics and Astronomy (DIFA), Via
Gobetti 93/2, I-40129, Bologna, Italy
7 INAF - Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125,
Firenze, Italy
8 INAF - Astronomical Observatory of Trieste, via G.B. Tiepolo 11, I-34143
Trieste, Italy
9 IFPU - Institute for Fundamental Physics of the Universe, via Beirut 2,
34151, Trieste, Italy
10 INAF - IASF Milano, Via A. Corti 12, I-20133, Milano, Italy
11 Space Telescope Science Institute, 3700 San Martin Dr., Baltimore, MD
21218, USA E-mail<EMAIL_ADDRESS>Universities Physics Alliance

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present a rest-frame UV-optical ($\lambda=2500{-}6400$ Å) stacked spectrum
representative of massive quiescent galaxies at $1.0<z<1.3$ with
log$(M_{*}/\mathrm{M_{\odot}})>10.8$. The stack is constructed using VANDELS
survey data, combined with new KMOS observations. We apply two independent
full-spectral-fitting approaches, obtaining consistent stellar ages and
metallicities. We measure a total metallicity, [Z/H] = $-0.13\pm 0.08$, and an
iron abundance, [Fe/H] = $-0.18\pm 0.08$, representing falls of $\sim 0.3$ dex
and $\sim 0.15$ dex respectively compared with the local Universe. We also
measure the alpha enhancement via the magnesium abundance, obtaining [Mg/Fe] =
0.23$\pm$0.12, consistent with similar-mass galaxies in the local Universe,
indicating no evolution in the average alpha enhancement of
log$(M_{*}/\mathrm{M_{\odot}})\sim 11$ quiescent galaxies over the last $\sim
8$ Gyr. This suggests the very high alpha enhancements recently reported for
several very bright $z\sim 1-2$ quiescent galaxies are due to their extreme
masses, log$(M_{*}/\mathrm{M_{\odot}})\gtrsim 11.5$, in accordance with the
well-known downsizing trend, rather than being typical of the $z\gtrsim 1$
population. The metallicity evolution we observe with redshift (falling [Z/H],
[Fe/H], but constant [Mg/Fe]) is consistent with other recent studies. We
recover a mean stellar age of $2.5^{+0.6}_{-0.4}$ Gyr, corresponding to a
formation redshift, $z_{\mathrm{form}}=2.4^{+0.6}_{-0.3}$. Recent studies have
obtained varying average formation redshifts for $z\gtrsim 1$ massive
quiescent galaxies, and, as these studies report consistent metallicities, we
identify different star-formation-history models as the most likely cause.
Larger spectroscopic samples from upcoming ground-based instruments will
provide precise constraints on ages and metallicities at $z\gtrsim 1$.
Combining these results with precise $z>2$ quiescent-galaxy stellar-mass
functions from James Webb Space Telescope data will provide an independent
test of formation redshifts derived from spectral fitting.

###### keywords:

keyword1 – keyword2 – keyword3

††pubyear: 2021††pagerange: The stellar metallicities of massive quiescent
galaxies at $\mathbf{1.0<z<1.3}$ from KMOS+VANDELS–References

##  1 Introduction

In the local, present-day Universe, the massive galaxy population
($\mathrm{log}_{10}(M_{*}\ /\ \mathrm{M_{\odot}})\gtrsim 10.5$) is dominated
by quiescent galaxies, which have shut down (quenched) their star-formation
activity (e.g. McLeod et al. 2021). The formation and quenching processes
leading to the rise of this dominant population, across at least the last 12
Gyr since the first known quiescent galaxies at redshift, $z\simeq 4$, are
therefore of central importance to our understanding of galaxy evolution.

Studying the massive quiescent galaxy population in the local Universe
presents several key challenges. Firstly, as the rate of change in stellar
population spectra is roughly logarithmic with age (e.g. Ocvirk et al. 2006),
constraints on ages via spectral fitting become steadily less precise for
older stellar populations. This means that formation redshifts for local
quiescent galaxies, with stellar population ages $\gtrsim 10$ Gyr, are highly
uncertain. In addition, massive galaxies gradually accrete new stellar
populations via merger events, which have the potential to change their
physical sizes, as well as the average ages and metallicities of their stellar
populations. This obscures the signatures of the dominant physical processes
that acted on these galaxies during their main epoch of formation.

To understand the rise of massive quiescent galaxies, it is therefore
necessary to conduct detailed observational studies across the whole history
of the Universe, from the first $1{-}2$ Gyr to the present day. By measuring
how the distributions of key physical parameters evolve as a function of
redshift, it should be possible to disentangle the degenerate effects of
different processes acting at different times, providing strong constraints on
the key physical ingredients required to produce the local massive galaxy
population.

As a practical first step towards achieving this goal, much attention is
currently focused on measuring the number densities of quiescent galaxies, as
well as the distributions of their physical sizes, mean stellar ages and
metallicities, as a function of stellar mass and redshift. Armed with this
information across a large fraction of cosmic history, we may then aspire to
build a self-consistent model for the assembly of the quiescent population,
including the physical processes that influence massive galaxies during their
main epoch of formation and quenching, the subsequent growth of individual
quiescent galaxies through merger events, and the growth of the quiescent
population as a whole via new galaxies quenching their star-formation
activity.

The largest single factor currently limiting progress in measuring the
evolution of these four key properties (number density, size, stellar age and
metallicity) across cosmic time is the availability of high-quality
observational data for faint, red, massive quiescent galaxies in the high-
redshift Universe. The best-constrained parameters for the quiescent
population currently are number densities and physical sizes. These are now
widely studied out to $z\sim 3$, as they can be reliably constrained via high-
quality photometric data (e.g. McLure et al. 2013; van der Wel et al. 2014;
Straatman et al. 2014, 2016; Cecchi et al. 2019; Mowla et al. 2019a; Mowla et
al. 2019b; Suess et al. 2019a, b; Girelli et al. 2019; Merlin et al. 2019;
Sherman et al. 2020; Carnall et al. 2020; Marsan et al. 2020; Santini et al.
2021).

Stellar ages (more generally star-formation histories; SFHs) have also been
studied using photometric data (e.g. Pacifici et al. 2016; Carnall et al.
2018). However, the age-metallicity-dust degeneracy in galaxy spectral shapes
results in relatively weak constraints, meaning the applied priors
significantly impact the results obtained (Carnall et al. 2019b; Leja et al.
2019a). Recently, the increasing availability of medium to high signal-to-
noise ratio (SNR) spectroscopic data for quiescent galaxies at intermediate
redshifts (e.g. van der Wel et al. 2016; McLure et al. 2018 \- hereafter M18),
combined with sophisticated full-spectral-fitting approaches, have produced
the first strong constraints on the SFHs of representative samples out to
$z\gtrsim 1$ (e.g. Wu et al. 2018a, b; Belli et al. 2019; Carnall et al.
2019a; Estrada-Carpenter et al. 2019, 2020; Wild et al. 2020; Tacchella et al.
2021).

However, systematic differences still exist in the results of these studies,
likely due to a combination of different assumed SFH models, and the fact that
stellar metallicities (which are strongly degenerate with ages) are relatively
weakly constrained by intermediate-SNR data, meaning the applied priors still
play a significant role.

Measuring the stellar metallicities of quiescent galaxies represents another
step-change in observational difficulty, with strong constraints ($\sim 0.1$
dex) only available at continuum SNR $\gtrsim 10{-}15\,$Å-1 in the rest-frame
optical (e.g. Gallazzi et al. 2005; Pacifici et al. 2012). This is compounded
at $z\gtrsim 1$ by the key features being shifted into the near-IR, where
stronger atmospheric absorption and emission make continuum observations from
the ground far more challenging.

Despite these challenges, stellar metallicity measurements are highly
valuable, as they are strongly constraining on galaxy-formation models, being
intimately linked to the physics of star formation and gas recycling (e.g.
Maiolino & Mannucci 2019). Whilst substantial progress has been made out to
$z\sim 3$ with the brighter, bluer continua of star-forming galaxies (e.g.
Steidel et al. 2016; Cullen et al. 2019, 2021), to date, studies of quiescent-
galaxy stellar metallicities at $z>1$ are rare, and often restricted to either
individual bright objects, or low-resolution grism spectra (e.g. Whitaker et
al. 2013; Onodera et al. 2015; Lonoce et al. 2015; Kriek et al. 2016).

Whilst the upcoming James Webb Space Telescope (JWST) will produce
exceptionally high-quality data for limited numbers of objects at $z>1$, truly
statistical studies are still several years away, awaiting the advent of the
next generation of ground-based multi-object spectrographs (e.g. Cirasuolo et
al. 2020).

In this work, we present the first determination of the stellar metallicities
of a mass-selected sample of quiescent galaxies at $z>1$. We combine rest-
frame near-UV data from the VANDELS survey (M18; Pentericci et al. 2018;
Garilli et al. 2021) with new rest-frame optical KMOS $YJ-$band spectroscopy,
to produce a stacked spectrum covering rest-frame $2500{-}6400$ Å for UVJ-
selected galaxies at $1.0<z<1.3$ with
$\mathrm{log_{10}(}M_{*}/\mathrm{M_{\odot})}>10.8$. We fit our stack using
both Bagpipes (Carnall et al., 2018) and Alf (Conroy & van Dokkum, 2012;
Conroy et al., 2018), obtaining consistent ages and metallicities.

At $1.0<z<1.3$, the VANDELS spectra span rest-frame wavelengths, $\lambda\sim
2500{-}4500$ Å. This means they include the Balmer/4000 Å break region
critical for precise age determination. However, generally, they do not
include the key Fe and Mg absorption features at $4500{-}5500$ Å most commonly
used to measure stellar metallicities in the local Universe (e.g. Gallazzi et
al. 2005).

Stellar metallicity was fitted as a free parameter in the VANDELS full-
spectral-fitting analysis of Carnall et al. (2019a). However, due to the lack
of these key Fe and Mg features, combined with the fact that empirical
stellar-population models, generally still accepted as more accurate than
those based on theoretical spectra (e.g. Coelho et al. 2020), are only
available at $\gtrsim 3500$ Å, we chose not to report our derived stellar
metallicities in that work. Instead, we obtained further observations with
KMOS to gain access to these key rest-frame optical features, to ensure
reliable results, and to facilitate direct comparisons with previous studies
in the local Universe.

The structure of this paper is as follows. In Sections 2 and 3, we introduce
our VANDELS and KMOS spectroscopic datasets respectively. In Section 4, we
describe the selection of our mass-complete sample, as well as the process of
constructing and fitting our representative stacked spectrum with Bagpipes and
Alf. Our results are presented in Section 5, and discussed in Section 6. We
present our conclusions in Section 7. All magnitudes are quoted in the AB
system. For cosmological calculations, we adopt $\Omega_{M}=0.3$,
$\Omega_{\Lambda}=0.7$ and $H_{0}$ = 70 $\mathrm{km\ s^{-1}\ Mpc^{-1}}$. All
times, $t$, are measured forwards from the beginning of the Universe. We
assume a Kroupa (2001) initial mass function. We assume the Solar abundances
of Asplund et al. (2009), such that $\mathrm{Z_{\odot}}=0.0142$.

##  2 VANDELS data and sample selection

VANDELS (M18; Pentericci et al. 2018; Garilli et al. 2021) is a large ESO
public spectroscopic survey of the high-redshift Universe, using the VIsible
Multi Object Spectrograph (VIMOS) instrument on the Very Large Telescope
(VLT). The primary aim of VANDELS is to detect continuum emission at high SNR
for high-redshift galaxies, moving beyond redshift acquisition to study galaxy
physical properties within the first 5 billion years prior to $z=1$.

###  2.1 The parent photometric sample

The parent photometric sample for this study consists of 812 massive quiescent
galaxies, selected by M18 as potential targets for the VANDELS survey. These
objects were selected from four photometric catalogues: the CANDELS GOODS
South and UDS catalogues of Guo et al. (2013) and Galametz et al. (2013), and
two further ground-based photometric catalogues, purpose-built for VANDELS.
These cover the regions immediately surrounding the CANDELS footprints. This
approach was necessary as the VIMOS spectrograph field of view covers a larger
area than the CANDELS HST imaging.

The parent sample was selected from these four photometric catalogues by the
following process (described in full detail in section 4 of M18). Objects were
first required to meet the following apparent magnitude and photometric
redshift, $z_{\mathrm{phot}}$, criteria

  * •

$H\leq 22.5$

  * •

$i\leq 25$

  * •

$1.0\leq z_{\mathrm{phot}}\leq 2.5$.

For the CANDELS catalogues, the photometric redshifts used were those
published by the CANDELS team (Dahlen et al., 2013). For the ground-based
catalogues, photometric redshifts were generated through a similar process by
the VANDELS team, taking the median of results obtained using a variety of
public codes. The $H$-band magnitude cut limits the sample to objects with
$\mathrm{log_{10}(}M_{*}/\mathrm{M_{\odot})}\gtrsim 10$. The $i$-band cut was
implemented to make sure the faintest objects would be detected in the VANDELS
spectra, and is not relevant to this study, which is focused on a brighter,
mass-selected sub-sample.

To select only quiescent galaxies, the following, permissive, rest-frame UVJ
magnitude selection criteria were then applied:

  * •

$U-V>0.88(V-J)+0.49$

  * •

$V-J<1.6$

  * •

$U-V>1.2$

From the 812 galaxies selected, 64 per cent have $z_{\mathrm{phot}}<1.3$, with
the sample being mass complete down to
$\mathrm{log_{10}(}M_{*}/\mathrm{M_{\odot})}=10.3$ at $1.0<z<1.3$ (Carnall et
al., 2019a). These 812 objects are referred to as the parent sample throughout
the rest of this work.

###  2.2 VANDELS spectroscopic observations

From the parent sample described in the previous sub-section, objects were
assigned to slits at random to be observed as part of VANDELS. Spectra were
obtained for 281 massive quiescent galaxies, roughly one third of the parent
sample. These represent roughly 13 per cent of the VANDELS survey, which
includes $\sim 2000$ objects in total. The remaining 87 per cent is composed
of star-forming galaxies at $z>2.4$. Objects were observed for 20, 40 or 80
hours, depending on their $i$-band magnitudes, with a mean integration of
$\sim 45$ hours.

The VANDELS spectroscopic data (described in full in Pentericci et al. 2018)
covers observed-frame wavelengths from $4800{-}10300$ Å, with spectral
resolving power, $R=\lambda/\Delta\lambda\sim 600$. The VANDELS data have high
SNR for these redshifts. The spectra of quiescent galaxies within the mass-
complete section of the parent sample discussed in Section 2.1 (those with
$1.0<z<1.3$ and $\mathrm{log_{10}(}M_{*}/\mathrm{M_{\odot})}>10.3$) have a
median continuum SNR $\sim$ 11 Å-1 at 7500 Å.

Spectroscopic redshifts were measured by the VANDELS team using the pandora.ez
software (Garilli et al., 2010). Spectra were assigned redshift quality flags
following Le Fèvre et al. (2013), with all but 12 of the 281 observed massive
quiescent galaxies assigned quality flags 3 or 4, corresponding to 95 and 100
per cent probabilities of correct identification respectively (in fact,
Garilli et al. 2021 demonstrate $\sim 99$ per cent reliability for the VANDELS
flag 3 and 4 redshifts combined). The 12 spectra with lower quality flags, as
well as one low-redshift interloper, were excluded from this work, leaving 268
spectra. Three objects with spectroscopic redshifts in the range $0.96<z<1.0$
were retained. All VANDELS spectroscopic data used in this work comes from the
final public data release, DR4, described in Garilli et al. (2021).

##  3 KMOS Data and Sample Selection

We observed four KMOS pointings during ESO P104, from October 2019 $-$ January
2020, under programme ID 0104.B$-$0885(A). Two pointings targeted each of the
two VANDELS fields (UDS and GOODS South). Observations were in the $YJ-$band,
providing wavelength coverage from $\sim 1-1.35\,\mathrm{\mu m}$ at $R\sim
3600$. During the period in which our data were taken, 23 KMOS arms were
functional, meaning we obtained spectra for a total of 92 objects. Each
pointing was observed for 8 hours on source, with 300-second exposures, and an
ABAB nodding pattern between object and sky positions. The mean seeing for our
observations was $0.6^{\prime\prime}$.

###  3.1 KMOS sample selection

Targets for our KMOS observations were drawn almost exclusively from the
parent sample of 812 objects described in Section 2.1. To maximise the utility
of our $YJ-$band observations, objects were prioritised by their $J-$band
magnitude. Objects were assigned to KMOS IFUs using the KMOS ArM Allocator
(KARMA) tool, with priority classes being defined as follows

  1. 1.

Parent sample with $J<21.5$

  2. 2.

Parent sample with $J>21.5$

  3. 3.

Post-starburst galaxies from Wilkinson et al. (2021).

The priority 3 targets, obtained via private communication, are unrelated to
this study, and were added only as fillers to make sure all IFUs were
assigned. The positions of the four KMOS pointings were optimised to target
the maximum possible number of priority class 1 objects, as well as to
maximise the overlap with VANDELS spectroscopy. In total, we allocated IFUs to
62 objects with priority 1, 29 objects with priority 2, and 1 object with
priority 3.

From the parent sample of 812 objects, a total of 273 have $J<21.5$, meaning
objects of priority class 1 were $\sim 4$ times more likely to be observed
than those of priority class 2. From the 91 parent sample objects observed
with KMOS, 51 objects were also observed by VANDELS, as described in Section
2.

###  3.2 KMOS data reduction

Our KMOS data were reduced using a combination of the standard Esorex pipeline
recipes and custom code, optimised for these data. The KMOS pipeline v4.0.0
was used to produce flux-calibrated, sky-subtracted cubes for each individual
pair of 300 second (object, sky) exposures. All our targets with $J\lesssim
22$ (approximately 75 out of 92 objects, including all those used in this
paper) are clearly detected in each single-exposure, wavelength-collapsed
cube.

Upon inspection of the wavelength-collapsed cubes, $1{-}2$ spaxel
($0.2{-}0.4^{\prime\prime}$) shifts in object centroids were noted with
respect to the WCS coordinates of each object. New centroid positions were
therefore measured for each frame by selecting the brightest pixel within a
1${}^{\prime\prime}-$diameter circular aperture centred on the WCS position.
These new centroids were verified by manual inspection of each cube.

Also noted were substantial sky-line residuals still present in the data
cubes, along with varying systematic shifts in pixel values away from zero
across the whole wavelength axis. To address these issues, we implement a
further, custom, sky-subtraction step. We first mask all pixels within a
1${}^{\prime\prime}-$diameter circular aperture centred on the new object
centroid pixel, as well as all pixels bordering the edge of the detector. We
then subtract the median of the remaining pixels from the cube slice at each
wavelength. This approach is only possible in this specific instance, as all
our targets, being $z>1$ quiescent galaxies, are extremely compact, with
effective radii, $r_{\mathrm{e}}<<1^{\prime\prime}$ (e.g. McLure et al. 2013;
van der Wel et al. 2014).

Exposures were then aligned according to their updated centroid positions and
median stacked, with uncertainties calculated via the robust median absolute
deviation (MAD) indicator. Finally, 1D spectra were extracted within
1${}^{\prime\prime}-$diameter circular apertures, using the Horne (1986)
optimal extraction algorithm.

Both the custom re-centroiding and additional sky-subtraction steps are
critical to recovering the expected Fe, Mg and Na absorption features in the
spectra of individual objects. Using the default, pipeline-combined cubes only
produces visible continuum features in the spectrum of our brightest target,
with $J=19.6$. However, with these additional steps, the strongest features
(Mg i 5170 Å, and Na i 5895 Å) are visible in almost all spectra for objects
with $J<21.5$.

###  3.3 KMOS redshift measurement

As described above, 51 out of 91 objects from the parent sample observed with
KMOS already have secure spectroscopic redshifts from VANDELS. For the
remaining 40 objects, we followed the same process as described in Section
2.2, using Pandora.ez, to measure redshifts. In Table 1, we report 25 new
spectroscopic redshifts, $z_{\mathrm{KMOS}}$, measured from our KMOS data,
along with their associated quality flags. The remaining 15 objects are all
significantly fainter ($J>21.8$), meaning that no reliable redshifts could be
measured.

Figure 1: The distribution of the sample in $J-$band magnitude versus redshift. Circles denote the parent sample from M18 (see Section 2.1). Objects included in our mass-complete photometric sample, selected as described in Section 4.1, are marked with squares. The colours denote the availability of spectroscopic data from VANDELS and KMOS for each object. The dashed horizontal line is the threshold for inclusion as a priority 1 target in our KMOS observations (see Section 3.1). Objects below were assigned priority 2. Table 1: Spectroscopic redshifts measured from our KMOS data (see Section 3.3). From the 91 parent sample objects observed with KMOS, 51 already have robust VANDELS spectroscopic redshifts. The 25 objects in this table are those of the remaining 40 for which redshifts could be measured. The remaining 15 objects are all fainter than $J=21.8$.  ID | RA | DEC | $J$ | $z_{\mathrm{KMOS}}$ | Flag
---|---|---|---|---|---
CDFS-017418 | 53.15497 |  $-$27.76891 | 19.86 | 1.0956 | 4
UDS-196414 | 34.48731 |  0$-$5.09687 | 20.23 | 1.0921 | 4
UDS-205452 | 34.50079 |  0$-$5.05545 | 20.59 | 1.0945 | 4
UDS-190420 | 34.51167 |  0$-$5.12379 | 20.60 | 1.0945 | 3
UDS-200825 | 34.52746 |  0$-$5.07677 | 20.74 | 1.0369 | 3
UDS-201280 | 34.51318 |  0$-$5.07575 | 20.78 | 1.4140 | 2
CDFS-004529 | 53.08042 |  $-$27.87204 | 20.78 | 1.0967 | 3
UDS-196179 | 34.50759 |  0$-$5.09889 | 20.91 | 1.2710 | 3
CDFS-016336 | 53.22897 |  $-$27.77253 | 20.95 | 1.0386 | 3
CDFS-020067 | 53.15878 |  $-$27.74239 | 20.97 | 1.2210 | 3
CDFS-022694 | 53.14845 |  $-$27.71946 | 20.99 | 1.2221 | 3
CDFS-003952 | 53.07276 |  $-$27.87632 | 21.03 | 1.1014 | 4
CDFS-004376 | 53.07153 |  $-$27.87246 | 21.05 | 1.0976 | 3
UDS-197616 | 34.43864 |  0$-$5.09219 | 21.07 | 1.6523 | 4
UDS-015126 | 34.31596 |  0$-$5.19366 | 21.09 | 1.3611 | 2
UDS-207822 | 34.48340 |  0$-$5.04458 | 21.18 | 1.0334 | 3
UDS-024934 | 34.52885 |  0$-$5.12719 | 21.22 | 1.0954 | 3
UDS-009642 | 34.33053 |  0$-$5.22374 | 21.26 | 1.0735 | 3
UDS-013785 | 34.32431 |  0$-$5.20136 | 21.27 | 1.0915 | 4
UDS-013519 | 34.33534 |  0$-$5.20168 | 21.40 | 1.5321 | 3
UDS-192952 | 34.48285 |  0$-$5.11423 | 21.43 | 1.3295 | 3
UDS-020224 | 34.34611 |  0$-$5.16692 | 21.44 | 1.0830 | 2
UDS-005970 | 34.31475 |  0$-$5.24326 | 21.45 | 1.7127 | 2
UDS-010643 | 34.38065 |  0$-$5.21789 | 21.46 | 1.2624 | 3
CDFS-014839 | 53.16516 |  $-$27.78587 | 21.53 | 1.3175 | 2
Figure 2: The stacked rest-frame UV-optical spectrum for quiescent galaxies at
$1.0<z<1.3$ with $\mathrm{log_{10}(}M_{*}/\mathrm{M_{\odot})}>10.8$. The
stacked spectrum is shown in blue in the top panel, with the error spectrum
shown in green. The SNR per Å and number of objects contributing to the stack
are shown in the middle and lower panels respectively. Key age and metallicity
sensitive features are marked in black. The red line overlaid on the spectrum
shows the posterior median model fitted to the stack with Bagpipes, as
described in Section 4.3. This includes a calibration polynomial, which makes
corrections of $\sim 5$ per cent. We derive a stellar metallicity of [Z/H] =
$-0.13\pm 0.08$ and a formation redshift of
$z_{\mathrm{form}}=2.4^{+0.6}_{-0.3}$. The gray vertical bands show regions
masked from fitting (see Section 4).

##  4 Stacking analysis

The KMOS spectra we obtained for objects with $J<21.5$ have a median SNR $\sim
4$ Å-1 at 12000 Å ($\sim 14$ per resolution element at $R=1000$). This is
unfortunately not sufficient to constrain the stellar metallicities of
individual objects. In this section, we therefore define a mass-complete
sample, for which we can construct a representative stacked spectrum from the
VANDELS + KMOS data described in Sections 2 and 3. We fit these data to
constrain the average stellar metallicities of massive quiescent galaxies at
$1.0<z<1.3$.

###  4.1 Sample selection

We begin by re-fitting the photometric data described in Section 2.1 for the
812 objects in the M18 parent sample with Bagpipes. We use the best available
redshifts, $z_{\mathrm{best}}$, defined by the following ranking:

  1. 1.

VANDELS DR4 redshift if flag 3 or 4; else

  2. 2.

KMOS redshift (Table 1) if flag 3 or 4; else

  3. 3.

Other spectroscopic redshift (see M18); else

  4. 4.

Photometric redshift, described in Section 2.1.

For consistency with other published work, we assume a Calzetti et al. (2000)
dust attenuation curve and Solar metallicity. We fit for five free parameters:
$V-$band attenuation ($A_{V}$), total stellar mass formed, and the three shape
parameters of a double-power-law SFH model (e.g. Carnall et al. 2018). The
priors assumed for these five parameters are the same as shown in Table 2,
which provides the full list of parameters and priors we use for our more-
sophisticated fits to our stacked spectrum plus photometry in Section 4.3.

Bagpipes uses the 2016 updated version of the Bruzual & Charlot (2003) stellar
population
models111https://www.bruzual.org/~gbruzual/bc03/Updated_version_2016, using
the MILES stellar spectral library (Falcón-Barroso et al. 2011) and updated
stellar evolutionary tracks of Bressan et al. (2012) and Marigo et al. (2013).

From these fits, we obtain stellar masses, $M_{*}$, and rest-frame UVJ
magnitudes. We now define a mass-complete sample for which we can construct a
representative stacked spectrum from our combined VANDELS + KMOS spectroscopic
datasets. We begin by imposing $z_{\mathrm{best}}<1.3$, leaving 512 objects.
We then follow up on the permissive UVJ criteria of M18 (see Section 2.1) by
requiring $(U-V)>0.88\times(V-J)+0.69$. This is the diagonal UVJ cut proposed
by Williams et al. (2009) for quiescent galaxy selection at $z<0.5$. We use
this criterion for our $z\gtrsim 1$ sample, as it has been shown by Carnall et
al. (2018); Carnall et al. (2019a) to consistently select objects with sSFR <
0.2 $t_{\mathrm{H}}^{-1}$ across a wide redshift range (where $t_{\mathrm{H}}$
is the age of the Universe as a function of redshift). This is a widely used
criterion for separating star-forming and quiescent galaxies (e.g. Pacifici et
al. 2016). This further reduces the sample to 409 objects.

As our KMOS observations targeted a more-limited, brighter sub-sample than
VANDELS, the selection criteria detailed in Section 3.1 are the most important
for defining our mass-completeness limit. We define this limit as the lowest
stellar mass for which 90 per cent of more massive galaxies have $J<21.5$,
which is the criterion for inclusion as a priority 1 target in our KMOS
observations. This lowest mass is approximately
$\mathrm{log_{10}(}M_{*}/\mathrm{M_{\odot})}=10.8$, and we therefore adopt
this as our mass-completeness limit. Imposing this stellar-mass criterion
returns 176 objects, with a median stellar mass of
$\mathrm{log_{10}(}M_{*}/\mathrm{M_{\odot})}=11.05$. VANDELS spectroscopy is
available for 77 of these, whereas KMOS spectroscopy is available for 37
objects. A total of 23 objects have both VANDELS and KMOS spectra. In both
cases, objects with spectroscopic data represent a random draw from our
176-object mass-complete sample. All VANDELS and KMOS objects have secure
(flag 3 or 4) spectroscopic redshifts.

Fig. 1 shows the distribution of galaxies in $J-$band magnitude versus
$z_{\mathrm{best}}$. The M18 parent sample is shown with gray circles, whereas
objects in our mass-complete sample are highlighted with open black squares.
The availability of VANDELS and KMOS spectroscopic data is indicated by
different coloured fills of these open black squares, as indicated in the
figure.

To summarise, our mass-complete sample is selected from the M18 parent sample,
introduced in Section 2.1, as follows

  * •

$1.0<z_{\mathrm{best}}<1.3$

  * •

$\mathrm{log_{10}(}M_{*}/\mathrm{M_{\odot})}>10.8$

  * •

$(U-V)>0.88\times(V-J)+0.69$

with $M_{*}$ and UVJ magnitudes determined by Bagpipes fitting.

###  4.2 Stacking procedure

We stack the 77 VANDELS and 37 KMOS spectra for our mass-complete sample to
produce a single representative stacked spectrum for massive quiescent
galaxies at 1.0 < z < 1.3 with
$\mathrm{log_{10}(}M_{*}/\mathrm{M_{\odot})}>10.8$. The individual VANDELS and
KMOS spectra were first shifted to the rest-frame, then flux-normalised, using
rest-frame wavelengths, $\lambda=2500{-}4750$ Å for VANDELS spectra and
$\lambda=5000{-}6000$ Å for KMOS spectra.

The spectra were then resampled to a common wavelength grid using SpectRes
(Carnall 2017). The VANDELS and KMOS spectra were then median stacked
separately to produce two stacks, with uncertainties calculated via the MAD
estimator. Pixels with strong sky line contamination were masked prior to
stacking. The two stacked spectra were then multiplied by the median of the
normalisation factors applied to their input spectra. The two stacks, which do
not overlap in wavelength, were then combined, and finally binned down to 5 Å
sampling.

The stack covers rest-frame wavelengths from $2500{-}6400$ Å, with the
transition from VANDELS to KMOS data at $\simeq 4750$ Å. Fewer than 10 objects
have wavelength coverage between $4700{-}4800$ Å, and we therefore mask this
region from the stack. We do not attempt to match the flux normalisations of
the two stacks at this stage. Instead, during the fitting procedure described
in Section 4.3, we fit a spectrophotometric calibration polynomial to both
sections of the stack separately, allowing their relative normalisations to be
fitted.

The combined stacked spectrum is shown in blue in the top panel of Fig. 2. The
error spectrum is shown in green in the same panel. The SNR per Å and the
number of objects contributing to the stack are shown in the lower two panels.
The wavelengths of key age and metallicity sensitive absorption features are
labelled in black.

We also generate stacked photometry for our mass-complete sample, by taking
the posterior median model fitted to each of the 176 objects in Section 4.1,
shifting this to the median redshift of our sample ($z=1.15$), then
calculating fluxes through a series of UV-IR filters. We use $UBVRiz$, HST
F606W, F814W, F125W and F160W, the HAWKI $K-$band, and $Spitzer-$IRAC Channels
1 and 2. We then produce stacked photometry following the same normalisation
and median stacking process as detailed above for the spectroscopic data.

###  4.3 BAGPIPES fitting of the stacked spectrum and photometry

To constrain the average stellar metallicity of our mass-complete sample, we
fit our stacked spectrum and photometry with
Bagpipes222https://bagpipes.readthedocs.io (Carnall et al. 2018). We apply the
fitting methodology developed in Carnall et al. (2019a), described in full
detail in section 4 of that work. We here provide a brief summary of the
method, including a description of the minor changes that have been made to
the fitted model for this work. A full list of the 19 free parameters of our
model, along with their associated priors, is given in Table 2.

As above, we fit a double-power-law SFH, this time allowing stellar
metallicity to vary with a logarithmic prior from $0.01{-}2.5$ Z⊙. Dust
attenuation is modelled using the form proposed by Salim et al. (2018), which
parameterises dust curve shape with a power-law deviation, $\delta$, from the
Calzetti et al. (2000) model. Emission lines are included in the fit, using a
method based on that of Byler et al. (2017) with the Cloudy photoionization
code (Ferland et al. 2017). The lifetime assumed for the stellar birth clouds
that give rise to nebular emission is 10 Myr, and $A_{V}$ is doubled for
emission from stars younger than this, as well as nebular line and continuum
emission.

Two separate second-order multiplicative Chebyshev polynomials are fitted to
the VANDELS and KMOS portions of the stacked spectrum to model any
imperfections in spectrophotometric calibration. We model the covariance
matrix for our spectroscopic data as follows. The diagonal terms are given by
the square of the green error spectrum plotted in Fig. 2, multiplied by a
factor, $a^{2}$, to allow for potential underestimation of uncertainties. We
fit $a$ with a logarithmic prior from $0.1{-}10$. The off-diagonal terms are
modelled with a Gaussian process, using an exponential-squared kernel.

The stacked spectrum is shifted to the median redshift of our sample
($z=1.15$) for fitting, and redshift is allowed to vary within a narrow range
about this value. A Gaussian prior is applied, with a standard deviation of
0.001, and a maximum deviation of 0.005. Velocity dispersion is modelled by
convolving the spectral model with a Gaussian kernel in velocity space, with
the standard deviation permitted to vary from $100{-}500$ km s-1 with a
logarithmic prior.

We exclude rest-frame wavelengths $<3550$ Å from the fit, as the MILES library
does not provide coverage bluer than this, with the Bruzual & Charlot (2003)
models instead employing a combination of theoretical stellar spectral
libraries. We mask the [O ii] line at 3727 Å, as it is currently unclear
whether [O ii] emission in quiescent galaxies originates from star-forming
regions, which is the only source of line emission in our Bagpipes model.
Finally, we also mask the Na D absorption feature at 5895 Å, as this has a
potential strong interstellar medium component (however, see Conroy et al.
2014).

To sample the posterior distribution for a model, Bagpipes uses the MultiNest
nested sampling algorithm (Skilling 2006; Feroz & Hobson 2008; Feroz et al.
2009; Feroz et al. 2019), via the PyMultiNest interface (Buchner et al. 2014).
For fitting our stacked data, MultiNest was run with 1000 live points,
requiring $\sim 500$ CPU hours.

Table 2: Parameters and priors for the Bagpipes model we fit to our stacked spectrum and photometry. The model is described in Section 4.3, and is adapted from the model presented in section 4 of Carnall et al. (2019b). The first 10 parameters describe our physical model, whereas the final 9 describe our systematic uncertainties model. For Gaussian priors, $\mu$ is the mean and $\sigma$ the standard deviation. Logarithmic priors are uniform in log base ten of the parameter. Component | Parameter | Symbol / Unit | Range | Prior | Hyperparameters
---|---|---|---|---|---
Global | Redshift | $z$ | (1.145, 1.155) | Gaussian | $\mu=1.15$ |  $\sigma$ = 0.001
| Velocity Dispersion |  $\sigma_{\mathrm{vel}}$ / km s-1 | (100, 500) | logarithmic |  |
SFH | Stellar mass formed | $M_{*}\ /\ \mathrm{M_{\odot}}$ | (1, $10^{13}$) | logarithmic |  |
| Metallicity | $Z\ /\ \mathrm{Z_{\odot}}$ | (0.01, 2.5) | logarithmic |  |
| Falling slope | $\alpha$ | (0.1, 1000) | logarithmic |  |
| Rising slope | $\beta$ | (0.1, 1000) | logarithmic |  |
| Peak time |  $\tau$ / Gyr | (0.1, $t_{\mathrm{H}}$) | uniform |  |
Dust |  $V-$band attenuation |  $A_{V}$ / mag | (0, 4) | uniform |  |
| Deviation from Calzetti et al. (2000) slope | $\delta$ | ($-0.3$, 0.3) | Gaussian | $\mu=0$ |  $\sigma$ = 0.1
| Strength of 2175 Å bump | $B$ | (0, 5) | uniform |  |
Calibration | VANDELS zero order | $P_{V0}$ | (0.5, 1.5) | Gaussian | $\mu=1$ |  $\sigma$ = 0.25
| VANDELS first order | $P_{V1}$ | ($-0.5$, 0.5) | Gaussian | $\mu=0$ |  $\sigma$ = 0.25
| VANDELS second order | $P_{V2}$ | ($-0.5$, 0.5) | Gaussian | $\mu=0$ |  $\sigma$ = 0.25
| KMOS zero order | $P_{K0}$ | (0.5, 1.5) | Gaussian | $\mu=1$ |  $\sigma$ = 0.25
| KMOS first order | $P_{K1}$ | ($-0.5$, 0.5) | Gaussian | $\mu=0$ |  $\sigma$ = 0.25
| KMOS second order | $P_{K2}$ | ($-0.5$, 0.5) | Gaussian | $\mu=0$ |  $\sigma$ = 0.25
Noise | White noise scaling | $a$ | (0.1, 10) | logarithmic |  |
| Correlated noise amplitude |  $b$ / $f_{\mathrm{max}}$ | (0.0001, 1) | logarithmic |  |
| Correlation length |  $l$ / $\Delta\lambda$ | (0.01, 1) | logarithmic |  |
Figure 3: Top: stellar metallicity measurements for massive quiescent galaxies
at $z\gtrsim 1$, compared with the local Universe ($z\lesssim 0.1$). The
metallicity we measure for our stack with Alf ([Z/H] = $-0.11\pm 0.12$) is not
shown, as it is almost identical to our Bagpipes result ([Z/H] = $-0.13\pm
0.08$). Our results, which are broadly consistent with previous work at
similar redshifts, suggest that [Z/H] was $\sim 0.3$ dex lower at $z\gtrsim 1$
compared with the present day. Lower left: Fe abundances for the subset of
studies, including ours, that use Alf. We find [Fe/H] $\sim 0.15$ dex lower
than the local Universe. Lower right: alpha enhancement, traced by [Mg/Fe]. We
find $\sim 0.2$ dex enhancement, suggesting no evolution across the last $\sim
8$ Gyr since $z\gtrsim 1$. The dashed horizontal lines show the Solar value.

###  4.4 ALF fitting of the stacked spectrum

In addition to fitting our stacked data with Bagpipes, we also carried out an
independent analysis with the Alf code333https://www.github.com/cconroy20/alf
(Conroy & van Dokkum 2012; Conroy et al. 2018). This decision was made,
firstly to provide a cross-check on our results via a more established method,
and secondly as the Bruzual & Charlot (2003) models assume scaled-Solar
abundances. By contrast, Alf allows individual element abundances to vary
separately, allowing us to constrain the level of alpha enhancement in our
target population, a valuable indicator of their formation timescales (e.g.
Thomas et al. 2003).

The Alf code is designed for fitting optical to near-IR continuum spectroscopy
for old ($\gtrsim$ 1 Gyr) stellar populations. Originally designed for
constraining the initial mass function, recent applications have focused on
stellar metallicities, with the code having been developed to fit individual
abundances for up to 19 elements.

The code also makes use of the MILES stellar spectral library, and includes
empirical stellar population models spanning $0.37-2.4\,\mu$m. We therefore
fit our stacked spectrum across a wavelength range similar to that described
in Section 4.3, this time also omitting rest-frame wavelengths from
$3550{-}3700$ Å. Alf continuum normalises input spectra using a high-order
polynomial, with one order per 100 Å of rest-frame spectral coverage, prior to
fitting. The code does not currently include the capability to fit photometric
data, so our stacked photometry was not used in this analysis. The code uses
the MCMC sampler Emcee (Foreman-Mackey et al. 2013).

We run Alf in simple mode, which includes 13 free parameters: redshift,
velocity dispersion, stellar age (a single burst SFH is assumed), total
stellar metallicity, [Z/H], and abundances for 9 individual elements,
including Fe and Mg.

##  5 Results

The posterior median model fitted to our stacked spectrum with Bagpipes is
shown in red in the top panel of Fig. 2. We measure a stellar metallicity of
[Z/H] = $-0.13\pm 0.08$. We also measure a mean stellar age of
$2.5^{+0.6}_{-0.4}$ Gyr, which, at $z=1.15$, corresponds to a mean formation
time, $t_{\mathrm{form}}=2.7^{+0.4}_{-0.6}$ Gyr after the Big Bang, or a
formation redshift of $z_{\mathrm{form}}=2.4^{+0.6}_{-0.3}$. We measure a dust
curve slope consistent with Calzetti et al. (2000), and an $A_{V}=0.48\pm
0.08$ magnitudes. This is in good agreement with Belli et al. (2019) and
Carnall et al. (2019a) (but see appendix B of van der Wel et al. 2021).

With Alf, we obtain element abundances [Fe/H] = $-0.18\pm 0.08$ and [Mg/H] =
$0.07\pm 0.09$. The enhancement of Mg relative to Fe can be used as a proxy
for alpha enhancement, and we obtain [Mg/Fe] = $0.23\pm 0.12$. We convert from
Fe and Mg abundances to total metallicity using [Z/H] = [Fe/H] + 0.94[Mg/H]
(e.g. Thomas et al. 2003), for consistency with other recent work (e.g. Kriek
et al. 2019). This yields a total metallicity of [Z/H] = $-0.11\pm 0.12$. The
metallicities measured for our stacked spectrum with Bagpipes and Alf are
therefore fully consistent. The stellar population ages are also consistent,
with Alf returning an age of $2.4\pm 0.3$ Gyr. Our stellar metallicity results
are summarised in Fig. 3, and compared to results from the literature in
Section 6.

##  6 Discussion

###  6.1 Stellar Metallicities

####  6.1.1 Comparisons with other results at $z\gtrsim 1$

Whilst this work presents the first measurement of the average stellar
metallicity of a mass-selected sample of massive quiescent galaxies at $z>1$,
several other studies have analysed individual bright objects, or magnitude-
selected samples, at similar redshifts. In Fig. 3, we compare our results for
[Z/H], [Fe/H] and [Mg/Fe] with other recent work at $1.0<z\lesssim 1.5$, as
well as results from the local Universe, which will be discussed in Section
6.1.2.

Kriek et al. (2019) report stellar metallicities for three massive quiescent
galaxies at $z\sim 1.4$ via full spectral fitting of Keck-LRIS and MOSFIRE
spectroscopy, also with the Alf code. These results are shown in Fig. 3 with
orange squares, and are in good agreement with the measurements we derive from
our stacked spectrum. Their most massive object is more Mg-enhanced, leading
to a higher total metallicity, consistent with the positive correlation
between mass and metallicity observed in the local Universe.

Lonoce et al. (2015); Lonoce et al. (2020) also report stellar metallicities
for two galaxies at $z\sim 1.4$ via full spectral fitting, which are shown in
the top-left panel of Fig. 3 with red triangles. These results are also
consistent with our findings, though it should be noted that Lonoce et al.
(2020) obtain strongly contrasting results via an alternative, spectral-index-
fitting analysis: [Z/H] $\sim-0.7$ and $+0.6$ for their lower and higher-mass
galaxies, respectively.

Onodera et al. (2015) report a stellar metallicity of [Z/H] =
$0.24^{+0.20}_{-0.14}$ via spectral index fitting, for a stacked spectrum
constructed from a magnitude-selected sample ($K<21.5$) of 24 quiescent
galaxies at $1.2\lesssim z\lesssim 2.0$. This is considerably higher than our
result at $1.0<z<1.3$, however, their relatively large reported uncertainties
mean the two results are not strongly in tension. It should also be noted that
the average mass of galaxies in the sample of Onodera et al. (2015) is 0.35
dex higher than our sample.

Finally, we show, with gray error bars, the stellar metallicities we obtained
for 53 individual massive quiescent galaxies in Carnall et al. (2019a), via
Bagpipes full spectral fitting of intermediate-SNR ($\sim 10$ Å-1) VANDELS
rest-frame near-UV spectroscopy. The Carnall et al. (2019a) results are in
very good agreement with the new analysis we present in this work, which
differs from our previous work in three key respects: the use of stacking, the
inclusion of additional rest-frame optical KMOS data, and the exclusion of
wavelengths, $\lambda<3550$ Å, which were fitted with theoretical stellar
templates in Carnall et al. (2019a). Our new result validates the
metallicities measured by our previous analysis, and hence the ages/SFHs that
were the main result of that work.

The individual-object results of Carnall et al. (2019a) shown in Fig. 3 do not
suggest a strong stellar mass-metallicity relation at $z\gtrsim 1$ above
log$(M_{*}/\mathrm{M_{\odot}})=10.5$, and are consistent with a flat
relationship. However, a relatively weak mass-metallicity relation, as seen in
the local Universe, is by no means ruled out. We have investigated the
possibility of splitting the mass-complete sample in this work into two mass
bins to further probe the mass-metallicity relation at $z\gtrsim 1$, however
our KMOS dataset does not provide sufficient SNR to obtain meaningful
constraints when splitting our sample into two or more bins.

We conclude that our new results are generally consistent with previous work
at $z\gtrsim 1$. In particular, multiple full spectral fitting analyses agree
on average stellar metallicities, [Z/H] $\sim-0.1$ for massive quiescent
galaxies with log${}_{10}(M_{*}/\mathrm{M}_{\odot})\sim 11$ at $z\gtrsim 1$.
Our measurement of [Fe/H] $\sim-0.2$ is in good agreement with Kriek et al.
(2019), as is our measurement of $\sim 0.2$ dex of Mg enhancement over Fe.
Kriek et al. (2016, 2019) show that this level of alpha enhancement suggests
formation timescales of $0.2{-}1.0$ Gyr, again in good agreement with the SFHs
reported by Carnall et al. (2019a).

####  6.1.2 Comparisons with the Local Universe

In the top-left panel of Fig. 3, we also show three determinations of the
stellar mass-metallicity relation in the local Universe ($z\lesssim 0.1$).
These three sets of results are all in reasonable agreement, however there are
key differences in methodology and sample selection that should be considered
when making comparisons. All three studies use data from the Sloan Digital Sky
Survey (SDSS; York et al. 2000).

Gallazzi et al. (2005) and Panter et al. (2008) report mean stellar
metallicity as a function of stellar mass for all galaxies, including both
star-forming and quiescent objects. However, the dominance of quiescent
galaxies at high masses in the local Universe means these results are still
reasonably comparable with the $z\gtrsim 1$ results shown. In addition,
Gallazzi et al. (2005) report $r-$band light-weighted metallicities, rather
than the mass-weighted quantities reported by the other studies. This may
explain the higher metallicities they report, as their results are weighted
more heavily towards younger stars.

The results of Conroy et al. (2014) are the most comparable with our results
at $z\gtrsim 1$, both in terms of methodology and sample selection. They
select quiescent galaxies by requiring no H$\alpha$ or [O ii] emission, and
require that objects lie on the Fundamental Plane, using the central slice
defined by Graves et al. (2010). At log$(M_{*}/\mathrm{M_{\odot}})\sim 11$,
they report [Z/H] = 0.21. Our result therefore implies, at fixed stellar mass,
$\sim 0.3$ dex evolution in the average [Z/H] for massive quiescent galaxies
across the $\sim 8$ Gyr from $z\sim 1.15$ to the present. From the lower-right
panel of Fig. 3, it can be seen that, surprisingly, this 0.3 dex evolution in
[Z/H] is not accompanied by any change in alpha (Mg) enhancement, with [Mg/Fe]
remaining at $\sim 0.2$.

Alpha elements are primarily produced by core-collapse supernovae (CCSNe),
which are the end-point in the evolution of massive stars. Enrichment with
alpha elements therefore rapidly follows the onset of star formation, within a
few Myr. By contrast, Fe-peak elements are produced by both CCSNe and Type Ia
supernovae (SNIa), which occur on much longer ($\sim$ Gyr) timescales. This
means that both Fe abundance and alpha enhancement in stars are linked to the
formation timescale of the galaxy, with lower Fe, more highly alpha enhanced
stellar populations expected in galaxies that formed on short timescales (e.g.
Weinberg et al. 2017; Kobayashi et al. 2020).

As the number density of massive quiescent galaxies steadily rises over cosmic
time, those that exist at $z\gtrsim 1$ are expected to be a biased sub-sample
of the local quiescent galaxy population: the earliest and fastest formed.
This is often referred to as progenitor bias. We would therefore expect to
see, on average, both lower Fe abundance and greater alpha enhancement for
$z\gtrsim 1$ quiescent galaxies, compared with the local population.

Beverage et al. (2021) have recently reported similar evolution of the
quiescent galaxy population from the local Universe to $z\sim 0.7$: falling
[Z/H] and [Fe/H], but constant [Mg/Fe] with increasing redshift (see also Choi
et al. 2014; Leethochawalit et al. 2019). They suggest the expected increase
in [Mg/Fe] might begin to become apparent at higher redshifts, pointing to
several individual-object studies that find very high [Mg/Fe] for individual
ultra-massive galaxies at $z\sim 2$ (Kriek et al. 2016; Jafariyazani et al.
2020).

Our results demonstrate the expected increase in [Mg/Fe] is still not apparent
at $z\gtrsim 1$ for typical log$(M_{*}/\mathrm{M_{\odot}})\sim 11$ quiescent
galaxies. This is in agreement with Kriek et al. (2019) at $z\sim 1.4$
(bottom-right panel of Fig. 3), who report higher alpha enhancement than the
local Universe for only a single, extremely massive galaxy with
log$(M_{*}/\mathrm{M_{\odot}})\sim 11.7$. Our results suggest the very high
[Mg/Fe] values (implying ultra-short, $\sim 100$ Myr, formation timescales)
recently reported by Kriek et al. (2016, 2019) and Jafariyazani et al. (2020)
for ultra-massive galaxies at $z\sim 1-2$ are more closely associated with
their exceptionally high masses, in accordance with the well-known downsizing
trend, rather than their high redshifts.

A detailed understanding of these issues will require a dataset of similar
quality to SDSS at $z\gtrsim 1$. Several planned surveys with upcoming
instruments, such as the Multi-Object Optical and Near-infrared Spectrograph
(MOONS) GTO survey MOONRISE (Cirasuolo et al. 2020; Maiolino et al. 2020),
will provide much larger numbers of high-SNR spectra at $z\gtrsim 1$. This
will allow individual metallicity determinations for large numbers of bright
targets, and high-SNR stacking experiments to be performed with fine
resolution in stellar mass and other key parameters, such as physical size.

###  6.2 Stellar Ages

As discussed in Section 6.1.1, the stellar metallicity we derive for our
stacked spectrum is in good agreement with those we obtained for individual
galaxies in Carnall et al. (2019a). In that work, we derived the following
relationship between formation time, $t_{\mathrm{form}}$ (which we define as
the time corresponding to the mean stellar age), and stellar mass for
quiescent galaxies at $z\gtrsim 1$:

$\bigg{(}\dfrac{t_{\mathrm{form}}}{\mathrm{Gyr}}\bigg{)}\ =\
2.56^{+0.12}_{-0.10}\ -\ 1.48^{+0.34}_{-0.39}\
\mathrm{log}_{10}\bigg{(}\dfrac{M_{*}}{10^{11}\mathrm{M_{\odot}}}\bigg{)}.$
(1)

At the median stellar mass of
$\mathrm{log_{10}(}M_{*}/\mathrm{M_{\odot})}=11.05$ for the mass-complete
sample from which we construct our stack, this relationship predicts a
formation time of 2.5$\pm$0.1 Gyr after the Big Bang. We derive consistent
results from both the Bagpipes and Alf analyses of our stack:
$t_{\mathrm{form}}=2.7^{+0.4}_{-0.6}$ Gyr and 2.4$\pm$0.3 Gyr respectively.

Several other recent studies have reported ages for massive quiescent galaxies
at $z\gtrsim 1$ using similar methodologies (Belli et al. 2019; Estrada-
Carpenter et al. 2020; Tacchella et al. 2021). These results are summarised in
fig. 10 of Tacchella et al. (2021), demonstrating that the formation times
they derive for their $z\sim 0.8$ sample are earlier on average than those we
derive at $z\gtrsim 1$. This is the opposite to the effect that would be
expected due to progenitor bias.

Due to the strong age-metallicity degeneracy in galaxy spectral shapes,
discrepant metallicities is an obvious potential cause for the unexpected
differences in derived ages between these studies. However, as can be seen
from fig. 19 of Tacchella et al. (2021), the average metallicity they derive
is just below Solar444The average metallicity derived for the UVJ-quiescent
sample in Tacchella et al. (2021) is [Z/H] = $-0.13\pm 0.02$ (private
communication)., as would be expected considering the comparison between local
Universe and $z\gtrsim 1$ metallicities shown in the top panel of Fig. 3. The
metallicities we derive in this work and Carnall et al. (2019a) are therefore
fully consistent with Tacchella et al. (2021), and are ruled out as the cause
of the differences in derived ages.

The remaining likely cause for these differences lies in the SFH models used
by the respective studies. In this work and Carnall et al. (2019a), we have
used the double-power-law SFH model introduced by Carnall et al. (2018). This
is a parametric model, and similar to the models used by Belli et al. (2019),
who report results consistent with Carnall et al. (2019a). By contrast,
Tacchella et al. (2021) use the non-parametric, continuity SFH prior
introduced by Leja et al. (2019a), which is known to return older stellar
populations than more traditional parametric models (Leja et al., 2019b).
Interestingly, Estrada-Carpenter et al. (2020) also use the non-parametric
continuity prior, finding results consistent with Tacchella et al. (2021).

It seems likely, therefore, that the underlying cause for the discrepant
results in the recent literature for the formation redshifts of massive
quiescent galaxies at $z\gtrsim 1$ is the use of different SFH models, rather
than conflicting derived metallicities. These SFH models constitute a set of
prior beliefs about when and how galaxies form, introducing an unavoidable
subjectivity into results obtained in the absence of strongly constraining
data (e.g. Ocvirk et al. 2006; Carnall et al. 2019b; Johnson et al. 2021).

To evaluate the success of these subjective prior choices, one should perform
cross-validation checks, testing whether results obtained via spectral fitting
are consistent with other, better established results in the literature (e.g.
Wuyts et al. 2011; Carnall et al. 2019b). In this case, we would wish to use
the SFHs derived for $z\gtrsim 1$ massive quiescent galaxies to predict the
stellar-mass function for quiescent galaxies at $z\sim 2-4$, then to compare
these predictions with direct observations of $z>2$ mass functions and/or
number densities.

To perform this challenging experiment, we are currently missing two vital
components. Firstly, the number of massive quiescent galaxies at $z\gtrsim 1$
with high-SNR rest-frame UV-optical spectroscopy is currently no more than a
few hundred. This, in combination with the significant uncertainties on
individual recovered SFHs, leads to huge statistical uncertainties in
predicted mass functions at $z>2$. To make useful predictions, vastly larger
samples will be required, numbering tens to hundreds of thousands of galaxies.
A number of upcoming surveys, such as MOONRISE (Maiolino et al. 2020), will
provide the necessary datasets over the next $\sim 5$ years.

The second missing component is a precise determination of the observed
stellar-mass function for massive quiescent galaxies at $z>2$, with current
studies suffering from considerable uncertainties in sample selection (e.g.
Schreiber et al. 2018; Carnall et al. 2020). Large imaging surveys with JWST
(e.g. PRIMER; Dunlop et al. 2021) will provide the extremely deep, high-
resolution infra-red imaging necessary to select mass-complete samples of
quiescent galaxies with confidence out to the highest redshifts, resulting in
precise constraints on quiescent galaxy stellar-mass functions at $z>2$.

Once these datasets are in place, it should be possible to clearly distinguish
between the parametric and non-parametric SFH models that currently return
substantially different formation redshifts for massive quiescent galaxies at
$z\gtrsim 1$.

##  7 Conclusion

In this work, we have combined data from the VANDELS survey with new KMOS
observations to construct a representative stacked spectrum for quiescent
galaxies at $1.0<z<1.3$ with $\mathrm{log}_{10}(M_{*}\ /\
\mathrm{M_{\odot}})>10.8$, covering rest-frame wavelengths,
$\lambda=2500{-}6400$ Å. The stacked spectrum is shown in Fig. 2. We also
report 25 new spectroscopic redshifts (Table 1) from our KMOS data.

We fit our stacked spectrum with Bagpipes, obtaining a stellar metallicity,
[Z/H] = $-0.13\pm 0.08$. We also obtain a formation time,
$t_{\mathrm{form}}=2.7^{+0.4}_{-0.6}$ Gyr after the Big Bang, corresponding to
a formation redshift of $z_{\mathrm{form}}=2.4^{+0.6}_{-0.3}$. Both of these
results are consistent with the results we presented in Carnall et al.
(2019a), which were obtained by fitting only VANDELS data, and working with
individual spectra rather than stacking. We also fit our stacked spectrum with
the Alf code, obtaining a consistent result for [Z/H], as well as an iron
abundance, [Fe/H] = $-0.18\pm 0.08$, and alpha enhancement, [Mg/Fe] = $0.23\pm
0.12$ (see Fig. 3).

By comparing our results at $z\gtrsim 1$ with results from the local Universe
($z\lesssim 0.1$), we demonstrate that the average [Z/H] for
$\mathrm{log}_{10}(M_{*}\ /\ \mathrm{M_{\odot}})\sim 11$ quiescent galaxies
has risen by $\sim 0.3$ dex across the $\sim 8$ Gyr since $z\gtrsim 1$,
whereas [Fe/H] has risen by $\sim 0.15$ dex. However, the alpha enhancement,
[Mg/Fe], we measure is the same as found by Conroy et al. (2014) in the local
Universe, implying no evolution at fixed stellar mass across at least the last
$\sim 8$ Gyr.

Given that $z\gtrsim 1$ massive quiescent galaxies are a biased sub-sample of
the local quiescent population (those that formed fastest and earliest), we
would expect them to have lower [Fe/H], as observed, but we would also expect
greater alpha enhancement. Our finding of no redshift evolution in [Mg/Fe] at
$\mathrm{log}_{10}(M_{*}\ /\ \mathrm{M_{\odot}})\sim 11$ from $z\gtrsim 1$ to
the present is in agreement with the results of Kriek et al. (2019) at $z\sim
1.4$ and Beverage et al. (2021) at $z\sim 0.7$.

This result suggests that the highly alpha enhanced, ultra-massive
($\mathrm{log}_{10}(M_{*}\ /\ \mathrm{M_{\odot}})\gtrsim 11.5$) galaxies
recently reported at $z\sim 1-2$ by some authors (e.g. Kriek et al. 2016,
2019; Leethochawalit et al. 2019) are highly alpha enhanced due to their
extreme masses, in accordance with the well-known downsizing trend, rather
than being typical of the $z\gtrsim 1$ quiescent population.

Recently, Tacchella et al. (2021) have highlighted differences in the average
formation redshifts measured for $z\gtrsim 1$ massive quiescent galaxies by
several recent studies. In particular, they, in agreement with Estrada-
Carpenter et al. (2020), find earlier formation than Belli et al. (2019) and
Carnall et al. (2019a). Given the latter two studies are at higher redshifts,
this is the opposite to the effect that would be expected due to progenitor
bias.

We demonstrate that the metallicities recovered by Tacchella et al. (2021) at
$z\sim 0.8$ are consistent with the results of this work and Carnall et al.
(2019a). This means discrepancies in recovered metallicities are ruled out as
the cause of the differences in recovered ages. We therefore conclude that the
differences in formation redshifts obtained by different authors are likely
due to the use of different star-formation-history models, with Estrada-
Carpenter et al. (2020) and Tacchella et al. (2021) using the Leja et al.
(2019a) continuity non-parametric SFH model, and Belli et al. (2019) and
Carnall et al. (2019a) using more-traditional parametric models.

To determine which of these models more accurately represents the SFHs of
$z\gtrsim 1$ massive quiescent galaxies, we will require much larger
statistical samples at $z\gtrsim 1$ with high-SNR rest-frame UV-optical
spectroscopy, as will be provided by upcoming instruments such as MOONS
(Cirasuolo et al., 2020). This will allow us to make firm predictions for
$z>2$ quiescent-galaxy stellar-mass functions via spectral fitting with
different SFH models, which can be compared with the precise $z>2$ stellar-
mass functions that will be provided by data from the upcoming James Webb
Space Telescope.

## Acknowledgements

The authors would like to thank Corentin Schreiber, David Maltby and Omar
Almaini for valuable advice about KMOS observing strategies and SNRs. We would
also like to thank Michael Hilker at the ESO helpdesk for assistance with KMOS
data reduction. A. C. Carnall would like to thank the Leverhulme Trust for
their support via the Leverhulme Early Career Fellowship scheme. A. Cimatti
acknowledges support from the grant PRIN MIUR 2017$-$20173ML3WW_001. B.
Garilli acknowledges support from the grants Premiale MITIC 2017 and INAF PRIN
"Mainstream 2019". Based on observations made with ESO Telescopes at the La
Silla or Paranal Observatories under programme ID(s) 194.A-2003(A-T) and
0104.B-0885(A).

## Data Availability

The VANDELS survey is a European Southern Observatory Public Spectroscopic
Survey. The full spectroscopic dataset, together with photometric catalogues
and derived quantities are available from http://vandels.inaf.it, as well as
from the ESO archive https://www.eso.org/qi. The KMOS data used in this work
is available from the ESO archive under programme ID 0104.B-0885(A). Reduced
KMOS data products, our stacked spectrum and model posteriors may be made
available upon reasonable request.

## References

  * Asplund et al. (2009) Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481
  * Belli et al. (2019) Belli S., Newman A. B., Ellis R. S., 2019, ApJ, 874, 17
  * Beverage et al. (2021) Beverage A. G., Kriek M., Conroy C., Bezanson R., Franx M., van der Wel A., 2021, arXiv e-prints, arXiv:2105.12750
  * Bressan et al. (2012) Bressan A., Marigo P., Girardi L., Salasnich B., Dal Cero C., Rubele S., Nanni A., 2012, MNRAS, 427, 127
  * Bruzual & Charlot (2003) Bruzual G., Charlot S., 2003, MNRAS, 344, 1000
  * Buchner et al. (2014) Buchner J., et al., 2014, A&A, 564, A125
  * Byler et al. (2017) Byler N., Dalcanton J. J., Conroy C., Johnson B. D., 2017, ApJ, 840, 44
  * Calzetti et al. (2000) Calzetti D., Armus L., Bohlin R. C., Kinney A. L., Koornneef J., Storchi-Bergmann T., 2000, ApJ, 533, 682
  * Carnall (2017) Carnall A. C., 2017, arXiv e-prints, arXiv:1705.05165
  * Carnall et al. (2018) Carnall A. C., McLure R. J., Dunlop J. S., Davé R., 2018, MNRAS, 480, 4379
  * Carnall et al. (2019a) Carnall A. C., et al., 2019a, MNRAS, 490, 417
  * Carnall et al. (2019b) Carnall A. C., Leja J., Johnson B. D., McLure R. J., Dunlop J. S., Conroy C., 2019b, ApJ, 873, 44
  * Carnall et al. (2020) Carnall A. C., et al., 2020, MNRAS, 496, 695
  * Cecchi et al. (2019) Cecchi R., Bolzonella M., Cimatti A., Girelli G., 2019, ApJ, 880, L14
  * Choi et al. (2014) Choi J., Conroy C., Moustakas J., Graves G. J., Holden B. P., Brodwin M., Brown M. J. I., van Dokkum P. G., 2014, ApJ, 792, 95
  * Cirasuolo et al. (2020) Cirasuolo M., et al., 2020, The Messenger, 180, 10
  * Coelho et al. (2020) Coelho P. R. T., Bruzual G., Charlot S., 2020, MNRAS, 491, 2025
  * Conroy & van Dokkum (2012) Conroy C., van Dokkum P. G., 2012, ApJ, 760, 71
  * Conroy et al. (2014) Conroy C., Graves G. J., van Dokkum P. G., 2014, ApJ, 780, 33
  * Conroy et al. (2018) Conroy C., Villaume A., van Dokkum P. G., Lind K., 2018, ApJ, 854, 139
  * Cullen et al. (2019) Cullen F., et al., 2019, MNRAS, p. 1344
  * Cullen et al. (2021) Cullen F., et al., 2021, MNRAS,
  * Dahlen et al. (2013) Dahlen T., et al., 2013, ApJ, 775, 93
  * Dunlop et al. (2021) Dunlop J. S., et al., 2021, PRIMER: Public Release IMaging for Extragalactic Research, JWST Proposal. Cycle 1
  * Estrada-Carpenter et al. (2019) Estrada-Carpenter V., et al., 2019, ApJ, 870, 133
  * Estrada-Carpenter et al. (2020) Estrada-Carpenter V., et al., 2020, ApJ, 898, 171
  * Falcón-Barroso et al. (2011) Falcón-Barroso J., Sánchez-Blázquez P., Vazdekis A., Ricciardelli E., Cardiel N., Cenarro A. J., Gorgas J., Peletier R. F., 2011, A&A, 532, A95
  * Ferland et al. (2017) Ferland G. J., et al., 2017, Rev. Mex. Astron. Astrofis., 53, 385
  * Feroz & Hobson (2008) Feroz F., Hobson M. P., 2008, MNRAS, 384, 449
  * Feroz et al. (2009) Feroz F., Hobson M. P., Bridges M., 2009, MNRAS, 398, 1601
  * Feroz et al. (2019) Feroz F., Hobson M. P., Cameron E., Pettitt A. N., 2019, The Open Journal of Astrophysics, 2, 10
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Galametz et al. (2013) Galametz A., et al., 2013, ApJS, 206, 10
  * Gallazzi et al. (2005) Gallazzi A., Charlot S., Brinchmann J., White S. D. M., Tremonti C. A., 2005, MNRAS, 362, 41
  * Garilli et al. (2010) Garilli B., Fumana M., Franzetti P., Paioro L., Scodeggio M., Le Fèvre O., Paltani S., Scaramella R., 2010, PASP, 122, 827
  * Garilli et al. (2021) Garilli B., et al., 2021, A&A, 647, A150
  * Girelli et al. (2019) Girelli G., Bolzonella M., Cimatti A., 2019, A&A, 632, A80
  * Graves et al. (2010) Graves G. J., Faber S. M., Schiavon R. P., 2010, ApJ, 721, 278
  * Guo et al. (2013) Guo Y., et al., 2013, ApJS, 207, 24
  * Horne (1986) Horne K., 1986, PASP, 98, 609
  * Jafariyazani et al. (2020) Jafariyazani M., Newman A. B., Mobasher B., Belli S., Ellis R. S., Patel S. G., 2020, ApJ, 897, L42
  * Johnson et al. (2021) Johnson B. D., Leja J., Conroy C., Speagle J. S., 2021, ApJS, 254, 22
  * Kobayashi et al. (2020) Kobayashi C., Karakas A. I., Lugaro M., 2020, ApJ, 900, 179
  * Kriek et al. (2016) Kriek M., et al., 2016, Nature, 540, 248
  * Kriek et al. (2019) Kriek M., et al., 2019, ApJ, 880, L31
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Le Fèvre et al. (2013) Le Fèvre O., et al., 2013, A&A, 559, A14
  * Leethochawalit et al. (2019) Leethochawalit N., Kirby E. N., Ellis R. S., Moran S. M., Treu T., 2019, ApJ, 885, 100
  * Leja et al. (2019a) Leja J., Carnall A. C., Johnson B. D., Conroy C., Speagle J. S., 2019a, ApJ, 876, 3
  * Leja et al. (2019b) Leja J., et al., 2019b, ApJ, 877, 140
  * Lonoce et al. (2015) Lonoce I., et al., 2015, MNRAS, 454, 3912
  * Lonoce et al. (2020) Lonoce I., Maraston C., Thomas D., Longhetti M., Parikh T., Guarnieri P., Comparat J., 2020, MNRAS, 492, 326
  * Maiolino & Mannucci (2019) Maiolino R., Mannucci F., 2019, A&ARv, 27, 3
  * Maiolino et al. (2020) Maiolino R., et al., 2020, The Messenger, 180, 24
  * Marigo et al. (2013) Marigo P., Bressan A., Nanni A., Girardi L., Pumo M. L., 2013, MNRAS, 434, 488
  * Marsan et al. (2020) Marsan Z. C., et al., 2020, arXiv e-prints, arXiv:2010.04725
  * McLeod et al. (2021) McLeod D. J., McLure R. J., Dunlop J. S., Cullen F., Carnall A. C., Duncan K., 2021, MNRAS, 503, 4413
  * McLure et al. (2013) McLure R. J., et al., 2013, MNRAS, 428, 1088
  * McLure et al. (2018) McLure R. J., et al., 2018, MNRAS, 479, 25
  * Merlin et al. (2019) Merlin E., et al., 2019, MNRAS, 490, 3309
  * Mowla et al. (2019a) Mowla L., van der Wel A., van Dokkum P., Miller T. B., 2019a, ApJ, 872, L13
  * Mowla et al. (2019b) Mowla L. A., et al., 2019b, ApJ, 880, 57
  * Ocvirk et al. (2006) Ocvirk P., Pichon C., Lançon A., Thiébaut E., 2006, MNRAS, 365, 46
  * Onodera et al. (2015) Onodera M., et al., 2015, ApJ, 808, 161
  * Pacifici et al. (2012) Pacifici C., Charlot S., Blaizot J., Brinchmann J., 2012, MNRAS, 421, 2002
  * Pacifici et al. (2016) Pacifici C., et al., 2016, ApJ, 832, 79
  * Panter et al. (2008) Panter B., Jimenez R., Heavens A. F., Charlot S., 2008, MNRAS, 391, 1117
  * Pentericci et al. (2018) Pentericci L., et al., 2018, A&A, 616, A174
  * Salim et al. (2018) Salim S., Boquien M., Lee J. C., 2018, ApJ, 859, 11
  * Santini et al. (2021) Santini P., et al., 2021, A&A, 652, A30
  * Schreiber et al. (2018) Schreiber C., et al., 2018, A&A, 618, A85
  * Sherman et al. (2020) Sherman S., et al., 2020, MNRAS, 491, 3318
  * Skilling (2006) Skilling J., 2006, Bayesian Anal., 1, 833
  * Steidel et al. (2016) Steidel C. C., Strom A. L., Pettini M., Rudie G. C., Reddy N. A., Trainor R. F., 2016, ApJ, 826, 159
  * Straatman et al. (2014) Straatman C. M. S., et al., 2014, ApJ, 783, L14
  * Straatman et al. (2016) Straatman C. M. S., et al., 2016, ApJ, 830, 51
  * Suess et al. (2019a) Suess K. A., Kriek M., Price S. H., Barro G., 2019a, ApJ, 877, 103
  * Suess et al. (2019b) Suess K. A., Kriek M., Price S. H., Barro G., 2019b, ApJ, 885, L22
  * Tacchella et al. (2021) Tacchella S., Conroy C., Faber S., Hernquist L., Eisenstein D., Johnson B., Leja J., 2021, arXiv e-prints, arXiv:2102.12494
  * Thomas et al. (2003) Thomas D., Maraston C., Bender R., 2003, MNRAS, 339, 897
  * Weinberg et al. (2017) Weinberg D. H., Andrews B. H., Freudenburg J., 2017, ApJ, 837, 183
  * Whitaker et al. (2013) Whitaker K. E., et al., 2013, ApJ, 770, L39
  * Wild et al. (2020) Wild V., et al., 2020, MNRAS, 494, 529
  * Wilkinson et al. (2021) Wilkinson A., Almaini O., Wild V., Maltby D., Hartley W. G., Simpson C., Rowlands K., 2021, MNRAS, 504, 4533
  * Williams et al. (2009) Williams R. J., Quadri R. F., Franx M., van Dokkum P., Labbé I., 2009, ApJ, 691, 1879
  * Wu et al. (2018a) Wu P.-F., et al., 2018a, ApJ, 855, 85
  * Wu et al. (2018b) Wu P.-F., et al., 2018b, ApJ, 868, 37
  * Wuyts et al. (2011) Wuyts S., et al., 2011, ApJ, 738, 106
  * York et al. (2000) York D. G., et al., 2000, AJ, 120, 1579
  * van der Wel et al. (2014) van der Wel A., et al., 2014, ApJ, 788, 28
  * van der Wel et al. (2016) van der Wel A., et al., 2016, ApJS, 223, 29
  * van der Wel et al. (2021) van der Wel A., et al., 2021, arXiv e-prints, arXiv:2108.00744