An Exquisitely Deep View of Quenching Galaxies through the Gravitational Lens: Stellar Population, Morphology, and Ionized Gas

2.1.1. Source e1341

2.1.2. Source m2129

2.1.3. Source m0451

2.1.4. Source m1423

3.1.1. HST Photometry

3.1.2. IRAC and WISE Photometry

We use detailed mass models of each cluster to derive the intrinsic properties (magnification and source morphology) of each target. These parametric models are constrained by the location of strongly lensed background systems using the software Lenstool (Jullo et al. 2007). For each cluster, the optimization of the mass model is performed through a Markov chain Monte Carlo (MCMC) minimizing the root mean square (rms) between the observed and the predicted locations of multiple systems. The magnification and error bars for the targets provided in Table 1 are derived using both the best model (providing the best reproduction of the multiple images) as well as the full set of MCMC realizations sampling the posterior distribution of model parameters. Details on these mass models are presented in the reference papers. Specifically, eMACSJ1341.9-2442 has one multiple image system with five images (Ebeling et al. 2018). The lensed target e1341 is multiply imaged. MACSJ2129.4-0741 has two systems with nine images (Toft et al. 2017). The lensed target m2129 is just outside the region of constraints, located 11'' away from the multiply imaged region. MS0451.6-0305 has 17 systems with 49 images (Jauzac et al. 2021). The lensed target m0451 is multiply imaged. MACSJ1423.8-2404 has 3 systems with 12 images (Limousin et al. 2010). The lensed target m1423 is 20'' away from the multiply imaged region. Table 1 lists the rms of the observed locations of multiple systems.

The same mass models are used to reconstruct the intrinsic source morphology of each target. We directly ray-trace the HST/WFC3 F160W image of each image onto a regular grid on the source plane, and perform the same reconstruction for an HST/WFC3 PSF model to derive our spatial PSF on the source plane. These reconstructions were obtained for the full set of MCMC realizations, to estimate the statistical errors in the source reconstructions.

The exception is for target e1341: we used the closest HST/WFC3 band in wavelength (F140W), as F160W imaging data are not available. The lensing arc, im1, of e1341 is additionally galaxy–galaxy lensed so that the reconstructed image does not cover the entire source once reconstructed. Therefore, only im2 and im3 are reconstructed for e1341.

We used the panchromatic VLT/X-SHOOTER Echelle spectrograph (Vernet et al. 2011) to obtain stellar continuum spectra for our targets. All spectra were taken with slit widths of 10 in the UVB arm, and 09 for the VIS and NIR arms. This corresponds to nominal spectral resolutions ¹¹ of R ≈ 5400, 8900, and 5600 for the UVB, VIS, and NIR arms, respectively. The observation log is provided in Table 2.

Table 2. VLT/X-SHOOTER Observation Log

	e1341 (im1)	m2129	m0451 (im1)	m1423
R.A. (J2000)	13h42m02374	21h29m2234	04h54m1265	14h23m4697
Decl. (J2000)	− 24d41m5211	− 07d41m3106	− 03d01m1650	+ 24d05m2852
Integration time (h) for UVB/VIS/NIR	4.14/4.43/4.80	2.29/1.74/2.67	4.12/4.36/7.87	11.64/12.04/13.33
S/N (Å−1)	11.6	6.3	15.5	9.4
Slit angle (north to east)	10°	−10°	35°	−80°, −10°

Notes. Positions are based on F160W, or F140W if the former is unavailable. The listed integration times only include useful data and not discarded data. The S/N per Å (rest-frame) is estimated from the rest-frame wavelength region between 4110 and 4210 Å just redward of the Hδ feature, with clipping against strong skylines and assuming the same correction for correlated noise following the procedures as done in Toft et al. (2017).

Download table as:  ASCIITypeset image

The reduction procedures for m2129 are described in Toft et al. (2017). The other three targets were reduced using the X-SHOOTER pipeline (version 2.6.8; Modigliani et al. 2010), supplemented by our customized scripts. The reduction steps were as follows. We used standard pipeline recipes up to the creation of the rectified Echelle orders (ORDER2D), with a few exceptions: First, we removed in the UVB frames fixed pattern noise with a Fourier analysis, similar to the method described by Schönebeck et al. (2014). Second, we removed a temporarily and spatially changing readout pattern from the NIR frames, which manifests itself as a large-scale parabola-like pattern in the rectified frames. This pattern can be efficiently removed by subtracting the median of unexposed pixels in each detector row. ¹² Third, the pipeline version being used creates small-scale ring-like features around bad pixels during rectification. To avoid problems with these rings, we increased the masks around detected cosmic ray hits, a modification that we directly implemented in the pipeline (Zabl et al. 2015).

Instead of processing all exposures simultaneously, we created ORDER2D files from only a single pair of nodding positions at a time. We stacked all nodding pairs for each order with our scripts, and subsequently merged the order stacks to a single 2D spectrum. In both steps, we weighted during the combining at each wavelength with the spatial median of the inverse pipeline variance. This weighting is different from the pipeline in two ways: First, the pipeline does not apply weights in the combining of the exposures. Second, during the order merging, the pipeline uses a pixel-based inverse variance weight, which we found to be somewhat problematic, as the pipeline variance estimate is a biased estimator.

Before stacking the nodding pairs, we corrected each of the 2D frames to a heliocentric standard and a vacuum wavelength scale. Flux calibration was done based on response functions determined from standard stars taken during the same nights as the science observations, or if no standard was taken during the same night, from the nearest available night. Telluric corrections were derived for each science exposure using observations of telluric standard stars observed close in time and airmass. In practice, we first determined an atmospheric model using molecfit (version 1.5.1; Kausch et al. 2015; Smette et al. 2015) from the telluric standard stars, and we subsequently evaluated this model at the exact airmass of the science observations. Finally, a substantial part of our data were observed without a working atmospheric dispersion corrector. Therefore, we implemented a software-based correction, which applies for the relevant exposures the appropriate wavelength-dependent spatial offset to the 2D frames and corrects approximately for the increased slit loss. The rectified 2D spectra are shown in Figure 2.

The 1D spectra, as shown in Figure 3, were extracted directly from the stacked ORDER2D files and were subsequently merged with inverse variance weighting. We did the extraction with a simple aperture extending from −05 to +05 around the pointing center.

We model the X-SHOOTER spectra using synthetic models from the Bagpipes population synthesis fitting code (Carnall et al. 2018), which uses the updated versions of the Bruzual & Charlot (2003) population synthesis models. Prior to the fitting, we correct the spectra and photometry for galactic reddening following the approach of Schlafly & Finkbeiner (2011), using the Fitzpatrick (1999) reddening law and R_V = 3.1 appropriate for high Galactic latitudes.

We adopt the delayed exponentially declining parameterization of the star formation history (SFH),

$$\begin{eqnarray}&&\mathrm{SFR}(t)\propto (t-{t}_{0}){e}^{-(t-{t}_{0})/\tau }.\end{eqnarray} \tag{ 1 }$$

To facilitate comparisons with literature measurements, we also provide the median stellar age and star formation duration calculated from the SFH following Pacifici et al. (2016), which are shown to be more robust than the age since the onset of star formation in the delayed-τ model (Belli et al. 2019). The median stellar age, t₅₀, is defined as the lookback time from the spectroscopic redshift z_spec at which 50% of the M_⋆ has been assembled. Here, t₁₀ − t₉₀ is defined as the duration between 10% and 90% of the M_⋆ being assembled. The provided uncertainties of t₅₀ and t₁₀ − t₉₀ are derived from posterior distributions of the SFH parameters.

In addition to the SFH parameters, we fit for the systemic redshift, stellar velocity dispersion (σ_⋆), stellar mass normalization (Kroupa 2001 initial mass function), stellar metallicity (Z_⋆/Z_⊙) ¹³ , and dust attenuation (A_V ) of a dust screen following the Calzetti et al. (2000) reddening law. The adopted priors on these fitting parameters are indicated in Table 3. Bagpipes includes nebular emission lines calculated with the Cloudy photoionization code (Ferland et al. 2017), where we assume the line-emitting gas has the same metallicity and dust attenuation as the stars.

Table 3. Stellar Population Properties

	zspec	log(μM⋆/M⊙)	log(sSFR/yr−1)	AV	log(Z⋆/Z⊙)	σ⋆	t0	τ	t50	t10 − t90
						(km s−1)	(Gyr)	(Gyr)	(Gyr)	(Gyr)
e1341	1.5954 ± 0.0001	11.62 ± 0.07	−11.20 ± 0.46	0.4 ± 0.1	−0.12 ± 0.07	160 ± 9	${1.43}_{-0.08}^{+0.06}$	${0.14}_{-0.03}^{+0.01}$	1.19 ± 0.05	0.47 ± 0.07
m2129	2.1487 ± 0.0002	11.86 ± 0.05	−12.77 ± 1.83	1.1 ± 0.1	−0.10 ± 0.11	318 ± 28	0.68 ± 0.06	${0.04}_{-0.01}^{+0.02}$	0.61 ± 0.05	0.14 ± 0.05
m0451	2.9223 ± 0.0001	11.65 ± 0.06	−9.77 ± 0.08	0.9 ± 0.1	−0.30 ± 0.07	170 ± 12	0.42 ± 0.03	0.051 ± 0.005	0.33 ± 0.02	0.17 ± 0.01
m1423	3.2092 ± 0.0002	11.01 ± 0.07	−8.43 ± 0.10	0.8 ± 0.1	−0.01 ± 0.10	277 ± 19	0.15 ± 0.01	${0.021}_{-0.001}^{+0.002}$	0.12 ± 0.01	0.07 ± 0.01
Prior	${ \mathcal N }({z}_{0},0.001(1+z))$	${ \mathcal N }(11.3,0.5)$	⋯	Uniform(0,2)	$Z/{Z}_{\odot }:{ \mathcal N }(1,0.5)$	${ \mathcal N }({\sigma }_{0},50)$	${ \mathcal N }(0.8,0.3)$	${ \mathcal N }(0.1,0.2)$	⋯	⋯
Limits	⋯	(10, 12.8)	⋯	⋯	Z/Z⊙: (0.2, 2)	(100, 450)	(0.05, tuniv)	(0.02, 0.3)	⋯	⋯

Notes. We report the median values and quote the average difference of 16th and 84th percentiles of the posterior distributions from the median values as uncertainties. See Section 3.4 for details on the derivation of these parameters. Note that the uncertainties are propagated from photometric and spectral flux errors and do not include systematic uncertainties inherent to stellar population synthesis. The solar metallicity, Z_⊙, is defined as 0.02 following the definition of Bruzual & Charlot (2003). A comparison of the derived parameters of m2129 with values published in Toft et al. (2017) and Newman et al. (2018a) is presented in Table 10. The priors used in the population synthesis fits are indicated in the bottom two rows, where ${ \mathcal N }(\mu ,\sigma )$ indicates a Gaussian distribution with mean μ and standard deviation σ. Numbers z₀ and σ₀ are initial guesses based on basic template fits to the individual X-SHOOTER spectra, and the broad prior width is adopted there to allow substantial flexibility around those values. There are implicit priors on sSFR, t₅₀, and t₁₀ − t₉₀ arising from the priors on SFH parameters t₀ and τ as specified.

Download table as:  ASCIITypeset image

The X-SHOOTER spectra of the four lensed targets at different redshifts are masked such that all targets have similar rest-frame spectral coverage from ≈2500 Å to just beyond [O iii] λ5007, though the targets have different gaps within this range due to the redshifted NIR atmospheric windows. To fit the combined observed data of the multiple X-SHOOTER channels and broadband photometry, we allow for an additional low-order "calibration" function multiplied to the spectra to reduce the effect of normalization differences between the spectra and photometry (Carnall et al. 2019b). Therefore, the photometry constrains the overall shape of the spectral energy distribution (SED) and the spectra constrain features on wavelength scales smaller than the calibration function (i.e., stellar absorption lines). Finally, we allow a final parameter that scales the X-SHOOTER spectrum uncertainties and applies the appropriate normalization penalty to the likelihood (Carnall et al. 2019b).

The posterior probability distribution function of the fit parameters is sampled within Bagpipes using the MultiNest algorithm (Feroz et al. 2019). The maximum a posteriori parameters of the chain are presented in Table 3, and the parameter covariances are shown in Figure 15. The best-fit models are shown in Figures 16–19. The delensed stellar masses and SFR are presented in Table 4.

Table 4. Delensed Stellar Population Properties

	log(M⋆)	log(SFR100)	M⋆/${{ \mathcal M }}^{* }$ a	M⋆/${{ \mathcal M }}^{* }$ q
e1341	10.15 ± 0.14	$-{1.0}_{-0.7}^{+0.3}$	0.2 ± 0.1	0.3 ± 0.1
m2129	11.20 ± 0.05	$-{1.6}_{-2.5}^{+1.4}$	2.7 ± 0.5	${2.3}_{-0.7}^{+1.0}$
m0451	10.61 ± 0.10	0.8 ± 0.1	0.4 ± 0.1	${0.7}_{-0.5}^{+0.4}$
m1423	10.58 ± 0.08	2.2 ± 0.1	0.1 ± 0.1	${0.6}_{-1.0}^{+0.6}$

Notes. Stellar population parameters for our sample after correcting for lensing magnification. The errors are propagated from the spectral fitting and the lens model. Masses and star formation rates are quoted in logarithmic units of M_⊙ and M_⊙ yr⁻¹, respectively. We use the stellar mass function presented in Muzzin et al. (2013) to estimate the characteristic masses of all galaxies (${{ \mathcal M }}^{* }$ _a) and those of quiescent galaxies only (${{ \mathcal M }}^{* }$ _q) at their respective redshifts.

Download table as:  ASCIITypeset image

Our sample probes galaxies near the characteristic mass of the stellar mass function of quenched galaxies. Table 4 shows the intrinsic stellar masses of the sample, and how they compare with the characteristic masses of galaxies using the stellar mass function presented in Muzzin et al. (2013). Intrinsic stellar masses are derived by dividing the lensed stellar mass by their magnification factors, propagating the respective uncertainties. Three out of four galaxies in the sample are less massive than the characteristic masses of galaxies at their respective redshifts. Thanks to lensing, we reach for the first time less extreme quenched galaxies at z ≳ 1.5 that are more representative of the overall galaxy population, as illustrated in Figure 4.

Zoom In Zoom Out Reset image size

The spectral fitting results presented in Section 3.4 indicate that our sample of quiescent galaxies have undergone rapid star formation histories. The reported uncertainties are propagated from photometric and spectral flux errors, and they do not include systematic uncertainties such as variation in IMF, SFH, etc. Their star formation timescales, t₁₀ − t₉₀, are at most a few hundred Myr and are short compared to their median stellar population ages and the ages of the universe at the observed redshifts. As shown in Figure 5, their specific star formation rates (sSFR) are lower than typical star-forming galaxies by at least an order of magnitude, except for m1423. Its SFR from spectral fitting is higher than that derived from [O ii], suggesting that the discrepancy is due to recent rapid quenching. We argue in Section 4.4.1 that m1423 is likely caught in transition from being star-forming to quiescent.

Zoom In Zoom Out Reset image size

Our spectral fitting assumes, like many other works, a parameterized SFH. Specifically, our work assumes the commonly used delayed-τ parameterization. Although these assumptions are quite standard, and are likely reasonable assumptions for high-redshift galaxies, it may still be insufficient to capture the individual SFH of some galaxy types, as argued by some works (e.g., Maraston et al. 2010; Behroozi et al. 2013; Gladders et al. 2013; Abramson et al. 2016; Pacifici et al. 2016; Ciesla et al. 2017). We note, however, that most of the concerns raised about the use of (delayed)-τ models do not apply to high-redshift quiescent galaxies like our sample. The short time between their formation and observation implies that there is little effect for model degeneracies, at least compared to low-redshift galaxies that are much more evolved. Possible alternatives are log-normal SFHs (Gladders et al. 2013) or nonparametric models (Iyer et al. 2019; Leja et al. 2019; Akhshik et al. 2020). The diversity of SFH models used in different studies makes a direct comparison of resulting parameters difficult, and certainly has a non-negligible effect on scaling relations derived (Carnall et al. 2019a). With Bagpipes, we find similar results for the stellar masses, metallicities, and timescales when adopting "double power-law" and "log-normal" SFHs as defined by Carnall et al. (2018). The only parameter with a significant dependence on the SFH model is the recent star formation rate, which is significantly lower for these other parameterizations that fall off more sharply than the delayed-τ model at late times for a given t₅₀ and t₁₀ − t₉₀. Despite this systematic uncertainty on the absolute value of the recent star formation rate, all of the parameterized fits agree that the specific star formation rate for this sample is well below that of typical star-forming galaxies at similar redshifts. Investigating the impact of assuming a nonparametric SFH on derived properties will be explored in a future work.

The metal content of galaxies is a powerful probe of their enrichment histories. Stellar metallicity correlates with stellar mass in quiescent, early-type galaxies at z = 0.1–0.7 (Gallazzi et al. 2005, 2014). As it requires deep spectra to accurately measure stellar metallicity, only a handful of quiescent galaxies at z ≳ 1 have robust stellar metallicity measurements thus far. Lensing magnification offers a unique opportunity to obtain such a measurement in galaxies with stellar mass below the characteristic mass out to z > 3.

In Figure 6, we show the stellar mass–metallicity relation (MZR) and the stellar mass–age relation. Only literature measurements with errors smaller than Δlog(Z_⋆/Z_⊙) < 0.3 are shown in Figure 6. The only works presenting stellar metallicity measurements of quiescent galaxies at z > 1.9 are Kriek et al. (2016) ¹⁴ , Newman et al. (2018a), and Jafariyazani et al. (2020). Only one galaxy, MRG-M0150, in the sample of Newman et al. (2018a) is plotted, because we exclude those without lensing models, which are needed to constrain the intrinsic M_⋆. As for MRG-M0138, we obtain a metallicity ¹⁵ based on the abundance measurement presented in Jafariyazani et al. (2020), which provided a more detailed analysis than that in Newman et al. (2018a). As m2129 is part of the sample in this work as well as in Newman et al. (2018a), to avoid duplication we show only our measurements here. ¹⁶ We also include measurements of individual quiescent or early-type galaxies at z = 1–2 from Lonoce et al. (2015) and Kriek et al. (2019), as well as composite spectra (Onodera et al. 2015; Estrada-Carpenter et al. 2019; Saracco et al. 2019). ¹⁷ We only consider literature measurements of Z_⋆ when they are based on detection of metal absorption features and determined to within 0.3 dex accuracy. We take the MZR for quiescent galaxies from Gallazzi et al. (2014) for z ≈ 0.1 (SDSS) and z ≈ 0.7, as well as the best-fit relation presented in Ferreras et al. (2019b) for 19 quiescent galaxies at z = 0.6–1.3 (median z = 0.85).

Zoom In Zoom Out Reset image size

Our sample provides new insights regarding the stellar MZR relation, because they probe the most distant and least massive, quenching galaxies among all z > 1 literature studies. Sources e1341 and m1423 have stellar metallicities on par with literature measurements extrapolated to higher redshifts and lower stellar masses, while m0451 and m2129 lie lower. The stellar metallicities of our sample span from ≈0.5 to 1 solar value, somewhat lower than those of more massive, quiescent galaxies at comparable redshifts (Lonoce et al. 2015; Kriek et al. 2016; Newman et al. 2018a; Jafariyazani et al. 2020), but are still higher those of star-forming galaxies at z = 2.5–5 (Halliday et al. 2008; Cullen et al. 2019). While these differences may be due to diverse chemical enrichment histories (and thus dependence on redshift and stellar masses), variations in metallicity calibrations and stellar populations probed are also plausible explanations. Indeed, existing stellar metallicity measurements of m2129 range from subsolar (log(Z_⋆/Z_⊙) = −0.6 ± 0.5; Toft et al. 2017) to mildly supersolar (log(Z_⋆/Z_⊙) = 0.16 ± 0.13; Newman et al. 2018a), as detailed in Appendix C. This work provides an intermediate measurement of log(Z_⋆/Z_⊙) = −0.10 ± 0.11. Moreover, while the works obtained using Bruzual & Charlot (2003) stellar population models assume a solar metallicity value of Z_⊙, more recent works using the alf code (Conroy & van Dokkum 2012; Conroy et al. 2018) assume more recent measurements of Z_⊙ = 0.0142 (Asplund et al. 2009). Systematic effects in deriving stellar metallicities are thus a significant source of uncertainty that needs to be better quantified in future works. We will discuss the implications of these results in Section 5.

Singly ionized magnesium Mg ii λλ 2796,2803 (hereafter Mg ii) absorption in the rest-frame UV spectrum can arise from warm atomic gas (T ∼ 10⁴ K) or stellar photospheres (late B-, A-, F-, and G-type stars; Fanelli et al. 1992; Snow et al. 1994). The adjacent Mg i λ2853 traces neutral magnesium and has lower equivalent widths than Mg ii in both the interstellar medium and stars (Maraston et al. 2009; Coil et al. 2011; Zhu et al. 2015). When present in stars, Mg i absorption is stronger in later-type stars (F- and G-types) and absent in O- and B-type stars (Fanelli et al. 1992; Rodríguez-Merino et al. 2005).

To facilitate spectral line measurements, we model the local stellar continua as a straight line using the model spectra within the line-free bandpasses specified in Table 1 of Maraston et al. (2009). Equivalent widths (EW) and line fluxes of Mg ii and Mg i are then measured by summation of the X-SHOOTER spectra normalized by the local continua over the central bandpasses as shown in Figure 7. We caution that emission line infilling might lead to underestimated absorption line EW and fluxes. Errors are propagated from the noise on the X-SHOOTER spectra. The spectral line measurements are provided in Table 5.

Zoom In Zoom Out Reset image size

Mg ii absorption is ubiquitous in our sample with different profile shapes, sometimes accompanied with redshifted emission, as shown in Figure 7. Mg i, being a weaker feature than Mg ii, is detected in some of the galaxies as well. Source e1341 has deep Mg ii and Mg i absorption at the systemic velocity, as well as a broad blueshifted Mg ii absorption wing that extends to ≈−2000 km s⁻¹ that is not apparent in Mg i. Although m2129 has the lowest S/N spectrum within the sample, a comparison of the Mg ii and Mg i profiles shows tentative absorption at systemic velocity accompanied by redshifted emission. Source m0451 shows a pair of slightly blueshifted Mg ii doublet absorption near systemic velocity accompanied by redshifted emission. Its Mg ii λ 2796 and Mg ii λ 2803 absorptions are shifted by ≈−200 km s⁻¹ from the systemic redshift. In addition to gas associated with m0451, there might also be contribution from its close companion, tidal interaction, and/or intragroup environment: the quiescent galaxy has a star-forming companion galaxy located ≈23 kpc away on the source plane (gal 2; MacKenzie et al. 2014) with the same redshift as measured from its CO J = 3 → 2 emission (z_{CO,gal 2} = 2.921 ± 0.001; Shen et al. 2021). This pair of galaxies is embedded in a z = 2.9 group as discussed in Section 2.1. Source m1423 has a highly blueshifted Mg ii absorption by up to ≈−1500 km s⁻¹ from the systemic velocity. In all cases, the Mg ii absorption at its trough is consistent with having zero flux, i.e., black. This suggests that the absorbing gas has high column density and covering fraction.

How do we interpret the origins of the Mg features in our sample? One way to discern between these origins is to examine the line profile. Photospheric absorption takes place at the systemic velocity. In fact, deep UV spectroscopy provided the first confirmations of evolved stellar populations dominated by stars older than ≳0.5 Gyr at z ≈ 1.5–2 with optical spectrographs (Dunlop et al. 1996; Spinrad et al. 1997; Cimatti et al. 2004, 2008; McCarthy et al. 2004; Daddi et al. 2005). On the other hand, interstellar or circumgalactic medium absorption is not restricted to the systemic velocity, and can be blueshifted (outflow) or redshifted (infall). While part of the absorption can be attributed to stellar photospheres, given the stellar metallicities of our sample (Section 4.2), gas absorption is required to account for the absorption depth as well as the blueshifted absorption and redshifted emission components. The Mg ii profiles are unlikely to be dominated by circumstellar winds, as Mg ii in stellar winds usually have modest velocities, of order ≈100 km s⁻¹ (Praderie et al. 1980; Snow et al. 1994) except for supergiants and giants. Furthermore, the Mg ii absorption in stars is not accompanied by redshifted Mg ii emission (Snow et al. 1994).

Galactic outflow is a likely explanation for the Mg ii profiles in m0451 and m1423. Gas outflowing into the line of sight obscures light from the stellar continuum creating the blueshifted absorption. Gas outflowing away from us scatters photons back into our line of sight, to create the redshifted emission (Weiner et al. 2009; Rubin et al. 2011). Galactic outflows are commonly seen in massive post-starburst galaxies at z < 1.5 (Tremonti et al. 2007; Coil et al. 2011; Sell et al. 2014; Maltby et al. 2019). In a galactic outflow, Mg i shows an absorption profile that is weaker than but similar to that of Mg ii, and it has an EW of a quarter or less than that of Mg ii (Coil et al. 2011). The similarities of the Mg i and Mg ii profiles in m1423, m0451, and perhaps m2129 suggest that outflows are present, as Mg i can only have photospheric absorption but not emission in stars (Fanelli et al. 1992). It is unlikely that the gas outflows reported in the youngest targets here are due to circumnuclear outflow. The rest-frame UV continua are well-modeled by the stellar populations without the need to invoke an AGN continuum (Figures 16–19). The Mg ii profiles of our targets are distinct from those of quasar winds in broad absorption line quasars (e.g., Trump et al. 2006). In Section 4.4, we rule out the presence of type 1 AGN in our sample. Thus, the Mg features in our sample provide unambiguous evidence for galactic-scale gas outflow in addition to evolved stellar populations.

The Mg ii profile encapsulates information about the structure of a galactic outflow, i.e., covering fraction, opacity, orientation, etc. Detailed modeling of the Mg ii profile is beyond the scope of this paper and will be deferred to a future work. Observations of rest-frame optical emission lines could further constrain the warm gas in quiescent galaxies, as we shall explore in Section 4.4. The implications of galactic-scale outflows are discussed in a broader context in Section 5.4.

Warm (T ∼ 10⁴ K) gas in galaxies can be photoionized by AGN or hot stars (both young and old), or excited by shocks. The strengths and ratios of emission lines can inform us of their production mechanism (Kewley et al. 2019, and references therein). Locally, massive quiescent galaxies sometimes have low-ionization emission regions (LIER). ¹⁸ Characterizing the gas conditions of quenching galaxies will provide clues to understand how they become quiescent.

The brightest emission lines within our spectral coverage are [O ii] λλ3726,3729 (hereafter [O ii]), [O iii] λλ4959,5007, and Hα. The oxygen emission lines are shown in velocity space along with the normalized Mg ii spectra in Figure 8, and are highlighted in orange in Figure 3. Emission line EWs and line fluxes are provided in Table 5. The measurements are computed as the excess emission to the best-fitting stellar continua (Section 4.1) using the observed spectra without binning or smoothing, measured within the passband range specified in Table 1 of Westfall et al. (2019).

Zoom In Zoom Out Reset image size

4.4.1. SFR Inferred from [O ii]

4.4.2. [O iii]-to-[O ii] Flux Ratio

4.4.3. [O iii] as a Tracer for Nuclear Activity

To constrain the intrinsic morphology of our galaxies on their source planes, we fit Sérsic profiles (Sérsic 1963; Sersic 1968) to the lensing-reconstructed images. Our procedures follow the approach described in Toft et al. (2017). In summary, to propagate the uncertainty introduced by the lensing model, we generate 100 realizations of each reconstructed image, and obtain the best-fitting morphological parameters to each realization using the GALFIT code ¹⁹ (Peng et al. 2002, 2010, 2011). The best-fitting parameters are reported in Table 7. The quoted uncertainties are computed from the standard deviation of 100 realizations (i.e., error due to the lensing model) and GALFIT errors added in quadrature. The reconstructed images, best-fitting Sérsic models, and the residual images are shown in Figure 9. Multiple-image systems (e1341 and m0451) provide additional constraints on the systematic uncertainties associated with the image reconstruction. Reassuringly, the Sérsic index n and effective radii r_eff are in good agreement across the multiple images. The axis ratio q, on the other hand, appears less constrained. While the axis ratio measurements of e1341 are consistent across the two multiple images, those of m0451 are more discrepant. The axis ratio is q ≈ 0.4 for the most magnified image (im1), compared to q ≈ 0.6–0.7 for the other two less magnified images. Although m0451 has a robust lensing model based on 17 multiply imaged systems (Section 3.2), there is possibly a systematic uncertainty in the reconstruction of the most magnified source, as the position angles vary across the three images (Table 7). The morphological parameters of the least magnified image (im3) should be least affected by lensing systematic uncertainties as long as it is resolved. Overall, multiple measurements enable us to assess the lensing systematic uncertainties in the morphological parameters. The r_eff and n are more robust to lens modeling uncertainties than q.

Zoom In Zoom Out Reset image size

We examine how the stellar morphologies of lensed quenched galaxies compare to those of unlensed ones. Figure 10 shows a comparison of the lensed quenched galaxies in this work as well as those presented in Newman et al. (2018a), compared to other unlensed spectroscopic samples of quiescent galaxies (Belli et al. 2017; Bezanson et al. 2018; Stockmann et al. 2020). The Sérsic indices of our sample are low with n < 2. The apparent axis ratios have intermediate values of q ≈ 0.4–0.7, below the median redshift relation of the apparent axis ratio presented in Figure 4 of Hill et al. (2019). These findings suggest that the lensed quenched galaxies of this work clearly have lower Sérsic indices and axis ratios compared to unlensed ones. All but one target (e1341) lies within the 1σ scatter of the mass–size relation of early-type galaxies at their respective redshifts (van der Wel et al. 2014), as shown in Figure 11. Source e1341 has r_eff ≈ 2.7 kpc (average of the estimate from its two reconstructed images), roughly four times larger than the mass–size relation of early-type galaxies at its redshift and M_⋆, corresponding to ≈4.3× the scatter of the relation. This places e1341 above the late-type galaxy mass–size relation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is unclear why lensed quenched galaxies of our sample are more disk-like than unlensed ones. The shapes of galaxies evolve with stellar mass and redshift. In the local universe, massive galaxies tend to have rounder stellar light profiles (higher n and q) than lower-mass galaxies (Krajnović et al. 2013). The trend is less pronounced at higher redshift, as stellar disks are more prevalent (van der Wel et al. 2011; Chang et al. 2013a, 2013b; Hill et al. 2019). In Figure 20, we show the stellar mass dependence of n and q. The targets of our lensed sample have lower n than unlensed quiescent galaxies at z = 0.6–1.0 (Bezanson et al. 2018) and z = 1.5–2.5 (Belli et al. 2017). The lack of an unlensed comparison sample at matching z and M_⋆ precludes us from concluding whether the disky profiles of our sample is due to their high redshifts and/or low stellar masses. Furthermore, lensed arcs that have flatter light profiles are preferentially selected as spectroscopic targets in our survey for spatially resolved studies. Another possible reason for the discrepancy is lens model uncertainty. If the magnifications were underestimated, for example due to additional substructure missed in the model, this could lead to artificially flatter light profiles. Finally, HST cannot adequately resolve the light profiles of distant compact quiescent galaxies, such that unlensed quiescent galaxies might appear rounder than they actually are. Higher spatial resolution imaging for a large sample of distant quiescent galaxies with JWST will help address the cause of this discrepancy.

We discuss the implications of these findings in Section 5. A caveat is that different rest-frame wavelengths are traced for the four galaxies in our sample (Table 7), due to their different redshifts and filters used. van der Wel et al. (2014, Equation (2)) presented a scaling relation to correct for the wavelength dependence of the effective radii. However, this involves the use of the average size gradient, ${\rm{\Delta }}\mathrm{log}{R}_{\mathrm{eff}}/{\rm{\Delta }}\mathrm{log}\lambda$, which is not well-characterized for early-type galaxies at our redshift and mass range. By using their equation, we derive a correction factor of ≲5% in effective radii. This caveat is thus expected to be insignificant and we thus did not correct for the wavelength dependence of the sizes.

Source m1423 is at the border between quiescent and star-forming (or early- and late-type) galaxies according to its position on the UVJ diagram (Figure 14). We note that its small effective radius suggests that it is nearly as compact as early-type galaxies, if we extrapolate the 3DHST+CANDELS structural relation at z = 2.5–3.0 (van der Wel et al. 2014) up to its redshift of 3.21.

To figure out what our sample of lensed quenched galaxies would evolve into by z ≈ 0, we estimate their halo mass (M_halo) evolution in the following way. Taking their delensed stellar masses, we compute their number densities using the galaxy stellar mass function of the COSMOS field (Muzzin et al. 2013). We use the Number Density Evolution Redshift Code ²⁰ (nd-redshift; Behroozi et al. 2013) to calculate the halo mass evolution for a given number density of a galaxy population. The code tracks the number density evolution due to mass accretion and mergers. In this procedure, we use a Monte Carlo simulation to propagate errors on stellar masses, the stellar mass function (due to Poisson uncertainties, photometric redshift errors, and cosmic variance), and the number density evolution. Figure 12 illustrates the expected halo mass evolution. To infer possible z ≈ 0 descendants, we overplot the halo mass distribution of the MASSIVE sample, as well as several well-known galaxy clusters including Coma, Perseus, Fornax I, and the Local Group (Crook et al. 2007; Li & White 2008; Veale et al. 2018, and references therein).

Zoom In Zoom Out Reset image size

Our calculation suggests that the embedding halos of our sample of lensed quenched galaxies would evolve into intermediate-sized galaxy groups/clusters with log(M_halo) ≈ 12–14, more like the Local Group than the Coma cluster. Their halo masses at z = 0 are below the median of those of the MASSIVE sample, and more than an order of magnitude less than those of the most massive galaxy clusters like Coma and Perseus. The halo mass evolution of our sample lies below those of various cluster surveys like CLASH (Postman et al. 2012) or GOGREEN (Balogh et al. 2017). The findings are in line with our expectations, given that our sample is more representative and numerous than the most massive galaxies at their epoch. An implicit assumption of this calculation is that the lensed quenched galaxies reside in field environments, such that the COSMOS stellar mass function provides a reliable estimate of their number densities. It is known that m0451 resides in an overdense environment that resembles a compact group (Section 2.1; MacKenzie et al. 2014; Shen et al. 2021). As for the other targets, there are no indications thus far that they reside in overdensities—although it cannot be ruled out, because of the inherent challenge in quantifying the environment of a lensed volume. Thus, for m0451, the halo mass by z = 0 could be higher than estimated. Another limitation is that the two highest-redshift targets, m0451 and m1423, have stellar masses below the 95% completeness limit of the COSMOS survey, so their number densities are less certain. The various conversions used to estimate halo masses could add further uncertainty to this comparison. At any rate, this exercise serves to provide an order-of-magnitude estimate of the z = 0 descendant halo mass. It is justified to conclude that halos containing the lensed galaxies studied in this work would not evolve into Coma-like clusters by z = 0, unless they are highly clustered.

Our sample of galaxies had already assembled a significant mass of stars when the universe was young. Our analysis in Section 4.1 indicates that their star formation histories are rapid, having formed 80% of stellar mass within 70–470 Myr. These timescales are comparable to or shorter than the median ages of the stellar populations, and much shorter than the age of the universe at their respective redshifts (≈2–4 Gyr).

How do our results compare with studies of other distant quenched galaxies? A meaningful comparison can only be made if the star formation timescale is measured with the same methods (full spectral fitting versus line indices, parameterization of star formation history) and defined in the same way. Given the variety of methods adopted in the literature, here we only attempt to conduct an order-of-magnitude comparison in order to get an impression of how the derived values compare. As near-infrared absorption line spectroscopy is time-consuming, constraints on the star formation history of quenched galaxies are limited to the most luminous ones that are more massive than our sample. Bearing these differences in mind, the rapid star formation duration of our sample is in qualitative agreement with those of the massive, quiescent galaxies presented in van de Sande et al. (2013; τ = 10–80 Myr, z = 1.5–2.1), Newman et al. (2018a; τ < 180 Myr, z = 1.9–2.6), and Valentino et al. (2020; τ = 10–16 Myr, z = 3.8–4.0), and are similar to or shorter than the quiescent galaxies at z = 1.4–2.6 (τ = 0.2–3.2 Gyr; Zick et al. 2018; see also Kriek et al. 2019).

Altogether, these findings provide evidence of a rapid buildup of stars through accelerated growth. The mere fact that massive, quiescent galaxies exist in a young universe requires that they formed stars at a higher rate before (see also Pacifici et al. 2016). As inferred from the median SFH, the peak SFR of our sample ranges from 1800 to 6000 M_⊙ yr⁻¹. The most massive (log(M_⋆/M_⊙) > 11), quiescent galaxies are shown to be consistent with having z = 3–6 submillimeter galaxies (SMG; with mean duty cycle of ≈40 Myr) as their progenitors, by means of number density and size comparison (Toft et al. 2014). Our sample probing a lower stellar mass regime appears consistent with having experienced an SMG phase, although at lower stellar masses than those presented in Toft et al. (2017). It is worth noting that m0451, together with all its companion galaxies within the lensed group, forms a submillimeter arc with total SFR = (450 ± 50) M_⊙ yr⁻¹ (MacKenzie et al. 2014). The group members span two orders of magnitude in SFR, and those with CO detections are expected to deplete their molecular gas within <0.14–1.0 Gyr (Shen et al. 2021). The lensed group in which m0451 resides may be a example of accelerated growth in dense environments. Compact star-forming galaxies at z ≳ 2–3 have been proposed as another progenitor for z < 2 quiescent galaxies (van Dokkum et al. 2015; Barro et al. 2017; Gómez-Guijarro et al. 2019), given their similarly compact stellar sizes, number densities, and modest SFR (∼100 M_⊙ yr⁻¹). While it is plausible that our quenched galaxies went through such a phase of compact star formation, their star-forming progenitors are likely at z ≈ 2−4 as inferred from their best-fit SFHs. Current work on compact star-forming galaxies is limited to z < 3, as HST cannot probe the rest-frame optical sizes of galaxies beyond z = 3. JWST will soon enable a census of z > 3 galaxies in order to identify the progenitors of these early-quenched galaxies.

Aside from number density and size comparisons, one could also gain insights into possible progenitors of our quenched galaxies using timescales. The star formation duration, t₁₀ − t₉₀, of our quenched galaxies is shorter than the molecular gas depletion time of "main-sequence" star-forming galaxies at similar z and stellar mass (0.5–0.7 Gyr; Tacconi et al. 2018) by a factor of a few to more than an order of magnitude. These "main-sequence" star-forming galaxies are thus unlikely to be the immediate progenitor to our quenched galaxy sample. The only distant galaxies with depletion times as short as <200 Myr are submillimeter galaxies (e.g., Bothwell et al. 2013), compact star-forming galaxies (e.g., Spilker et al. 2016; Popping et al. 2017), and starbursting radio galaxies (Man et al. 2019). These galaxies have star formation rates above or within the "main sequence" of star-forming galaxies, but what sets them apart from the "main sequence" is their high star formation rates for their molecular gas mass, i.e., they have high star formation efficiency and an equivalently short gas depletion timescale (Elbaz et al. 2018). The star formation histories of our quenching galaxies are certainly compatible with having experienced such a rapid, efficient phase of star formation prior to quenching.

Having discussed the possible progenitors of our sample in Section 5.2, we now turn to their subsequent evolution in order to fully explore their evolutionary scenario over the next 9–12 Gyr until the present day. The intermediate stellar masses of our sample imply that they would evolve into relatively massive galaxies in the local universe. Identifying their low-redshift descendants is like finding needles in haystacks, as massive galaxies become more numerous and more of them quench their star formation over time (Muzzin et al. 2013; Davidzon et al. 2017). Fortunately, the robust constraints on the stellar populations and morphologies provide clues to infer their evolution.

Archeological studies of nearby massive, early-type galaxies suggest that the bulk of their stars formed early by z ≈ 1−2 (Thomas et al. 2005; McDermid et al. 2015). ²¹ Are the early-quenched galaxies presented in this work the precursors of the metal-rich cores of nearby early-type galaxies? We compare the spatially integrated stellar metallicities of our sample to the resolved measurements of nearby massive early-type galaxies. Numerous integral field spectroscopic surveys provide stellar metallicity gradient measurements. Most such surveys of nearby massive galaxies resolve well within once or twice the effective radii (e.g., SAURON, ATLAS3D, MANGA, SAMI; Kuntschner et al. 2010; González Delgado et al. 2015; Martín-Navarro et al. 2018; Bernardi et al. 2019; Ferreras et al. 2019a; Oyarzún et al. 2019; Krajnović et al. 2020), yet few probe beyond the cores, because of the limited field of view. Thus, we compare our integrated stellar metallicities with results from the MASSIVE and CALIFA surveys that provide gradient measurements to the largest spatial extent, out to ≈2.5 r_eff. We use the gradient fits of the MASSIVE sample (Greene et al. 2015, Table 2), in units of kpc, to compute the stellar metallicity as [Z/H] = [Fe/H] + 0.94 [Mg/Fe] (Thomas et al. 2003, Equation (4)). As for the CALIFA survey, we show the median stellar metallicity gradient of early-type galaxies presented in Zibetti et al. (2020), in bins of semimajor axis in kpc. We also include the stellar metallicity gradient of NGC 4889 for comparison. NGC 4889 is one of the two brightest cluster galaxies in the Coma Cluster (Coccato et al. 2010). Last, we compare our results against the stellar metallicity gradients of the "relic galaxies," i.e., local massive, quiescent galaxies that are compact in size (Trujillo et al. 2014; Ferré-Mateu et al. 2017).

The comparison is shown in Figure 13, where we label our sample at the source-plane effective radii of the spectroscopic images. The stellar metallicities of our sample are lower than those found in the cores of local early-type galaxies. Further chemical enrichment, perhaps by gas-rich mergers and/or star formation rejuvenation, needs to take place if they are to evolve into metal-rich cores of local massive, early-type galaxies. Kriek et al. (2016) reported that a dry minor merging scenario could explain the chemical evolution of a massive, quiescent, Mg-enhanced galaxy into its local counterparts. If the targets presented in this work would merge with more metal-rich quiescent galaxies (Figure 6) in their subsequent evolution, the metal content of the merger remnant would be more in line with z < 1 mass-metallicity relations. There are two potential caveats with regard to the comparison shown in Figure 13. First, it is not straightforward to compare a spatially integrated measurement with a gradient. A better comparison can be made by deriving a spatially integrated measurement from the gradient fit and the luminosity or mass profile. Second, the derived metallicities depend on the calibration (absorption line indices or full spectral fitting) as well as the assumptions involved (e.g., star formation history, metal yield of various stellar types). Resolved elemental abundance analysis is required to address these issues (see Jafariyazani et al. 2020). A detailed comparison addressing these caveats is beyond the scope of this paper.

Zoom In Zoom Out Reset image size

The stellar mass range of our sample suggests that they are more likely to be fast rotators rather than slow rotators, if we take the stellar kinematics and mass distribution of z ≈ 0 as a reference (Emsellem et al. 2011; Veale et al. 2017a). Fast rotators are expected to be more than an order of magnitude more numerous than slow rotators at z ≈ 2 (Khochfar et al. 2011), as dry merging is a dominant mechanism to reduce the spin of galaxies (Naab et al. 2006; Lagos et al. 2018a, 2018b) and should be more prevalent at later cosmic times than wet mergers (Hopkins et al. 2010). Like fast rotators, the apparent axis ratios of our sample span a wide range, as shown in Figure 10. Their low Sérsic index (n < 2) lends further support to them being fast rotators, as all slow rotators have n > 2 (Krajnović et al. 2013). Resolved absorption line spectroscopy is needed to confirm the kinematic nature of our sample. Indeed, m2129 is shown to have rotation-dominated stellar kinematics (Toft et al. 2017; Newman et al. 2018b). JWST will enable us to obtain resolved kinematics for these compact galaxies in the near future.

Future evolution of early-quenched galaxies depends on their ability to rejuvenate star formation, if molecular gas is made available for star formation again, e.g., through gas accretion, mergers. Studies of the star formation histories of z < 1 quenched galaxies do reveal that a small fraction had evidence for rejuvenated star formation in the recent past (Chauke et al. 2018). A recent analysis of a lensed quiescent galaxy at z = 1.9, MRG-S0851, reveals evidence for star formation rejuvenation in the inner kpc within the past ∼100 Myr. If representative of galaxies having similar spectral energy distribution, the abundance implies that ≈1% of massive quiescent galaxy at z = 1–2 are potentially experiencing star formation rejuvenation (Akhshik et al. 2021). In two future works, we will report on the molecular gas content of quenching galaxies including targets of this work (A. Man et al. 2021, in preparation; Whitaker et al. 2021).

Gravitational lensing and deep spectroscopy have provided us an exquisitely deep view into the properties of ordinary galaxies as they quench their star formation. What insights can we gain regarding the mechanism responsible for their decline in star formation activity? Gas needs to be brought into galaxies and sufficiently cool and settle in order to form stars. Star formation quenches, temporarily or permanently, if one or more of these necessary conditions is lacking. In this subsection, we explore how our observations enable us to constrain the cause of quenching.

An important discriminant of star formation quenching mechanisms is the timescale over which they operate. The stellar population analysis in Section 4.1 suggests that our sample has experienced a rapid star formation history, with τ < 0.2 Gyr. An alternative illustration is provided by examining their position on the UVJ diagram, a diagnostic first developed to separate star-forming galaxies from quiescent ones (Wuyts et al. 2007; Williams et al. 2009; Muzzin et al. 2013; Whitaker et al. 2013). Galaxies evolve on the UVJ diagram as they quench their star formation and become old (e.g., Barro et al. 2014; Merlin et al. 2018). In Figure 14, we compare the rest-frame (U – V) and (V – J) colors of our sample with the fast and slow quenching models. The (U – V) and (V – J) colors listed in Table 9 are computed from the SED models. ²² The model tracks shown in Figure 14 are from Belli et al. (2019), assuming two τ models for fast quenching (τ = 0.1 Gyr) and slow quenching (τ = 1 Gyr), respectively. The colors are dust-corrected by assuming an evolving A_V that declines with the SFR over time, starting from 2 to 0.4, and assuming R_V = 4.05 and the Calzetti et al. (2000) dust attenuation law. It is apparent that the slow quenching model is incompatible with our sample: a galaxy that follows the slow quenching track with τ ≳ 1 Gyr would stay star-forming in the first 3 Gyr since the onset of its star formation, i.e., within three times of the e-folding timescale. This corroborates our findings in Section 4.1, which indicate rapid star formation timescales of τ ≲ 0.2 Gyr and t₁₀ − t₉₀ = 0.07 – 0.47 Gyr. This conclusion stands even for higher initial values of A_V or adopting R_V ranging from 3 to 6. The lower (V – J) values of m1423 and m0451 can be attributed to their shorter star formation histories (τ < 0.1 Gyr), stronger burst strength, or possibly different dust correction (see also Barro et al. 2014; Merlin et al. 2018; Wu et al. 2020). Overall, the UVJ diagram provides an independent illustration of a fast star formation quenching scenario, in qualitative agreement with spectroscopic investigations of z ≈ 2 quiescent galaxies (Kriek et al. 2016; Newman et al. 2018a; Zick et al. 2018; Belli et al. 2019).

Zoom In Zoom Out Reset image size

Stellar morphology provides another clue to understanding quenching. The lensing-reconstructed HST images show that their stellar structures are disk-like (Section 4.5). This suggests that quenching does not necessarily require or coincide with bulge formation. Similar findings have been reported in morphological analyses of photometrically selected massive, early-type galaxies (Bruce et al. 2012, z = 1–3; Chang et al. 2013a, z = 0.6–1.8), as well as in the analysis of stellar kinematics of massive, quenched galaxies at z = 2–2.6 (Toft et al. 2017; Newman et al. 2018b). This work confirms that the lack of coincidence between quenching and morphological transformation extends to intermediate-mass galaxies (log(M_⋆/M_⊙) = 10.2–10.6) at z = 1.6–3.2. Their disky morphologies are in line with the prevalence of flat, oblate morphology among photometrically selected quiescent galaxies at z ≈ 2.5 being higher than that of their counterparts in the present day (Stockton et al. 2008; van der Wel et al. 2011; Chang et al. 2013a, 2013b; Hill et al. 2019). In fact, massive UVJ-selected quiescent galaxies at z = 2–3.5 are as flat as star-forming galaxies (Hill et al. 2019). In terms of stellar morphology, they resemble local fast rotators rather than slow rotators. Quenched galaxies at z ≈ 2 may thus have some level of rotation in their stellar velocity fields, as already confirmed in spatially resolved studies (Toft et al. 2017; Newman et al. 2018b). The compact sizes and thus high stellar mass densities, as well as their disky kinematics, all point toward a dissipative formation process such as a gas-rich merger (Cox et al. 2006; Naab et al. 2006, 2014; Wuyts et al. 2010; Lagos et al. 2018b). Although the merger fraction of massive galaxies at z ≳ 2 is less than 10% (Man et al. 2016a), there is emerging evidence that mergers may be more prevalent among quenched galaxies (Glazebrook et al. 2017; Stockmann et al. 2020). The absence of a prominent stellar bulge or bar disfavors morphological quenching (Martig et al. 2009) or bar quenching (Khoperskov et al. 2018) as plausible quenching mechanisms. These processes, however, may help to maintain the quiescence of star formation in the long run. If the targets presented in this work are to evolve into present-day slow rotators, they would need to substantially grow their bulges later on in a process unrelated to quenching, e.g., dry major mergers.

Another clue comes from the gas conditions of quenching galaxies. The faint emission lines in our sample suggests the presence of warm (T ∼ 10⁴ K), low-ionization gas as discussed in Section 4.4. Furthermore, the youngest quenching galaxies show evidence for outflowing warm gas in their Mg ii profiles (Sections 4.3, 4.4). Outflowing warm gas is detected in post-starburst galaxies at z < 1.4 (Tremonti et al. 2007; Coil et al. 2011; Sell et al. 2014; Baron et al. 2017, 2018). Maltby et al. (2019) reported tentative hints of faster outflows as seen in the Mg ii blueshifted absorption in younger post-starburst galaxies. Our findings are in good agreement with the expectation that Mg ii absorption is more likely to be detected in more recently quenched galaxies (≲500 Myr as inferred by their light-weighted ages) than older ones, as shown in a study of galactic wind in K+A galaxies (Coil et al. 2011). Spatially extended, redshifted Lyα emission is present in m0451 (Jauzac et al. 2021, z_Lyα = 2.93) across the three multiple images as detected with the VLT/MUSE spectrograph. The extended Lyα emission in quenched galaxies may originate from photon scattering, galactic winds, or unresolved satellite galaxies (Taniguchi et al. 2015).

While these findings do not readily imply that outflows cause star formation quenching, understanding the relation between the two phenomena could better constrain how quenching took place. Although AGN are commonly thought to drive fast outflows ejecting gas beyond the gravitational potential of massive galaxies (Tremonti et al. 2007; Baron et al. 2017, 2018; Förster Schreiber et al. 2019), compact starbursts may also be capable of doing so (Sell et al. 2014; Rupke et al. 2019). Does the fast outflow become less prevalent as the stellar populations age because their energy sources (young stars and/or AGN) have faded, or is it because previous outflows have cleared them of warm gas? The latter is unlikely to be the case: the Mg ii absorption of the entire sample is deeper than expected from photospheric absorption alone, suggesting that warm gas is present after quenching and not completely evacuated from the galaxies (see Jafariyazani et al. 2020 for a related discussion on NaD absorption). On the other hand, SMBH accretion rates can vary over much shorter timescales than galaxy star formation (Hickox et al. 2014). It is plausible that an AGN outflow persists for ∼100 Myr after the AGN has turned off (King et al. 2011; Zubovas & King 2014). There is also ample evidence for the peak of AGN accretion to occur ≳100–250 Myr after the peak of star formation (Schawinski et al. 2009; Wild et al. 2010; Volonteri et al. 2015a, 2015b; Baron et al. 2018; Falkendal et al. 2019). The answers to these questions are fundamental to understanding how star formation quenching proceeds. Observations with the JWST and other multiwavelength facilities will tackle these questions.

What do our results imply for how galaxies quench their star formation? A different perspective can be gained by asking whether a mechanism is needed to stop galaxies from forming stars, or if galaxies simply have an accelerated star formation episode and stay quiescent thereafter (Abramson et al. 2016). It is plausible that the mechanism that quenches star formation of galaxies is not the one responsible for maintaining their quiescence, given the very different timescales considered. Absorption line spectroscopy of z ≳ 2 massive, quenching galaxies including our sample provides solid indications that they have experienced a rapid star formation history in the first few billion years after the Big Bang. If quenching is simply a rapid decline of star formation without the need for an external actor, the question then becomes what causes a rapid gas conversion to stars ("starburst"). Compressive motions and effective angular momentum loss of gas can facilitate starbursts, and so can gas-rich major mergers, both within a Gyr or so (Barnes & Hernquist 1991; Hopkins et al. 2008). Gas outflows driven by AGN and/or starbursts can act in concert with rapid gas consumption during starburst to rapidly clear galaxies of cold gas (Man et al. 2019, and references therein). Simulations by Su et al. (2019) demonstrate that stellar feedback, morphological feedback, magnetic field, thermal conduction, and stellar cosmic rays do not effectively quench the star formation of massive galaxies. This supports the idea that other processes like AGN feedback are required, although there is a broad variety of AGN feedback implementations in simulations (e.g., Sijacki et al. 2007; Booth & Schaye 2009; Dubois et al. 2012; Su et al. 2020). This work presents tentative evidence that, as star formation quenches, the mode of SMBH growth transitions from radiative-mode to jet-mode (Section 4.4.3). This is in good agreement with low-z studies (Best & Heckman 2012), and corroborates findings of faint AGN in massive quiescent galaxies at z ≳ 1 (Olsen et al. 2013; Man et al. 2016b; Barišić et al. 2017; Aird et al. 2019). More observations are needed to inform whether and how AGN affect star formation of host galaxies, particularly faint AGN that are more ubiquitous than the quasar phase. On the other hand, figuring out their future evolution until the present day would require the knowledge of whether molecular gas is present or can be made available again to form stars. In a future work, we will quantify the amount of molecular gas of our sample using ALMA observations. To identify whether and how gas cools, we need measurements of the gas temperature, density, and metallicity of various gas phases to constrain the cooling functions. These measurements are crucial to properly characterize the gas conditions in order to understand how and why massive galaxies experience a decline in their star formation. Ultimately, we need deep, multiwavelength spectroscopic observations to constrain the relevant timescales and trace any past and ongoing star formation activity, gas accretion and outflows, and AGN activity. Only then can we begin to disentangle the intricate relations between stars, gas, black holes, and how they shape the evolution of galaxies.