Fast, Slow, Early, Late: Quenching Massive Galaxies at z ∼ 0.8

The spectroscopic data have been taken as part of the HALO7D program, a survey conducted in CANDELS fields (Grogin et al. 2011; Koekemoer et al. 2011) with the Keck II/DEIMOS instrument (Faber et al. 2003). HALO7D is a multisemester program with the main goal of surveying faint halo stars with Hubble Space Telescope (HST) measured proper motions in order to measure their line-of-sight velocities and chemical abundances, giving 6D phase-space information and chemical abundances for hundreds of remote Milky Way halo stars (Cunningham et al. 2019a, 2019b). The targeted fields, together with the deep exposures necessary to reach the faintest stars in the Milky Way halo, are an opportunity for a novel synergy of extragalactic and Galactic science. In addition to the primary halo star targets—which only occupy about a quarter of slitlets on a given DEIMOS mask—spectra for extragalactic targets have been taken. These data have been used to study galactic winds in z ∼ 1 (Yesuf et al. 2017; W. Wang et al. 2021, in preparation), internal galaxy kinematics (G. Barro et al. 2021, in preparation), and dwarf galaxies (Y. Guo et al. 2021, in preparation). Here we focus on the highest-priority filler sample of galaxies, i.e., massive star-forming and quiescent galaxies at z ∼ 0.8.

2.1.1. Sample Selection

2.1.2. Observations and Data Reduction

The observations and data reductions are described in detail in Cunningham et al. (2019a). In summary, the HALO7D observations use the 600 line mm^–1 grating on DEIMOS centered around 7200 Å with the GG455 order-blocking filter. This setup gives a nominal wavelength coverage of 4600–9500 Å at a resolution (FWHM) of ∼ 3.5 Å for a 1'' slit width and 0.65 Å pixel^–1 dispersion. The slit position angles are set to within ±30% of the parallactic angle to minimize light loss in the blue owing to atmospheric dispersion. In total, 232 galaxies were observed. Useful data (continuum signal-to-noise ratio (S/N) of at least 5 Å^–1, no major artifacts in the data, and no active galactic nucleus (AGN) with point sources or with broad lines) have been collected for 161 galaxies. For the remainder of this paper, we focus on those galaxies. The exposure times of these galaxies range between 4.0 and 48.6 hr, with an average of 11.1 hr. The HALO7D observations were reduced using the automated DEEP2/DEIMOS spec2d pipeline developed by the DEEP2 team (Cooper et al. 2012; Newman et al. 2013), which, among other things, performs the sky subtraction. Wavelength regions that are heavily affected by skylines are masked, making up about 15% of all the pixels. Calibrations were done using a quartz lamp for flat-fielding and red NeKrArXe lamps for wavelength calibration. We do not perform any flux calibration of the spectra since we directly model the spectroscopic flux calibration during fitting (Section 4.1). The instrumental line-spread function (LSF) has been measured from these arc lamps and the night skylines, as described in Appendix A. Importantly, our analyses of both photometry and spectroscopy in this work make the simplifying assumption that all galaxies are spatially uniform. In future work, we will account for and exploit spatial variations in colors and stellar populations that do exist (e.g., Szomoru et al. 2013; Tacchella et al. 2015b; Mosleh et al. 2017; Suess et al. 2019).

Figure 1 shows the stellar mass (M_⋆) as a function of observed redshift (z_obs) for our sample of 161 galaxies. The stellar masses are obtained from our SED modeling, as described in Section 4. The color scaling of the points indicates the S/N of the spectra, measured in an observed-frame wavelength window of λ_obs = 7000–9200 Å. Our sample spans a wide range in redshift with z_obs ≈ 0.4–1.2 and about two orders in stellar mass (M_⋆ ≈ 10¹⁰–10¹² M_⊙). The core of the sample lies in the redshift interval 0.6–1.0, with an average redshift of 〈z_obs〉 = 0.8. The median stellar mass is $\mathrm{log}{M}_{\star }/{M}_{\odot }=11.0$. The galaxies at z_obs < 0.6 have the highest S/N. In the core redshift range of our sample (z_obs = 0.6–1.0), no clear trend of redshift with S/N exists. The quality of our spectra is comparable to that of the LEGA-C survey (van der Wel et al. 2016), but our sample is smaller while probing a larger redshift range. Importantly, although the redshift range probed by our galaxies is rather large, thanks to the broad wavelength coverage (λ_obs = 4600–9500 Å), key absorption features are covered by all galaxies. Specifically, the core sample of our galaxies probes the hydrogen absorption lines from H10 (found at 3799 Å) to Hβ (at 4863 Å), the calcium H and K lines (at 3934 and 3969 Å), the CN line (at 4160 Å), the Mg i b triplet (at 5176 Å), and several other Mg (at 5530 Å), Ca (including at 4227 and 4455 Å), and Fe lines (including at 4383, 4531, 4668, and 5270 Å). Only galaxies probing the highest redshifts (z_obs > 1) do not have coverage of the MgIb triplet and the Hβ line.

Zoom In Zoom Out Reset image size

We match the 161 HALO7D galaxies to photometric catalogs. As described in the previous section, most of our galaxies (82.8% of the sample) lie in the CANDELS survey footprint. Specifically, 19.0%, 20.9%, and 42.9% of the sample lies in COSMOS, EGS, and GOODS-N, respectively. We match those galaxies with 3D-HST photometric catalogs (Skelton et al. 2014). Specifically, our galaxies are covered by between 17 (the EGS field) and 44 (the COSMOS field) photometric bands spanning a range of 0.3–8 μm in the observed frame. The photometry is supplemented by Spitzer/MIPS 24 μm fluxes from Whitaker et al. (2014). The MIPS 24 μm coverage is important because the rest-frame MIR wavelengths are dominated by warm dust emission, a key empirical proxy for obscured star formation (Kennicutt 1998). The other 28 galaxies (17.2% of the sample) are lying not in the 3D-HST footprint but in the extended EGS region. We use the UV–IR photometry in the Rainbow database ¹⁹ published by Barro et al. (2011a) for those objects.

The 3D-HST team self-consistently rederives zero-points for each instrument and filter in order to bring data from different instruments onto a common flux scale. Details are described in Skelton et al. (2014). Since this process is imperfect, we adopt the procedure by Leja et al. (2019b) and add the zero-point correction for each band of photometry to the flux errors in quadrature. This effect varies from 0% to 28% of the total flux, depending on the photometric band. Additionally, a 5% minimum error is enforced for each band of photometry to allow for systematic errors in the physical models for stellar, gas, and dust emission.

In this section, we compare our galaxy sample to the underlying galaxy population at z = 0.6–1.0, the core redshift range of our sample. Figure 2 shows from top left to bottom right the planes of SFR − M_⋆, rest-frame UVJ colors, Σ₁–M_⋆, and R_e–M_⋆. The circles indicate our sample, while the background hexbins show the whole galaxy population of CANDELS/3D-HST. Specifically, the stellar population parameters of the CANDELS/3D-HST comparison sample (M_⋆, SFR, and UVJ rest-frame colors) have been taken from Leja et al. (2019b), while H-band half-light size R _e and the central stellar mass surface density within 1 kpc Σ₁ are obtained from van der Wel et al. (2014). In particular, we estimate Σ₁ (Cheung et al. 2012; Saracco et al. 2012; Fang et al. 2013; van Dokkum et al. 2014; Tacchella et al. 2015a; Barro et al. 2017; Tacchella et al. 2017; Woo et al. 2017) by computing the fraction of the total luminosity in the H band within 1 kpc from the single Sérsic fits by van der Wel et al. (2014), assuming a constant mass-to-light ratio throughout the galaxy. The blue line in the top left panel shows the SFMS (Leja et al. 2021), the black line in the top right panel shows the UVJ-quiescent box as defined in Whitaker et al. (2012a), the red line in the bottom left panel marks the Σ₁ − M⋆ relation for quiescent galaxies as measured in Barro et al. (2017), and the blue and red lines show the R_e–M_⋆ relations for star-forming and quiescent galaxies from van der Wel et al. (2014), respectively. The color-coding in all panels corresponds to the doubling number ${\mathscr{D}}$, a measure of the number of times the stellar mass would double over the age of the universe at the current sSFR (see Equation (2)). As mentioned in the Introduction, we classify galaxies as star-forming, transitioning, and quiescent if ${\mathscr{D}}(z)\gt 1/3$, $1/20\lt {\mathscr{D}}(z)\lt 1/3$, and ${\mathscr{D}}(z)\lt 1/20$, respectively. These cuts correspond to blue, green, and red coloring in Figure 2, verifying that these cuts are meaningful.

Zoom In Zoom Out Reset image size

As shown in Figure 2, the galaxies of our sample follow the overall trends of the massive galaxy population. At a given M_⋆, star-forming galaxies have lower Σ₁ and larger R _e than quiescent galaxies. By selection (Section 2.1.1), our sample is unbiased above 10¹¹ M_⊙, while below this limit it is skewed toward quiescent galaxies (low SFRs and UVJ-quiescent colors). Importantly, our quiescent galaxies span the whole parameter space of Σ₁–M_⋆ and R_e–M_⋆ of quiescent galaxies. The star-forming galaxies in our sample are typically massive with M_⋆ > 10¹¹ M_⊙. Several star-forming galaxies (based on the doubling number ${\mathscr{D}}$) in our sample lie in the UVJ-quiescent box, consistent with the expected contamination rate of roughly 20% (Leja et al. 2019c). We take a more detailed look at those objects and the UVJ color–color diagram in Appendix E.

Throughout this paper, we compare our observational measurements with similar measurements from the cosmological, hydrodynamical simulation IllustrisTNG (Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2018; Pillepich et al. 2018; Springel et al. 2018; Nelson et al. 2019b). A detailed comparison between IllustrisTNG and our observations will allow us to assess the validity of adopted subgrid models in IllustrisTNG. In particular, we will focus on processes that are related to shaping SFHs, such as quenching, which in IllustrisTNG is mainly driven by black hole feedback: galaxies in IllustrisTNG quench once the energy from black hole kinetic winds at low accretion rates becomes larger than the gravitational binding energy of gas within the galaxy stellar radius (e.g., Terrazas et al. 2020). This occurs at a particular black hole mass threshold.

Throughout this paper, we focus on the intermediate-sized box TNG100-1 (TNG100 for the remainder of the paper), which combines a moderate resolution with a large volume to allow us to track the evolution of massive quiescent galaxies at z ≈ 0.5−2. The baryonic mass resolution of this box is 1.4 × 10⁶ M_⊙, and the box size is 110.7³ Mpc³. We randomly select a sample of 1340 out of the 5703 galaxies with M_⋆ > 10¹⁰ M_⊙ from the snapshot at z = 0.7. For each galaxy in TNG100, the stellar properties are extracted by considering all the bound particles. About a third (475 galaxies) of those galaxies are quiescent, consistent with observational estimates (Donnari et al. 2019). In order to perform a fair comparison between our observational measurements and the measurements from TNG100, we "project" the theoretical quantities into the observational space (spectra and photometry). Specifically, we predict the DEIMOS spectra and photometry by using the stellar particles' ages and metallicities and by drawing the other parameters of the galaxy SED, such as dust attenuation, dust emission, and AGN emission, from the prior (see Section 3). For this, we use the same stellar population synthesis models as in the fitting (Section 3.1). We add noise to the spectra and photometry by drawing from the noise distribution of our observational data (see Figure 1). After predicting realistic spectroscopic and photometric data from TNG100 galaxies, we run the same analysis on them as on our observations (Sections 3 and 4) and apply the same UVJ/sSFR selection as in the observational sample when comparing to the observations in the individual figures.

Although we perform this "apples-to-apples" comparison between observations and IllustrisTNG, the comparison still has its limitations. First, even though the observations and simulations probe the same comoving volume within a factor of 2, cosmic variance and sample selection (both observational and simulated samples are not mass complete) can affect the interpretation of the rarest objects. Second, because of the finite mass resolution, the intermediate-sized box TNG100 is not able to resolve the first-forming galaxies (e.g., Vogelsberger et al. 2020). The high-resolution box TNG50 (Nelson et al. 2019a; Pillepich et al. 2019) has an order of magnitude higher baryonic mass resolution, but its volume is a factor ∼8 smaller. Hence, massive quiescent galaxies are not well probed, in particular at higher redshifts. Third, aperture effects might play a role: even though our spectra include most of the light (at least 80% for most of our galaxies as estimated from the slit geometry and the HST H-band morphology), this may still be different from including all the bound stellar particles in IllustrisTNG. This should be kept in mind when we interpret the comparison of the observations with IllustrisTNG.

Additionally, these synthetic mock observations of the IllustrisTNG galaxies allow us to assess systematic uncertainties in the estimated parameters compared to the true values. A preliminary examination reveals that our inferred mass-weighted ages and star formation timescales are on average overestimated by ∼20% and ∼10%, respectively. We do not find any bias for the quenching epoch or quenching timescale over the whole sample, though there are a fraction of galaxies (about 20%) that quench rapidly in the simulation (τ_quench < 200 Myr) for which we overestimate the quenching timescale (τ_quench ≈ 0.5–1.0 Gyr). These galaxies typically quench early (more than 1 Gyr before the epoch of observation), making it difficult for us to pick up the signal of "fast quenching" in the data (in addition to having SFH bins on scales of Gyr on long look-back times; see the next section). Although these limitations do not alter our conclusions, they are of course important to consider in more detail. We postpone a more thorough analysis of this to an upcoming publication.

For stellar population synthesis, the Flexible Stellar Population Synthesis (FSPS) package ²⁰ is used (Conroy et al. 2009). In this work we use the MIST stellar evolutionary tracks and isochrones (Choi et al. 2016; Dotter 2016) with the MILES stellar spectral library (Falcón-Barroso et al. 2011). The MIST models are based on MESA, an open-source stellar evolution package (Paxton et al. 2011, 2013, 2015, 2018).

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 prior that is uniform in $\mathrm{log}({Z}_{\star }/{Z}_{\odot })$ between −1.0 and 0.19, where Z_⊙ = 0.0142. The upper limit of the stellar metallicity prior is given by the MILES stellar templates. This upper limit might lead to an underestimation of the stellar metallicity in the fitting, which itself would result in an overestimation of the stellar ages. Although we cannot rule this out completely, we find that our mass-weighted ages are on average overestimated by only ∼ 20% (0.05 dex) using the mock spectra of IllustrisTNG (Section 2.4). We intend to constrain the abundance pattern in more detail in an upcoming publication. Finally, a Chabrier (2003) initial mass function is assumed throughout this work.

In our fiducial setup, we assume a "nonparametric," piecewise constant SFH (Cid Fernandes et al. 2005; Ocvirk et al. 2006; Tojeiro et al. 2007). "Nonparametric" here means that no particular shape for the SFH is assumed and that an arbitrary function in SFH space can be reasonably approximated. Lower et al. (2020) showed that this flexible nonparametric approach outperforms traditional parametric forms (such as exponentially declining or lognormal SFHs) in capturing variations in galaxy SFHs, leading to significantly improved stellar masses in SED fitting. We assume that the SFH can be described by N_SFH time bins, where the SFR within each bin is constant. We fix N_SFH = 10 for the science analysis; we explore varying the number of time bins in Appendices C and D. Increasing the N_SFH to 14 does not affect our results (see also Leja et al. 2019a), but the fits of 30 galaxies (19% of the sample) do not converge within a reasonable amount of time (i.e., 14 days on a single CPU). There are approaches that determine the appropriate number of time bins on the fly, such as adaptively binning in time (Tojeiro et al. 2007) or using evidence comparison to determine the optimal number of bins (Dye et al. 2008; Iyer et al. 2019). As discussed in Leja et al. (2019a), we use a piecewise model with a fixed number of bins because it is more scalable in a sampling framework.

The N_SFH time bins are specified in look-back time. Throughout the paper, we define the look-back time to be the time prior to the epoch of observation. Four bins are fixed at 0–30 Myr, 30–100 Myr, 100–300 Myr, and 300–1000 Myr to capture variation in the recent SFH of galaxies. We model a maximally old population with a fifth bin at (0.95t_H − t_H), where t_H is the age of the universe at the observed redshift. The remaining N_SFH − 5 bins are spaced equally in logarithmic time between 1 Gyr and 0.95t_H.

Since we include "more bins than the data warrant" and let the sampler fully map the interbin covariances allowed by the prior and the data, a potential failure mode for this is overfitting, which is caused by an excess of model flexibility and results in overestimated uncertainties. This is in contrast to the classic dangers of "underfitting," whereby model parameters are overly constrained when too little parameter space is permitted. The danger posed by overly flexible models can be alleviated by choosing a prior that weights for physically plausible SFH forms. This is a complex problem that has been explored in detail in Leja et al. (2019a; see also Carnall et al. 2019a). We follow this analysis by adopting a continuity prior, which enforces smoothness by weighting against sharp transitions in SFR(t), similar to the regularization schemes from Ocvirk et al. (2006) and Tojeiro et al. (2007). The prior is tuned to allow similar transitions in SFR to those of galaxies in the Illustris hydrodynamical simulations (Vogelsberger et al. 2014; Diemer et al. 2017), though it is deliberately designed to include broader behavior than seen in these simulations since we do not want to assume these models to be the truth. The resulting prior probability densities for SFR(t), sSFR, mass-weighted age, star formation timescale, and quenching timescale are shown in Figure 3.

Zoom In Zoom Out Reset image size

In addition to this nonparametric approach, we also consider parametric SFHs in order to be able to compare to previous literature and to explore possible biases. In the parametric approach, we assume that the SFH follows a delayed τ-model:

$\begin{eqnarray}&&\mathrm{SFR}(t)=(t-{t}_{a}){e}^{-(t-{t}_{a})/\tau }.\end{eqnarray} \tag{ 3 }$

The parameter τ is varied with a logarithmic prior between $-1.0\lt \mathrm{log}(\tau )\lt 10.0$, and the parameter t _a is varied with a uniform prior between 0 and the age of the universe at galaxies z_obs (t_H(z_obs)). The results for the parametric SFHs are shown and discussed in Appendices C and D. Briefly, parametric SFHs typically lead to younger ages than nonparametric SFHs, which biases the stellar mass low and the sSFR high.

We model dust attenuation using a two-component dust attenuation model with a flexible attenuation curve. Specifically, we use the two-component Charlot & Fall (2000) dust attenuation model, which postulates separate birth-cloud and diffuse dust screens. The birth-cloud component (τ_dust,1(λ)) in our model attenuates nebular emission and stellar emission only from stars formed in the past 10 Myr:

$\begin{eqnarray}&&{\tau }_{\mathrm{dust},1}(\lambda )={\hat{\tau }}_{\mathrm{dust},1}{\left(\displaystyle \frac{\lambda }{5500\,{\rm{\mathring{\rm A} }}}\right)}^{-1}.\end{eqnarray} \tag{ 4 }$

The diffuse component (τ_dust,2(λ)) has a variable attenuation curve and attenuates all stellar and nebular emission from the galaxy. We use the prescription from Noll et al. (2009):

$\begin{eqnarray}&&{\tau }_{\mathrm{dust},2}(\lambda )=\displaystyle \frac{{\hat{\tau }}_{\mathrm{dust},2}}{4.05}(k^{\prime} (\lambda )+D(\lambda )){\left(\displaystyle \frac{\lambda }{5500\,{\rm{\mathring{\rm A} }}}\right)}^{n}.\end{eqnarray} \tag{ 5 }$

${\hat{\tau }}_{\mathrm{dust},2}$ controls the normalization of the diffuse dust, n is the diffuse dust attenuation index, $k^{\prime} (\lambda )$ is the Calzetti et al. (2000) attenuation curve, and D(λ) is a Lorentzian-like Drude profile describing the UV dust bump. We tie the strength of the UV dust absorption bump to the best-fit diffuse dust attenuation index, following the results of Kriek & Conroy (2013). The free parameters in Equation (5) are therefore ${\hat{\tau }}_{\mathrm{dust},2}$ and n. We adopt a flat prior for ${\hat{\tau }}_{\mathrm{dust},2}$ ($0\lt {\hat{\tau }}_{\mathrm{dust},2}\lt 4.0$) and n (−1.0 < n < 0.4). The upper limit on n is chosen to disallow a flat attenuation curve, which would cause ${\hat{\tau }}_{\mathrm{dust},2}$ to be nearly fully degenerate with the normalization of the SED.

Although ${\hat{\tau }}_{\mathrm{dust},1}$ and ${\hat{\tau }}_{\mathrm{dust},2}$ have a similar effect on the SED and are often degenerate, it is important to distinguish between these parameters to properly predict emission lines, in particular the line equivalent widths. The total optical depth toward nebular emission lines is roughly twice that of the stellar continuum (Calzetti et al. 1994; Kashino et al. 2013; Price et al. 2014). In our dust attenuation model, this means ${\hat{\tau }}_{\mathrm{dust},1}\sim {\hat{\tau }}_{\mathrm{dust},2}$, since ${\hat{\tau }}_{\mathrm{dust},2}$ is applied to the entire emission from the galaxy. We adopt a joint prior on the ratio of the two in order to allow for some reasonable variation around the fiducial results in the literature: a clipped normal centered on 1 with a width of 0.3 in the range of $0\lt {\hat{\tau }}_{\mathrm{dust},1}/{\hat{\tau }}_{\mathrm{dust},2}\lt 2.0$.

We assume energy balance, i.e., all the energy attenuated by dust is reemitted in the IR (da Cunha et al. 2008). Thanks to this assumption, the MIR photometry delivers additional constraints on the total amount of dust attenuation and on the dust-free stellar SED. However, in order to apply energy balance and to compute L_IR from the UV–MIR SED, we need to make some assumptions about the shape of the IR SED.

We use the Draine & Li (2007) dust emission templates to describe the shape of the IR SED, which are based on the silicate-graphite-PAH model of interstellar dust (Mathis et al. 1977; Draine & Lee 1984). These templates have three free parameters controlling the shape of the IR SED: ${U}_{\min }$, γ _e , and q_PAH. ${U}_{\min }$ and γ _e together control the shape and location of the thermal dust emission bump in the IR SED, while q_PAH describes the fraction of total dust mass that is in PAHs. This last parameter is particularly important because a substantial fraction, or even a majority, of the MIR emission comes from strong PAH emission features. Since our photometry only includes bands up to MIPS 24 μm, we use informative priors for ${U}_{\min }$ and γ _e , while assuming a flat prior for q_PAH (see Table 1). The adopted priors are consistent with both the SINGS sample (Draine et al. 2007) and the Brown et al. (2014) galaxies with Herschel photometry and lead to a minimal amount of bias in L_IR, SFR, and dust attenuation in galaxies without far-IR photometry. This is discussed in detail in Appendix C of Leja et al. (2017).

Building on Leja et al. (2018), we adopt the AGN templates from Nenkova et al. (2008a) and Nenkova et al. (2008b). The CLUMPY AGN templates are incorporated in FSPS (Conroy et al. 2009), and a detailed description is given in the FSPS documentation. Only the dust emission from the central torus is included in this model; it is assumed that the UV and optical emission from the central engine is fully obscured by the AGN dust torus. This is a viable assumption since we discarded AGNs with point sources or with broad lines. Our AGN model has two free parameters: f_AGN, the ratio of the bolometric luminosity between the galaxy and the AGN, and τ_AGN, the optical depth of an individual dust clump at 5500 Å. A log-uniform prior is adopted for f_AGN, with an allowed range of 10⁻⁵ < f_AGN < 3. A log-uniform prior describes the observed power-law distribution of black hole accretion rates (Aird et al. 2017; Georgakakis et al. 2017; Caplar et al. 2018). A log-uniform prior on τ_AGN is adopted between 5 < τ_AGN < 150, as the SED response to logarithmic changes in τ_AGN is approximately linear (see Figure 1 in Leja et al. 2018).

We adopt the standard approach to generating nebular emission in FSPS, whereby the ionizing continuum from the model stellar populations is assumed to be fully absorbed by the gas and emitted as both line and continuum emission. The nebular line and continuum emission is generated using a CLOUDY (Ferland et al. 1998, 2013) grid within FSPS, as described in Byler et al. (2017). We assume for the gas-phase metallicity a uniform prior between $-2.0\lt \mathrm{log}({Z}_{\mathrm{gas}}/{Z}_{\odot }\,)\,\lt 0.5$ and for the ionization parameter U a uniform prior between $-4.0\lt \mathrm{log}(U)\lt -1.0$. Furthermore, we assume a flat prior for the gas-phase velocity dispersion (30 < σ_gas/(km s⁻¹) <300). In addition, motivated by the complexity of the physics that produce nebular emission lines (Kewley et al. 2019 and references therein), we take a flexible approach to model the nebular line amplitudes (Section 4.3).

As described in Section 2.1.2, the spectra are not flux-calibrated. At each likelihood call, we match the model spectrum to the normalization of the spectroscopic data by fitting a polynomial in wavelength to their ratio. Our approach implements the Chebyshev polynomial calibration model, computed at each likelihood call from a simple least-squares maximum likelihood fit to the ratio of the data to the calibrated model spectrum, excluding the regions where emission lines may be present. The order of the polynomial m is determined by $m=({\lambda }_{\max }-{\lambda }_{\min })/100{\rm{\mathring{\rm A} }}$ within each wavelength interval (see also Kelson et al. 2000 and Conroy et al. 2018), with typically m ≈ 20. This order is chosen to account for broad continuum mismatch issues but is not so flexible that it could overfit broad absorption features. We have experimented with this approach by changing the order of the polynomial and find that the results are generally insensitive to this choice.

Using the maximum likelihood fit for the calibration has the advantage of computational speed. However, ideally one would marginalize the likelihood of the data over all possible calibration polynomials for each model call. Naively, this can be done at the cost of introducing m additional model parameters describing the polynomial coefficients, but this is computationally prohibitive at present. It is possible to analytically marginalize the likelihood over all possible coefficient values, but this has not yet been implemented (but see Carnall et al. 2019b).

The net effect of our approach is that the large-scale continuum shape and normalization of the model are set by the photometry. By fitting a moderate-order polynomial to the ratio between the observed and physical model spectrum at each likelihood call, the spectroscopic calibration model basically removes all information content from the continuum shape of the spectroscopic data. This means that the continuum shape of the observed spectrum does not inform any of the galaxy's physical parameters. Instead, information about physical parameters that affect the continuum shape derives from the photometry, which does not include any multiplicative calibration model. Therefore, there is no degeneracy between the spectroscopic calibration model and the galaxy's parameters, such as the dust content or the SFH.

We find that the standard fitting procedure is sensitive to outliers, that is, spectroscopic data points that are not well described by our model, because of inaccurate uncertainties or limitations of the model itself. We mitigate this problem by becoming insensitive to "bad" spectral data points. Specifically, we use a mixture model to describe outliers, following the approach described in Hogg et al. (2010; see also Sharma 2017; Press 1997, pp. 49–60). These outlier pixels do not include the masked pixels (Section 2.1.2).

This model alters the likelihood by assuming some possibility that any given spectral pixel is an outlier, f_out. The likelihood is calculated by marginalizing over f_out for each pixel; thus, no individual pixels are uniquely identified as outliers (see Equation (7) below). It is assumed that outlier pixels have their uncertainties inflated by a factor of 50. f_out is a free parameter in the model, and typically 0.1% of pixels in each fit are outliers.

In addition, to account for possible under- or overestimates of the spectroscopic uncertainties, we introduce a parameter that multiplies the spectroscopic uncertainties by some constant factor (j_spec) before calculating the likelihood. This is a free parameter in the model; in general, we find values very close to 1, indicating that the spectroscopic noise is not broadly under- or overestimated.

As mentioned in Section 3.6, nebular emission lines are included in the model spectrum. Each emission line is modeled as a Gaussian with a variable width and amplitude. We fit for the velocity dispersion of the gas (σ_gas), while the emission-line amplitudes are marginalized over in each fitting step, as described in Johnson et al. (2021). In our fiducial fitting, we use the maximum likelihood amplitude for each emission line, which means that we totally decouple the emission lines from the SFH since the emission lines for most of our quiescent galaxies are believed to be emitted from low-ionization emission-line regions, shocks, and/or AGNs.

We describe here the fitting methodology. We fit the physical model described in Section 3 to the observational data (spectroscopy with photometry) presented in Section 2, together with the models for systematic effects described in Sections 4.1 and 4.2. We denote the parameters of the model with Θ.

We assume that the uncertainties on the photometric fluxes σ_i are Gaussian and independent. In this case, the scatter of observed fluxes f_i about their true values follows a χ² distribution. The log-likelihood function for the photometry then follows:

$\begin{eqnarray}&&\mathrm{ln}({{ \mathcal L }}_{\mathrm{phot}})=-0.5\sum _{i}^{{N}_{\mathrm{band}}}{\left(\displaystyle \frac{{f}_{i}-{m}_{i}({\rm{\Theta }})}{{\sigma }_{i}}\right)}^{2},\end{eqnarray} \tag{ 6 }$

where m_i (Θ) is the model prediction for the observed flux f_i and N_band is the number of photometric bands.

We make the same assumptions about the distribution of uncertainties when fitting the spectrum. However, the likelihood equation is modified by the outlier model such that

$\begin{eqnarray}&&{{ \mathcal L }}_{\mathrm{spec}}=(1-{f}_{\mathrm{out}}){ \mathcal L }(f,m,\sigma )+{f}_{\mathrm{out}}{ \mathcal L }(f,m,{\sigma }_{\mathrm{out}}),\end{eqnarray} \tag{ 7 }$

where ${ \mathcal L }$ is calculated following Equation (6) (instead of summing over the photometric bands, we sum over spectral pixels in wavelength) and σ_out = 50σ as described in Section 4.2.

The total likelihood is thus

$\begin{eqnarray}&&\mathrm{ln}{{ \mathcal L }}_{\mathrm{total}}=\mathrm{ln}{{ \mathcal L }}_{\mathrm{phot}}+\mathrm{ln}{{ \mathcal L }}_{\mathrm{spec}}+\mathrm{ln}{{ \mathcal L }}_{\mathrm{eline}},\end{eqnarray} \tag{ 8 }$

where ${{ \mathcal L }}_{\mathrm{eline}}$ represents a penalty term that takes into account the prior on emission-line amplitudes. We set this term to 0 in our fidicual analysis.

Our fiducial model assumes 10 SFHs bins, leading to a total of 27 free parameters (Table 1). Sampling our posterior distribution with the dynamic nested sampling algorithm dynesty (Speagle 2020) therefore requires several million evaluations of our log-likelihood function. Each likelihood call takes about 50 ms. Fitting each galaxy therefore requires roughly 100 CPU hours.

We verified the fits to the photometry and spectroscopy for each individual galaxy. We checked that none of the posteriors pile up at the edges of the priors. This is particularly true for the stellar metallicity, which alleviates the concern that our stellar age estimates are biased high. The median reduced χ² calculated with the model drawn from the posteriors for the photometry and for the spectroscopy are ${1.08}_{-0.52}^{+1.02}$ and ${0.83}_{-0.17}^{+0.22}$, respectively. The χ² for the spectroscopy is slightly below 1, which is consistent with our obtained jitter terms j_spec to be clustered at 1 (we do not allow for j_spec < 1.0). We have also inspected the stacked residuals of the photometry, finding that the model does reproduce the data well. We found a weak but significant trend related to the K-band and IRAC photometry. Specifically, the average K band has a χ of about 1 (i.e., model underestimated the K-band flux), while the IRAC bands show a gradient so that the model K-to-IRAC colors are too red with respect to observations. A possible cause for this trend could be thermally pulsing AGB (TP-AGB) stars. In principle, FSPS allows us to choose the normalization of the TP-AGB stars, though it is currently unfeasible to marginalize over this within Prospector. Nevertheless, this feature should be investigated in the future. We also investigated the stacked residual of the observed- and rest-frame spectra, finding no significant trend, which shows that the removing of skylines and modeling of emission lines overall worked well.

We show the observational data of an example galaxy along with the fitted model in Figure 4. The galaxy has a redshift of z_obs = 0.852. The spectroscopic data have S/N = 23.8 Å^–1 with an exposure time of 9.3 hr. The model fits the photometric and spectroscopic data well. The residuals are distributed around 0. The reduced χ² for the photometry and spectroscopy are 0.5 and 0.7, respectively.

Zoom In Zoom Out Reset image size

The resulting posteriors of some key quantities of this fit are shown in Figure 5. This joint posterior plot shows the stellar mass (M_⋆), sSFR, dust opacity in the V band (${\hat{\tau }}_{\mathrm{dust},2}$), stellar age (t_age), and stellar metallicity (Z_⋆). We find for this galaxy a stellar mass of $\mathrm{log}{M}_{\star }/{M}_{\odot }={10.87}_{-0.03}^{+0.03}$, a low sSFR with $\mathrm{logsSFR}/{\mathrm{yr}}^{-1}=-{11.59}_{-0.73}^{+0.38}$, a dust opacity of ${\hat{\tau }}_{\mathrm{dust},2}\,={0.16}_{-0.04}^{+0.04}$, an age of ${t}_{\mathrm{age}}/\mathrm{Gyr}={2.51}_{-0.36}^{+0.45}$, and roughly solar metallicity with $\mathrm{log}{Z}_{\star }/{Z}_{\odot }={0.04}_{-0.07}^{+0.06}$. Here and throughout the paper the age t_age corresponds to the mass-weighted age. From this fit, we find that this galaxy is quiescent with a doubling number of ${\mathscr{D}}=0.02$, i.e., it takes this galaxy about 50 Hubble times (age of the universe at z_obs) to double its mass with its current sSFR. Although this galaxy is not actively forming stars, the galaxy is overall young, consistent with an SFH that rose through most of cosmic time and only declined in the past ∼1 Gyr.

Zoom In Zoom Out Reset image size

In addition, we discuss in Appendix B the gain in fitting both photometry and spectroscopy. In summary, the spectroscopy constrains the metallicity, while the photometry constrains the dust attenuation. However, both the photometry and the spectroscopy are needed to break the dust–age–metallicity degeneracy and derive SFHs.

Figure 6 presents the SFH of all galaxies in our sample as a function of observed redshift (z_obs; from left to right) and of stellar mass (M_⋆; from top to bottom). We plot the SFHs as sSFR versus look-back time (SFR vs. look-back time is shown in Figure 7). The solid lines show individual galaxies, while the shaded region shows the 16th−84th percentiles. The coloring of the lines corresponds to whether galaxies at the time of observations (z_obs) are star-forming (blue), transitioning (green), or quiescent (red). Rejuvenating galaxies, i.e., galaxies that were quiescent in the past and at the epoch of observation in the transition or star-forming region, are shown as orange lines. Plotting in sSFR(t) gives the advantage of being able to directly read off from the figure whether a galaxy is still actively growing owing to star formation. In order to help guide the eye, the dashed and solid black lines indicate the boundary between the star-forming and transition regions (sSFR = 1/[3t_H(z)]) and between the transition and quiescent regions (sSFR = 1/[20t_H(z)]), meaning that galaxies that at least double their mass within three times the age of the universe are considered star-forming, while galaxies that do not double their mass within 20 times the age of the universe are considered quiescent (see Equation (2)). The green shaded region marks the transition region.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The key result presented in Figure 6 is the diversity of SFHs in our sample. At a given observed redshift and stellar mass, we find a large diversity of pathways through the sSFR−time space. This is consistent with other studies that find that massive galaxies form a highly diverse population at z ≳ 1 (e.g., van Dokkum et al. 2011). By definition, galaxies start with high sSFRs at early times with sSFR ≈ 10⁻⁸ − 10⁻⁹ yr⁻¹, which means that these galaxies double their mass every few hundred million years. This phase of star formation lasts only a few hundred Myr in some galaxies, while it lasts for several Gyr in other galaxies. There is a large diversity in both when galaxies start transitioning to the quiescent region (i.e., when they start quenching) and how long this transition takes (quenching timescale). We quantify both the duration of the star-forming phase (i.e., star formation timescale) and the quenching timescale in more detail in the upcoming sections.

Similarly, Figure 7 shows the SFR versus look-back time tracks for our sample, again highlighting the diversity of pathways. In particular, this figure highlights the range of different SFRs at early cosmic times. Although there is the overall trend that star-forming galaxies form later than quiescent galaxies (i.e., they peak at later cosmic times), there are several outliers that do not follow this trend. Furthermore, even at fixed M_⋆ and z_obs, quiescent galaxies themselves show a large diversity in early SFRs and the peak times, consistent with many pathways to quiescence. This is also highlighted separately in Figure 8, which shows the sSFR tracks for massive ($\mathrm{log}({M}_{\star }/{M}_{\odot })\gt 11.0$) quiescent galaxies with z_obs = 0.6–1.0. Quiescent galaxies quench fast, slow, early, and late—and not all galaxies that quench early quench fast, nor do all galaxies that quench late quench slowly.

Zoom In Zoom Out Reset image size

Despite this diversity, there is some overall coherence, in particular regarding the sSFR tracks in the star-forming phase. Therefore, it is worth studying the median SFHs for different samples. Figure 9 shows the median SFH (sSFR vs. cosmic time) for our sample, focusing now on galaxies with z_obs = 0.6–1.0. We show the median SFH for quiescent galaxies in the mass bins $\mathrm{log}({M}_{\star }/{M}_{\odot })=10.5\mbox{--}11.0$ and $\mathrm{log}({M}_{\star }/{M}_{\odot })=11.0\mbox{--}11.5$ as dotted and solid lines, respectively. The median SFH of star-forming galaxies with $\mathrm{log}({M}_{\star }/{M}_{\odot })=11.0\mbox{--}11.5$ is shown as a blue line. We do not show the low-mass star-forming galaxies because only a few galaxies are in this bin (see Figure 2). The gray band indicates the rescaled specific dark matter accretion rate with ${\dot{M}}_{{\rm{h}}}/{M}_{{\rm{h}}}\propto {\left(1+z\right)}^{2.5}$ (Wechsler et al. 2002; Neistein & Dekel 2008; Dekel et al. 2013).

Zoom In Zoom Out Reset image size

Figure 9 shows that during this early phase of star formation (within 2 Gyr of the big bang, i.e., z > 3) star-forming and quiescent galaxies have similar sSFRs and are consistent with direct measurements of SFR and M_⋆ (i.e., SFMS) at high z by Salmon et al. (2015). This shows that one can in principle use this archaeological approach to estimate the SFR and M_⋆ of the galaxy population in observationally inaccessible parameter space (i.e., high z and low M_⋆; e.g., Iyer et al. 2018), though caution must be exercised since these SFHs include all stellar mass ever accreted, i.e., it is difficult to correct for the effects of merging (see Section 6.2). The median sSFR track of star-forming galaxies follows the independent estimates of the SFMS also at lower redshifts (Leja et al. 2021), which is an important consistency and validation check of our obtained SFHs. Furthermore, the median SFH of star-forming galaxies lies within the gray band at all times, indicating that the sSFR evolution is consistent with the evolution of the specific mass accretion rate of dark matter halos.

Following this early phase, the median SFH of star-forming galaxies (blue line) tracks well the simple relation of the dark matter accretion rate, and it is also consistent with lower-z SFR and M_⋆ measurements by Leja et al. (2021). We find that the quiescent galaxies at z_obs ≈ 0.8 decouple from the SFH of star-forming galaxies at z ∼ 2–3 and the average quenching timescale is roughly 1.5–2.0 Gyr. Interestingly, higher-mass galaxies ($\mathrm{log}({M}_{\star }/{M}_{\odot })=11.0\mbox{--}11.5$) transition on average slightly earlier than the lower-mass galaxies ($\mathrm{log}({M}_{\star }/{M}_{\odot })=10.5\mbox{--}11.0$), though the average quenching timescale seems not to depend strongly on stellar mass. We discuss this further in Section 5.4.

The SFHs of the galaxies in our sample show a large diversity. Nevertheless, there are some common features. SFHs can be described in three parts: a star-forming phase, a transition phase, and the quiescent phase. We now use a few simple parameters to describe the overall shape, including the mass-weighted age, the formation time, the star formation timescale, and the quenching timescale (Table 2).

First, we focus on the mass-weighted age t_age, which can be directly computed from our SFHs. Importantly, t_age is defined as look-back time from the epoch of observation. Therefore, t_age depends on the epoch of observation (z_obs); Figure 10 shows t_age as a function of z_obs for massive galaxies ($\mathrm{log}{M}_{\star }/{M}_{\odot }\gt 10.5$). We confirm the large diversity of SFHs: at a given z_obs, we find a large variety of mass-weighted ages. At z_obs ≈ 0.8, some galaxies have relatively young ages (1 – 2 Gyr), while other galaxies are rather old with ages of >6 Gyr. The colors of the points correspond to the mass doubling ${\mathscr{D}}$ (Equation (2)), indicating that the quiescent galaxies (red points) are typically older than star-forming galaxies (blue points). However, there is a large overlap between these two kinds of galaxies, indicating that the star formation activity at the epoch of observation cannot tell the full story about the past SFH.

Zoom In Zoom Out Reset image size

The thin gray to thick black lines show the passive evolutionary paths of simple stellar populations (SSPs) with different formation redshifts, ranging from z_form = 1 to z_form = 10. Since our SFHs are typically not well described with SSPs, this is only meant to guide the eye, highlighting that our sample spans a wide range of formation redshifts (see the next figure). We find that several of our galaxies are close to being maximally old (z_form ≳ 10), i.e., a few galaxies are close to the upper boundary of allowed ages given by the age of the universe.

Figure 10 also compares our measurements of the mass-weighted age with other estimates in the recent literature. At redshift 0, we indicate with the red shaded region the current age constraints of massive quiescent galaxies from a range of literature (Gallazzi et al. 2005, 2021; Spolaor et al. 2010; Trussler et al. 2020). Carnall et al. (2019b), Belli et al. (2019), Estrada-Carpenter et al. (2020), and Stockmann et al. (2020) use a similar approach to that presented here, where extended (parametric and nonparametric) SFHs are fit to individual spectra. Choi et al. (2014) use full-spectrum stellar population synthesis modeling on a stack of quiescent galaxies quoting SSP-equivalent ages, while Onodera et al. (2015) use the Lick absorption-line indices to infer the light-weighted age on a stack of 24 massive quiescent galaxies. Similarly, Gallazzi et al. (2014) inferred light-weighted ages of individual galaxies at z ∼ 0.7 via Lick indices (we only plot the median of the 10¹¹ M_⊙ bin). Finally, Kriek et al. (2016) perform full-spectrum fitting on an individual quiescent galaxy at z = 2.1. Although this literature list is extensive, it is not complete, i.e., there are several other studies that constrain the ages of massive galaxies at intermediate redshifts that we have not plotted here (e.g., Ferreras et al. 2017; Jørgensen et al. 2017). In summary, our measurements are overall consistent with those measurements, which also indicate a large variety of formation redshifts. However, we acknowledge that the definition of "age" adopted in the different studies also leads to some spread and that ages—in particular, SSP-equivalent ages—are not straightforwardly translated in formation redshifts.

The difficulty with interpreting Figure 10 is that t_age at two different epochs cannot directly be compared with each other since t_age evolves with epoch for a quiescently evolving system. It is therefore more informative to compute the epoch that corresponds to that age: the formation redshift z _f . We plot z _f as a function of M_⋆ in Figure 11. The left panel shows all of the galaxies in our sample, while the middle and right panels only plot the UVJ-quiescent objects. The points are colored according to their star formation timescale τ_SF, which is defined by the time between t₂₀ and t₈₀ (time when 20% and 80% of the stellar mass was formed, respectively).

Zoom In Zoom Out Reset image size

Figure 11 shows that galaxies span a wide range of z _f at fixed M_⋆. As expected, star-forming galaxies have typically more recent formation redshifts. The most massive, quiescent galaxies ($\mathrm{log}({M}_{\star }/{M}_{\odot })\gtrsim 11.5$) formed early with z _f ≳ 4.0, with some having a formation of z _f ≈ 10. We find only a weak trend with stellar mass: more massive galaxies have formed slightly earlier. Fitting galaxies with ($\mathrm{log}({M}_{\star }/{M}_{\odot })\gt 10.5$) leads to ${t}_{\mathrm{form}}/\mathrm{Gyr}=(2.1\pm 0.1)-(1.3\pm 0.4)\cdot \mathrm{log}({M}_{\star }/({10}^{11}\,{M}_{\odot }))$, where t_form is the cosmic time that corresponds to z _f . Beyond this trend, our low-mass sample of quiescent galaxies shows that several of those objects interestingly formed as early as the high-mass galaxies. Finally, Figure 11 also shows a trend with the star formation timescale τ_SF: galaxies that formed early formed their stars in a shorter amount of time than galaxies that formed later. We discuss this further in the next section.

We compare the formation redshift of our sample with literature values in the middle panel of Figure 11, focusing only on UVJ-quiescent galaxies. The measurements of Estrada-Carpenter et al. (2020) span a similar range in z _f to ours, though they are probing a narrower mass range (i.e., miss the highest- and lowest-mass galaxies). Belli et al. (2019) span a smaller range in z _f , reporting no objects with z _f > 6 but finding a stronger trend with stellar mass. Carnall et al. (2019b) describe a similar mass trend to that found in our sample (see below), but their formation redshifts are systematically lower (see also Siudek et al. 2017). The Stockmann et al. (2020) measurements, probing the most massive galaxies, lie between our measurements and the measurements of Carnall et al. (2019b).

The blue line in the right panel of Figure 11 shows the best-fit relation for our sample: ${t}_{\mathrm{form}}/\mathrm{Gyr}=2.6\,-1.5\,\cdot \mathrm{log}({M}_{\star }/({10}^{11}\,{M}_{\odot }))$. As also reported in Carnall et al. (2019b), z = 0.7 galaxies in IllustrisTNG (TNG100; Nelson et al. 2018; Pillepich et al. 2018) seem to be inconsistent with this mass trend: in TNG100, the formation redshifts of quiescent galaxies seem not to depend on stellar mass. As described in Section 2.4, we go one step further than Carnall et al. (2019b) by projecting the TNG100 quantities into the observational space and perform the same measurements on the predicted spectra and photometry as in our observations (Section 2.4). The TNG100 results are shown in the right panel of Figure 11 as gray crosses. We confirm the negligible dependence of z _f on stellar mass: quiescent TNG galaxies at z = 0.7 seem to have a formation redshift of z _f ≈ 2 − 4, independent of M_⋆. Furthermore, TNG100 produces a similar wide distribution of z _f to that of our observations, but it seems to underproduce the earliest-forming galaxies with z _f ≳ 6. This also holds when looking at the true SFHs and ages, i.e., without processing the simulated galaxies through the observational pipeline. This could point to missing physics in the IllustrisTNG model. However, the lack of early-forming galaxies could also be explained by resolution effects, i.e., the TNG100 box might not resolve the generation of the first galaxies. The higher-resolution TNG50 box could (at least partially) alleviate this problem, but the volume probed is limited (TNG50 is roughly an order of magnitude smaller than TNG100), making cosmic variance effects more severe.

We focus further on the relation between the star formation timescale τ_SF, the formation redshift z _f , and the stellar mass M_⋆. Figure 12 plots τ_SF as a function of M_⋆, colored by z _f , focusing on UVJ-quiescent galaxies. At M_⋆ ≈ 10¹¹ M_⊙, τ_SF ranges from a few hundred Myr to 4 Gyr, highlighting again the diversity of pathways that achieve the same final stellar mass. As already indicated in Figure 11, there is a clear trend that early-forming galaxies have shorter star formation timescales. While this is consistent with observational studies at z = 0 that are based on the average ages and element abundance ratios of elliptical galaxies (e.g., Thomas et al. 2005; Graves et al. 2010), our model assumptions (i.e., our SFH prior) allow for a different outcome (see also Appendix D); this means that the data unequivocally prefer a correlation between early formation and shorter SFH timescales. As shown in Figure 19 and verified on the whole data set, the degeneracy in the fitting is along "older age (higher z _f )—larger τ_SF." This can be understood by the following. The SFH is typically more constrained at recent times than at early times, implying that if one reduces the age, one needs to reduce the SFH at early times. This then leads to a shortening of τ_SF. This is the opposite trend shown in Figure 13, where older galaxies (with larger z _f ) have typically a shorter τ_SF.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Furthermore, our measurements show that galaxies with recent formation redshifts (z _f < 3) all have rather long star formation timescales (τ_SF > 2 Gyr), highlighting that galaxies that formed most of their mass in just 2 Gyr are uncommon at lower redshifts ("late bloomers"; Section 6.3). Additionally, there is a weak trend that more massive galaxies have a shorter τ_SF. In summary, this 3D plane can be described by

$\begin{eqnarray}\begin{array}{rcl}{\tau }_{\mathrm{SF}}/\mathrm{Gyr} & = & (4.1\pm 0.2)-(0.3\pm 0.2)\cdot \mathrm{log}\left(\displaystyle \frac{{M}_{\star }}{{10}^{11}\,{M}_{\odot }}\right)\\ & & -(3.1\pm 0.2)\cdot \mathrm{log}(1+{z}_{{}_{f}}),\end{array}\end{eqnarray} \tag{ 9 }$

confirming that the mass dependence is rather weak ($\propto {M}_{\star }^{-0.3}$), while the scaling with epoch of formation is strong ($\propto {\left(1+{z}_{{\rm{f}}}\right)}^{-3.1}$). This fit is shown as solid lines in Figure 12 for z _f = 2, 4, 6, and 10, suggesting that this fit reproduces our measurements well.

Having described when and over which timescales star formation occurred in galaxies, we now move on to characterizing over which timescales galaxies cease their star formation. The analysis and figures in this section include all galaxies that quenched at some point in their past, i.e., mostly UVJ-quiescent galaxies, but also galaxies that went through a rejuvenation event recently. We quantify this by the quenching timescale τ_quench, which measures the time spent in the transition region between 1/[20t_H(z)] < sSFR < 1/[3t_H(z)] (between the dashed and solid black lines in Figure 6). Furthermore, we measure the quenching epoch z_quench (Table 2), which represents when the galaxy is halfway through the transition region. Typically, the epoch of quenching is more tightly constrained than the quenching timescale (Appendix D). Changing the definition of z_quench to the epoch when quenching starts has a negligible impact on our results in this section and the conclusion of this work.

Figure 13 shows τ_quench as a function of z_quench, color-coded by M_⋆. The forbidden region, where galaxies quench on longer timescales than the age of the universe, is marked in gray. The left panel of Figure 13 shows our measurements, while the right panel shows the data from IllustrisTNG (TNG100). We find that the observed galaxies span a wide range of τ_quench and z_quench: consistent with Figure 8, some galaxies transition early (z_quench > 3), while some galaxies quench recently (z_quench ≈ 0.8 − 1.0). Furthermore, τ_quench spans basically the full allowed range at all z_quench. We find a median quenching timescale of ${\tau }_{\mathrm{quench}}={1.1}_{-0.7}^{+1.2}$ Gyr and ${\tau }_{\mathrm{quench}}/{\tau }_{{\rm{H}}}({z}_{\mathrm{quench}})\,={0.23}_{-0.16}^{+0.24}$ when normalized by the age of the universe.

Our observational measurements are compared to our measurements from IllustrisTNG, processed in the same way (Section 2.4). The quiescent galaxies from TNG100 have a median ${\tau }_{\mathrm{quench}}={0.7}_{-0.4}^{+0.7}\,\mathrm{Gyr}$ and typically quench recently with z_quench ≈ 1–2—consistent with the typical measurements of our observations. However, TNG100 galaxies span a narrower range in both τ_quench and z_quench. There are no galaxies that quenched early (i.e., z_quench > 3), consistent with the absence of early-forming galaxies (z _f > 6) shown in Figure 11. This may indicate that the quenching pathways in TNG100 are too monotonic and more diversity is needed. On the other hand, although the volumes probed by our observations and simulations are comparable, cosmic variance could still play a role concerning the abundance of rare objects. For example, Merlin et al. (2019) and Valentino et al. (2020) find quiescent objects at z > 3 in the larger TNG300 box. However, those studies directly probe these high-redshift snapshots, while we perform (consistent with our observations) an archaeological stellar population approach (we expect that mergers lower z_quench).

A median quenching timescale in TNG100 of ${\tau }_{\mathrm{quench}}\,={0.7}_{-0.4}^{+0.7}\,\mathrm{Gyr}$ for galaxies at z = 0.7 is consistent with the timescales quoted in Nelson et al. (2018), who find that the transition timescale through the green valley for massive z ∼ 0 galaxies is roughly 1 Gyr. They find that lower-mass z ∼ 0 galaxies transition on longer timescales (∼2 Gyr). In our analysis at z = 0.7, this trend is weaker, i.e., galaxies at all masses have similar quenching timescales: in our observations, we find ${\tau }_{\mathrm{quench}}={1.2}_{-0.8}^{+1.1}\,\mathrm{Gyr}$ and ${\tau }_{\mathrm{quench}}={1.0}_{-0.9}^{+0.6}\,\mathrm{Gyr}$ for galaxies with $\mathrm{log}({M}_{\star }/{M}_{\odot })=10.5$ and $\mathrm{log}({M}_{\star }/{M}_{\odot })=11.5$, respectively. This could point toward a difference due to epoch of observations (environmental quenching setting in at lower redshifts; e.g., Webb et al. 2020) or due to method (following the SFH evolution of the progenitor vs. an archaeological stellar population approach)—something to look into in more detail in the future.

We investigate in Figure 14 how τ_quench depends on stellar mass (M_⋆; left panel) and halo mass (M _h ; right panel). The halo masses are obtained for a subsample (only those that are plotted in Figure 14) of our galaxies from Fossati et al. (2017), who leverage the spectroscopic and grism redshifts from the 3D-HST survey to derive densities in fixed apertures to characterize the environment of galaxies brighter than JH₁₄₀ < 24 mag in the redshift range 0.5 < z < 3.0. The uncertainties on these halo masses are large. We compared these halo masses with the central stellar mass density (Σ₁)—a proxy for the central velocity dispersion, which has been suggested to be correlated with the halo mass (Schechter 2015; Zahid et al. 2016; Utsumi et al. 2020). We find indeed a correlation between Σ₁ and M _h and that galaxies follow the relation by Zahid et al. (2016), although with a large scatter in M _h of about 0.5 dex at fixed Σ₁. The Σ₁ − M _h connection has been extensively discussed in Chen et al. (2020): smaller galaxies have smaller halos and lower-density centers, but there is a large scatter for star-forming galaxies. The scatter of the Σ₁ − M _h relation gets even larger after quenching, when additional, post-quenching physical processes set in. This is something we will explore more in the future when studying the connection between the SFHs and the morphologies of the galaxies.

Zoom In Zoom Out Reset image size

Although the uncertainties are large, Figure 14 shows that galaxies in massive halos with $\mathrm{log}({M}_{{\rm{h}}}/{M}_{\odot })\gt 13.0$ all have short quenching timescales (τ_quench < 1 Gyr), i.e., there is an absence of galaxies with long quenching timescales. The points in Figure 14 are color-coded by z _f , highlighting that these massive halos also host galaxies that formed early (z _f > 5). At lower halo masses ($\mathrm{log}({M}_{{\rm{h}}}/{M}_{\odot })\approx 12.0-13.0$), τ_quench is also correlated with z _f > 5, which is not surprising, since quiescent galaxies that formed recently need also to quench quickly, as described above. A similar, but weaker, trend can be found for τ_quench − M_⋆.

One of the key results of our analysis is the large diversity of SFHs. Specifically, galaxies at a given epoch and stellar mass show a broad range of different paths in the sSFR versus look-back time plane (Figure 6). In all bins, even in the most massive bin, we find galaxies of all types: star-forming, transitioning, quiescent, and rejuvenating. Furthermore, galaxies transition from the star-forming to the quiescent region at different epochs and over a range of timescales (Figures 8 and 13). Consequently, there is a wide range of formation redshifts (Figure 11) and star formation timescales (Figure 12) at fixed stellar mass.

Despite this diversity, there are general trends. Unsurprisingly, there is a relation between the formation redshift z _f and the star formation timescale τ_SF: early-forming galaxies form their stars on shorter timescales (Figure 12 and Equation (9)). This relation only depends weakly on stellar mass. Looking at the median evolution of the sSFR with cosmic time (Figure 9), we find that the sSFRs of star-forming galaxies follow $\propto {\left(1+z\right)}^{2.5}$, consistent with the specific accretion rate of dark matter halos and the SFMS measured independently at different epochs. This is already a first indication that the galaxies' SFMS describes an evolutionary track along which individual galaxies evolve on average.

The similarity between the sSFRs of star-forming galaxies and the specific accretion rate of dark matter halos is expected from theoretical studies. Based on numerical simulations, Dekel et al. (2013) show that specific gas inflow rate and sSFR of galaxies scale with the specific dark matter accretion rate of their halos, as expected from analytical considerations using the extended Press−Schechter approximation (Bond et al. 1991). Furthermore, simple empirical models that link the growth of galaxies to their dark matter halos are able to successfully describe galaxies at low and high z (Behroozi et al. 2013; Moster et al. 2013; Rodríguez-Puebla et al. 2017; Tacchella et al. 2018). Importantly, this rule for star-forming galaxies ($\mathrm{sSFR}\propto {\left(1+z\right)}^{2.5}$) does not directly translate into a simple stellar-to-halo mass relation at earlier cosmic time because the star formation efficiency (i.e., the slope of the stellar-to-halo mass relation) possibly depends itself on halo mass and redshifts (see, e.g., Figure 3 in Rodríguez-Puebla et al. 2016).

Focusing on quiescent galaxies, we find that they decouple on average from the SFMS early on (z > 3) and quench around z_quench ≈ 1–2 (Figure 9). This is consistent with recent findings of a quiescent galaxy population at z > 3 (Gobat et al. 2012; Hill et al. 2017; Kubo et al. 2018; Schreiber et al. 2018; Girelli et al. 2019; Tanaka et al. 2019; Carnall et al. 2020; D'Eugenio et al. 2020; Forrest et al. 2020a, 2020b; Saracco et al. 2020; Shahidi et al. 2020; Valentino et al. 2020; Santini et al. 2021). In particular, Saracco et al. (2020) report a star formation timescale of a few hundred Myr and a short quenching timescale (≲150 Myr) for a quiescent galaxy at z = 3.352, which is consistent with our estimates. We find that the majority of those objects evolve passively to lower redshifts, but there are a few cases where rejuvenation plays a role (Section 6.5). Consistently, we find that some galaxies are maximally old, meaning that they have formation redshifts of z _f ∼ 10 (Figure 11). These superold systems are the most massive galaxies in our sample ($\mathrm{log}({M}_{\star }/{M}_{\odot })\gt 11.2$) and also live in the highest-mass halos ($\mathrm{log}({M}_{{\rm{h}}}/{M}_{\odot })\gt 13.0$)—progenitors of today's slow rotators in the centers of clusters (e.g., Cappellari 2016). Clearly, these early formation redshifts are, of course, exciting news for the James Webb Space Telescope (JWST), which will be able to shed new light on the formation of those objects. These objects should be easily detected and characterized in detail, because we estimate them to have SFRs between 10 and 1000 M_⊙ yr⁻¹ (Figure 7). However, it is also possible that these early-forming galaxies were still in pieces at these high redshifts, making it harder for JWST and other high-z surveys to discover them.

We find that the formation redshift depends on stellar mass, but only weakly (Figure 11), which is consistent with previous results for galaxies at this and earlier epochs (Wu et al. 2018b; Carnall et al. 2019b; Estrada-Carpenter et al. 2019; Ferreras et al. 2019; Morishita et al. 2019). As mentioned in the Introduction, local galaxies show a tighter relation between the stellar age and the stellar velocity dispersion than between the stellar age and stellar mass (Graves et al. 2010). As previously noted by Carnall et al. (2019b), TNG100 is not consistent with the stellar mass trend and, overall, predicts lower formation redshifts (younger ages) than seen. It seems that TNG100 is missing early-forming quiescent objects. This is of interest since it could point toward diverse pathways to quiescence—more diverse than what currently is implemented in TNG100, where the black hole mass is the main determinant of whether a galaxy quenches (e.g., Terrazas et al. 2020). On the other hand, as outlined in Section 2.4, effects of cosmic variance, sample selection, and resolution could also contribute to this absence of such galaxies.

Although all the high-mass systems have formed early, the diversity increases toward lower masses. Interestingly, several of our low-mass, quiescent galaxies ($\mathrm{log}({M}_{\star }/{M}_{\odot })\,=10.0-10.5$), which have not been probed by the literature previously, have formation redshifts of z _f ≈ 3−8, similar to what we find in TNG100.

As highlighted in the Introduction, a fundamental question concerning the SFMS is how galaxies evolve about this scaling relation. Are galaxies moving about this relation on short or on long timescales relative to the age of the universe? In the most extreme case, where galaxies actually do not evolve about the SFMS but follow independent (lognormal) SFH trajectories, the SFMS itself has no deeper implication and equal-mass galaxies on the SFMS need not form an evolutionary cohort (Gladders et al. 2013; Abramson et al. 2016). In such a scenario, we would expect that the scatter of the SFMS shows a strong age gradient (see Figures 9 and 10 in Abramson et al. 2016), because galaxies at the upper envelope of the SFMS can only have formed recently. On the other hand, if galaxies fluctuate about the SFMS ridgeline on short timescales, similarly massed galaxies have similar SFHs and the evolution of the SFMS can be used to track the SFH of galaxies.

We first focus on tracking how galaxies evolve in the SFR − M_⋆ plane. Figure 15 plots SFH tracks in the space of SFR − M_⋆ (top panels) and Δ_MS − M_⋆. Δ_MS is the log-distance from the SFMS, i.e., ${{\rm{\Delta }}}_{\mathrm{MS}}=\mathrm{log}(\mathrm{sSFR}/{\mathrm{sSFR}}_{\mathrm{MS}})$, where sSFR_MS is the sSFR of the SFMS at a given look-back time/redshift and stellar mass of the galaxy' history. We use the SFMS prescription of Leja et al. (2021), which is based on the 3D-HST photometry and stellar population analysis presented in Leja et al. (2019b), meaning that the SFMS is estimated from direct measurements of SFR and M_⋆ of galaxies at z = 0.5 − 3.0. The left, middle, and right panels show the tracks of quiescent, star-forming, and rejuvenating galaxies, respectively. The classification is based on Equation (2). Rejuvenating galaxies are galaxies that were quiescent in the past and have ${\mathscr{D}}\gt 1/20$ at the epoch of observation. Each track consists of a range of look-back times t_lb, ranging back to t_lb = 4.5 Gyr, which corresponds to z ≈ 3. It is important to keep in mind when interpreting these tracks that they do not necessarily reflect the SFH of a single galaxy, since galaxies can merge and the resulting stellar population mixes in situ and ex situ formed stars. We discuss this further in Section 6.6.

Zoom In Zoom Out Reset image size

Figure 15 shows that quiescent galaxies have similar tracks. This is not too surprising since when galaxies significantly reduce their SFR they will not move further to the right in this figure. Obviously, mergers can increase M_⋆ of quiescent galaxies, but our SFH measurements are insensitive to this: we are only able to infer the combined SFH of the main progenitor and all of its previously merged subcomponents. We find that by z ∼ 3 most of the quiescent galaxies are on the SFMS, consistent with what we found in Figure 9. We will discuss these tracks and quenching more generally in Section 6.4. Similarly, we will follow up on the rejuvenating galaxies in Section 6.5.

We now focus on star-forming galaxies (middle panel of Figure 15), which unfortunately only includes a limited number of galaxies in our sample (31 galaxies, out of which 25 are at z_obs = 0.6–1.0 and are displayed in the figure). The key point to notice is that the star-forming galaxies show significant movement relative to the SFMS in recent times. Some galaxies moved more than 1 dex in the past 100 Myr. Galaxies on average cross the SFMS twice, while only one galaxy never crosses the SFMS and three galaxies cross the SFMS four times. This also holds when considering the uncertainties on the inferred SFHs. Galaxies move from above to below the SFMS on a timescale of ${t}_{\mathrm{cross}}={0.33}_{-0.28}^{+1.74}\,\mathrm{Gyr}$, while they move from below to above on a timescale of ${t}_{\mathrm{cross}}\,={4.02}_{-3.58}^{+1.08}\,\mathrm{Gyr}$. This only considers the first crossing after the observations. Later oscillations are more difficult to estimate because of our time bins, which leads to oscillation timescales of 1–2 Gyr. We also notice that galaxies slightly lie above the SFMS at large look-back times (see also quiescent galaxy panels on the left), which could be related to a systematic difference between the SFMS prescription and this work, or the accretion of smaller galaxies with time, leading to a bias high in SFR relative to the SFMS. ²¹

To further study this, we look explicitly for an age gradient across the SFMS to understand whether there is some long-term correlation of whether galaxies are above or below the SFMS ridgeline. Figure 16 shows the formation redshift z _f (top panel), the mass-weighted age t_age (middle panel), and the ratio of the SFR measured over the previous 10 Myr and 1 Gyr (bottom panel) as a function of the distance from the SFMS Δ_MS. The vertical black line marks the SFMS, while the gray band indicates its scatter of ∼0.3 dex. We plot star-forming (blue circles), transitioning (green triangles), quiescent (red squares), and rejuvenating (orange stars) galaxies, classified by our mass-doubling criterion specified in Equation (2).

Zoom In Zoom Out Reset image size

The top panel of Figure 16 shows that there is no significant gradient in z _f with Δ_MS. We notice that the galaxies above the SFMS (Δ_MS > 0.5 dex) have all formed recently (z _f < 1.5), which can be naturally explained by a burst of star formation adding significantly stellar mass, leading to a decrease of z _f . However, recent formation redshifts z _f are also not uncommon on the SFMS, which shows a large diversity of z _f . This is also true when we study this trend at fixed stellar mass.

The middle panel of Figure 16 indicates a weak, not significant trend of t_age with Δ_MS. This holds also when studying this trend at fixed stellar mass. Again, galaxies above the SFMS are the youngest galaxies in our sample, indicating that this recent burst of star formation significantly adds stellar mass. These results point toward galaxy oscillations about the SFMS on timescales significantly shorter than the age of the universe.

The bottom panel of Figure 16 shows the logarithm of the ratio of the SFR measured over short (10 Myr) and over long (1 Gyr) timescales, i.e., a positive number indicates that the SFR has increased, while a negative number indicates that the SFR has decreased over the past 1 Gyr. There is no significant trend across the SFMS for star-forming galaxies; only galaxies above the SFMS have recently increased their SFRs significantly. As expected, transition and rejuvenating galaxies—roughly at the same distance from the SFMS—have significantly different ratios, consistent with the idea that transition galaxies are reducing the SFRs, while rejuvenating galaxies are increasing their SFRs. Quiescent galaxies typically have a negative ratio, indicating decreasing SFHs.

In summary, Figure 15 (together with Figure 9) indicates that galaxies naturally move along the SFMS. Importantly, we find clear indications of short-term variations in the SFH of star-forming galaxies, indicating that at least some of these galaxies move on rather short timescales (few hundred Myr) about the SFMS. Figure 16 shows no strong trends of the formation redshift or galaxy age across the SFMS; only "star-bursting" galaxies show younger ages. It is important to stress that we focus here on z ≈ 0.8. We expect that these gradients in z _f and t_age across the SFMS are expected to increase in strength toward lower redshifts, because environmental processes will play a larger role in regulating star formation (e.g., Poggianti et al. 2006; Sánchez-Blázquez et al. 2009; Jablonka et al. 2013; Cantale et al. 2016; Webb et al. 2020; van der Burg et al. 2020). An indication for this is, for example, the dependence of the clustering strength on the distance from the SFMS (Berti et al. 2021).

These results are consistent with numerical simulations, which show that galaxies oscillate about the SFMS (Tacchella et al. 2016; Matthee & Schaye 2019). On which timescales galaxies oscillate is still uncertain observationally, but a fruitful avenue to characterize this is the use of the temporal power spectrum density (Caplar & Tacchella 2019). This framework allows us to assess the importance of SFR fluctuations on a wide range of timescales, and it can tell us about different physical processes regulating star formation that act on different timescales, ranging from molecular clouds to gas accretion (Tacchella et al. 2020). Interestingly, the predictions for the power spectrum density from galaxy formation models (numerical as well as semianalytical models) vary widely, highlighting that this is an interesting space to constrain those models (Iyer et al. 2020). Future, more detailed measurements of the power spectrum density via different SFR and stellar age indicators will be important to assess how galaxies evolve about the SFMS (e.g., Guo et al. 2016; Broussard et al. 2019; Caplar & Tacchella 2019; Emami et al. 2019; Faisst et al. 2019; Wang & Lilly 2020).

It is challenging for the "grow-and-quench" framework to explain a large population of galaxies that are extremely young. To this end, Dressler et al. (2016, 2018) searched for "late bloomers," which are massive ($\mathrm{log}({M}_{\star }/{M}_{\odot })\gt 10.0$) galaxies at z < 1 that formed the majority of their stars within ∼2 Gyr of the epoch of observation. Dressler et al. (2018) found that late bloomers account for ∼ 20% of z ∼ 0.6 galaxies with masses of the modern (z ≈ 0) Milky Way, with a moderate dependence on mass. They conclude that this strongly contradicts paradigms in which galaxies are thought to simply grow along the SFMS and quench. Mamon et al. (2020), using SDSS spectra, have also looked for such late-bloomer galaxies in the local universe (but called them very young galaxies (VYGs)), finding that they are common (∼50%) at $\mathrm{log}({M}_{\star }/{M}_{\odot })\approx 8.0$ but rare at high stellar masses (∼0.1% at $\mathrm{log}({M}_{\star }/{M}_{\odot })\approx 11.5$).

Figure 17 shows the fraction of stellar mass formed in the previous 2 Gyr as a function of M_⋆. We focus on $\mathrm{log}({M}_{\star }/{M}_{\odot })\gt 11$ since this is the mass range where our sample is unbiased (see Section 2.1.1). Our measurements are indicated with circles (color-coded by Δ_MS), while the green crosses show the TNG100 data processed in the same way as our data. We find the expected correlation: galaxies above the SFMS have a higher fraction of stellar mass formed in the previous 2 Gyr.

Zoom In Zoom Out Reset image size

Only a few of our galaxies are late bloomers (are above the horizontal dashed black line), which is consistent with our mass-weighted age estimates (Figure 10). We find an observed late-bloomer fraction of ${f}_{\mathrm{LB}}={3.6}_{-1.8}^{+8.9} \%$ in the mass range of $\mathrm{log}({M}_{\star }/{M}_{\odot })=11.0-11.5$. In order to understand the nature of these late bloomers, we correlate t₂₀ (look-back time by when 20% of the stellar mass was formed) and τ_SF with the fraction of mass formed in the previous 2 Gyr (y-axis of Figure 17). We find that t₂₀ smoothly decreases with this mass fraction, which indicates that late bloomers seem to have started their star formation more recently than nonlate bloomers. Interestingly, late bloomers seem to show some diversity in τ_SF, i.e., there are objects with τ_SF ∼ 1.5 Gyr as well as objects with τ_SF > 3 Gyr. This indicates that late bloomers can form through both a recent burst of star formation and a more continuously increasing SFH. This is also confirmed when looking at the SFHs themselves (i.e., SFR vs. time; Figure 7).

In comparison with TNG100, we do not find a significantly different distribution in Figure 17, but there are no late bloomers at all in TNG. The fraction is consistent with 0, i.e., ${f}_{\mathrm{LB}}={0.0}_{-0.0}^{+1.3} \%$. Dressler et al. (2018) quote at 〈z〉 = 0.6 f_LB ≈ 5%−10% for galaxies with $\mathrm{log}({M}_{\star }/{M}_{\odot })=11.0$, which is slightly higher than what we find in our observations and TNG100. The weak disagreement between Dressler et al. (2018) and our measurement could point toward an overestimate of the late-bloomer fraction in Dressler et al. (2018)—in particular toward lower M_⋆—caused by the assumed prior in their work. They fit for the mass formed in N fixed time bins, which is not agnostic about the shape of SFR(t) but instead prefers rising SFHs and high instantaneous sSFRs (Leja et al. 2019a). This can be understood intuitively by imagining a series of random draws for each of the N time bins that add up to a fixed total mass by construction (typically, this total mass is loosely constrained by the normalization of the SED). Most of the SFHs thus constructed will have significant mass in only one or two time bins, which naturally results in "bursty" SFHs and in SFHs for star-forming galaxies that are rising. These effects could lead to an overestimation of the late-bloomer fraction. A separate issue is the uncertainty in the assumed dust attenuation law. In order to estimate the stellar mass formed in the previous 2 Gyr (in particular, the age range of 0.5−2.0 Gyr), accurate estimates of the far-UV and near-UV luminosities are needed (see discussion in Appendix B). This means performing an accurate dust correction. We know that the dust attenuation law varies from galaxy to galaxy (Kriek & Conroy 2013; Reddy et al. 2015; Salim & Narayanan 2020; Shivaei et al. 2020), as expected from theoretical models (Smith et al. 2015; Narayanan et al. 2018), and there is a need for a flexible attenuation law when reproducing the physical properties of a large variety of objects (Lo Faro et al. 2017). We marginalize over different slopes of the dust attenuation law, while Dressler et al. (2018) assume a fixed Calzetti et al. (2000) law but allow for different dust amounts in each SFH time bin. These different approaches could also lead to a difference in f_LB.

Quiescent galaxies decouple from the median SFH of star-forming galaxies and the typical accretion rate of dark matter (Figure 9), thereby falling significantly below the SFMS (>2 dex). We show that galaxies transition over a wide range of timescales, from a few tens of Myr to a few Gyr (Figure 13). The typical quenching timescale in our sample is ${1.1}_{-0.7}^{+1.2}\,\mathrm{Gyr}$. About 25% of our galaxies quench on short timescales with τ_quench < 500 Myr. This is consistent with complementary studies of post-starburst galaxies, which infer that such fast-quenching galaxies are responsible for about 20%–50% at 0.5 < z < 1.5 (Belli et al. 2019; Wild et al. 2020).

Focusing on the left panel of Figure 15, we find that galaxies quench over a wide range of stellar masses. The lowest-mass systems quench at $\mathrm{log}({M}_{\star }/{M}_{\odot })\approx 10.3$, while the highest-mass systems have $\mathrm{log}({M}_{\star }/{M}_{\odot })\approx 11.6$. As mentioned before, there is a high probability, in particular for the most massive systems, that these galaxies have merged and consist of multiple progenitors. This means that the quenching mass range in reality might be narrower. For example, in TNG100, galaxies typically quench in the range $\mathrm{log}({M}_{\star }/{M}_{\odot })\approx 10.5\mbox{--}11.0$. Above this mass range, the ex situ mass fraction increases (Tacchella et al. 2019). Furthermore, these tracks are nearly vertical, indicating little mass growth during quenching. Specifically, we measure a relative stellar mass growth (in comparison with final stellar mass) in the "green valley" of ${f}_{{\rm{\Delta }}{M}_{\star }}={7.9}_{-7.6}^{+25.0} \%$, which, by definition, is correlated with the quenching timescale.

In the literature, several studies highlight the need for two quenching paths—slow and fast—in both observations (e.g., Yesuf et al. 2014; Barro et al. 2016; Wu et al. 2018a; Belli et al. 2019; Carnall et al. 2019b) and simulations (e.g., Rodríguez Montero et al. 2019). We do not find a bimodal distribution of quenching timescales in Figure 13, which could be related to the rather large uncertainties of τ_quench. At face value, our quenching timescale distribution and the wide quenching mass range can be explained by a combination of a quenching mechanism with a "gate keeper" (Tacchella et al. 2016): while the galaxy grows in stellar mass on the SFMS, the galaxy tries to quench and depart the SFMS via stellar and black hole feedback, which itself can be activated through mergers and violent gas instabilities (Hernquist 1989; Zolotov et al. 2015). However, as long as the replenishment time of fresh gas is short enough, i.e., there is enough fuel available for sustaining star formation, the galaxy is forced to keep forming stars and stay on the SFMS. Only when the halo is hot enough (the halo mass plays the role of the aforementioned "gate keeper") can fresh gas not be supplied after another quenching attempt, and the galaxy is able to quench. A result of this combination of a quenching mechanism with a "gate keeper" is a wide range of quenching timescales and quenching masses since the exact timing of when the quenching mechanism is activated (e.g., a merger takes place) and the timing of when the critical halo mass is reached differ from galaxy to galaxy. This is also consistent with Figure 14, where the most massive halos (presumably the halos with the most shock-heated gas) host galaxies that were able to quench rapidly.

Consistently, Chen et al. (2020) argue that the wide quenching mass range can be explained by varying black halo masses that beat against halo masses that also vary across this mass range. This picture accounts nicely for the width and stable mean mass of the quenching mass range, and it naturally accommodates a "fast" quenching channel precipitated by gas-rich mergers. However, this work has only little to say about how long it takes to quench, since it focuses mainly on preventive feedback, which has a natural timescale of the order of the halo dynamical timescale. Quenching timescales should be shorter at early times since halos are denser, but we do not find a such trend here. Ejective feedback might lead to more variation in the quenching timescale.

An additional observational avenue to constrain the quenching mechanism(s) is to study the spatial progression of quenching within galaxies. The first steps in this direction have been done in high-redshift galaxies, both for star-forming galaxies by studying the distribution of star formation (Wuyts et al. 2013; Tacchella et al. 2015a; Nelson et al. 2016; Abdurro'uf 2018; Morselli et al. 2019) and for quiescent galaxies by studying the stellar populations on spatially resolved scales (Akhshik et al. 2020; Jafariyazani et al. 2020).

As mentioned above, rejuvenation in our sample is uncommon. We use a rather stringent criterion, where galaxies are rejuvenating if they were quiescent in the past and are either in the star-forming or the transition region at the time of observation (based on the criterion in Equation (2)). Specifically, only two galaxies in our sample of 161 galaxies rejuvenate back to be star-forming, corresponding to 1.2%. There are seven other galaxies (4.3%) that rejuvenate, but only to the transition region. Interestingly, as highlighted in Figure 16, these rejuvenating galaxies all have high formation redshifts (z _f ≈ 5−8) and are rather old (half-mass age of t_age > 5 Gyr), which indicates that only little mass has formed in these recent rejuvenation events (<10%). Also, we do not see any indication of AGNs in these systems (f_AGN < 1%).

Consistent with high formation redshift and old age, we find that rejuvenated galaxies prefer to lie in massive halos. In particular, we find that the fraction of rejuvenating galaxies at $\mathrm{log}({M}_{{\rm{h}}}/{M}_{\odot })\lt 13.0$ is only f ≈ 4%, while at $\mathrm{log}({M}_{{\rm{h}}}/{M}_{\odot })\gt 13.0$ it is f ≈ 30%. This could signal a cooling flow failure as seen in local massive ellipticals such as Perseus A (NGC 1275; e.g., Canning et al. 2014) and Abell 1795 (e.g., Ehlert et al. 2015), which are massive central cluster galaxies but are hosting star formation. Galaxy–galaxy merger, which can bring in new fuel for star formation, might also play a role, but because of the high relative velocity of galaxies in these high-density environments (i.e., (proto)clusters), mergers should be less common.

The rarity of rejuvenating galaxies is consistent with theoretical expectations (Trayford et al. 2016; Pandya et al. 2017) and other observational measurements (Belli et al. 2017; Chauke et al. 2019; Akhshik et al. 2021). TNG100 predicts a rejuvenation fraction above $\mathrm{log}({M}_{\star }/{M}_{\odot })\gt 10.5$ of 1.6%. In Chauke et al. (2019), the authors infer that 16% of the z ∼ 0.8 quiescent population has experienced rejuvenation events in the redshift range 0.7 < z < 1.5 after reaching quiescence at some earlier time. This is a different criterion than the one we are using. With their criterion, a double-peaked SFH (separated by an SFR ≈ 0) is considered rejuvenated (see Figure 4 of their work). This raises, of course, the interesting question of the exact definition of rejuvenation. Since it is challenging to recover a short burst of star formation several Gyr back in time and our prior weights against such a behavior (Leja et al. 2019a), our definition of rejuvenation is a conservative one and focuses on the current epoch only.

Furthermore, we believe that the quoted 16% in Chauke et al. (2019) is an overestimation of rejuvenation events. First, a double-peaked SFH (with quiescence in between) does not imply rejuvenation. It could also be that two galaxies, each with a different formation redshift, have merged. The resulting merger remnant will have a double-peaked SFH, but no actual rejuvenation took place. Second, Chauke et al. (2019) fit for the mass formed in N fixed time bins, which leads to a prior in SFHs that tends to form the majority of the mass in one or two time bins (see Section 6.3). Nevertheless, we agree on their main conclusion that these rejuvenation events only contribute little to the total stellar mass and that they are not an important evolutionary channel when considering the growth of the quiescent galaxy population.

Clearly, it would be of great interest to study the morphology of these rejuvenating galaxies, which we will do in the future. These galaxies might be consistent with objects such as presented in Mancini et al. (2019), where the SFHs have been measured for bulge and disk components separately for 10 galaxies at 0.45 < z < 1. These authors have argued that the bending of the SFMS is due significantly to quiescent galaxies acquiring star-forming disks, which rejuvenate the systems.

A limitation in our study is the sample size of 161 galaxies, which makes it difficult to split the sample by physical properties, such as mass, star formation, or formation redshift. Large samples of both star-forming and quiescent galaxies will allow us, for example, to shed more light on the correlation between quenching timescale and other galaxy properties. Furthermore, our sample only includes a few star-forming galaxies; increasing the sample size of star-forming galaxies will help us to understand how galaxies evolve about the SFMS in more detail.

In the future we are planning to expand the sample size by including galaxies from other spectroscopic surveys, such as LEGA-C (van der Wel et al. 2016) and VANDELS (McLure et al. 2018; Pentericci et al. 2018). In addition, for increasing the number of objects, the inclusion of these surveys will help to span a wider range of epochs of observation. This will be particularly useful for understanding how galaxies evolve post-quenching, including both major and minor mergers. As highlighted above, a major limitation when studying SFHs is the uncertainty related to mixing several galaxies into a single galaxy and misinterpreting the measured SFH as the evolutionary path of one single object. Comparing the SFH of in situ and ex situ stars in TNG100 does not actually show a significant difference, which might alleviate some of these concerns. It seems that galaxies, which will later be accreted by the main progenitor, have similar SFHs to the main progenitor itself, i.e., there is some overarching coherence scale for the star formation in galaxies in a certain environment. Nevertheless, a more detailed investigation is needed.

Understanding how representative a measured SFH is in comparison to the evolutionary path of the main progenitor is an open question. Probing galaxies at different epochs could help with this (but see Abramson et al. 2016). Furthermore, connecting the measured SFHs to the morphology of the galaxies is also of great interest (Belli et al. 2015; Fagioli et al. 2016; Almaini et al. 2017; Faisst et al. 2017; Williams et al. 2017; Maltby et al. 2018; Wu et al. 2018a; Estrada-Carpenter et al. 2020; Suess et al. 2020)—something we will focus on in an upcoming publication. This is of interest not only for understanding the evolution of galaxy morphology with epoch but also because the morphology of galaxies can give us clues about the occurrence of mergers.

As described in Section 3, the high-resolution empirical stellar library MILES has a comparable resolution to our observed galaxy spectra. In particular, the observed spectra of our low-z galaxies have a lower resolution than the MILES templates. Therefore, the MILES templates need to be matched to the LSF of our observations by convolving them with a kernel. This instrumental LSF of the observation describes the observed shape of an intrinsically very narrow spectral line at the wavelength λ, due to purely instrumental broadening, including the effect of integrating over the detector pixels.

We construct the LSF for each HALO7D/CANDELS field separately. Specifically, we measure the Gaussian width σ_λ as a function of wavelength from fitting Gaussians to arc lamps and night skylines for all the slits that fell in each field. We then fit a quadratic function to all the width measurements as a function of wavelength. These LSFs are shown in Figure 18. The solid lines show the quadratic LSFs for each individual field. We measure a typical width of σ_λ ≈ 1.5−2.0 Å and find little difference from field to field.

Zoom In Zoom Out Reset image size

Since we assume that the LSF is Gaussian, the kernel with which the templates need to be convolved is also Gaussian with a width σ_diff. Specifically, we determine the broadening σ_diff of the rest-frame MILES library spectra required to match the observed LSF by

$\begin{eqnarray}&&{\sigma }_{\mathrm{diff}}^{2}={\sigma }_{{\rm{v}}}^{2}-{\sigma }_{{\rm{v}},\mathrm{MILES}}^{2},\end{eqnarray} \tag{ A1 }$

where σ_v = c × σ_λ /λ_obs and σ_v,MILES = c × FWHM_MILES/2.355/λ_rest with FWHM_MILES = 2.54 Å. As noted in the main text, σ_diff becomes negligibly small for high-redshift galaxies in our sample. This means that the observed LSF is dominated by the physical stellar velocity dispersion—which is fitted—and the "uncorrectable" difference between the MILES and instrumental LSF is not important to the final total LSF. Finally, the effects of the LSF on the inferred stellar population parameters are negligible.

As mentioned in the main text, we assume a nonparametric SFH for fitting the photometric and spectroscopic data. We call it nonparametric, because we do not assume any particular parametric form for the shape of the SFH. Nevertheless, the number of time bins (N_SFH) describing the SFH is a parameter in our prescription. Our fiducial assumption is N_SFH = 10. We have explored fewer and more bins. Obviously, more bins increase the time resolution, which is important for determining, among other quantities, the star formation and quenching timescales.

Figure 22 compares the result of some of our key quantities (stellar mass, sSFR, dust attenuation, age, and star formation timescale τ_SF) for our fiducial N_SFH = 10 run with our N_SFH = 14 run. The 1-to-1 line is indicated as a solid black line. This figure shows that these two runs deliver consistent results, which is comforting. Furthermore, this is consistent with the analysis by Leja et al. (2019a), who explored varying the number of time bins between N_SFH = 4 and 14 and showed that the results of the mock analysis are largely insensitive to the number of bins as long as N_SFH ≳ 4. We choose our fiducial number of time bins to be 10 because the N_SFH = 14 run leads to nonconvergence of 30 galaxies within a reasonable amount of time (i.e., 14 CPU days).

Zoom In Zoom Out Reset image size

We highlight here some of the differences that arise when assuming a parametric SFH instead of our fiducial nonparametric SFH. We assume a delayed-tau model for the parametric SFH as defined in Equation (3). As shown in Appendix D, assuming a parametric SFH imposes strong priors on several key quantities, including age (formation redshift), star formation timescale, and quenching timescale (Figures 25 and 26).

Figure 23 shows the posteriors for key quantities (M_⋆, sSFR, ${\hat{\tau }}_{\mathrm{dust},2}$, t_age, and τ_SF) assuming a nonparametric (red) and a parametric (blue) SFH. These are the posteriors of the same galaxy (ID 16490) as shown in Figure 5. The posterior SFHs are significantly different in two cases: while the nonparametric and parametric SFHs are similar in the last ∼1.5 Gyr (where the S/N is the highest), they diverge at earlier times. In particular, the parametric SFH results in a much higher peak SFR of >80 M_⋆ yr⁻¹ than the nonparametric one (peak SFR of ∼30 M_⋆ yr⁻¹). Given by its parametric form, the parametric SFH declines rapidly at earlier times and is consistent with 0 at a look-back time beyond 2 Gyr. This leads to a short star formation timescale, consistent with Figure 25, which shows that long star formation timescales are omitted by this prior. On the other hand, the nonparametric SFH peaks and declines more slowly toward earlier cosmic times. This results in older ages, larger M_⋆, and lower sSFR in the nonparametric SFH case than in the parametric SFH case.

Zoom In Zoom Out Reset image size

Figure 24 compares the key quantities of the nonparametric SFH run (x-axis) with the ones of the parametric SFH run (y-axis). Consistent with Figure 23, we find a significant difference for M_⋆, sSFR, and t_age. There are no significant differences for the dust opacity ${\hat{\tau }}_{\mathrm{dust},2}$. This is overall consistent with Carnall et al. (2019a), who showed that parametric SFH models (including delayed-tau models) bias measurements such as SFR and stellar masses by imposing a strong prior preference for young stellar populations. They conclude—and we support it with our analysis here—that this makes it challenging for parametric SFH models to understand mass assembly in galaxies. Nonparametric SFH models are better suited because they are less biased and return more accurate errors than parametric SFHs (Leja et al. 2019a).

Zoom In Zoom Out Reset image size