The VANDELS survey: the star-formation histories of massive quiescent galaxies at 1.0 < z < 1.3
A. C. Carnall R. J. McLure J. S. Dunlop F. Cullen D. J. McLeod V. Wild B. D. Johnson S. Appleby R. Davé R. Amorin M. Bolzonella M. Castellano A. Cimatti O. Cucciati A. Gargiulo B. Garilli F. Marchi L. Pentericci L. Pozzetti C. Schreiber M. Talia G. Zamorani


firstpage–lastpage
1ES 1927+654: An AGN Caught Changing Look on a Timescale of Months
    [
    December 30, 2023
==================================================================




We present a Bayesian full-spectral-fitting analysis of 75 massive (M_* > 10^10.3M_⊙) UVJ-selected galaxies at redshifts of 1.0 < z < 1.3, combining extremely deep rest-frame ultraviolet spectroscopy from VANDELS with multi-wavelength photometry. By the use of a sophisticated physical plus systematic uncertainties model, constructed within the  code, we place strong constraints on the star-formation histories (SFHs) of individual objects. We firstly constrain the stellar mass vs stellar age relationship, finding a steep trend towards earlier average formation time with increasing stellar mass (downsizing) of 1.48^+0.34_-0.39 Gyr per decade in mass, although this shows signs of flattening at M_* > 10^11M_⊙. We show that this is consistent with other spectroscopic studies from 0 < z < 2. This relationship places strong constraints on the AGN-feedback models used in cosmological simulations. We demonstrate that, although the relationships predicted by Simba and IllustrisTNG agree well with observations at z=0.1, they are too shallow at z=1, predicting an evolution of ≲0.5 Gyr per decade in mass. Secondly, we consider the connections between green-valley, post-starburst and quiescent galaxies, using our inferred SFH shapes and the distributions of galaxy physical properties on the UVJ diagram. The majority of our lowest-mass galaxies (M_* ∼ 10^10.5M_⊙) are consistent with formation in recent (z<2), intense starburst events, with timescales of ≲500 Myr. A second class of objects experience extended star-formation epochs before rapidly quenching, passing through both green-valley and post-starburst phases. The most massive galaxies in our sample are extreme systems: already old by z=1, they formed at z∼5 and quenched by z=3. However, we find evidence for their continued evolution through both AGN and rejuvenated star-formation activity.


galaxies: evolution – galaxies: star formation – methods: statistical 




§ INTRODUCTION


Understanding the origin of the colour bimodality in the local galaxy population remains one of the most important goals of extragalactic astronomy. The emerging picture is one in which feedback of energy from baryonic processes plays a central role in the quenching of star-formation activity, with supernovae and active galactic nuclei (AGN) thought to dominate at lower and higher masses respectively. A range of other factors are also thought to play important roles, for example mergers and environmental effects such as ram-pressure stripping. However, we still lack a detailed physical description of these and other relevant processes.

The inherent complexity of galaxy formation physics means that large-scale numerical simulations are required to connect theoretical models with observable properties. The extreme computational expense of such simulations precludes, for the present, a Bayesian statistical approach to parameter estimation and model selection in this context. Instead, predictions from individual simulations are made for a number of observable properties, such as the distributions of galaxy colours, stellar masses and star-formation rates (SFRs), which can be compared with observational results. These comparisons provide qualitative insights as to how the underlying physical model should be refined (e.g. ; ; ; ; ; ). 

However, numerical simulations have now reached a level of complexity such that a range of models with varying physical ingredients can produce good approximations of the standard set of well-constrained observable properties. The challenge for new observational studies therefore is to provide precise measurements of a wider range of physical parameters which are highly constraining on galaxy formation models, such as stellar ages, star-formation histories (SFHs), stellar and nebular metallicities and levels of dust attenuation. Crucially, these studies must be performed on large and representative samples of galaxies. They must also be extended to high redshift, in order to constrain the evolution of galaxy properties across cosmic time.

An area of particular theoretical interest in recent years has been the AGN-feedback processes which quench star-formation in the most massive galaxies (e.g. ; ; ; ). These processes should leave strong imprints, not just on the bright end of the galaxy luminosity function, but also on the star-formation histories of massive quiescent galaxies (e.g. ). Additional constraints are also available from the physical properties of galaxies transitioning between the star-forming and quiescent populations. If strong constraints can be placed on these more subtle indicators, it will be possible to begin ruling out models which are capable of matching simpler observables.

A huge literature exists on the stellar ages of massive quiescent galaxies (e.g. ; ;  ; ; ; ; ; ; ; ; ; ; ; ; ; ). However, these measurements are challenging to make for several reasons.

As stellar populations age their luminosities fall rapidly, meaning that the evidence of earlier star-formation episodes can easily be lost in the glare of younger stellar populations. Additionally, galaxy spectral energy distributions (SEDs) suffer from strong, interrelated degeneracies between different physical properties, such as the age-metallicity-dust degeneracy (e.g. ; ; ). Because of these issues, photometric data often fail to strongly constrain galaxy physical parameters, meaning the applied priors can significantly impact the results obtained (e.g. ; ). Finally, significant systematic uncertainties exist in the empirical models used to interpret observational data (e.g. ).

The situation has been improved in the local Universe by the advent of large systematic surveys designed to obtain high signal-to-noise ratio (SNR) continuum spectroscopy. These data are more strongly constraining on subtle galaxy physical parameters (e.g. ; ), however, until recently such data has been extremely scarce at higher redshifts. Additionally, interpreting spectroscopic data in a way that makes full use of the available information is challenging, both in terms of the complexity of the required models and the computational expense of fitting these models to data (see Section <ref>). The use of simplifying assumptions to reduce the complexity of the problem, for example by fixing nuisance parameters to fiducial values, typically leads to biases in derived physical parameter values and to underestimated uncertainties (e.g. ; ).

Despite these challenges, a consensus has emerged around several important results. Firstly, at fixed observed redshift, less-massive galaxies are found to have younger stellar populations than their more massive counterparts (e.g. ; ; ). This is often referred to as downsizing, or mass-accelerated evolution. Secondly, at fixed stellar mass, a trend towards lower average formation redshift is found with decreasing observed redshift. A combination of factors contribute to this, including new galaxies joining the red sequence (e.g. ; ; ), mergers (e.g. ; ), and periods of rejuvenated star-formation activity (e.g. ). 

Finally, there is considerable evidence for at least two distinct quenching mechanisms with different timescales, which change in relative importance with observed redshift (e.g. ; ; ; ; ; ). Rapid quenching, often associated with post-starburst galaxies (e.g. ), is thought to dominate at high redshift (z≳1), whereas slower quenching, associated with green-valley galaxies, is thought to dominate at lower redshift (z≲1). Whilst these three fundamental results have gained broad acceptance, precise quantitative measurements are still lacking, in particular at high redshift.

Within the last year, two new, large, high-redshift spectroscopic surveys have been completed: Lega-C (; ) and VANDELS (; ). These surveys have greatly expanded the availability of high-SNR continuum spectroscopy within the first eight billion years of cosmic history, providing new opportunities for placing strong constraints on subtle galaxy physical parameters (e.g. ; ; ).

In parallel, a new generation of spectral modelling and fitting tools has been developed (e.g. ; ; ; Johnson et al. in prep). These codes include complex, flexible physical models of the kind necessary to reproduce the properties of observed spectroscopic data, and make use of modern computational and statistical methods to fit these models to data within a fully Bayesian framework. This allows the recovery of full posterior distributions for physical parameters, meaning realistic uncertainties can be obtained, including an understanding of complex, multi-parameter degeneracies.

In this work we present the first analysis of a sample of extremely deep rest-frame near-ultraviolet spectra from VANDELS. Our targets are 75 UVJ-selected galaxies with stellar masses of log_10(M_*/M_⊙) > 10.3 at redshifts of 1.0 < z < 1.3. We analyse our spectra, in parallel with multi-wavelength photometry, using Bayesian Analysis of Galaxies for Physical Inference and Parameter EStimation (; ). The  code is used to fit a complex physical plus systematic uncertainties model, allowing us to obtain strong yet realistic constraints on the physical parameters of our target galaxies (see Section <ref>). In particular we will discuss their SFHs, quantifying the downsizing trend at z∼1 and comparing our results to predictions from cosmological simulations. We will also consider the properties of the post-starburst and green-valley galaxies in our sample, in an attempt to understand the evolutionary pathways of galaxies towards the red sequence.

The structure of this work is as follows. In Section <ref>, we introduce VANDELS, and give details of the selection of our sample. Then, in Section <ref>, we give details of the physical model we construct within  to describe our targets. In Section <ref> we discuss spectral fitting approaches, and introduce our spectroscopic plus photometric fitting methodology. We present our results in Section <ref>, discuss these results in Section <ref>, and present our conclusions in Section <ref>. All magnitudes are quoted in the AB system. For cosmological calculations we adopt Ω_M = 0.3, Ω_Λ = 0.7 and H_0 = 70 km s^-1 Mpc^-1. All times, t, are measured forwards from the beginning of the Universe. For posterior distributions we quote 50th percentile values and 16th-84th percentile ranges. We assume a <cit.> initial mass function.








§ VANDELS DATA AND SAMPLE SELECTION


VANDELS <cit.> is a large, recently completed ESO Public Spectroscopic Survey using the VIMOS instrument on the VLT at Paranal Observatory. The survey targeted 2106 high-redshift galaxies in the UKIDSS Ultra-Deep Survey (UDS) field and Chandra Deep Field South (CDFS). Whilst 87 per cent of the VANDELS targets are star-forming galaxies at z>2.4, the final 13 per cent are massive, UVJ-selected passive galaxies at 1.0 < z < 2.5. In this section we describe our data, as well as the selection of our mass complete sample of 75 objects from the VANDELS DR2 public release[<https://www.eso.org/sci/publications/announcements/sciann17139.html>]




 §.§ Photometric catalogues and parent sample


The VANDELS photometric catalogues and sample selection procedure are both described in full  in <cit.>. Here we present a brief summary of the key points relevant to this work. Both sets of VANDELS pointings in UDS and  CDFS are centred on the CANDELS fields <cit.>. Because the VIMOS field of view is larger than the areas imaged by CANDELS, the VANDELS photometric catalogues were supplemented with a variety of ground-based public imaging data. Because of this, the VANDELS sample is drawn from four different photometric catalogues, each spanning a UV-NIR observed wavelength range from ∼0.3-5 μm.

Each of these catalogues was subjected to an extensive SED fitting campaign to construct derived-parameter catalogues, including robust photometric redshifts, z_phot, stellar masses and rest-frame magnitudes. The initial VANDELS sample selection was performed using a 2015 version of these catalogues, which has since been supplemented by deeper data. The photometric data used in this work, as well as the stellar masses used in Section <ref>, come from the v1.0 internal catalogues, the final versions of which will be made public as part of the final VANDELS data release.

The VANDELS passive sample was selected by the following process. Firstly, objects were required to have 1.0 < z_phot < 2.5. Objects were then selected to have observed H-band magnitudes of H < 22.5, corresponding to stellar masses of log(M_*/ M_⊙) ≳ 10. Next, objects were selected by rest-frame UVJ colours (e.g. ; ). In order to ensure that all targets would be detected with sufficient SNR in the VIMOS spectra, a final selection criterion was applied, requiring an observed i-band magnitude of i < 25. This slightly biases the full VANDELS passive sample against the faintest and reddest objects. This process, summarised below, results in a parent sample of 812 objects.



  * 1.0 < z_phot < 2.5

  * H < 22.5

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

  * U - V > 1.2

  * V - J < 1.6

  * i < 25.





 §.§ VANDELS spectroscopic observations


The VANDELS observations are described in full in <cit.>, and so we again provide here only a brief summary of the relevant points. Of the 812 objects selected by the process detailed in Section <ref>, a random sample of 268 were assigned slits and observed. All observations were conducted using the MR grism, providing R≃600 spectroscopy spanning an observed wavelength range from λ = 4800-10000Å. Objects were observed for either 20, 40 or 80 hours, depending on their i-band magnitudes, to obtain SNRs of 15-20 per resolution element (∼10Å) in the i-band. Spectroscopic redshifts, z_spec, were measured and verified by the VANDELS team, as described in section 5 of <cit.>.

A known issue with the VANDELS spectra is a systematic drop in flux at the blue end (λ≲ 5600Å). This region typically has a low SNR in the passive spectra, but an empirical correction was derived and implemented based on the bluer star-forming sample (see section 4.1 of ). This median correction is applied to all of the VANDELS spectra, however, object-to-object variations persist at levels of up to ∼30 per cent. This calibration uncertainty is fitted as part of our model, as discussed in Section <ref>.

Because of the rapid build-up of the red sequence across our target redshift range, the passive sample is heavily weighted towards lower redshifts, with 88 per cent having z_phot < 1.5. This means that, for the vast majority of the sample, these spectra contain a full suite of rest-frame UV-optical absorption features, including MgUV, Ca H and K, the 2640Å, 2900Å and 4000Å breaks and Hδ and higher order Balmer lines, as well as the [O ii] 3727Å emission line.




 §.§ The 1.0 < z < 1.3 mass complete sample


As discussed in Section <ref>, the VANDELS passive sample is not mass complete across the whole redshift range from 1.0 < z < 2.5. Furthermore, the full suite of rest-frame UV spectral features described in Section <ref> is only available at the lower end of this redshift range. Based on these considerations, it was decided to impose additional redshift and stellar-mass limits on the sample to define a mass complete sample for which all of these features are available.

We first apply a limit of z_spec < 1.3, such that the 4000Å break falls blue-ward of the strong sky-line contamination long-ward of 9250Å. We then return to the v1.0 photometry and derived parameter catalogues and re-apply the initial VANDELS passive sample selection criteria, excepting the i-band magnitude limit (see Section <ref>). For our reduced redshift range, 98 per cent of objects with stellar masses of log(M_*/ M_⊙) > 10.3 meet the i-band limit imposed in Section <ref>. 

We therefore impose this mass limit, meaning our final sample is a random draw from a 98 per cent mass complete sample. We finally require that objects have received > 90 per cent of their final exposure time in VANDELS DR2, and have a spectroscopic redshift quality flag of 3 or 4 (corresponding to > 95 per cent probability of being correct), resulting in a sample of 75 objects. The final sample has a near-uniform distribution in redshift within our chosen limits. Due to the VANDELS observing strategy, our spectra all have similar SNRs. The distribution of 7500Å SNRs is approximately Gaussian, with a mean value of 12.5Å^-1 and standard deviation of 5.3Å^-1. Key information for each of the objects in our sample is provided in Table <ref>.








§ PHYSICAL MODEL


This section describes the physical model we construct within  to describe our targets.  provides a Python framework for self-consistently modelling the stellar, nebular, dust and intergalactic medium properties of galaxy spectra. Nebular emission is modelled by post-processing stellar templates through the Cloudy photoionization code <cit.>, using a method based on that of <cit.>.  also provides several optional dust-attenuation prescriptions, and assumes energy balance, with attenuated light re-radiated in the infrared assuming the models of <cit.>. Full details are provided in section 3.1 of <cit.>.


For all of our objects, we vary the observed redshift within a narrow range centred on the spectroscopic redshift, measured as described in section 5 of <cit.>. We impose a Gaussian prior centred on z_spec, with standard deviation, σ=0.05, and allow deviations of up to 3σ in either direction. Velocity dispersion is modelled by convolution of the spectral model with a Gaussian kernel in velocity space. We apply a logarithmic prior to velocity dispersion, σ_vel, between 40 and 400 km s^-1. A summary of the parameters and priors of our model is provided in Table <ref>. The fitting of  this model to our data is described in Section <ref>.




 §.§ Stellar population model


In this work we use the default  stellar population models, which are the 2016 updated version of the <cit.> models[<https://www.bruzual.org/ gbruzual/bc03/Updated_version_2016>] (see ). These models have been updated to include the MILES stellar spectral library <cit.>, providing ∼2.5Å resolution from 3525-7500Å, and updated stellar evolutionary tracks (; ).

We model the chemical-enrichment histories of our galaxies with a delta-function, assuming that all stars within the galaxy have the same metal content with scaled-Solar abundances. This single metallicity is varied with a logarithmic prior between -2 < log_10(Z_*/Z_⊙) < 0.4 (we define Solar metallicity, Z_⊙ = 0.02).

We parameterise the SFHs of our galaxies using the double-power-law model described by <cit.>. The SFH is described by 

    SFR(t) ∝ [(t/τ)^α + (t/τ)^-β]^-1


where α is the falling slope, β is the rising slope and τ is related to the peak time. In  <cit.> we demonstrate that this model produces unbiased estimates of the redshifts of formation and quenching for quiescent galaxies from the Mufasa simulation across a wide redshift range.



 §.§ Dust attenuation model


Recent studies at high redshift have favoured an average attenuation curve slope of n ≃ 0.7, where A_λ∝λ^-n (e.g. ; ), similar to that found by <cit.> for local galaxies. However, several studies suggest significant object-by-object variation (e.g. ; ). We therefore model dust attenuation with the modified <cit.> model described in section 3.1.4 of <cit.>. We place a Gaussian prior on n, with mean, μ=0.7 and standard deviation, σ=0.3. We set permissive lower and upper limits of 0.3 and 1.5 respectively. We likewise allow a wide range of V-band attenuations, A_V, for stellar continuum emission, adopting a uniform prior from 0 < A_V < 8. We adopt a fixed value of ϵ=2 for the ratio of attenuation between stellar birth clouds and the wider interstellar medium (ISM). We also adopt a value of 10 Myr for t_BC, the lifetime of stellar birth clouds, meaning that A_V is doubled for emission from stars formed in the last 10 Myr.




 §.§ Nebular emission model


Our UVJ selection is designed to identify galaxies with low levels of ongoing star-formation, and hence nebular emission. However, it is important for our model to be capable of reproducing the spectra of dusty star-forming galaxies which can contaminate UVJ-selected samples. We therefore implement the nebular emission model described in section 3.1.3 of <cit.>, with a fixed ionization parameter of log_10(U) = -3. Light from stars which formed more recently than the lifetime we assume for stellar birth clouds (t_BC = 10 Myr) is processed through our nebular model. The resulting nebular continuum and line emission is attenuated by twice the ISM A_V, as described in Section <ref>.




§ COMBINING OUR SPECTROSCOPIC AND PHOTOMETRIC DATA


This section describes the fitting of the physical model described in Section <ref> to the combined datasets described in Section <ref>. We begin in Section <ref> by reviewing the literature on galaxy spectral fitting. Then, in Section <ref>, we show that a simple approach, which does not allow for systematic uncertainties, fails to describe our data. In Section <ref> we construct a model for these systematics, and in Section <ref> we fit our combined physical plus systematics model to our joint datasets.




 §.§ Historical approaches to spectral fitting


Historically, galaxy spectral fitting techniques have been applied to photometric data (e.g. ), and this has remained a popular approach for several reasons. Photometric data are widely available across a wide range of wavelengths, and both random and systematic uncertainties are relatively simple to characterise (e.g. ; ). It is also far simpler to construct models for the broad-band colours of galaxies than for detailed spectral features (e.g. ; ), and this has been shown to be sufficient for estimating basic physical properties such as stellar masses (e.g. ).

Such analyses typically assume that the uncertainties on photometric fluxes are well-determined, Gaussian distributed and independent. In this case, the scatter of observed fluxes f_i about their true values follows a chi-squared distribution (where i runs over a number, N_bands, of photometric bandpasses). A physical model which is a function of some parameters, Θ, can therefore be fitted to these observations using the log-likelihood function

    ln(ℒ_phot) =  K -0.5∑_i^N_bands(f_i - m_i(Θ)σ_i)^2


where K is a constant, m_i(Θ) is the model prediction for the observed flux f_i, and σ_i is the corresponding uncertainty.

Significant failures of the above assumptions are fairly simple to identify by assessing the quality of fit, typically using the minimum reduced chi-squared value, and are rare enough that the affected objects can simply be excluded from the analysis. More-subtle failures can be modelled, for example by applying variable zero-point offsets to each band (e.g. ), or by asserting that the uncertainty be greater than some fixed fraction of the observed flux (typically 5 per cent; e.g. ; ). This prevents uncertainties from being underestimated in the high-SNR regime, where the precision of the photometric calibration dominates the error budget.

However, as described in Section <ref>, photometric data are limited in their ability to constrain more-subtle galaxy physical parameters. It has been shown that spectroscopic observations have the potential to improve this situation, however accurate spectrophotometric calibration is notoriously challenging, owing to the need to correct for a range of atmospheric and instrumental effects, such as differential atmospheric refraction, telluric contamination and characterisation of the sensitivity function of the detector. Even for comparatively well-calibrated spectra, wavelength-dependent uncertainties are known to exist, typically at levels of ∼10 per cent (e.g. ; ; ). Historically therefore, analyses of spectroscopic data have been limited to individual spectral features (commonly Lick indices), such that results are independent of spectrophotometric calibration (e.g. ; ; ). Whilst these analyses have produced many extremely valuable results, they do not make use of the full information content of spectroscopic data (e.g. ). 

More recently, with the advent of large, well calibrated spectroscopic surveys such as the Sloan Digital Sky Survey (SDSS; ), attention has shifted towards full-spectral-fitting methods, which attempt to model and fit the whole information content of spectroscopic data (e.g. ; ; ; ; ). The simplest approach to full spectral fitting is to make the same assumption of well-determined, independent Gaussian uncertainties on each spectral pixel flux f_j (where j runs over the number of pixels in the spectrum, N_pix). In this case, the log-likelihood function can again be written as

    ln(ℒ_spec) =  K -0.5∑_j^N_pix(y_j - m_j(Θ)σ_j)^2


where K is a constant, m_j(Θ) is the model prediction for the pixel flux f_j, and σ_j is the corresponding uncertainty. 

This approach has however been demonstrated to be less successful in describing spectroscopic data than photometry. <cit.> note that the quality of the fits they obtain are typically poor, and attribute this to both inadequacies in the models they fit to their data, and the difficulty of obtaining reliable uncertainties. Furthermore, in <cit.>, the authors report that improvements to the SDSS spectrophotometric calibration have significantly changed their inferred SFHs. These and other authors also note the challenges which exist in the exploration of the higher-dimensional parameter spaces of the more complex models required to fit spectroscopic data. 

As the availability of high-quality panchromatic photometric data has increased, interest has grown in the use of photometry to complement spectroscopic analyses. These joint analyses promise improvements in our understanding by taking advantage of both the broad wavelength coverage and excellent calibration of photometry, and the strong constraints on subtle physical parameters offered by spectroscopy (e.g. ; ).






 §.§ A simple approach to fitting our joint datasets


When jointly analysing several datasets, the likelihood function is the product of the separate likelihoods. The log-likelihood function for the joint analysis of spectroscopy and photometry is therefore ln(ℒ) =  ln(ℒ_phot) + ln(ℒ_spec). In this section we take a simple approach, by attempting to jointly fit our spectroscopic and photometric datasets using the physical model described in Section <ref>, and the ansatz for ln(ℒ_phot) and ln(ℒ_spec) given in Equations <ref> and <ref> respectively. This approach will henceforth be described as our simple model. The nine parameters of this model are summarised at the start of Table <ref>. The other six parameters below are only included in the final model, which is described in Sections <ref> and <ref>.

For our photometric data we employ the common methods discussed in Section <ref> for dealing with systematic calibration uncertainties. We first apply the photometric zero-point offsets calculated for the VANDELS photometric catalogues by <cit.>. We also assert that the uncertainty for each band must be ≥ 5 per cent of the observed flux, except for the two IRAC channels where a threshold of 10 per cent is used. The VANDELS DR2 spectra were binned by a factor of two to a sampling of 5Å. Fitting was carried out using the wavelength range from 5200Å < λ < 9250Å where the detector sensitivity is high and sky-line contamination is minimal compared to longer wavelengths. Sampling of the posterior was carried out using the MultiNest algorithm <cit.>.

The results of this method applied to an example object from our sample are shown in the middle panels of Fig. <ref>. It can be seen that the overall shape of the posterior distribution matches the spectroscopic data well. However, on closer inspection (and comparison to the above panel), it can be seen that the depths of individual absorption features are poorly reproduced. The most obvious failure however is in reproducing the observed photometry from 1-2 μm, where our model significantly overestimates the observed fluxes. The object shown in Fig. <ref> is typical of the sample, and these issues were observed for the majority of our galaxies.

The fact that the joint fit is incapable of matching both datasets simultaneously means that the two are inconsistent under the assumptions we have made. The joint fit adopts a region of parameter space which best describes the spectroscopic data at the expense of the photometry because there are a larger number of spectral pixels than photometric bands, and hence more terms in the log-likelihood function which depend on the spectroscopy than the photometry. Even though the joint fit is dominated by the spectroscopy, the quality of fit to the spectroscopic data is still poor, in accordance with the results of <cit.>.

Issues of this nature have been commonly observed in similar analyses, leading to suggestions that the spectroscopic data should be somehow down-weighted in the likelihood function, in order to give “equal consideration” to both datasets. However, apart from being statistically unjustified, this cannot solve the underlying issues of inconsistency between the datasets and poor quality of fit to the spectroscopic data. Instead it is necessary to understand the causes of these issues, so that these effects can be included in the model we fit to our joint datasets.

One possible reason for inconsistencies between our datasets would be aperture bias, as our photometry is measured within 2^''-diameter apertures, whilst our spectroscopic observations use a 1^'' slit. However, similar studies at lower redshifts have found aperture effects to be small (e.g. ), and we would expect a smaller effect still, due to the larger angular diameter distance to our target redshift range and the smaller physical sizes of quiescent galaxies at higher redshifts (e.g. ). As a check, we take effective radius measurements, r_e, from <cit.> for the 26 objects which have CANDELS imaging. The mean r_e is 0.35^'', and no correlation exists between  r_e and the implied inconsistency between our datasets. Aperture bias also cannot explain the poor quality of fit to the spectroscopic data.

As described in Section <ref>, there is a known issue with the VANDELS spectrophotometric calibration, which is typically low at the blue end, reddening the spectra with respect to the observed photometry. As the joint fits are dominated by the spectroscopic data, this would be expected to cause the fitted models to be redder than the observed photometry, which is consistent with the overestimation of the 1-2 μm photometry observed in the middle panels of Fig. <ref>. Incorrect spectroscopic calibration is also likely to lead to a poor quality of fit, as it is unlikely that models will exist within the parameter space being explored which can simultaneously match the perturbed spectral shape and absorption features present. We therefore conclude that systematic calibration uncertainties in our spectroscopic data are the most likely cause of the issues we identify, and move on to develop a model for these systematic uncertainties.




 §.§ Modelling spectroscopic systematic uncertainties


In general, as discussed in Section <ref>, it is extremely complex to construct a physical model for the atmospheric and instrumental effects to which both our spectroscopic data and our empirical stellar-population models are subjected. We therefore take a different approach, by constructing a flexible empirical model for systematic perturbations of our spectroscopic data about our physical model. 

We will split these perturbations by the general form they take: either additive or multiplicative (e.g. ). We will refer to these as noise and calibration offsets respectively. We then construct flexible models for these offsets by introducing nuisance parameters, Φ, into the spectroscopic log-likelihood function. We can then later marginalise these nuisance parameters out of our posterior distribution, in order to obtain a posterior for our physical parameters which includes uncertainties due to systematic effects.

We modify the log-likelihood function for spectroscopy presented in Equation <ref> as follows. Firstly we generalise our model, m(Θ) to  m(Θ, Φ) by dividing through by a multiplicative polynomial calibration model, P_j(Φ), such that 

    m_j(Θ, Φ) = m_j(Θ)P_j(Φ).


This model will be discussed in Section <ref>.

Secondly, we drop the assumption that the uncertainties on our observed spectroscopic fluxes are independent, allowing additive correlated noise between our spectral pixels. We hence replace the second term on the right of Equation <ref> with a matrix equation, in which the inverse of the covariance matrix is multiplied on both sides by the vector of residuals between our observed and model fluxes. Our covariance matrix, 𝐂(Φ), will be drawn from a Gaussian process model, which is described in Section <ref>. Our spectroscopic log-likelihood function is now

    ln(ℒ_spec) =  K - ln(|𝐂(Φ)|) - Δ^T 𝐂(Φ)^-1 Δ


where Δ = y_j - m_j(Θ, Φ) is the vector of residuals between our observed and model fluxes and K is a constant. Equation <ref> is simply a generalisation of Equation <ref>, and reduces back to Equation <ref> for the case in which P_j(Φ) = 1 and 𝐂(Φ) is diagonal with elements σ_j^2.





  §.§.§ Polynomial calibration model


A method recently adopted by several authors for addressing spectrophotometric calibration uncertainties is to perturb the spectroscopic data by a polynomial function of wavelength. One approach is to set the polynomial coefficients before fitting a physical model to the data, by comparing synthetic photometry derived from the spectrum to observed photometry in the same wavelength range (e.g. ). This  is computationally simple, however it requires a significant number of photometric observations in the spectroscopic wavelength range. Another approach is to fit the polynomial coefficients at the same time as the parameters of the physical model, which has the advantage of incorporating calibration uncertainties into the uncertainties on physical parameters, but is more computationally expensive (e.g. ; ).

We take the latter approach, parameterising P_j(Φ) in Equation <ref> with a second-order Chebyshev polynomial. Because our physical model is divided by P_j(Φ), the polynomial posterior can be thought of as the multiplicative offset which would need to be applied to our data to correct its calibration. We caution however that any issues with the spectrophotometric calibration of our model spectra will also be incorporated into this polynomial. 

The choice of a second order polynomial was made firstly due to computational constraints, and secondly because the calibration issues with the VANDELS spectra are known to be well approximated by a quadratic function of wavelength (see fig. 4 of ). A more flexible choice for this model would be a multiplicative Gaussian process model (e.g. Johnson et al. in prep), and we intend to explore this option in future work.

We apply Gaussian priors to all three polynomial coefficients with standard deviations of σ = 0.25. The prior means are μ=1 for the zero order, and μ=0 for the first and second order. This means that the prior mean and median for P_j(Φ) are equal to 1 for all wavelengths, equivalent to no change in the calibration. The maximum deviation allowed for any polynomial order is 2σ from the prior mean.




  §.§.§ Gaussian process noise model


As discussed in Section <ref>, it is common to assume that uncertainties on observational data are independently Gaussian distributed with well-known variances. There are many good reasons to suspect that these assumptions do not hold in the case of spectroscopic data. 

Firstly, it is standard practice to oversample the resolution element of the optical system by at least a factor of two, leading to local covariances between pixels, although this can be mitigated to some extent by binning adjacent pixels. Secondly, there is also good evidence that the error spectra determined from typical data reduction pipelines are underestimates of the true pixel variances (e.g. ; ). A common approach is to expand the variances for all spectral pixels by the median residual determined from an initial round of fitting (e.g. ). Finally, a range of effects from template mismatch to poor sky subtraction have the potential for introducing correlated additive offsets between the data and models being fitted (e.g. ). 

By modifying the log-likelihood function presented in Equation <ref> to that of Equation <ref>, we have relaxed the assumption of independence in our spectroscopic uncertainties. We now parameterise our covariance matrix, 𝐂(Φ), in terms of both independent (white) noise and covariant noise between pixels. We will fit these parameters alongside those of our physical and polynomial calibration models. The form we assume for our covariance matrix is

    𝐂_jk(Φ) = a^2σ_j σ_k δ_jk + b^2exp(-(λ_j - λ_k)^22l^2)


where σ_j, k are the uncertainties on our pixel fluxes, λ_j, k are the central wavelengths of our pixels, δ_jk is the Kronecker delta function, and a, b and l are free parameters to be fitted.

The first term in Equation <ref> deals with the uncorrelated noise on our data. As we suspect that our uncertainties may be underestimated, we allow their magnitude to vary by a^2, where a is assigned a logarithmic prior between 0.1 and 10 (e.g. see section 6 of ). This is similar to the iterative approaches of other authors, however the fact that this parameter is allowed to vary during fitting means that its uncertainty is propagated into the uncertainties on our physical parameters.

The second term in Equation <ref> is drawn from a Gaussian process model, and allows us to model covariant noise between our spectral pixels. Gaussian process regression is implemented in  using the George Python package (). We adopt an exponential-squared kernel and fit the normalisation b and correlation length l. We assign logarithmic priors to both of these quantities. We define b in units of the maximum flux in the observed spectrum, f_max, and allow values from 10^-4 to 1. The maximum flux is used as the unit of b such that the same range of prior values can be used for each spectrum. The mean or median flux value could also have been used. 

Similarly, we define l in units of the wavelength range covered by the spectral data, Δλ (in this case 4050Å), and allow values from 0.01 to 1. The minimum correlation length (∼40Å) was chosen to prevent the Gaussian process model from reproducing individual absorption and emission features in our spectra. Our Gaussian process model is intended to model poor sky subtraction and template mismatch between our models and data, as demonstrated with an additive polynomial by <cit.>. As currently implemented, it cannot model covariances between adjacent spectral fluxes due to oversampling of the resolution element, or resampling from an initial non-uniform wavelength sampling. A term in Equation <ref> to account for this is a possible extension to our model (e.g. ).




 §.§ Final fitting of our joint datasets
 

In this section we describe our final fitting methodology, from which all of the results presented in Section <ref> are derived. We again fit the physical model described in Section <ref>, however we now also fit the models for systematic effects introduced in Section <ref>, by exchanging Equation <ref> for Equation <ref> in our log-likelihood function. 

Our photometric data are treated in the same way as described in Section <ref>, and we use the same wavelength range and binning for our spectroscopy. In addition, we also mask several spectral regions which experience strong telluric contamination, leading to residuals such as those visible in the top left panel of Fig. <ref> at  ∼7600Å. The regions masked are 6860-6920Å, 7150-7340Å, and 7575-7725Å. Finally, we mask the rest-frame region from 3702-3752Å, containing the [O ii] emission line. This is because the excitation mechanism for low-level line emission in quiescent galaxies is still controversial, with AGN and ionization from old stars both thought to contribute (e.g. ; ; ; ). By contrast, the only mechanism which can excite [O ii] emission in our  physical model is ionization from young stars, meaning that our inferred SFRs could be biased by [O ii] emission excited by other processes. We will compare our observed [O ii] equivalent widths to our inferred specific star-formation rates (sSFRs) in Section <ref>.

Our final model has 15 free parameters, summarised in Table <ref>. Sampling our posterior distribution with MultiNest therefore requires several million evaluations of our log-likelihood function, each of which is relatively computationally expensive, in particular the inversion of the covariance matrix. Fitting each galaxy therefore requires several tens of CPU hours, limiting the scalability of this method.

The posterior distribution for our final model fitted to the object discussed in Section <ref> is shown in the bottom panels of Fig. <ref>. Both the spectroscopic and photometric data can now be seen to be well matched by our posterior distribution. In order to provide a quantitative comparison, we also show in Fig. <ref> minimum reduced chi-squared values, χ_ν^2, for both of the fits shown. In order for the comparison to be direct, the χ_ν^2 for the fit with the simple model has been divided by the maximum likelihood value of a^2 from the fit with the full model (see Equation <ref>). This means that the same variances are assumed for the spectroscopic data in both cases. When comparing to a chi-squared distribution with 794 degrees of freedom, the full model fit is statistically acceptable (24 per cent chance of a worse fit), whereas the simple model fit is strongly excluded by the data. The χ_ν^2 values shown in Fig. <ref> are typical of galaxies in our sample.

An expanded view of the bottom left panel of Fig. <ref> is shown in Fig. <ref>. Additionally, the posterior distributions for our polynomial calibration and Gaussian process noise models are shown below the observed spectrum. The posterior distribution for our polynomial calibration model can be seen to follow the expected form, with little change across most of the spectral range, and an increase of ∼30 per cent at the blue end. This is typical of the polynomial corrections we recover for objects in our sample, although the degree of correction at the blue end varies by an additional ∼30 per cent from object to object. These calibration corrections follow the expected form found by <cit.>, and eliminate the apparent inconsistencies between spectroscopic and photometric datasets discussed in Section <ref>. The corrections introduced by our Gaussian process noise model can be seen to be small, however the extra flexibility in continuum shape allows the absorption features present in the observed spectrum to be well fitted by our model.

Finally, the black line in the main panel of Fig. <ref> shows the posterior median for the physical model we fit to our observed spectrum. Assuming that systematics on our spectral models are negligible, this can be thought of as the best fit to our observational data corrected for systematic effects. For clarity, the black line in the top panel divided by the polynomial in the bottom panel then added to the Gaussian process model in the centre gives the median of the orange posterior distribution shown in the top panel.




§ RESULTS


In this section we present the results obtained by fitting our full model, as described in Sections <ref> and <ref>, to the sample described in Section <ref>.




 §.§ Quiescent and green valley sub-samples


Two main methods have been used to define samples of quiescent galaxies: selection by sSFR and selection by rest-frame UVJ colours (typically evolving with observed redshift; e.g. ). Several recent studies define quiescence by a time-evolving criterion of sSFR < 0.2/t_obs, where t_obs is the age of the Universe when the galaxy is observed (e.g. ; ). This was demonstrated by <cit.> to produce good agreement with UVJ selection using a non-evolving colour criterion of U - V > 0.88(V - J) + 0.69 from 0.25 < z < 3.75.

As detailed in Section <ref>, the VANDELS UVJ selection uses the more permissive 1.0 < z < 2.0 colour criterion from <cit.> of U - V > 0.88(V - J) + 0.49. We therefore apply a further, sSFR-based selection to our sample in order to facilitate comparisons with other recent work. We apply the slightly modified method of <cit.>, which uses the normalised SFR, or nSFR. This is defined as the SFR averaged over the most recent 100 Myr, SFR_100, as a fraction of the average SFR over the whole history of the galaxy. This can be written as

    nSFR = SFR_100t_obs∫_0^t_obsSFR(t) dt.


In fig. 5 of <cit.> we demonstrate that a selection criterion of nSFR < 0.1 produces good agreement with sSFR < 0.2/t_obs at all redshifts. This requires the ongoing SFR of the galaxy at the redshift of observation to be less than 10 per cent of its historical average. By the application of this criterion we separate our sample into 53 quiescent and 22 green-valley galaxies. This approach will be compared to UVJ-based selection criteria in Section <ref>.










 §.§ Stellar mass vs formation redshift


From the star-formation histories we infer for our quiescent sub-sample we measure several quantities. Firstly, we calculate the mean time (measured forwards from the Big Bang) at which the stars in our galaxies formed, t_form, given by

    t_form = ∫_0^t_obs t SFR(t) dt/∫_0^t_obsSFR(t) dt.


This corresponds to the mean stellar age, or mass-weighted age, and is similar to the median formation time, t_50, used in some studies. We then calculate the redshift corresponding to t_form, which we call the formation redshift.

We also calculate the history of the nSFR parameter over our inferred SFHs and extract the age of the Universe at which nSFR first falls below 0.1. This is the time at which the galaxy would first enter our quiescent sub-sample, and we therefore refer to this as the time of quenching, t_quench.

Times of formation and quenching inferred using the double-power-law SFH model described in Section <ref> were extensively validated using mock photometric observations of simulated quiescent galaxies by <cit.>. However, subsequent evidence (; Iyer et al. in prep.) suggests that this SFH model under-predicts the ages of star-forming galaxies. We therefore do not report inferred ages for our green-valley sub-sample.

The times of formation and quenching we infer for our quiescent sub-sample are shown in Fig. <ref> as a function of our inferred stellar masses. It can be seen that almost all of our lower-mass galaxies (M_* ≲ 10^11M_⊙) are found to have formed the bulk of their stellar populations at z < 3. By contrast, more-massive galaxies display a greater spread in formation redshift, with the oldest objects found to have formed their stars at z ∼ 5 and quenched by z=3.

In order to understand the average properties of our sample as a function of stellar mass, we split our objects into four mass bins. A 0.25 dex spacing was used between log_10(M_*/M_⊙) = 10.75 and 11.25, with additional bins for objects with M_* < 10^10.75M_⊙ and M_* > 10^11.25M_⊙. The mean formation redshift in each bin is shown with black errorbars in the top panel of Fig. <ref>. A clear trend is visible at lower masses, which shows some signs of flattening at M_* > 10^11M_⊙, as found at z=0.7 by <cit.>.

Following <cit.> we fit a linear relationship to describe the average t_form for our galaxies as a function of stellar mass, including an intrinsic scatter to account for effects unrelated to the stellar masses of our objects (e.g. galaxy environment). For the mean relationship we find

    (t_formGyr) = 2.56^+0.12_-0.10 - 1.48^+0.34_-0.39 log_10(M_*10^11M_⊙)


with an intrinsic scatter of 0.58^+0.09_-0.08 Gyr. The posterior for the mean relationship is shown in blue in the top panel of Fig. <ref>. This result will be discussed in Section <ref>. 

It is interesting to consider the distribution of SFRs for the objects shown in Fig. <ref>. For all objects with nSFR < 0.01 (those coloured dark red) the lower bound on the current level of star-formation is zero: we hence describe these objects as having no detectable star-formation. Galaxies with lighter colours have detectable star-formation at levels low enough that we still describe them as quiescent. It can be seen that no trend exists between t_form and nSFR in Fig. <ref>. This supports the findings of <cit.>, who attribute star-formation in quenched galaxies to stochastic processes such as minor mergers and rejuvenation events. 

A trend is visible between t_quench and nSFR, with more recently quenched galaxies having higher nSFR. However this result is not robust, as our double-power-law SFH model cannot reproduce rejuvenation events. Once star-formation drops to near zero it cannot rise again under this model, meaning quenching must be delayed until recent times to match any level of ongoing star formation. We hence conclude that we cannot reliably measure t_quench with our double-power-law model from UV spectroscopy for galaxies with detectable ongoing star-formation. A more advanced approach will be needed to model the details of these SFHs (e.g. ; ; ).







 §.§ Distribution in Dn4000 vs Hδ


Historically, a common method for inferring galaxy ages and sSFRs from UV-optical spectroscopy has been to 
measure the strengths of the 4000Å break (Dn4000) and the Balmer delta (Hδ) absorption feature (e.g. ; ). We therefore report these spectral indices for our galaxies, both to show the distribution of these parameters within the quiescent population at 1.0 < z < 1.3, and to check that the results we infer from our full-spectral-fitting method are in agreement with the expected relationships between these parameters. 

We measure Dn4000 from our spectra as the ratio of average fluxes between 3850-3950Å and 4000-4100Å, whilst masking out pixels which deviate by more than 3σ from the posterior median model fitted in Section <ref>, which are assumed to experience significant sky-line contamination (31 objects have pixels masked, with the average fraction masked being 2 per cent). We measure the rest-frame Hδ equivalent width, EW_Hδ, by fitting a first-order polynomial plus Gaussian model to the 100Å spectral region centred on Hδ. We then correct for nebular emission using the Hδ flux predicted by our fitted  model (see Section <ref>). These corrections are small (≲1Å) for our green valley sub-sample, and typically negligible for our quiescent sub-sample. For 8 of our 75 objects, Hδ falls within the observed wavelength ranges masked in Section <ref>, and is therefore strongly contaminated by sky lines, such that no measurement could be made. For these objects we apply the same fitting methodology to the posterior prediction for this spectral region from our fitted  model.

Our sample is shown on the EW_Hδ vs Dn4000 plane in the left panel of Fig. <ref>. Our sub-samples are significantly offset, with our green valley objects having stronger Hδ absorption and a weaker 4000Å break. There is a clear transition at 1.3 < Dn4000 < 1.4 and EW_Hδ∼4Å, with almost all of our quiescent sub-sample at higher Dn4000 and lower EW_Hδ. Four significant outliers are visible: three quiescent objects with low Dn4000 and strong Hδ absorption, and one green-valley object with a strong 4000Å break (Dn4000 ∼ 1.55). These objects will be discussed in Section <ref>.

In the right panel of Fig. <ref> our quiescent sub-sample is shown coloured by our inferred mass-weighted ages (see Section <ref>). A trend in age with Dn4000 is visible as expected, with our oldest objects having Dn4000 ≃ 1.7. Contours are plotted showing the distributions of quiescent galaxies at lower observed redshifts. Both samples are taken from <cit.> with the additional imposition of our nSFR < 0.1 criterion. The SDSS sample, shown in gray, is at 0.04 < z < 0.14; the Lega-C sample, shown in blue, is at 0.6 < z < 1.0. A similar evolution  of ∼0.2 in the average Dn4000 value can be seen from z∼0.1 to z∼0.8 and from z∼0.8 to our sample at an average redshift of z∼1.15. The cosmic time interval between SDSS and Lega-C is approximately four times that between Lega-C and VANDELS, demonstrating that the distribution of quiescent galaxies moves towards lower Dn4000 at an accelerating pace with increasing lookback time as expected.





 §.§ Trends with rest-frame UVJ colours


In this section we consider the positions of our galaxies on the UVJ diagram, and physical parameter trends with UVJ colours. Our sample is shown on the UVJ diagram in Fig. <ref>. 



  §.§.§ Trends with star-formation rate


The top-left panel of Fig. <ref> shows a direct comparison of our nSFR-based selection to both sets of UVJ selection criteria discussed in Section <ref>. Good agreement can be seen between the dashed criterion of U - V > 0.88(V - J) + 0.69 and the nSFR-based selection introduced by <cit.>. None of our quiescent sub-sample fall below the dashed line, whereas seven green valley objects fall above. All seven are very close to our nSFR threshold, with the majority having 16th posterior percentiles at nSFR < 0.1.

The top-right panel of Fig. <ref> shows our sample coloured by nSFR. At an observed redshift of z∼1, a nSFR value of 0.1 is equivalent to log_10(sSFR/yr^-1) ≃ -10.5. We observe a clear trend in nSFR (and hence sSFR) perpendicular to the red sequence, in agreement with recent results (e.g. ; ). It is worth noting that the solid box representing the more permissive UVJ selection criteria includes objects which are forming stars at up to ∼50-100 per cent of their historical average SFRs, meaning that these criteria alone should not be used to select high-redshift quiescent galaxies.

We therefore argue that the nSFR < 0.1 criterion is the most robust method for selecting quiescent galaxy samples. This method is truly redshift-independent, as it selects galaxies which are forming stars below a fixed fraction of their historical average SFRs. Both a fixed UVJ selection of U - V > 0.88(V - J) + 0.69 and sSFR selection proportional to t^-1 produce results similar to nSFR < 0.1. By contrast, we argue that both the original redshift-dependent UVJ selection criteria, and selection using a fixed sSFR threshold are less appropriate, as they include galaxies which are proportionally more highly star-forming  at higher redshifts.



  §.§.§ Trends with dust attenuation


Our sample is shown coloured by A_V in the bottom-left panel of Fig. <ref>. It can be seen that we find a strong trend in dust attenuation across our UVJ box, which follows the trend we observe in sSFR, perpendicular to the dust reddening vector. These results are in accordance with those of <cit.>, who observe a drop in the dust content of their star-forming galaxies as they move closer to the red sequence. We find that the majority of objects within the solid UVJ selection box which have A_V > 1 are those we identify as in the green valley, rather than truly quiescent. This confirms that, within the UVJ box, dust attenuation is more strongly related to sSFR than position along the dust-reddening vector. Our quiescent sub-sample is typically found to be less dusty, however there is still a noticeable trend with distance from the edge of the UVJ selection box, as was found by <cit.>.



  §.§.§ Trends with stellar age


The mass-weighted ages of our quiescent sub-sample (see Section <ref>) are shown on the UVJ diagram in the bottom-right panel of Fig. <ref>. Also plotted are ages predicted by the relationship derived by <cit.>. The lines which run perpendicular to the dashed UVJ selection are lines of constant age, which are shown from 0.5-4.5 Gyr in intervals of 0.5 Gyr. The ages we derive for our objects are in good agreement with <cit.>, despite several methodological differences (median vs mean stellar age, different SFH models), demonstrating that the ages of quiescent galaxies are less model-dependent than those of star-forming objects.

It is remarkable that such a clear trend in stellar age can exist parallel to the dust reddening vector on the UVJ diagram. Galaxies along the top-left edge of our distribution can be seen to follow a pure age sequence, with no evolution in dust attenuation. A population of totally quenched objects with younger stellar populations which still retained significantly more dust would disrupt this, as well as the age trend found by <cit.>.  This implies that quenching galaxies must lose most of their dust before their sSFRs drop to the extent where they can join the red sequence. These ideas will be explored further in Section <ref>.




 §.§ [O ii] emission properties


As described in Section <ref>, the [O ii] 3727Å emission line was masked during our spectral fitting analysis due to uncertainties as to the excitation mechanism in quiescent galaxies. We now consider the distribution of [O ii] emission in our sample. We first measure the rest-frame equivalent width of the line, EW_[O ii], using the same method as was applied to the Hδ feature in Section <ref>. 

A comparison of our measured [O ii] equivalent widths with our inferred sSFRs is shown in the left panel of Fig. <ref>. Our green-valley galaxies typically exhibit stronger [O ii] emission, with 77 per cent having EW_[O ii] < -5Å. By contrast, only 26 per cent of our quiescent sub-sample has [O ii] emission stronger than this threshold. This is despite the dustier nature of our green valley sub-sample (see Section <ref>). Whilst our green-valley sub-sample clearly exhibits stronger [O ii] than our quiescent sub-sample, it can be seen that at fixed EW_[O ii] our inferred sSFRs span a range as large as ∼2 dex. This is in agreement with previous results which identify [O ii] as a poor predictor of ongoing star-formation in quiescent galaxies (e.g. ).

Without rest-frame optical spectroscopy it is challenging to constrain possible AGN contributions to our [O ii] fluxes. However, the majority of local quiescent galaxies with detectable [O ii] emission have been shown to exhibit high [O ii]/Hα ratios (e.g. ) which are inconsistent with excitation by ongoing star-formation <cit.>. We therefore identify galaxies in our sample with strong [O ii] emission as likely hosts of low-level AGN activity, in particular the quiescent galaxies for which we find EW_[O ii] < -10Å. We will further consider the possibility of AGN activity in our galaxies in Section <ref>.

Ionization by hot low-mass stars has also been postulated as an explanation for line emission in quiescent galaxies (e.g. ). Recently, <cit.> reported a marginal difference in the stellar ages of local quiescent galaxies with and without visible emission lines. In principle, this measurement should be easier to make at z∼1, as the stellar populations of quiescent galaxies are considerably younger. Our quiescent sub-sample is shown coloured by stellar age in the right panel of Fig. <ref>. No clear correlation between [O ii] emission and stellar age is visible, however our sample is considerably smaller than those available in the local Universe, which may preclude the detection of this subtle effect. Future large high-redshift spectroscopic surveys will be a valuable tool for addressing this issue.



 §.§ Post-starburst and rejuvenated galaxies


Post-starburst galaxies (PSBs) are widely identified as one of two major transitional states between the star-forming population and red sequence (see Section <ref>). A variety of methods have been used to identify samples of PSBs, ranging from spectroscopic selection based on strong Hδ or Hβ absorption and a lack of emission lines (e.g. ), to principal component analyses (e.g. ), to selection by rest-frame UVJ magnitudes (e.g. ). Comparisons between different methods have found significant overlap (e.g. ), however a fully self-consistent set of criteria which returns objects with the desired properties is still to be agreed upon.

A detailed discussion of the physical properties of VANDELS PSBs will be presented by Wild et al. in prep. In this section we briefly discuss the distribution of PSBs within our sample, and the degree of consistency between different selection methods. The right panel of Fig. <ref> shows the EW_[O ii] vs EW_Hδ parameter space often used to spectroscopically select PSBs (e.g. ; ). We identify three spectroscopic PSBs, all of which are members of our quiescent sub-sample. These objects are marked with stars in Figs <ref>, <ref>, <ref> and <ref>.

Our spectroscopically identified PSBs are shown on the UVJ diagram in the top-right panel of Fig. <ref>. It can be seen that our PSBs occupy the region predicted by <cit.>, towards the bottom-left of the UVJ selection box. The larger of the two blue boxes is the PSB selection used by <cit.>, which assumes that PSB features are visible for quiescent galaxies with median stellar ages of 300-800 Myr. Both of the objects we find within this box are part of our green-valley sub-sample and, whilst both are close to our PSB selection box in Fig. <ref>, neither fulfils our PSB selection criteria. The object close to the top of the box has strong Hδ absorption, however it has EW_[O ii] = -6.8±0.3Å. By contrast, the object at the bottom of the box fulfils our [O ii] criterion but has EW_Hδ∼3Å. 

The objects we identify spectroscopically as PSBs occupy a region slightly further up the red sequence, within the region for which <cit.> predict ages of 800-1200 Myr. This region is highlighted with a smaller blue box above the one used by <cit.>. Our results suggest this is the region in which spectroscopic features usually associated with PSBs are strongest. The fourth object within our extended box, as well as the closest object above this box, have [O ii] consistent with our selection criteria but slightly weaker Hδ absorption. These results are in good agreement with <cit.>, who find that ∼50 per cent of galaxies identified by the photometric selection proposed by <cit.>, upon which the <cit.> UVJ selection is based, exhibit [O ii] and Hδ equivalent widths consistent with the spectroscopic criteria shown in Fig. <ref>.

Our findings suggest that the contribution of the PSB quenching channel to the growth of the red sequence is towards the upper end of the range found by <cit.>. However, as will be discussed in Section <ref>, the timescale over which PSB features are visible is not necessarily determined by the time taken for a galaxy to traverse the blue box on Fig. <ref> by passive evolution of its stellar population.

We finally note, based upon Fig. <ref>, the presence of one galaxy in our green valley sub-sample which displays Dn4000 > 1.5, but which has a blue continuum below this wavelength, and strong [O ii] emission. For this object we find a stellar mass of log_10(M_*/M_⊙) = 11.25±0.15 and a mass-weighted age of 3.6±0.5 Gyr, amongst the oldest and most massive in our sample. With currently available data it is challenging to discriminate between rejuvenated star-formation and AGN activity (though we find no evidence of AGN activity in either X-ray or radio datasets in Section <ref>). However, in either case, objects such as this are clearly of significant interest for assisting our understanding of continuing mass assembly in the oldest galaxies (e.g. ; ). Assuming no AGN contribution we infer a SFR of 8.9^+4.9_-3.2 M_⊙ yr^-1, meaning the stellar mass of this object would increase by ∼2-5 per cent over 100 Myr.









 §.§ Evidence of AGN activity


As discussed in Section <ref>, it is challenging to determine whether or not galaxies in our sample host an AGN, given only our UV-NIR photometry and rest-frame UV spectroscopy. Whilst line ratios from rest-frame optical spectroscopy would be the ideal tool for studying this, we here consider ancillary datasets in the X-ray, mid-infrared and radio to attempt to constrain any AGN contributions.

We first consider X-ray data from the Chandra seven megasecond source catalogue <cit.> and Subaru/XMM-Newton deep survey <cit.>. We find two matches within our sample, both from the Chandra seven megasecond source catalogue. These objects are shown as hexagons in the left panel of Fig. <ref>, and can be seen to have two of the four strongest [O ii] lines in our sample (the second highest is in the UDS). This suggests that strong [O ii] emission is a good predictor of X-ray AGN activity in quiescent galaxies.

We then consider radio data at 1.4GHz from the Very Large Array, which is available for both of our fields (; ). We find that three of our galaxies are detected, none of which are in common with our X-ray detected objects. All three are massive, with log_10(M_*/M_⊙) > 11, and all three sit close to our nSFR threshold, with log_10(sSFR/yr^-1) ≃ -10.5. There is no indication from their rest-frame UV continua or [O ii] lines that these galaxies host AGN, as has typically been found for radio galaxies at high redshift (e.g. ).

We finally consider publicly available Spitzer MIPS 24μm imaging, which is available for 57 of our galaxies (; ). For each of these objects we extract fluxes within 14^''-diameter apertures and manually inspect the images, flagging objects as isolated (33 objects) or potentially confused (24 objects). Of our isolated sources, 21 are members of our quiescent sub-sample, and 12 are members of our green valley sub-sample.

Only one of the two objects identified above as X-ray sources has MIPS 24μm coverage: the galaxy in our green valley sub-sample. This object has the strongest isolated detection in our sample by approximately a factor of two. Of our 21 isolated quiescent galaxies, only one has a detection at >3σ, whereas 7 out of 12 isolated green valley objects are detected above this threshold. For our eight objects with robust, isolated detections we use the calibration of <cit.> to convert our aperture-corrected fluxes into SFRs. We find that the X-ray detected source is the only object with significantly more 24μm-flux than expected, based on the SFRs we derive from our spectral fitting analysis.

The quiescent object which is detected at 24μm is one of the three identified as radio AGN above. This object is the most massive in our sample, with log_10(M_*/M_⊙) = 11.66±0.14, and one of the oldest, with a mass-weighted age of 3.5±0.6 Gyr. However, we find no evidence of an AGN contribution to its 24μm flux: our inferred SFR is 13.2^+5.6_-3.8 M_⊙ yr^-1, whereas its 24μm-predicted SFR is 14.3±2.5 M_⊙ yr^-1. This object is similar to the rejuvenated object discussed in Section <ref>, and seems consistent with the scenario proposed by <cit.>, in which cooling of gas in the hot halo begins to provide fuel for jet-mode AGN activity and rejuvenated star-formation ≳2 Gyr after quenching has occurred. 

The two 1.4GHz catalogues we consider jointly provide coverage down to 100μJy sensitivity for all our objects. This corresponds to a rest-frame 1.4GHz luminosity of 10^23.6 W Hz^-1 at z=1, or 10^23.9 W Hz^-1 at z=1.3. We find that 2 of our 4 objects with log_10(M_*/M_⊙) > 11.5 are detected above this level, whereas only 1 of our 39 objects  with 11.0 < log_10(M_*/M_⊙) < 11.5 is detected. This is broadly consistent with the local relationship between radio-loud AGN fraction and stellar mass (e.g. ; ), in agreement with other studies which find little redshift evolution (e.g. ; ).




 §.§ Stacking analyses


We finally perform a stacking analysis to demonstrate the average spectral properties of the galaxies in our two sub-samples. We de-redshift our spectra and normalise over the rest-frame wavelength range from 3200-3600Å, which is in the centre of our observed spectral range and contains no strong features. We then resample our spectra to a common wavelength grid using Spectres <cit.>. Median stacked green valley and quiescent spectra are shown in the top and bottom panels of Fig. <ref> respectively. The stacks have total exposure times of 745 and 1971 hours respectively.

A clear transition from Balmer to 4000Å break can be seen between the two spectra, demonstrating the older stellar populations of our quiescent galaxies. As discussed in Section <ref>, the median [O ii] flux is significantly higher in our green valley sub-sample. The continuum below 4000Å can also be seen to be bluer, with considerably more flux at <3000Å, indicating higher levels of ongoing star-formation.

Median stacked 24μm images are also shown in Fig. <ref>, using all 57 objects for which MIPS coverage is available (see Section <ref>). Our green valley stack shows a clear detection, with no corresponding detection in our quiescent stack.




§ DISCUSSION


In this section we further discuss our results. In Section <ref> we compare our stellar mass vs age relationship to results from the literature and recent simulations. In Section <ref>, we consider the connections between our quiescent, green-valley and post-starburst galaxies. In Section <ref> we discuss the shapes of the SFHs we infer for our quiescent sub-sample.



 §.§ The stellar mass vs stellar age relationship


As discussed in Section <ref>, the epoch of formation as a function of galaxy stellar mass is a key observable property which is strongly constraining on AGN-feedback models. In this section we compare our results to other observational studies, as well as making comparisons with simulations. 



  §.§.§ Stellar mass vs age from spectroscopic studies


The stellar mass vs stellar age relationship we derive for our quiescent sub-sample is compared to results from the literature in Fig. <ref>. In the left panel we show the mean relationship we derive between stellar mass and mass-weighted formation time for our quiescent sub-sample in Section <ref>, along with the positions of our individual quiescent galaxies.

In order to facilitate comparisons with earlier work, we also calculate the same mean relationship using r-band light-weighted formation times, t_r-band (e.g. ). For this quantity we obtain a mean relationship of

    (t_r-bandGyr) = 2.91^+0.08_-0.09 - 1.24^+0.27_-0.30 log_10(M_*10^11M_⊙)


with an intrinsic scatter of 0.51^+0.08_-0.07 Gyr. The slope of this relationship is slightly shallower than that which we obtain for t_form (Equation <ref>), and is offset towards later formation times by ∼350 Myr at log_10(M_*/M_⊙) = 11.

In the left panel of Fig. <ref>, we also show the 23 galaxies studied by <cit.> at observed redshifts of 1.5 < z < 2.5. We have shown in Section <ref> that our inferred ages agree well with the predictions of <cit.>, meaning that no significant systematic offsets should exist between the two samples. A similar range of formation redshifts can be seen for both samples, with the <cit.> objects appearing slightly offset towards earlier formation redshifts on average. The slope of the relationship can be seen to be similar in both studies.

In the right panel our r-band light-weighted formation times are compared to those found by <cit.> for quiescent galaxies at z∼0.1 and z∼0.7 respectively. Again, good agreement can be seen between the slope we infer for our relationship and those found by <cit.>. We thus conclude that, at log_10(M_*/M_⊙) > 10.3, an evolution of ∼1.5 Gyr in formation time per decade in stellar mass is a robust result which remains consistent from the local Universe to at least z∼2.

However, the relative vertical positions of our relationship and the <cit.> relationship at z∼0.7 do not follow the expected trend of earlier formation with increasing observed redshift at fixed stellar mass. This is likely to be due to the significant methodological differences between the two studies. <cit.> fit indices derived from rest-frame optical spectra and use an exponentially declining SFH model, whereas we apply a full-spectral-fitting approach to rest-frame UV spectra using a double-power-law SFH model. The use of different SFH models and priors is known to affect derived ages, which could plausibly give rise to this effect (e.g. ; ). It is additionally likely that the optical spectra of quiescent galaxies are dominated by an older population than dominates in the UV, meaning that fitting different spectral regions could result in different SFHs. 






  §.§.§ Stellar mass vs age in cosmological simulations


Given our conclusion in Section <ref> that the steep stellar-mass vs stellar age relationship observed for quiescent galaxies at low redshift is already in place by z∼2, it is interesting to consider whether this trend is reproduced by modern cosmological simulations. Historically, this relationship has been challenging to match, even in the local Universe, for both quiescent and star-forming galaxies (e.g. ; ).

We consider predictions from the 100 h^-1 Mpc box length runs of both Simba <cit.> and IllustrisTNG (e.g. ), using snapshots at z=0.1 and 1.0. In order to match these predictions as closely as possible with the observational studies discussed in Section <ref>, we apply apertures to the simulated galaxies as follows. For the z=0.1 snapshots, we apply 3^''-diameter circular apertures to each galaxy, for consistency with the SDSS observations used by <cit.>. For the z=1 snapshots, we apply 1^'' square apertures, to match the 1^'' slits used by <cit.>, <cit.> and this work (assuming a ∼1^'' region along the direction of the slit is extracted). We select all galaxies within these snapshots which meet our nSFR < 0.1 criterion (see Section <ref>), then calculate mass-weighted and r-band light-weighted formation times for the simulated galaxies. We use  to predict the r-band flux from each star particle.

In Fig. <ref> we show the average mass-weighted and r-band light-weighted formation times as a function of stellar mass for each snapshot. At each point along the horizontal axis, the median formation time of galaxies within a 0.25 dex mass range centred on that point is shown. The two simulations can be seen to be in reasonably good agreement, with discrepancies confined to levels of ≲250 Myr. The z=0.1 relationships predicted by these simulations have slopes consistent with the ∼1.5 Gyr per decade in mass evolution found by the observational studies shown. Additionally, the normalisations of the z=0.1 relationships in both simulations are in good agreement with the results of <cit.>, as found for the general z<0.2 galaxy population in IllustrisTNG by <cit.>.

However, at z=1 both simulations predict significantly weaker average stellar mass vs age relationships. Using mass-weighted formation times, Simba predicts an evolution of 0.38 Gyr per decade in stellar mass across the interval shown in Fig. <ref>, whereas IllustrisTNG predicts an evolution of 0.51 Gyr per decade. This suggests that these simulations do not accurately reproduce the detailed physical properties of massive quiescent galaxies at z>1. This is particularly interesting in the context of the results of <cit.>, who find that the precursors to these simulations (Mufasa and Illustris) significantly under-predict the number density of quiescent galaxies at 3<z<4.



  §.§.§ Stellar mass vs age from photometric studies


Several recent studies also attempt to probe the stellar mass vs age relationship using photometric data. This has the advantage of providing better statistics, as larger samples are available. However, as discussed in Section <ref>, the age-metallicity-dust degeneracy leads to larger uncertainties on individual age measurements. In this section we compare two recent photometric studies to the spectroscopic studies discussed in Section <ref>. 

<cit.> consider a sample of 845 objects with multi-band photometry from CANDELS, whereas <cit.> consider a sample of 9289 galaxies from UltraVISTA <cit.>. The results of these two studies are summarised in Fig. <ref>. We derive t_form for the stacked SFHs of <cit.> using the best-fitting double-power-law parameters reported in their table A1. Despite their very different methodologies, these studies can be seen to produce similar results. However, the slopes reported are considerably shallower than those shown in Fig. <ref>, at ≲0.5 Gyr per decade in stellar mass.

We attribute this inconsistency to the increased uncertainties on individual object SFHs when considering photometric data. In the presence of large uncertainties, a population of very massive, very old objects will be preferentially scattered towards later formation times.  This is due to the constraint imposed by the age of the Universe. Similarly, a population of younger, less massive quiescent galaxies will be preferentially scattered towards earlier formation times, as the time-evolution of galaxy stellar populations is much more rapid at younger ages, meaning later formation times will be strongly inconsistent with the observed spectrum.

Both of these effects act to flatten the trends observed in Fig. <ref> when compared to the spectroscopic analyses in Fig. <ref>. This finding highlights the importance of upcoming large spectroscopic surveys in providing stronger constraints on stellar ages than are available from photometric data alone (e.g. ; ).







 §.§ Connecting green-valley, post-starburst and quiescent galaxies


Much debate exists as to how galaxies evolve away from an initial position within the star-forming population, potentially through green-valley or post-starburst phases, to eventually join the red sequence. Several recent studies have considered this problem in the context of evolutionary tracks across the UVJ diagram (e.g. ; ). One of the most challenging aspects of this approach is modelling the time-evolution of dust attenuation, which has a significant impact on these evolutionary tracks.

<cit.> consider the number of galaxies which pass through the larger blue PSB selection box drawn on the top-right panel of Fig. <ref>, arguing for separate fast and slow quenching mechanisms which do and do not pass through the box respectively. They find that fast quenching plays a more important role at high redshift, in accordance with previous work (e.g. ; ; ; ). A toy model to describe both fast and slow quenching routes is presented in their fig. 12, assuming that A_V∝ SFR. For their exponentially declining SFH models this means galaxies lose most of their dust early on, and the tracks therefore approach the UVJ selection box from the bottom-left part of the diagram. 

However, for our green-valley sub-sample we find relatively high dust attenuation (A_V∼1-2) even very close to the dashed UVJ boundary in Fig. <ref>. This implies that these objects evolve differently to the scenarios proposed by <cit.>, approaching the UVJ selection box from further upwards and to the right on the UVJ diagram. This is expected, given their high masses and the fact that A_V is more strongly correlated with stellar mass than SFR in star-forming galaxies (e.g. ; ). In this section we discuss quenching scenarios for green-valley galaxies.





  §.§.§ Modelling the evolution of UVJ colours


We begin by considering the time-evolution of dust attenuation. As discussed in Section <ref>, we observe strong trends in both nSFR and A_V perpendicular to the dashed UVJ selection in Fig. <ref>, meaning nSFR and A_V are correlated. This relationship is shown in Fig. <ref> for galaxies with detectable levels of star-formation (nSFR > 0.01; see Section <ref>). We fit a linear relationship as described in Section <ref>, obtaining

    A_V = 0.52^+0.12_-0.11 log_10(nSFR) + 1.61^+0.15_-0.13


with an intrinsic scatter of 0.37^+0.06_-0.05 magnitudes. At lower nSFR we find that A_V reaches a minimum value of ∼ 0.5 as discussed in Section <ref>. At higher nSFR we do not have galaxies in our sample to test whether this relationship holds. In the discussion which follows we use a simple empirical ansatz for A_V, consistent with Equation <ref>. We assume that A_V = 0.5 log_10(nSFR) + 1.5 for 0.01 < nSFR < 1, then that A_V remains constant at 0.5 for lower nSFR (shown as a black line in Fig. <ref>). For simplicity we assume the <cit.> dust attenuation law, with attenuation doubled towards stars formed in the last 10 Myr, as in Section <ref>.

Armed with this relationship, we construct four representative double-power-law SFHs (see Equation <ref>) to describe a range of quenching scenarios. We consider “early quenching”, which occurs at z∼2, and “late quenching” at z∼1. For each of these two scenarios we consider fast and slow quenching paths, which have timescales of ∼100 Myr and ∼1 Gyr respectively, in accordance with <cit.>. The parameters of the four models are



  * Early fast quench: τ=3 Gyr, β=0.5, α=100

  * Early slow quench: τ=3 Gyr, β=0.5, α=10

  * Late fast quench: τ=5 Gyr, β=0.5, α=100

  * Late slow quench: τ=5 Gyr, β=0.5, α=10.


We model the evolution of the UVJ colours of these four mock galaxies forwards from the point at which their nSFR first falls below 1. this can be thought of as the time at which they enter the green valley, as this is roughly the highest nSFR we find in our green valley sub-sample (see Fig. <ref>).





  §.§.§ Relating the green-valley, post-starburst and quiescent populations



The tracks the four mock galaxies follow across the UVJ diagram are shown in Fig. <ref>. Their SFHs are shown in the inset panels. We highlight with green and red points the redshifts at which nSFR falls below 1 and 0.1 respectively, corresponding to the times at which these galaxies enter the green valley and then quench. It can be seen that each of these objects at first follows a curving path, the shape of which is determined by both quenching speed and the duration of star-formation activity before quenching begins (as well as depending strongly on the assumptions we have made regarding dust attenuation). At the end of this curving track, galaxies begin to passively age along a straight path towards the upper right of the UVJ box. The “late slow quench" model has not reached this point by z=0.5.

Galaxies which quench rapidly can be seen to briefly pass through the blue PSB selection box, both entering and leaving by the top-right edge. This picture is distinct from the PSB evolution channel of <cit.>, which enters the PSB box from the lower-left edge. This is a consequence of different assumptions regarding the evolution of dust: we assume that galaxies lose their dust rapidly during quenching, whereas in the <cit.> scenario dust attenuation is already low on approach to the green valley. How far our rapidly quenched galaxies enter into the PSB box depends on how extended in time star-formation activity is prior to quenching, as well as quenching speed. At high redshift, when star-formation cannot be very extended, galaxies fall further into the box, and spend more time in the PSB phase. At lower redshifts, the older stellar populations present in galaxies prevent their UVJ colours becoming blue enough to fall as deep into the PSB box.






We interpret these findings in the context of the results of <cit.>, <cit.> and <cit.>, who suggest a dual origin for the post-starburst population. High-redshift PSBs primarily follow a UVJ evolution similar to that suggested by <cit.>. They form in a brief, intense starburst, as suggested by <cit.>, but either never build up significant dust or shed their dust rapidly as star-formation begins to fall. This could be a consequence of strong AGN-driven outflows at high-redshift (e.g. ; ), which have been linked to z∼1 PSBs (). These galaxies enter the PSB box by the bottom-left (or possibly bottom-right) edge, and exit by the top-right. The morphological evidence presented by <cit.> and <cit.> supports this picture, with high-redshift PSBs found to be extremely compact, consistent with a major disruption event such as a merger.

By contrast, lower-redshift PSBs primarily form by the rapid quenching of normal star-forming galaxies <cit.>, which shed their dust more slowly, following a UVJ evolution similar to that shown in the bottom-left panel of Fig. <ref>. These objects briefly enter and leave the PSB box by the top-right edge. This is again consistent with <cit.>, who find that these low-redshift PSBs have less-concentrated structures, more similar to ordinary star-forming galaxies. At progressively lower redshifts, the progressively older stellar populations present in these objects prevent them from falling as far into the PSB box, explaining the decreasing number density of PSBs with redshift. This also explains the clustering of galaxies close to the top-right edge of the PSB box (e.g. see fig. 10 of ), which becomes more pronounced at lower redshifts.

A final piece of evidence for this scheme comes from the individual SFHs we infer for the three spectroscopic PSBs identified in Fig. <ref>. The posterior SFHs we infer for these objects are shown in Fig. <ref>. For two objects the SFHs we infer are very extended before rapidly quenching. These objects are consistent with the cusp of the track shown in the bottom-left panel of Fig. <ref>. Whilst rapidly quenched, these objects are not, in a literal sense, “post-starburst". By contrast, the third object is a literal post-starburst, with almost all its stellar mass formed within a ≲500 Myr period, approximately 1 Gyr before it is observed. This object has the youngest mass-weighted age of the galaxies in our quiescent sub-sample. It can be seen in Fig. <ref> to the extreme lower left of the quiescent population, and is the spectroscopic PSB closest to the <cit.> selection box. It also has the strongest Hδ absorption within our quiescent sub-sample (see Fig. <ref>). Further discussion of the SFHs and physical properties of VANDELS z∼1 PSBs will be presented by Wild et al. in prep.

Our model galaxies which quench more slowly can be seen to follow a more conventional path in Fig. <ref>, crossing the green valley approximately perpendicular to the edge of the UVJ selection box, and joining the red sequence higher up. These tracks are plausibly linked to  objects higher up the red sequence in Fig. <ref> which still retain higher A_V. However, most of their dust has already been lost by the time they cross the dashed line, preserving the strong age trend along the red sequence found by <cit.>.




 §.§ The star-formation histories of massive quiescent galaxies at 1.0 < z < 1.3


In light of the discussion of Section <ref>, it is interesting to consider the details of the SFH shapes we find for our quiescent sub-sample. The posterior median SFHs we infer are shown in Fig. <ref>, split into the same four mass bins used in Section <ref>. To demonstrate the magnitude of the typical uncertainties, in each bin a randomly selected object is highlighted in orange, and the 16^th-84^th percentiles of the posterior are shaded. The mean formation times for galaxies in each bin (corresponding to the black errorbars on the top panel of Fig. <ref>) are shown with black dashed vertical lines.

The average trend towards earlier formation with increasing stellar mass shown in Fig. <ref> is also visible in Fig. <ref>, however there is significant variation within each mass bin. In particular, in the highest mass bin, the SFH shapes are extremely diverse, which results in the average formation time being later than the next lowest mass bin. This demonstrates that, even at fixed stellar mass and observed redshift, the quiescent population contains galaxies with a wide range of formation histories.  The fact that knowledge of the stellar mass is not sufficient to make strong predictions about the SFH suggests again that a range of physical processes contribute to the quenching of star-formation.

In our lowest-mass bin the majority of objects formed recently and have very bursty SFHs with short ≲ 500 Myr formation timescales. A smaller number of similar objects are also present in our higher-mass bins, though they are a minority. The PSB shown in the bottom panel of Fig. <ref> is the youngest member of this group, though other members form later in cosmic time, as they have lower observed redshifts. Many or all of these objects are likely to have experienced a PSB phase at an earlier stage in their evolution. They are also plausibly linked with the population of highly star-forming submillimetre galaxies, the redshift distribution of which peaks at z∼2 (e.g. ). The number of such objects present in our VANDELS DR2 sample does not allow us to place strong statistical constraints on their number density as a function of redshift. However the full VANDELS sample will contain approximately twice as many objects consistent with our selection criteria.

The older galaxies in our sample display a broad range of quenching timescales. Those which quench rapidly are likely to have experienced both green-valley and PSB phases, consistent with our “early fast quench" model in Fig. <ref>. A fraction however also have SFHs consistent with our “early slow quench" model, with quenching timescales of ≳1 Gyr. The two oldest galaxies in our sample, both of which have log_10(M_*/M_⊙) > 11, quenched at z > 3. This corresponds to a comoving number density of 2±1×10^-5 Mpc^-3 for quiescent galaxies at z=3, consistent with , who find a number density of 1.4±0.3 × 10^-5 Mpc^-3 for quiescent galaxies with log_10(M_*/M_⊙) > 10.5 at 3 < z < 4.

The “early slow quench" objects may genuinely have experienced slow quenching at the highest redshifts, however two other scenarios are possible. Firstly, their SFHs may be composites of several systems which formed and quenched at different times, then merged to form the most massive quiescent galaxies. Secondly, it is possible that these systems were rapidly quenched, however, because their stellar populations are already old when observed, and hence slowly evolving, our observations cannot rule out slower quenching scenarios. Future instruments, such as NIRSpec on board the James Webb Space Telescope (JWST), which hold the promise of extending VANDELS-quality observations to quiescent galaxies at the highest redshifts, will provide invaluable constraints on the earlier evolution of these systems.




§ CONCLUSION


In this work we report SFHs for a sample of 75 UVJ-selected galaxies with stellar masses of log_10(M_*/M_⊙) > 10.3 at observed redshifts of 1.0 < z < 1.3. As described in Section <ref>, our data consist of deep rest-frame UV spectroscopy from VANDELS, as well as multi-wavelength photometry. Using the  code we fit our combined datasets with a sophisticated joint model for the physical properties of our galaxies and systematic uncertainties affecting our spectroscopic data, described in Sections <ref> and <ref>. The combination of extremely deep VANDELS spectroscopy with our sophisticated fitting methodology allows us to significantly improve upon previous analyses, obtaining strong, yet realistic constraints on the SFHs of our target galaxies.

We firstly quantify the average stellar mass vs stellar age relationship for massive quiescent galaxies at 1.0 < z < 1.3 (Section <ref>). We find a steep trend towards earlier mass-weighted formation times with increasing stellar mass (downsizing) of 1.48^+0.34_-0.39 Gyr per decade in mass (Equation <ref>). We observe some evidence for the flattening of this trend at the highest masses (M_* > 10^11M_⊙), as was reported by <cit.>. A slightly shallower trend of 1.24^+0.27_-0.30 Gyr per decade in mass is observed for r-band light-weighted formation times (Equation <ref>). As shown in Fig. <ref>, the slope of this relationship is in agreement with spectroscopic results at z∼0.1 and z∼0.7 from <cit.>, and at z∼1.75 from <cit.>. 

Recent photometric studies find weaker trends, of ≲0.5 Gyr per decade in mass (; ; see Fig. <ref>), which we attribute to larger uncertainties in individual age determinations due to the age-metallicity-dust degeneracy (see Section <ref>). We conclude that, at fixed observed redshift, an evolution in the stellar age vs stellar mass relationship of ∼1.5 Gyr per decade in mass is a robust result, which holds across the observed redshift interval from 0 < z < 2.

This result places strong constraints on the AGN-feedback models used in modern cosmological simulations. As such, in Fig. <ref>, we compare predictions from the Simba and IllustrisTNG simulations to our observational results. We conclude that the stellar mass vs stellar age relationships predicted by these simulations at z=0.1 are in good agreement with observations. However, at z=1, the relationships predicted are considerably flatter than our observational results, at ≲0.5 Gyr per decade in mass. This conclusion supports recent findings, which suggest that, whilst modern simulations now reproduce well the properties of local massive quiescent galaxies (e.g. ; ), agreement is increasingly poor with increasing observed redshift (e.g. ).

By considering the distributions we find for galaxy physical parameters on the UVJ diagram (Fig. <ref>), and the SFH shapes we infer for our sample (Fig. <ref>), we attempt to understand the connection between green-valley, post-starburst (PSB) and quiescent galaxies, and to constrain quenching mechanisms at z>1. We demonstrate that typical green-valley galaxies, if rapidly quenched, pass through a PSB phase en route to quiescence, and show that SFHs consistent with this evolution exist within our sample. 

We additionally identify a class of predominantly log_10(M_*/M_⊙) ∼ 10.5 galaxies which formed and quenched at z<2 in extreme starbursts with ≲500 Myr timescales. These objects are consistent with mergers and associated strong, AGN-driven outflows, and are plausibly related to submillimetre galaxies. These objects also pass through a PSB phase, supporting a dual origin for the PSB population (e.g. ; ; ). We finally find that some of our most massive galaxies appear to exhibit slow quenching at early times (z>2), though their quenching timescales are harder to constrain due to their older stellar populations. These objects are plausibly the result of mergers between galaxies with stellar populations formed at different times. To understand the earlier evolution of these systems in detail, deep continuum spectroscopy must be extended to the highest redshifts.

Our results demonstrate the power of large, high-redshift spectroscopic surveys for placing strong constraints on the evolution of subtle galaxy physical parameters across cosmic time. These results are important for furthering our understanding, as they are strongly constraining on models of galaxy formation. Upcoming instruments such as NIRSpec on JWST and MOONS at the VLT will greatly expand our high-redshift spectroscopic capabilities, and sophisticated fitting methodologies, such as presented in this work, will be key to realising their full potential.



§ ACKNOWLEDGEMENTS


The authors would like to thank Philip Best, Dylan Nelson and Po-Feng Wu for helpful discussions which have improved the quality of this work. A. C. Carnall, F. Cullen and S. Appleby acknowledge the support of the UK Science and Technology Facilities Council. A. Cimatti acknowledges the grants PRIN MIUR 2015, ASI n.I/023/12/0 and ASI n.2018-23-HH.0. This work is based on data products from observations made with ESO Telescopes at La Silla Paranal Observatory under ESO programme ID 194.A-2003(E-Q). This work is based in part on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a NASA contract. This research made use of Astropy, a community-developed core Python package for Astronomy <cit.>.

mnras
 


