The Typical Massive Quiescent Galaxy at z ∼ 3 is a Post-starburst

Timing the formation and quenching of the most massive (M_⋆ > 10¹¹M_⊙) passively evolving galaxies (PEGs) is the subject of intense debate. Several studies kept unveiling their existence at progressively increasing lookback times, placing their formation at z > 3–4 (Gobat et al. 2012; Glazebrook et al. 2017; Forrest et al. 2020; Valentino et al. 2020). Reproducing the observed number density of such galaxies at all epochs is a compelling concern of current galaxy evolution models in cosmological simulations, as the relative importance of the quenching mechanisms is not yet clear (Man & Belli 2018). To this end, several photometric samples of high-z PEGs exist (Straatman et al. 2014; Merlin et al. 2018), yet spectroscopic confirmation is still needed to precisely assess their redshifts and the degree of contamination by star-forming interlopers, especially at epochs when PEGs might be just starting to emerge. At z ∼ 3 the chances of collecting statistically meaningful samples of PEG spectra are hampered by their rarity and lack of prominent emission lines. Therefore, the assembly of sizeable samples of quiescent spectra is essential to improve upon the photometric age/dust-attenuation estimation for the bulk of their stellar populations, to draw relative comparisons of their spectral properties with respect to lower-redshift massive PEGs and to reveal any residual emission line that might be linked to star formation. The latter point provides useful information on the availability of unstable, cold molecular gas that is, itself, a poorly constrained quantity at this epoch. In this Letter we report on the average stellar population properties in a sizeable sample of nine spectroscopically confirmed 2.4 < z < 3.3 quiescent galaxies. We constrain their global residual star formation from the average [O ii] emission, showing implications on their average gas fraction. We assume a ΛCDM cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, ${{\rm{\Omega }}}_{M}=0.27,{{\rm{\Omega }}}_{{\rm{\Lambda }}}=0.73$, and a Salpeter initial mass function (IMF). Magnitudes are given in the AB photometric system.

We selected galaxies from a parent sample of more than 50 reliable PEGs candidates with K_tot < 22.5 in the McCracken et al. (2010) catalog with photometric redshifts between 2.5 ≤ z_phot ≤ 3.5 in the COSMOS field. A first selection was done using the BzK criterion (Daddi et al. 2004) selecting passive BzK galaxies plus objects formally classified as star-forming BzK but having a signal-to-noise ratio (S/N) < 5 in the B and z bands, to retain passive galaxies as they become fainter in such bands with increasing redshift and decreasing mass. We then selected from this sample UVJ quiescent galaxies (Williams et al. 2009), exploiting photometric redshifts calibrated for high-z PEGs (Onodera et al. 2012; Strazzullo et al. 2015). To minimize contamination from dusty star-forming galaxies, objects in McCracken et al. (2010) with Spitzer/MIPS 24 μm S/N > 4 were discarded, except for galaxies with high-confidence passive spectral energy distributions (SEDs). Objects having SED fits to optical–NIR broadband photometry consistent with dusty star-forming solutions and contaminated photometry in all images were also discarded. Among the most massive bona fide passive candidates, 10 galaxies were targeted for Hubble Space Telescope (HST) WFC3/IR G141 near-IR¹⁵ observations: 9 with H_AB < 22 (M_⋆ > 1.1 × 10¹¹M_⊙) plus 1 robust candidate with H_AB = 22.9 (M_⋆ = 8 × 10¹⁰M_⊙) among the highest-z objects. Scheduling from 1 to 3 orbits (17 in total) according to each source's flux provided low-resolution spectra (R = 130) with a mean S/N ∼ 4.5 per resolution element for each target and ancillary F160W imaging (HST GO 15229 program). The data were reduced and decontaminated from neighboring sources' overlapping spectra by means of the grizli software package (Brammer 2018).¹⁶ For each target a 2D spectrum was produced and drizzled to a scale of 0.8 times the native pixel size. The 1D spectra were then optimally extracted (Horne 1986) and fitted from 11000 to 16900 Å with Bruzual & Charlot (2003, hereafter BC03) composite stellar populations derived using different star formation histories (SFHs): allowing for a constant, exponentially declining, delayed exponentially declining, and a truncated SFH, where the star formation rate (SFR) drops to 0 after a timescale between τ = 0.001 − 1.0 Gyr. Ages were allowed to vary between t = 100 Myr and 10 Gyr, from redshift 5 to 0. The Calzetti et al. (2000) extinction law was adopted with A_V varying between 0 and 8 mag. BC03 templates were broadened to the grism resolution using a Gaussian kernel whose FWHM was derived via a Gaussian fit of the 2D light profile of each source.

The broadband photometry (Laigle et al. 2016) for each galaxy was also fitted imposing the grism-derived redshifts. The addition of the photometric information mostly leads to narrowing down the range of dust extinction values at any confidence level. The detailed analysis of individual galaxies will be presented in a forthcoming paper.

2.1. Quiescence of Individual Targets

To fully exploit the sample size we created an average spectrum. We scaled the wavelength vector of each galaxy to z_m and interpolated their fluxes and error spectra on a 5 Å rest-frame grid, increasing the error in each new pixel by the square root of the width ratio between the old and the new spectral bin, to account for the introduced noise correlation. Each error spectrum was interpolated in quadrature. Each spectrum was normalized to its average flux between λ_rest = 3800 and 3900 Å in order to select a wavelength range covered for all galaxies. We then created a stacked spectrum as the inverse variance weighted mean of the individual fluxes in each wavelength bin. The final error spectrum was computed as the error on the flux-weighted mean in each pixel. A jackknife resampling provides consistent results. The median z_spec uncertainty from our individual galaxies is dz ∼ 0.006 at 68% confidence. This does not affect the final spectral resolution, being three times smaller than the broadening from galaxy sizes. The average spectrum was rescaled to match the average broadband photometry of the sample at λ_rest = 3800–3900 Å, in turn obtained by stacking the individual galaxy best-fit SEDs (as done for the spectra). Figure 1 shows the average spectrum with clear Balmer absorption lines from Hη to Hγ, the [O ii] (unresolved) doublet at 3727 Å, and the Fe i absorption line. With a mean S/N ∼ 13 per pixel at 4000 Å this is, to date, one of the highest S/N spectra available for high-z PEGs and the first spanning a large spectral range around key spectral breaks.

Zoom In Zoom Out Reset image size

Table 1.  Main Properties of our Sample Derived from the Stacked Spectrum (See the Text)

〈z〉	M⋆	t50	AV	Dn4000	DB	EW[O ii]	F[O ii]	SFR
	(1011 M⊙)	(Gyr)	(mag)			(Å)	(erg s−1 cm−2)	(M⊙ yr−1)
2.808	1.8 ± 0.8	${0.30}_{-0.05}^{+0.20}$	${0.9}_{-0.5}^{+0.4}$	1.28 ± 0.05	1.88 ± 0.07	2.1 ± 0.6	3.1 ± 0.7 $\times {10}^{-18}$	7 ± 3

Note. Solar metallicity and a Salpeter IMF were assumed throughout.

Download table as:  ASCIITypeset image

We analyzed the stacked spectrum as in Section 2, this time limiting the age grid to the age of the universe at z_m = 2.808 (2.3 Gyr) and using the average FWHM derived for the individual targets as spectral broadening.¹⁷ The mean properties of our sample result in ${t}_{50}={0.30}_{-0.05}^{+0.20}$ Gyr and ${Av}\,={0.6}_{-0.4}^{+0.6}$ mag at fixed solar metallicity given the high 〈M_⋆〉 of our sample (Mancini et al. 2019; see Table 1). The 1σ uncertainties were estimated marginalizing over these two quantities, taking into account the simple parameterizations used for the SFHs.

We then measured the H_δA, D_n4000, and D_B spectral indices (Worthey & Ottaviani 1997; Balogh et al. 1999; Kriek et al. 2006). Since our spectral resolution is worse than in the original Lick system, the H_δA was measured on the best-fit BC03 model smoothed to the Lick resolution and adding the broadening due to a 200 km s⁻¹ stellar velocity dispersion corresponding to our 〈M_⋆〉. The D_n4000 and D_B instead are not affected much by resolution. The indices result in D_n4000 = 1.28 ± 0.05, D_B = 1.88 ± 0.07, and ${H}_{\delta A}\,={10.1}_{-0.6}^{+0.0}$ Å. The error on ${H}_{\delta A}$ was computed measuring the index maximum variation on the templates within the uncertainties on t₅₀ and A_V, accounting for the negative correlation between these two. The weakness of the D_n4000 compared the strength of the Balmer break (D_B) and of ${H}_{\delta A}$ classifies the spectrum as post-starburst dominated by A-type stars. We show the indices in Figure 2, where we also show ${H}_{\delta A}$ measured directly on the stacked spectrum to display the effect of low resolution. We checked whether our spectroscopic sample is representative of the parent photometric sample and quantified the impact of the selection bias introduced by the H-band selection at 2.5 ≤ z ≤ 3.5. A Kolmogorov–Smirnov test reveals that the probability of being drawn from the same M/L distribution is only 3%, confirming the selection bias. However, removing 30% of the highest M/L objects from the parent sample yields p = 0.5, leaving our HST sample consistent with the remaining 70%. The median M/L ratios of the two samples become identical once the highest M/L objects, 40% of the parent sample, are removed. The lowest M/L that a non-post-starburst galaxy (defined here as a stellar population characterized by D_B/D_n4000 < 1) would have relative to that of our average galaxy implies a cut (M/L_old ≥ M/L_HST + Δ M/L) that includes at most 30% of the unobserved sample. We thus conservatively conclude that at least 60%–70% of passive galaxies at z ∼ 3 are post-starburst, as the remaining "redder" galaxies could be as well post-starbursts with higher dust extinction, but also redshift interlopers, star-forming interlopers, or AGN (our single failed confirmation from spectroscopy is among the reddest and lowest M/L galaxies in the parent sample). Hence, the typical passive galaxy at z ∼ 3 is a post-starburst. This conclusion is corroborated by the global fast rise of the passive stellar mass function over 2.5 < z < 4 (Muzzin et al. 2013; Davidzon et al. 2017) consistent with a ∼1 dex increase of the number density, implying that a large proportion of the massive passive population in place at z ∼ 3 must have been very young and recently formed. In addition, from Figure 3 and from the strong Balmer lines in the published spectra, it appears that all z ≥ 3 spectroscopically confirmed passive galaxies from the literature are sharing our global spectral properties. Further efforts need to be taken to spectroscopically identify trulyold and passive galaxies (non-post-starburst) at high redshift, if they exist.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Given the tight anticorrelation between the age of a stellar population and the specific SFR (sSFR = SFR/〈M_⋆〉) of its host galaxy, expected in the case of a smooth intrinsic SFH (as opposed to a simple stellar population (SSP), or a complex SFH), we investigated whether the age decrease at fixed M_⋆ from z ∼ 1.5 to z ∼ 2.8 translates into a higher sSFR.

We used MPFIT to model the residuals of the fit using a single Gaussian centered and fixed at 3727 Å with a width matched to the data spectral resolution, measuring F_{[O ii]} = (3.1 ± 0.7) × 10⁻¹⁸ erg s⁻¹ cm⁻². The observed [O ii] flux was dereddened adopting the best-fit dust attenuation (A_V = 0.6 mag) using a value of f = 0.83 for the ratio between stellar and nebular extinction (Kashino et al. 2013) as in Gobat et al. (2017) in order to perform a consistent relative comparison. Following Kennicutt (1998) we then obtained an SFR = 7 ± 3 M_☉ yr⁻¹, which is a factor of ∼60 below the main-sequence level at 2.5 < z < 3.5 (Schreiber et al. 2015). This is in agreement with estimates for intermediate redshift quiescent galaxies from dereddened [O ii] luminosity (SFR_{[O ii]} = 4.5 M_⊙ yr⁻¹; Gobat et al. 2017), which is in turn consistent with their FIR luminosity (SFR_IR = ${4.8}_{-1.3}^{+1.8}$ M_⊙yr⁻¹; Gobat et al. 2018, hereafter G18). Therefore, the age evolution implied by the continuum appears to occur at constant sSFR. This is also suggested by the equivalent width (EW) of the line that is 2.1 ± 0.6 Å rest-frame. This value is almost half that obtained by Gobat et al. (2017) due to a higher UV continuum produced by ∼1 Gyr younger stars. Assuming that the [O ii] flux arises from a 100 Myr old stellar population, forming stars at a constant rate and with A_V = 0.6 mag, we subtracted the associated continuum component from the stacked spectrum to explore the effect on the age determination, deriving ${t}_{50}={0.30}_{-0.05}^{+0.20}$ Gyr and ${Av}={0.9}_{-0.5}^{+0.4}$ mag. Letting the age and extinction of the star-forming component vary (100–300 Myr and 0.6–1.5 mag, respectively) leaves t₅₀ unchanged while increasing A_V up to 1.0 mag. The flux at λ > 4400 Å varies around the 95% of the original flux, suggesting that the CSF component does not account for more than 1%–2% of 〈M_⋆〉 as also confirmed by the mass produced in a Δt ∼ 300 Myr at a constant SFR of 10 M_☉ yr⁻¹, namely, the upper limit inferred here. We thus conclude that the presence of the youngest stellar population does not impact our results significantly. The sSFR of our sample, (4.35 ± 2.47) × 10⁻¹¹ yr⁻¹, is consistent with the available estimates at z ∼ 1.5–1.8 (Sargent et al. 2015; Bezanson et al. 2019) and with the sSFR for the four quiescent galaxies with a spectral break in Schreiber et al. (2018) for which sSFR_{[O ii]} is available. However, the best-fit t₅₀ obtained here implies a formation redshift z_form ≳ 3.2, meaning that at z ∼ 4 our galaxies would be actively forming stars. To check that the average spectrum is not dominated by emission from unobscured AGN, we removed from the stack the two X-ray detected sources (${L}_{{X}_{2\mbox{--}10\mathrm{keV}}}={10}^{44.3}$ and 10^43.7 erg s⁻¹), getting ${F}_{{[{\rm{O}}{\rm{i}}{\rm{i}}]}_{\mathrm{noX}}}$ = (3.53 ± 0.90) × 10⁻¹⁸ erg s⁻¹cm⁻² consistent within the uncertainties. Stacking the available X-ray data in COSMOS for the undetected sources results in ${L}_{{X}_{2-10\mathrm{keV}}}\lt {10}^{43.1}$ erg s⁻¹. If part of the [O ii] luminosity were caused by low AGN activity or shocks, the intrinsic SFR would be even lower, strengthening the conclusion that the young spectral age does not directly map into a higher sSFR. Using the UV-extended Maraston et al. (2009) models the best-fitting template yields a similar t₅₀ ∼ 0.3 Gyr but produces a slightly worse χ² and a less solid [O ii] detection, with a flux that is 60% of the one measured with BC03 templates, at 2.9σ confidence.

It is interesting to place our z ∼ 3 results in the overall evolutionary context of passive galaxies. In order to be able to compare to the available literature we express, for convenience, the SFR constraints in terms of the available gas fraction through the relation f_mol = sSFR/SFE, where SFE is the star formation efficiency. We caution that this is just an alternative way to look at the SFR result, as we are using the same single constraint: we use both quantities interchangeably in the following and notably in Figure 4. We assume the same SFE derived in G18 (5 × 10⁻¹⁰ yr⁻¹), which is lower by ×2–3 than that of typical star-forming galaxies, noticing that such reduced SFE is also typical of post-starburst galaxies (Suess et al. 2017). Our SFR thus converts into M_mol = (1.5 ± 0.6) × 10¹⁰ M_⊙, hence f_mol ∼ 9±4%. We compare it to CO or dust-continuum-based gas fractions and upper limits (converted to Salpeter) for quiescent and post-starburst galaxies: Davis et al. (2014) and Saintonge et al. (2011) for local massive PEGs; Sargent et al. (2015), Bezanson et al. (2019), Spilker et al. (2018), Zavala et al. (2019), Rudnick et al. (2017), Suess et al. (2017), Spilker et al. (2018), Hayashi et al. (2018), and Gobat et al. (2018) for intermediate-z quiescent galaxies; Schreiber et al. (2018) and Valentino et al. (2020) for z ∼ 3–4 galaxies. Despite the uncertainties, our data at z ∼ 3 seem to disfavor the steep (1+z)^4–5 trend inferred from z = 0 to 1.5–2, suggesting a flattening in the M_mol/M_⋆ evolution (or, equivalently, of the sSFR).

Zoom In Zoom Out Reset image size

The published SFR_{[O ii]} for z > 3 galaxies with a clear spectral break in Schreiber et al. (2018) and Valentino et al. (2020) also seem to support this trend.

Our data thus suggest that a substantial fraction of the massive quiescent population at z ∼ 3 is approaching the quenching epoch, becoming intrinsically young and with spectra of post-starburst galaxies, but with a more or less constant f_mol (or, equivalently, sSFR) over 1.5 < z < 3–4. Various processes could be simultaneously acting: cosmological cold flows, AGN feedback, dust grain growth in the interstellar medium, satellite accretion or dust destruction by sputtering with a hot X-ray halo. It is beyond the scope of this Letter to investigate this further.

One could wonder if highly obscured star formation could be present and go unrecognized, given that their post-starburst nature might suggest that a highly obscured starburst might have been previously present. Stacking at 3 GHz (excluding two clear radio AGN detections at S_{3 GHz} 0.58 ± 0.03 mJy and 14 ± 4 μJy) results in a 3σ signal with S_{3 GHz} = 2.72 ± 0.93 μJy, which translates into a rest-frame L(1.4 GHz) ∼2 × 10²³ W Hz⁻¹, 4 times higher than in Gobat et al. (2018). This is probably just suggesting a higher radio AGN activity at fixed stellar mass closer to the quenching epoch, continuing the rapid evolution seen from z = 0 to 1.5–2. On the other hand, such a detection would also be formally consistent with 40–50 M_⊙ yr⁻¹ of obscured star formation (assuming the IR–radio correlation at z = 3; Delhaize et al. 2017). This would still place the typical object 8 times below the MS. While we tend to interpret strong Balmer absorption lines plus weak [O ii] emission as a sign of post-starburst galaxies with residual unobscured star formation, we cannot fully rule out stronger obscured star formation, which could be analogous to what seen locally for e(a) galaxies (Poggianti & Wu 2000). Only future Atacama Large Millimeter/submillimeter Array observations of several of these targets could conclusively solve this issue.

We thank the anonymous referee for their constructive report and Gabriel Brammer for his help with grizli. A.C. acknowledges grants PRIN MIUR 2015, PRIN 2017 (20173ML3WW_001), and ASI n.I/023/12/0. S.J. acknowledges MICIU grant AYA2017-84061-P, co-financed by FEDER. I.D. is supported by Marie Skłodowska-Curie grant agreement No. 788679.