<figure><img src="image_1.png" /></figure>

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

Allison W. S. Man<sup>1,2</sup>①, Johannes Zabl<sup>3,4</sup>①, Gabriel B. Brammer<sup>5,6</sup>①, Johan Richard<sup>3</sup>①, Sune Toft<sup>5,6</sup>①, Mikkel Stockmann<sup>5,6</sup>①, Anna R. Gallazzi<sup>7</sup>①, Stefano Zibetti<sup>7</sup>①, and Harald Ebeling<sup>8</sup>①

Received 2020 September 2; revised 2021 May 31; accepted 2021 June 10; published 2021 September 17

## Abstract

This work presents an in-depth analysis of four gravitationally lensed red galaxies at $z = 1.6–3.2$. The sources are magnified by factors of 2.7–30 by foreground clusters, enabling spectral and morphological measurements that are otherwise challenging. Our sample extends below the characteristic mass of the stellar mass function and is thus more representative of the quiescent galaxy population at $z > 1$ than previous spectroscopic studies. We analyze deep VLT/X-SHOOTER spectra and multiband Hubble Space Telescope photometry that cover the rest-frame UV-to-optical regime. The entire sample resembles stellar disks as inferred from lensing-reconstructed images. Through stellar population synthesis analysis, we infer that the targets are young (median age = 0.1–1.2 Gyr) and formed 80% of their stellar masses within 0.07–0.47 Gyr. Mg II $\lambda\lambda$ 2796, 2803 absorption is detected across the sample. Blueshifted absorption and/or redshifted emission of Mg II are found in the two youngest sources, indicative of a galactic-scale outflow of warm ($T \sim 10^4$ K) gas. The [O III] $\lambda$5007 luminosity is higher for the two young sources (median age less than 0.4 Gyr) than the two older ones, perhaps suggesting a decline in nuclear activity as quenching proceeds. Despite high-velocity ($v \approx 1500$ km s$^{-1}$) galactic-scale outflows seen in the most recently quenched galaxies, warm gas is still present to some extent long after quenching. Altogether, our results indicate that star formation quenching at high redshift must have been a rapid process ($<1$ Gyr) that does not synchronize with bulge formation or complete gas removal. Substantial bulge growth is required if they are to evolve into the metal-rich cores of present-day slow rotators.

*Unified Astronomy Thesaurus concepts: Galaxy evolution (594); Galaxy processes (614); Galaxy chemical evolution (580); Galaxy quenching (2040); Galaxy structure (622); Stellar populations (1622); Supermassive black holes (1663); Metallicity (1031); Warm ionized medium (1788); Active galactic nuclei (16); Galaxy ages (576); Galaxy stellar content (621)*

## 1. Introduction

Understanding how and when galaxies form stars is a cornerstone of galaxy evolution studies. Much effort has been dedicated to studying the most massive galaxies, as they are the brightest and thus easiest to observe. The most massive galaxies (with stellar masses a few times the characteristic mass of the stellar mass function) by nature have complex formation and accretion histories. The early, rapid growth of massive galaxies is driven by high gas accretion rates, fueling in situ star formation at $z \gtrsim 2$, followed by a later, longer phase of mass growth characterized by frequent merging with satellite galaxies (Oser et al. 2010; Lackner et al. 2012; Hirschmann et al. 2015; Rodriguez-Gomez et al. 2016; Wellons et al. 2016). This two-phase scenario provides a decent explanation for a broad set of observations of nearby, massive, early-type galaxies, including their stellar light profiles (D’Souza et al. 2014), element abundance ratios (Thomas et al. 2005), stellar metallicity and age gradients (Greene et al. 2013; Ferreras et al. 2019a; Zibetti et al. 2020), and diverse stellar velocity fields (Emsellem et al. 2011; Krajnović et al. 2013; Veale et al. 2017b).

The evolution of intermediate-mass galaxies (with stellar masses near the characteristic mass of the stellar mass function), on the other hand, is less constrained by observations than massive galaxies. The number density of quiescent galaxies declines steeply with stellar mass beyond the characteristic mass (e.g., Muzzin et al. 2013; Davidzon et al. 2017), i.e., intermediate-mass galaxies are orders of magnitude more numerous than the most massive ones. The relative dominance of *in situ* star formation and *ex situ* mass accretion in intermediate-mass galaxies is a matter of contention among various simulations (Oser et al. 2010; Lackner et al. 2012; Rodriguez-Gomez et al. 2016; Moster et al. 2020). Intermediate-mass quiescent galaxies at $z > 2$, observed shortly after these galaxies shut down their star formation, can thus be used to discern between different star formation and feedback models.

Related to understanding how galaxies form stars, an outstanding question in galaxy evolution studies is why galaxies experience a decline in star formation, a process loosely referred to as “quenching.” Spectroscopic confirmation of massive, quiescent galaxies at $z \gtrsim 2$ provides evidence for fast quenching of individual galaxies (e.g., Glazebrook et al. 2017; Forrest et al. 2020;

<footnote n="1"><sup>1</sup> Dunlap Institute for Astronomy and Astrophysics, University of Toronto, 50 St George Street, Toronto ON, M5S 3H4, Canada; aman@phas.ubc.ca, allisonmanws@gmail.com</footnote>

<footnote n="2"><sup>2</sup> Department of Physics & Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver BC, V6T 1Z1, Canada</footnote>

<footnote n="3"><sup>3</sup> Univ Lyon, Univ Lyon1, Ens de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230 Saint-Genis-Laval, France</footnote>

<footnote n="4"><sup>4</sup> Institute for Computational Astrophysics and Department of Astronomy & Physics, Saint Mary's University, 923 Robie Street, Halifax NS, B3H 3C3, Canada
<sup>5</sup> Cosmic Dawn Center (DAWN)</footnote>

<footnote n="6"><sup>6</sup> Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2100 Copenhagen, Denmark</footnote>

<footnote n="7"><sup>7</sup> INAF—Osservatorio Astrofisico di Arcetri, Largo Enrico Fermi 5, I-50125 Firenze, Italy</footnote>

<footnote n="8"><sup>8</sup> Institute for Astronomy, University of Hawai'i, Honolulu, HI 96822, USA</footnote>
Stockmann et al. 2020; Valentino et al. 2020). Their high stellar masses ($log(M_{*}/M_{\odot}) \gtrsim 11$) imply that they had a much higher star formation rate in the past, followed by a rapid decline in star formation rate that creates the absorption lines typical of post-starburst galaxies without emission of OB-type stars. Active galactic nucleus (AGN) feedback is a common explanation for fast quenching, at least in theory: supermassive black hole accretion creates winds that expel gas from galaxies (momentum injection) and/or heat gas (energy injection), leaving little cold gas available for further star formation (Alexander & Hickox 2012; Fabian 2012, and references therein). Because of the complex, multiscale nature of black hole-gas interaction, AGN feedback is among the most debated topics in galaxy studies, and consequently its implementations vary drastically across simulations (Di Matteo et al. 2005; Springel et al. 2005; Dubois et al. 2012; Sijacki et al. 2015; Croton et al. 2016; Bower et al. 2017).

On the other hand, a *population* of galaxies is said to have quenched as the fraction of red (passive or quiescent) galaxies becomes higher toward lower redshifts (e.g., Peng et al. 2010; Muzzin et al. 2013; Davidzon et al. 2017). While there is certainly a causal relation between the star formation histories and quenching of *individual* and a *population* of galaxies, they are not equivalent to each other, as the latter also involves population changes (e.g., formation of blue galaxies that cross a mass threshold, mergers). The two different methodologies result from the empirical nature of galaxy evolution studies: galaxies evolve on such long timescales that astronomers must infer their evolution by reconstructing their star formation and assembly histories through either spectral features or population changes. The former approach requires deep spectroscopy and is thus limited to small samples, while the latter approach is less accurate in separating star-forming galaxies from quiescent ones but bears the merit of larger samples that span a range of environments. The two are different yet complementary approaches to understanding how galaxies form stars and quench.

The question of what quenches galaxies can be posed another way: what explains the diversity of star formation histories among galaxies? Massive, quiescent galaxies formed the bulk of their stars earlier on than less massive ones, in general (Thomas et al. 2005, 2010; Gallazzi et al. 2006; Domínguez-Sánchez et al. 2016; Chauke et al. 2018; Wu et al. 2018). Star-forming galaxies have a more extended duration of star formation than quiescent galaxies at later times at a given mass (Ferreras et al. 2019b). It remains debated whether quenching is triggered by an event (e.g., AGN) or it is simply due to the exhaustion of star-forming fuel. Cosmological simulations in general require an energy source like AGN to quench star formation and to prevent galaxies from forming too many stars too early and efficiently (Su et al. 2019, and references therein). Some studies suggest that feedback can be provided by other sources like old stellar populations, i.e., Type Ia supernovae (Mathews 1990; Ciotti et al. 1991) or asymptotic giant branch stars (Conroy et al. 2015). Alternatively, gas can be prevented from cooling sufficiently by virial shocks (Birnboim & Dekel 2003) or gravitational heating (Khochfar & Ostriker 2008). Other models suggest that galaxy star formation can be made inefficient in the presence of a stellar bar or bulge (Martig et al. 2009; Khoperskov et al. 2018). A large number of quenching mechanisms have been proposed in the literature (see Man & Belli 2018, for an extensive list). The answer as to what quenches star formation in galaxies

likely depends on the timescale (early or late cosmic epochs, fast or slow) as well as the conditions of the galaxies (mass, environment, age of stellar populations). To distinguish between the many quenching mechanisms proposed, we need more constraints on quiescent galaxies beyond what is already known from correlations.

The most stringent constraints of distant quiescent galaxies come from deep spectroscopy. Absorption line strengths and the overall continuum shape provide valuable information on the age and metallicity of the integrated stellar population of galaxies, as well as their star formation histories. Several hours of integration time with eight-meter class telescopes are generally needed to secure detections of multiple absorption lines for constraining stellar populations. The advent of sensitive near-infrared spectrographs has enabled absorption line studies for several tens of $z > 1.5$ quiescent galaxies (Kriek et al. 2009, 2016; van Dokkum et al. 2009; Onodera et al. 2010, 2012; van de Sande et al. 2011, 2013; Toft et al. 2012; Bezanson et al. 2013; Belli et al. 2014, 2017; Glazebrook et al. 2017; Kado-Fong et al. 2017; Schreiber et al. 2018; Estrada-Carpenter et al. 2019; D'Eugenio et al. 2020; Forrest et al. 2020; Stockmann et al. 2020; Valentino et al. 2020). These results have confirmed the existence of distant quiescent galaxies with low star formation rates (through nondetections or weak detections of $H\alpha$ emission and [O II]), evolved stellar ages ($\approx 0.5-1$ Gyr), and in some cases high stellar metallicity. Despite the representative nature of intermediate-mass galaxies and their potential to constrain galaxy evolution models, no deep spectroscopic observations exist for quiescent galaxies at $z > 2$ around the characteristic masses to date.

Gravitational lensing is a promising way to deepen our understanding of the faint, compact quenching galaxies that are elusive. Foreground galaxy clusters can magnify background galaxies and amplify their fluxes, enabling us to obtain higher signal-to-noise (S/N) spectra within reasonable integration times for these otherwise faint galaxies. This allows us to probe intrinsically fainter galaxies that are more representative of the quiescent galaxy population than the most luminous ones (e.g., Stockmann et al. 2020). Lensed galaxies can also be magnified, so we can achieve a more resolved view of the starlight of these compact galaxies that are barely resolvable with the Hubble Space Telescope (HST). The increasing rarity of quiescent galaxies toward higher redshift implies that cluster-lensed quiescent galaxies are uncommon, typically less than one per galaxy cluster. Only nine such galaxies have been discovered to date (Muzzin et al. 2012; Geier et al. 2013; Newman et al. 2015, 2018a, 2018b; Hill et al. 2016; Toft et al. 2017; Ebeling et al. 2018; Akhshik et al. 2020). The most remarkable finding is the confirmation of a rotating stellar disk that has ceased star formation (Toft et al. 2017; Newman et al. 2018b).

We have systematically searched for lensed quiescent galaxies behind galaxy clusters for spectroscopic follow-up. This paper presents the analysis of a pilot sample and is structured as follows. Section 2 presents the target selection and provides a description of each source. Section 3.1 presents the details of the photometric measurements. Section 3.2 presents the details on the lens model and source reconstructions. In Section 3.3, we present the procedures for reducing the VLT/X-SHOOTER spectra. Section 3.4 describes the procedures for the stellar population synthesis fitting. In Section 4, we present the results of our analysis, while in Section 5, we discuss the
<figure><img src="image_2.png" /><figcaption>**Figure 1.** HST color images of the four spectroscopic targets and their respective foreground lensing clusters. The lensed quiescent galaxies are shown with white rectangles. The composite images are made using the F160W/F110W/F814W filter combination represented by red/green/blue colors on a logarithmic scale for m2129, m0451, and m1423. For e1341, the F140W/F110W/F814W filter combination is used instead as F160W observation is not available for that source. The images measure 90'' (e1341 and m0451) or 120'' (m2129 and m1423) on each side.</figcaption></figure>

broader implications of our findings. Section 6 summarizes our conclusions.

All magnitudes are quoted in the AB system (Oke & Gunn 1983), unless otherwise stated. The common logarithm with base 10 is used unless otherwise stated. We assume a Kroupa (2001) initial mass function. A cosmology of $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega_M = 0.3$, and $\Omega_\Lambda = 0.7$ is adopted throughout the paper.

## 2. Sample

The objective of our survey is to identify and to conduct spectroscopic follow-up of lensed, distant quiescent galaxies. We searched the HST archive for galaxy cluster imaging and identified red galaxies with broadband colors consistent with quiescent stellar populations at $z > 1$ that are gravitationally lensed by foreground clusters and bright enough for spectroscopic follow-up. Our spectroscopic observation campaign spanned over 12 semesters, and we refined our selection criteria over time. In general, we use infrared imaging to identify red

sources bright enough for spectroscopic follow-up and to construct reliable lens models. Source m2129 was selected using a distant red galaxy criterion ($J - K_S > 1.36$ (Franx et al. 2003) based on VLT/ISAAC imaging as described in Geier et al. (2013). Source m0451 was identified to be red in its HST/WFC3 color ($J_{F110W} - H_{F160W} = 1.35$) and known to be at $z \approx 3$ based on the lensing configuration. Source m1423 was selected because of its high photometric redshift ($z_{phot} > 3$) and low specific star formation rate (log(sSFR) $\approx -10$) based on spectral energy distribution fitting of the CLASH data, although these numbers have been revised with the spectroscopic analysis. Source e1341 was discovered from an HST Snapshot program of massive clusters (PI: Ebeling) as a spectacular candidate: the target is red, based on its Gemini/GMOS ($g'$, $r'$, and $i'$) and HST/WFC3 (F110W and F140W) images. Its spectral energy distribution is consistent with being a quiescent galaxy at $z \approx 1.35$ that is triply lensed, with one image close to the critical line, and has extremely high magnification (Ebeling et al. 2018). This paper presents the
<table><caption>**Table 1**
Gravitational Lensing Properties</caption><thead><tr><th></th><th>e1341</th><th>m2129</th><th>m0451</th><th>m1423</th></tr></thead><tbody><tr><td>z<sub>spec</sub></td><td>1.5954</td><td>2.1487</td><td>2.9223</td><td>3.2092</td></tr><tr><td>Lensing cluster</td><td>eMACSJ1341.9-2442</td><td>MACSJ2129.4-0741</td><td>MACSJ0454.1-0300<br>(a.k.a. MS0451.6-0305)</td><td>MACSJ1423.8-2404</td></tr><tr><td>z<sub>lens</sub></td><td>0.835</td><td>0.5889</td><td>0.5377</td><td>0.5431</td></tr><tr><td>Magnitude</td><td>18.5, 19.8, 20.4</td><td>20.1</td><td>20.9, 21.8, 22.0</td><td>22.1</td></tr><tr><td>μ</td><td>30 ± 8, 12.6 ± 1.5, 8.2 ± 0.9</td><td>4.6 ± 0.2</td><td>10.9 ± 2.1, 6.1 ± 0.6, 3.8 ± 0.3</td><td>2.7 ± 0.2</td></tr><tr><td>Mass model reference</td><td>Ebeling et al. (2018)</td><td>Toft et al. (2017)</td><td>Jauzac et al. (2021)</td><td>Limousin et al. (2010)</td></tr><tr><td>rms</td><td>0.45</td><td>0.60</td><td>0.73</td><td>0.82</td></tr></tbody></table>

Notes. The total AB magnitudes as measured in HST F160W (or F140W for e1341) are reported. The lensing magnification factors (μ) scale linearly with flux. The underlined values indicate the spectroscopic targets. The rms refers to the root mean square of the observed image positions of the multiple images for each lensing system.

analysis of our pilot sample, comprising four lensed quiescent galaxies with sufficiently high S/N spectra for absorption line analysis. The color images of the sample are shown in Figure 1. The targets were observed with the X-SHOOTER spectrograph (Vernet et al. 2011) at the Very Large Telescope. The observations and data reduction procedures are described in Section 3.3. Below, we briefly discuss the individual properties of each target, and we provide an overall summary of the lensing properties in Table 1.

## 2.1. Targets

### 2.1.1. Source e1341

Source e1341 is a triply imaged quiescent galaxy with the highest magnification factor known to date (Ebeling et al. 2018). The X-SHOOTER spectrum presented in this paper is based on observations conducted on 2017 July 24, 25, 26, 27, August 6, and 2018 March 21, 28 under ESO program ID 099. B-0912 (PI: M. Stockmann).

### 2.1.2. Source m2129

Source m2129 is a singly lensed massive galaxy at $z=2.1$ with quiescent SFR. Follow-up deep absorption line spectroscopy with VLT/X-SHOOTER was first presented in Geier et al. (2013). A more recent analysis of the stellar kinematics shows that m2129 is the first quiescent galaxy at $z>2$ confirmed to have rotation-dominated stellar kinematics (Toft et al. 2017; Newman et al. 2018b). The analysis presented here is based on the X-SHOOTER observations conducted on 2011 May 12, 14 and September 5, 14 under ESO program ID 087.B-0812 (PI: S. Toft), as also used in Geier et al. (2013) and Toft et al. (2017).

### 2.1.3. Source m0451

Source m0451 is situated in a compact group at $z \approx 2.9$ with at least eight members separated within <100 kpc in projection (MacKenzie et al. 2014; Shen et al. 2021). The quenched galaxy studied in this work is triply lensed, as are its star-forming companion galaxies, and together they form a giant submillimeter arc with a total unlensed SFR of $\approx(450 \pm 50) M_{\odot}yr^{-1}$. The system was originally identified by Takata et al. (2003), and the lensing model was presented in Zitrin et al. (2011) as multiple system 4. We adopt a magnification factor of $\mu = 10.9 \pm 2.1$ for our spectroscopic target as estimated by the latest lens model (Jauzac et al. 2021, ID G.2). The

X-SHOOTER spectrum analyzed in this work is based on the observations conducted on 2011 September 11, 13, 2015 December 5, 9, and 2016 January 9 under ESO programs ID 087.B-0812 and 096.B-0994 (PI: S. Toft).

### 2.1.4. Source m1423

Source m1423 is a singly lensed system. The X-SHOOTER spectrum presented here is based on observations conducted on 2014 April 18, May 9, 10, June 5, 8, 2017 March 26, April 12, July 1, 2, August 5, 6, and 2018 February 21, 22, 23, March 12, under ESO programs ID 093.B-0815 (PI: A. Man) and 097.B-1064 (PI: M. Stockmann).

## 3. Method

### 3.1. Photometry

#### 3.1.1. HST Photometry

All targets are selected from the clusters identified from the (extended) MAssive Cluster Survey (MACS and eMACS; Ebeling et al. 2001, 2007, 2010) targeting X-ray luminous, high-redshift clusters. Two of the clusters, MACS2129 and MACS1423, are part of the Cluster Lensing And Supernova survey with Hubble (CLASH; Postman et al. 2012). MS0451-03 has six-band HST imaging available through the Herschel Lensing Survey (Egami et al. 2010; Jauzac et al. 2021).

The HST imaging observations are processed with the grizli analysis software (Brammer 2019).<fnref n="9" /> Briefly, we first retrieve all available HST Advanced Camera for Surveys (ACS) and Wide-Field Camera 3 (WFC3) images described above that cover a given target. These observations are divided into associations defined as sets of exposures of a target in a given instrument and filter combination obtained at a single epoch (i.e., a “visit” in the Hubble observation planning parlance). We align the absolute astrometry of each visit to stars in the GAIA DR2 catalog (Gaia Collaboration et al. 2018), accounting for proper motions between the GAIA and HST observation epochs. Finally, we create mosaics in each instrument/filter combination using the AstroDrizzle software package, where we use the dithered exposures to identify additional cosmic rays and hot pixels with Astro-Drizzle parameters similar to those used by Momcheva et al. (2016). All filters of a given target are drizzled to a common 0.06 pixel grid.

<footnote n="9">⁹ https://github.com/gbrammer/grizli</footnote>
The total magnitude of the reddest HST filter is measured from the automatic aperture with the *SExtractor* code (Bertin & Arnouts 1996). Color magnitudes are measured from small apertures of diameter 1"/2 to optimize S/N, and then scaled to match the total magnitude of the reddest filter. We include the correlated noise in the noise budget introduced by *MultiDrizzle*, following the procedure described in Casertano et al. (2000). The photometry is listed in Table 8.

### 3.1.2. IRAC and WISE Photometry

All but one of our targets have Spitzer/IRAC 3.6 and 4.5 µm observations from the SURFSUP survey (Bradač et al. 2014; Huang et al. 2016). Given the significantly coarser spatial resolution of IRAC (point-spread function full width at half maximum PSF FWHM ≈ 1"/7 for channels 1 and 2) compared to HST, and that source blending is severe due to the high source density in cluster environments, aperture photometry is inadequate to measure IRAC fluxes. Here, we use the *TPHOT* code version 2.0.8¹⁰ (Merlin et al. 2015, 2016) to obtain deblended IRAC flux measurements. We use the reddest HST images, F160W, as priors for the morphology and brightness to deblend the IRAC emission. We first derive a convolution kernel to match the F160W images to the IRAC images based on their PSFs. PSFs for in-flight IRAC observations are used. The F160W images are then convolved with the respective kernels, which are then fitted to the observed IRAC emission maps. The best-fitting multiplicative scaling factor for each object is the deblended IRAC fluxes. The IRAC fluxes and errors are listed in Table 8. For e1341, we retrieve Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) images from the AllWISE data release. WISE has an angular resolution of 6"/1, 6"/4, 6"/5, and 12"/0 at 3.4, 4.6, 12, and 22 µm, with point-source sensitivity of 0.08, 0.11, 1, and 6 mJy in unconfused regions (5σ). We perform the deblending procedure with *TPHOT* as was done for IRAC, using the reddest image, F140W, as prior. The deblended WISE flux measurements are listed in Table 8. We convert the deblended WISE intensities to magnitudes and fluxes in Jy using the zero points derived for a flat spectrum (i.e., $F_ν$ is constant). In the AllWISE source catalog, W3 is shallow and shows only a weak detection (S/N = 3.3) detection at the position of im1 prior to deblending. After running *TPHOT*, the W3 fluxes of all three multiple images have S/N < 1.5 and are therefore individually undetected. W4 has no significant detection at all image positions, and the *TPHOT* fit did not converge. Thus, we consider W3 and W4 nondetections, given their shallow sensitivity.

### 3.2. Magnifications and Source Reconstructions

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.

### 3.3. Spectroscopic Observations and Data Reduction

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 1"/0 in the UVB arm, and 0"/9 for the VIS and NIR arms. This corresponds to nominal spectral resolutions<fnref n="11" /> of R ≈ 5400, 8900, and 5600 for the UVB, VIS, and NIR arms, respectively. The observation log is provided in Table 2.

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.<fnref n="12" /> Third, the pipeline version being used creates small-scale ring-like features around bad pixels during rectification. To avoid problems with these rings, we

<footnote n="11">¹¹ The quoted R values are strictly lower limits, as these values are derived for sources whose emission completely fill the slit. Our sources are generally compact, so R can be slightly higher than stated. The improvement in R is very small, compared to the intrinsic velocity dispersion of our targets as stated in Table 3, and does not affect our conclusions in any significant way.</footnote>

<footnote n="12">¹² Only detector rows covering the bluest NIR orders have enough unexposed pixels for this self-calibration to work. Fortunately, these are the orders most strongly affected by the issue.</footnote>

<footnote n="10">¹⁰ https://tphot.wordpress.com</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Rectified 2D X-SHOOTER spectra in rest-frame (top panel) and observed frames (bottom panel), ordered by descending redshifts from top to bottom. The top panel shows the spectra in linear grayscale for the rest-frame wavelength range from 0.35 to 0.53 µm (and for e1341 and m2129 also from 0.650 to 0.663 µm). The bottom panel shows the spectra for the observed wavelength range from 0.41 to 2.3 µm. The color scaling has three separate linear parts (black-blue/blue-gray/ gray-white) such that both high and low surface brightness parts are visualized. Important lines are annotated in the top panel and shown as vertical dashed orange lines. The ranges shown in the top panel are indicated in the bottom panel by the vertical yellow and green dotted lines. All spectra are shown for a spatial extent of ±2′′5. The outermost parts are affected by the nodding reduction and therefore do not fully represent the actual light profile there.</figcaption></figure>

<table><caption>**Table 2**
VLT/X-SHOOTER Observation Log</caption><thead><tr><th></th><th>e1341 (im1)</th><th>m2129</th><th>m0451 (im1)</th><th>m1423</th></tr></thead><tbody><tr><td>R.A. (J2000)</td><td>13h42m02s374</td><td>21h29m22s34</td><td>04h54m12s65</td><td>14h23m46s97</td></tr><tr><td>Decl. (J2000)</td><td>-24d41m52s11</td><td>-07d41m31s06</td><td>-03d01m16s50</td><td>+24d05m28s52</td></tr><tr><td>Integration time (h) for UVB/VIS/NIR</td><td>4.14/4.43/4.80</td><td>2.29/1.74/2.67</td><td>4.12/4.36/7.87</td><td>11.64/12.04/13.33</td></tr><tr><td>S/N (Å<sup>-1</sup>)</td><td>11.6</td><td>6.3</td><td>15.5</td><td>9.4</td></tr><tr><td>Slit angle (north to east)</td><td>10°</td><td>-10°</td><td>35°</td><td>-80°, -10°</td></tr></tbody></table>

**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).

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 -0′′5 to +0′′5 around the pointing center.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Extracted 1D X-SHOOTER spectra of our sample, sorted by descending redshifts and increasing stellar population ages from top to bottom. Relevant features are denoted by vertical dotted lines. The regions near the [O II] and [O III] emission lines are colored orange. The gray shades denote the flux errors. The spectra are shifted to their rest frame, and scaled to the flux at 4200 Å and arbitrarily shifted upward for visualization. The shown spectra are binned using a variance-weighted mean. The width of each bin is chosen so that the S/N per bin is similar over the full spectral range. The median bin sizes are rest-frame 2.9, 1.5, 1.9, and 2.3 Å for m1423, m0451, m2129, and e1341, respectively.</figcaption></figure>

### 3.4. Stellar Population Fitting

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),

$$SFR(t) \propto (t - t_0)e^{-(t - t_0)/\tau}. \quad (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-$\tau$ model (Belli et al. 2019). The median stellar age, $t_{50}$, is defined as the lookback time from the spectroscopic redshift $z_{spec}$ at which 50% of the $M_*$ has been assembled. Here, $t_{10} - t_{90}$ is defined as the duration between 10% and 90% of the $M_*$ being assembled. The provided uncertainties of $t_{50}$ and $t_{10} - t_{90}$ are derived from posterior distributions of the SFH parameters.

In addition to the SFH parameters, we fit for the systemic redshift, stellar velocity dispersion ($\sigma_*$), stellar mass normalization (Kroupa 2001 initial mass function), stellar metallicity ($Z_*/Z_\odot$)$^{13}$, 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.

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] $\lambda$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 MultiN-est 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.

## 4. Results

### 4.1. Rapid Star Formation and Old Ages

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

<footnote>$^{13}$ The solar metallicity, $Z_\odot$, is defined as 0.02 following the definition of Bruzual & Charlot (2003).</footnote>
<table><caption>**Table 3**
Stellar Population Properties</caption><thead><tr><th></th><th>z<sub>spec</sub></th><th>log(μM<sub>*</sub>/M<sub>☉</sub>)</th><th>log(sSFR/yr<sup>-1</sup>)</th><th>A<sub>V</sub></th><th>log(Z<sub>*</sub>/Z<sub>☉</sub>)</th><th>σ<sub>*</sub><br>(km s<sup>-1</sup>)</th><th>t<sub>0</sub><br>(Gyr)</th><th>τ<br>(Gyr)</th><th>t<sub>50</sub><br>(Gyr)</th><th>t<sub>10</sub> - t<sub>90</sub><br>(Gyr)</th></tr></thead><tbody><tr><td>e1341</td><td>1.5954 ± 0.0001</td><td>11.62 ± 0.07</td><td>-11.20 ± 0.46</td><td>0.4 ± 0.1</td><td>-0.12 ± 0.07</td><td>160 ± 9</td><td>1.43<sup>+0.06</sup><sub>-0.08</sub></td><td>0.14<sup>+0.01</sup><sub>-0.03</sub></td><td>1.19 ± 0.05</td><td>0.47 ± 0.07</td></tr><tr><td>m2129</td><td>2.1487 ± 0.0002</td><td>11.86 ± 0.05</td><td>-12.77 ± 1.83</td><td>1.1 ± 0.1</td><td>-0.10 ± 0.11</td><td>318 ± 28</td><td>0.68 ± 0.06</td><td>0.04<sup>+0.02</sup><sub>-0.01</sub></td><td>0.61 ± 0.05</td><td>0.14 ± 0.05</td></tr><tr><td>m0451</td><td>2.9223 ± 0.0001</td><td>11.65 ± 0.06</td><td>-9.77 ± 0.08</td><td>0.9 ± 0.1</td><td>-0.30 ± 0.07</td><td>170 ± 12</td><td>0.42 ± 0.03</td><td>0.051 ± 0.005</td><td>0.33 ± 0.02</td><td>0.17 ± 0.01</td></tr><tr><td>m1423</td><td>3.2092 ± 0.0002</td><td>11.01 ± 0.07</td><td>-8.43 ± 0.10</td><td>0.8 ± 0.1</td><td>-0.01 ± 0.10</td><td>277 ± 19</td><td>0.15 ± 0.01</td><td>0.021<sup>+0.002</sup><sub>-0.001</sub></td><td>0.12 ± 0.01</td><td>0.07 ± 0.01</td></tr><tr><td>Prior</td><td>N(z<sub>0</sub>, 0.001(1 + z))</td><td>N(11.3, 0.5)</td><td>...</td><td>Uniform(0,2)</td><td>Z/Z<sub>☉</sub>: N(1, 0.5)</td><td>N(σ<sub>0</sub>, 50)</td><td>N(0.8, 0.3)</td><td>N(0.1, 0.2)</td><td>...</td><td>...</td></tr><tr><td>Limits</td><td>...</td><td>(10, 12.8)</td><td>...</td><td>...</td><td>Z/Z<sub>☉</sub>: (0.2, 2)</td><td>(100, 450)</td><td>(0.05, t<sub>univ</sub>)</td><td>(0.02, 0.3)</td><td>...</td><td>...</td></tr></tbody></table>

**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<sub>☉</sub>, 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 N(μ, σ) indicates a Gaussian distribution with mean μ and standard deviation σ. Numbers z<sub>0</sub> and σ<sub>0</sub> 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<sub>50</sub>, and t<sub>10</sub> - t<sub>90</sub> arising from the priors on SFH parameters t<sub>0</sub> and τ as specified.
<table><caption>**Table 4**
Delensed Stellar Population Properties</caption><thead><tr><th></th><th>log(M<sub>*</sub>)</th><th>log(SFR<sub>100</sub>)</th><th>M<sub>*</sub>/M<sup>*</sup><sub>a</sub></th><th>M<sub>*</sub>/M<sup>*</sup><sub>q</sub></th></tr></thead><tbody><tr><td>e1341</td><td>10.15 ± 0.14</td><td>-1.0<sup>+0.3</sup><sub>-0.7</sub></td><td>0.2 ± 0.1</td><td>0.3 ± 0.1</td></tr><tr><td>m2129</td><td>11.20 ± 0.05</td><td>-1.6<sup>+1.4</sup><sub>-2.5</sub></td><td>2.7 ± 0.5</td><td>2.3<sup>+1.0</sup><sub>-0.7</sub></td></tr><tr><td>m0451</td><td>10.61 ± 0.10</td><td>0.8 ± 0.1</td><td>0.4 ± 0.1</td><td>0.7<sup>+0.4</sup><sub>-0.5</sub></td></tr><tr><td>m1423</td><td>10.58 ± 0.08</td><td>2.2 ± 0.1</td><td>0.1 ± 0.1</td><td>0.6<sup>+0.6</sup><sub>-1.0</sub></td></tr></tbody></table>

**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_{\odot}$ and $M_{\odot}$ yr$^{-1}$, respectively. We use the stellar mass function presented in Muzzin et al. (2013) to estimate the characteristic masses of all galaxies ($M_{\odot}^*$) and those of quiescent galaxies only ($M_{\odot}^*$) at their respective redshifts.

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 \gtrsim 1.5$ that are more representative of the overall galaxy population, as illustrated in Figure 4.

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_{10} - t_{90}$, 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.

Our spectral fitting assumes, like many other works, a parameterized SFH. Specifically, our work assumes the commonly used delayed-$\tau$ 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)-$\tau$ 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-

<figure><img src="image_5.png" /></figure>

**Figure 4.** Stellar masses and redshifts of our sample (large colored symbols), compared with samples of quiescent galaxies with high-quality spectra for absorption line analysis (Kriek et al. 2016, 2019; Bezanson et al. 2018; Newman et al. 2018a; Schreiber et al. 2018; Belli et al. 2019; Saracco et al. 2019; Stockmann et al. 2020). The stellar masses of lensed samples of this work and Newman et al. (2018a) are intrinsic, i.e., delensed.

<figure><img src="image_6.png" /></figure>

**Figure 5.** Star formation rate (SFR) and stellar mass ($M_*$) of our sample, split in three redshift bins. Corrections for lensing magnifications have been applied. Colored symbols denote our targets as labeled in the legend. Filled symbols show the SFR from stellar population fitting (Section 3.4). Open symbols show the SFR inferred from the [O II] emission line (Section 4.4.1), slightly shifted along the x-axis for clarity. The colored lines show the best-fitting SFR-$M_*$ relation presented in Speagle et al. (2014), calculated at the redshift of each target following the respective color scheme. The shades of color illustrate the observed scatter of 0.3 dex around these relations.

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-$\tau$ model at late times for a given $t_{50}$ and $t_{10} - t_{90}$. 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.

## 4.2. Stellar Mass-Metallicity Relation

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
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Stellar mass–metallicity relation (left) and stellar mass–age relation (right) of quiescent galaxies. Data points are colored by redshifts as indicated in the color bar. Filled symbols denote our sample, and $t_{50}$ is plotted as the y-axis on the right panel. Open symbols show literature measurements of galaxies with individual measurements: COSMOS-307881 (Lonoce et al. 2015), MRG-M0138 (Jafariyazani et al. 2020), COSMOS-11494 (Kriek et al. 2016), and MRG-M0150 (Newman et al. 2018a). We plot samples of massive quiescent galaxies (Kriek et al. 2019) and spheroidal galaxies in the cluster XLSSJ0223 (Saracco et al. 2019). Also included are measurements made on binned grism spectra from the CLEAR survey (Estrada-Carpenter et al. 2019) and stacked spectra (Onodera et al. 2015). The solid yellow and orange lines show $z \approx 0.1$ and $z \approx 0.7$ scaling relations to quiescent galaxies from Gallazzi et al. (2014), and the total scatter is shown as colored shades. The dotted line shows the best-fit relation to the quiescent galaxy sample of Ferreras et al. (2019b), at $\bar{z} = 0.9$. When error bars overlap, we apply a small shift of less than 0.03 along the x-axis to facilitate visualization.</figcaption></figure>

at $z \gtrsim 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 $\Delta\log(Z_*/Z_\odot) < 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)<fnref n="14" />, 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_*$.

0.3 dex accuracy. We take the MZR for quiescent galaxies from Gallazzi et al. (2014) for $z \approx 0.1$ (SDSS) and $z \approx 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$).

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 $\approx 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_\odot) = -0.6 \pm 0.5$; Toft et al. 2017) to mildly super-solar ($\log(Z_*/Z_\odot) = 0.16 \pm 0.13$; Newman et al. 2018a), as detailed in Appendix C. This work provides an intermediate measurement of $\log(Z_*/Z_\odot) = -0.10 \pm 0.11$. Moreover, while the works obtained using Bruzual & Charlot (2003) stellar population models assume a solar metallicity value of $Z_\odot$, more recent works using the alf code (Conroy & van Dokkum 2012; Conroy et al. 2018) assume more recent measurements of $Z_\odot = 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.

### 4.3. Magnesium Features at Rest-frame UV

Singly ionized magnesium Mg II λλ2796,2803 (hereafter Mg II) absorption in the rest-frame UV spectrum can arise from warm atomic gas ($T \sim 10^4$ K) or stellar photospheres (late B-, A-, F-, and G-type stars; Fanelli et al. 1992; Snow et al. 1994).

<footnote n="14">14 The plotted metallicity is computed from the [Mg/Fe] and [Fe/H] obtained from spectral fitting as provided in Kriek et al. (2016), i.e., $[Z/H] = [Fe/H] + 0.94 [\alpha/Fe] = (-0.25 \pm 0.11) + 0.94 (0.59 \pm 0.11) = (0.30 \pm 0.15)$.</footnote>

<footnote n="15">15 MRG-M0138 has [Mg/Fe] = 0.51 ± 0.05, which is unusually high given that the abundance of other α-elements follow closely the values of nearby early-type galaxies (Jafariyazani et al. 2020). For our purpose of estimating $[Z/H]$, we thus adopt the median [α/Fe] value of the cores of nearby massive early-type galaxies (Greene et al. 2019) rather than the measured value of [Mg/Fe], i.e., $[Z/H] = [Fe/H] + 0.94 [\alpha/Fe] = (0.26 \pm 0.04) + 0.94 \times (0.30 \pm 0.04) = (0.54 \pm 0.05)$.</footnote>

<footnote n="16">16 Table 10 provides a comparison of the properties of m2129 presented in this work, Toft et al. (2017), and Newman et al. (2018a).</footnote>

<footnote n="17">17 We multiplied the stellar mass presented in Onodera et al. (2015) by a factor of 0.67 to bring their Salpeter IMF to the Kroupa IMF assumed throughout this work, following the value used in Madau & Dickinson (2014). No correction was applied to stellar masses obtained using the Chabrier IMF, as the correction (0.04 dex) is negligible and does not affect our conclusions here.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** The rest-frame UV spectra zoomed in on Mg II λλ2796, 2803 and Mg I λ 2853 as indicated by the vertical dotted lines. The bottom x-axis refers to the velocity with respect to the Mg II λ 2803 at the systemic redshifts. The y-axis is in units of erg s⁻¹ cm⁻² Å⁻¹. The vertical red shades show the central bandpasses for Mg II and Mg I. The continua used for the EW and line flux measurements are shown as straight brown lines. The horizontal gray lines at the bottom of each panel show the continuum bandpasses used for fitting the continuum. The red continuum bandpass for Mg I λ 2853 is beyond the axis limit and is not shown in this figure. Refer to Section 4.3 for more details.</figcaption></figure>

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.

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<sub>CO,gal 2</sub> = 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 blue-shifted 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
<table><caption>**Table 5**
Spectral Line Measurements</caption><thead><tr><th></th><th>Mg II</th><th>Mg I</th><th>[O II]</th><th>[O III]λ4959</th><th>[O III]λ5007</th><th>Hα</th><th>O32</th><th>[O III]λ5007/λ4959</th><th>log(sSFR([O II])yr<sup>-1</sup>)</th></tr></thead><tbody><tr><td colspan="10" style="text-align:center;">Equivalent widths, W (Å)</td></tr><tr><td>e1341</td><td>16.4 ± 1.0</td><td>4.8 ± 0.6</td><td>-3.3 ± 0.4</td><td>0.2 ± 0.3</td><td>0.2 ± 0.3</td><td>0.4 ± 0.1</td><td>...</td><td>...</td><td>...</td></tr><tr><td>m2129</td><td>6.3 ± 4.0</td><td>6.6 ± 2.6</td><td>-1.5 ± 2.2</td><td>-0.2 ± 0.6</td><td>-3.3 ± 0.5</td><td>4.2 ± 0.8</td><td>...</td><td>...</td><td>...</td></tr><tr><td>m0451</td><td>5.9 ± 2.9</td><td>2.3 ± 3.9</td><td>-8.0 ± 1.4</td><td>-1.5 ± 0.5</td><td>-8.0 ± 0.7</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>m1423</td><td>2.6 ± 1.4</td><td>0.7 ± 1.8</td><td>-2.7 ± 1.5</td><td>-0.7 ± 0.9</td><td>-4.5 ± 0.8</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td colspan="10" style="text-align:center;">Line fluxes (×10<sup>-18</sup> erg s<sup>-1</sup> cm<sup>-2</sup>)</td></tr><tr><td>e1341</td><td>-69.1 ± 4.3</td><td>-24.4 ± 3.0</td><td>28.0 ± 3.4</td><td>-4.2 ± 5.2</td><td>-4.2 ± 5.7</td><td>6.3 ± 2.0</td><td>&lt;0.2</td><td>...</td><td>-11.0 ± 0.4</td></tr><tr><td>m2129</td><td>-17.5 ± 11.2</td><td>-20.7 ± 8.1</td><td>8.7 ± 11.9</td><td>3.3 ± 8.4</td><td>41.6 ± 5.8</td><td>38.5 ± 7.2</td><td>&gt;3.5</td><td>&gt;5.0</td><td>&lt;-11.3</td></tr><tr><td>m0451</td><td>-20.2 ± 9.9</td><td>-8.1 ± 14.2</td><td>37.1 ± 6.3</td><td>13.2 ± 4.4</td><td>67.6 ± 5.7</td><td>...</td><td>1.8 ± 0.3</td><td>5.1 ± 1.8</td><td>-10.2 ± 0.4</td></tr><tr><td>m1423</td><td>-6.8 ± 3.7</td><td>-1.7 ± 4.7</td><td>6.1 ± 3.3</td><td>1.8 ± 2.6</td><td>11.8 ± 2.2</td><td>...</td><td>1.9 ± 1.1</td><td>&gt;4.5</td><td>-10.3 ± 0.5</td></tr></tbody></table>

**Notes.** Rest-frame equivalent widths and observed frame line fluxes are provided. Emission lines have negative equivalent widths and positive line fluxes, and vice versa for absorption lines. No correction for lensing magnification or extinction has been applied. Note that some of the Mg II and Mg I measurements are likely affected by emission line infilling (Section 4.3). Here, [O II] refers to the blended [O II] λλ3726,3729 doublet. No correction has been applied to deblend Hα from the [N II] doublet that brackets it. O32 refers to the [O III] λ 5007/[O II] λλ3726,3729 flux ratio. sSFR([O II]) is the specific SFR inferred from the [O II] emission. Upper limits are 1σ.

<figure><img src="image_9.png" /><figcaption>**Figure 8.** The binned spectra are zoomed in on the [O II] and [O III] emission lines (middle and bottom panels) in velocity space, together with the Mg II absorption (top panel). The Bagpipes model continua are shown as gray lines.</figcaption></figure>

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.

### 4.4. Ionized Gas Emission

Warm ($T \sim 10^4$ 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).<fnref n="18" /> 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).

#### 4.4.1. SFR Inferred from [O II]

The [O II] luminosity is sometimes used as an SFR indicator in high-redshift galaxies via an empirical calibration to Hα luminosity. Its excitation is sensitive to the oxygen abundance and the ionization state of the gas (Kennicutt 1998; Kennicutt & Evans 2012, and references therein). The use of [O II] as an SFR indicator in post-starburst galaxies has been brought to question, as [O II] can arise from AGN, shocks, post-asymptotic giant branch stars, and cooling flows (Yan et al. 2006). We use the [O II] luminosity as a consistency check for the SFR provided by the stellar population fitting. We apply the

<footnote n="18">¹⁸ Spatially resolved observations show that emission lines in massive, quiescent galaxies are not necessarily confined to the nuclear region (Singh et al. 2013; Belfiore et al. 2016). Therefore, we adopt the acronym LIER, rather than LINER, without reference to the nuclear region as used in earlier studies.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Source-plane morphologies of our sample reconstructed from the lens models. The reddest HST filters available are used. For each image, the first column shows the reconstructed image on the source plane. The second column shows the best-fitting Sérsic model as determined with GALFIT. The third column shows the residual of the image subtracted by the model. The last image shows the effective point-spread function, i.e., what point sources (stars) on the image plane would look like on the source plane. The first three columns of each multiple image are always shown with the same color scale, although the color scale is not identical across the images.</figcaption></figure>

Kennicutt (1998, Equation (3)) conversion adjusted to the Kroupa (2001) IMF to infer the SFR from the [O II] line fluxes and provide them in Table 5. The inferred log(sSFR) range from −10.2 to < −11.3, on par with the low values obtained from the full spectral fitting (Table 3). Reassuringly, the SFRs inferred from [O II] are fully consistent with being quiescent stellar populations.

As H II regions may not be the sole supplier of ionizing photons, the SFRs inferred from [O II] are upper limits. Although the inferred SFRs are not corrected for extinction, we expect the correction to be minor in these galaxies if the nebular extinction is similar to that of the stellar continuum.

The most discrepant measurement is for m1423 (Figure 5). Source m1423 is most likely caught in a rapidly quenching phase, as emission line SFR indicators are sensitive to more recent SF than the UV continuum (Kennicutt & Evans 2012). This is in line with the transitional nature suggested by its *U*−*V* and *V*−*J* colors, as discussed in Section 5.4 and shown in Figure 14. Another possibility is that the [O II] emission is heavily dust-obscured in m1423. Its stellar continuum has $A_V = 0.8 \pm 0.1$, implying an attenuation of 1.2 mag at the wavelength of [O II], which is too small to fully account for the discrepant SFR estimate from spectral fitting. Unless the nebular emission is substantially more dust-obscured than the stellar continuum, this is unlikely to be the cause for discrepancy.

### 4.4.2. [O III]-to-[O II] Flux Ratio

The [O III] $\lambda$5007/[O II] $\lambda\lambda$3726,3729 flux ratio, hereafter O32, is commonly used as an ionization parameter diagnostic. Type 1 quasars have a distribution of O32 that peaks at $\approx 5$, while type 2 quasars peak at O32 $\approx 2$ (Zakamska et al. 2003, Figure 7). LIERs are expected to have O32 near unity (Netzer 1990). Nonactive galaxies with stellar populations aged $\ge 1$ Gyr have O32 $\lesssim 0.6$ (Johansson et al. 2016). The variation of O32 reflects the different ionization parameters as well as the oxygen abundance and ISM pressure across galaxies (Kewley & Dopita 2002; Kewley et al. 2019).

Our sample spans a range of O32 as shown in Table 5. O32 is higher than unity in the youngest three targets, m1423, m0451, and m2129, in line with the expectation that they host type 2 AGN as discussed in Section 4.4.3. Their O32 and [O III] are comparable to type 2 AGN at $z \approx 0$ (Kauffmann et al. 2003; Silverman et al. 2009). O32 is poorly constrained for the oldest two targets, e1341 and m2129, but their low values along with low [O III] luminosities suggest that they are unlikely to host type 1 AGN.

Discerning the source of ionizing photons provides additional clues to quenching mechanisms. If a past, luminous AGN episode was responsible for star formation quenching, one might expect to see signatures in the structure and ionization state of the gas (King et al. 2011; Zubovas & King 2014). We discuss the implications of these findings in Section 5.
<figure><img src="image_11.png" /></figure>

**Figure 10.** The Sérsic indices and axis ratios of lensed quiescent galaxies from this work and Newman et al. (2018a) are individually labeled. They are contrasted against unlensed quiescent galaxies from the LEGA-C survey (Bezanson et al. 2018, *z* = 0.6–1.0), from the MOSFIRE survey of the CANDELS field (Belli et al. 2017, *z* = 1.5–2.5), and from the X-SHOOTER survey of the COSMOS field (Stockmann et al. 2020, *z* = 2.0–2.7). The distribution of Sérsic indices and axis ratios of these three samples are plotted as histograms. The stellar mass dependence of the Sérsic indices and axis ratios are shown in Figure 20.

<figure><img src="image_12.png" /></figure>

**Figure 11.** The effective radii of our sample, normalized by their sizes as $r_{\text{eff.}}$, $n = r_{\text{eff.}}/(M_{*}/(5 \times 10^{10} M_{\odot})^{n})$, are denoted as filled colored symbols. We adopt $\alpha$ values for early-type galaxies at the nearest redshift bins (van der Wel et al. 2014, Table 1). The red and blue solid curves represent the size evolution of early- and late-type galaxies with $M_{*} = 5 \times 10^{10} M_{\odot}$, following the parameterized redshift evolution presented in van der Wel et al. (2014, Figure 6).

To unambiguously discern the source of ionizing photons, it is imperative to access more emission lines in order to place them on diagnostic plots. In particular, the [O I]-to-H$\alpha$ ratio is an important addition to characterize the hardness of the radiation field, to distinguish between H II regions, Seyferts and LIERs (Kewley et al. 2006; Johansson et al. 2016, Figure 5). Observations of the [S II] and [N II] doublets will help to quantify the importance of shocks (Yan & Blanton 2012).

<figure><img src="image_13.png" /></figure>

**Figure 12.** The colored lines represent the evolution of the halo mass of the lensed quenched galaxies as discussed in Section 5.1. The colored shades represent the $1\sigma$ uncertainties propagated from the number densities and their evolution. As for m1423, labeled as a dotted line, the shades only indicate the uncertainties for the number density evolution, as the lower number density error bounds reach below $10^{-5}$ Mpc$^{-3}$ such that the nd-redshift does not return reliable results. The boxplot shows the MASSIVE sample at $z \approx 0$ (Veale et al. 2018): the central line represents the median, the box represents the interquartile range, and the whiskers represent the full range. The boxplot has been offset in redshift slightly to improve visualization. Several well-known nearby galaxy clusters are individually labeled for comparison (Crook et al. 2007; Li & White 2008).

<figure><img src="image_14.png" /></figure>

**Figure 13.** Stellar metallicities of our sample, compared with the metallicity gradient of nearby galaxies. The green curves with crosses are median metallicity gradients of early-type galaxies in the CALIFA survey (Zibetti et al. 2020). The blue curves are fits of massive early-type galaxies in the MASSIVE survey (Greene et al. 2015). The dotted portions of the lines are extrapolations of the best fits to small radii. Both samples are shown in bins of stellar velocity dispersions. The gray curve represents the best fit to the metallicity gradient of NGC 4889, one of the two bright central galaxies (BCG) of the Coma cluster (Coccato et al. 2010). The orange shades represent the metallicity gradients of local massive, compact, quiescent galaxies (or “relic galaxies”; Ferré-Mateu et al. 2017).

Spatially resolved emission line ratios will help distinguish between different ionizing sources, as commonly done in nearby galaxies (Singh et al. 2013; Belfiore et al. 2016).

### 4.4.3. [O III] as a Tracer for Nuclear Activity

[O III] $\lambda$ 5007 luminosity, $L_{[O III]}$, is a proxy of the AGN accretion rate (Heckman & Best 2014, and references therein). The line is detected in our sample except for the oldest target e1341. The [O III] FWHM line width is <1000 km s$^{-1}$, ruling out the presence of Seyfert 1 AGN in the sample (Padovani et al. 2017, and references therein). While [Ne V] $\lambda$ 3426 is a reliable tracer for AGN activity (e.g., Feltre et al. 2016), it is
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Rest-frame U–V and V–J magnitudes of our sample are denoted as filled circles, color-coded by their $t_{50}$. The contours show the 1$\sigma$, 2$\sigma$, and 3$\sigma$ distributions of the colors propagated from the stellar population fitting. The dotted line shows the color criteria for separating quiescent (upper left corner) from star-forming galaxies as in Muzzin et al. (2013), while the solid black line shows the modified criteria presented in Belli et al. (2019). The colored curves show the evolution of fast and slow quenching model tracks presented in Belli et al. (2019), with $\tau = 0.1$ Gyr (blue) and $\tau = 1$ Gyr (green), respectively. The crosses mark the time intervals at 0.5, 1, 3, and 5 Gyr.</figcaption></figure>

typically an order of magnitude fainter than [O III] $\lambda$ 5007 in AGN (Zakamska et al. 2003), and is below our detection limit even if present. Only m0451 has both of the [O III] $\lambda\lambda$4959, 5007 emission lines detected, and the ratio is within 1.5$\sigma$ of the theoretical value of 2.98 (Storey & Zeippen 2000).

The intrinsic (delensed) $L_{[O III]}$ are shown in Table 6. These $L_{[O III]}$ are more than an order of magnitude below those of powerful radio galaxies at $z \approx 2$ during the quasar feedback phase (Nesvadba et al. 2017). In the absence of an alternative tracer of the AGN bolometric luminosity ($L_{AGN}$) such as the IR SED or X-ray observations, we estimate the $L_{AGN}$ from $L_{[O III]}$ (without dust extinction correction) using a mean bolometric correction of 3500 assuming the same factor applies to type 1 and type 2 AGN (Heckman et al. 2004). Attributing the entirety of the [O III] emission to AGN, this implies that $L_{AGN} = 1.1–1.6 \times 10^{45}$ erg s$^{-1}$ for the three targets with detected [O III] on par with or below those at the faint end of AGN samples at $z \approx 2$ (Circosta et al. 2018; Leung et al. 2019). These $L_{AGN}$ imply a mass outflow rate of $<10 M_{\odot}$ yr$^{-1}$ following the best-fit relation between the $L_{AGN}$ and mass outflow rate (Leung et al. 2019, Figure 9). The mass outflow rates are likely upper limits, as this derivation assumes that all the [O III] emission is photoionized by the AGN. The outflowing ionized mass are thus insignificant and unlikely to escape the deep gravitational potentials of these compact, massive galaxies.

The Eddington ratio, $L_{AGN}/L_{Edd}$, is a useful quantity to
examine the accretion mode of supermassive black holes

(SMBHs). Radiative-mode SMBHs accrete at 1%–10% of $L_{Edd}$, whereas jet-mode SMBHs accrete much less efficiently, at less than 1% of $L_{Edd}$ (Best & Heckman 2012). The luminosity of the classical Eddington limit for each target is derived as $L_{Edd} = 3.3 \times 10^4 M_{BH}$ (Heckman & Best 2014, Equation (4)). Black hole masses are inferred from the stellar velocity dispersions (Table 3) using the local $M_{BH}-\sigma_{\star}$ relation presented in McConnell & Ma (2013): $\log(M_{BH}) = 8.32 + 5.64 \log(\sigma_{\star}/200 \text{ km s}^{-1})$. The resulting Eddington ratios are listed in Table 6. Based on this calculation, the older two targets, m2129 and e1341, host jet-mode SMBH. Source m0451 hosts a radiative-mode SMBH, while m1423 straddles the threshold. Unavoidably, there are systematic uncertainties involved in the calculation of $L_{Edd}$, such as the fraction of $L_{[O III]}$ photoionized by the AGN, the AGN bolometric correction, the contribution of rotation in the measured stellar velocity dispersions, and the scatter in the $M_{BH}-\sigma_{\star}$ relation. Accurately constraining these quantities is beyond the scope of this paper. If our sample evolves along the $M_{BH}-\sigma_{\star}$ relation, their supermassive black holes must have already grown considerably. The weak [O III] emission suggests that our sample are well past their peak AGN episode and do not harbor copious amount of outflowing, line-emitting ionized gas. Altogether, our results suggest a transition of AGN accretion from radiative-mode to jet-mode throughout the star formation quenching process. We discuss the implications of these findings in Section 5.

4.5. Morphological Analysis

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 proce-
dures 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 recon-
structed image, and obtain the best-fitting morphological
parameters to each realization using the GALFIT code<fnref n="19" /> (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 reconstruc-
tion. Reassuringly, the Sérsic index *n* and effective radii $r_{\text{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 \approx 0.4$ for the most magnified image (im1), compared to
$q \approx 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

<footnote n="19">¹⁹ https://usersobs.carnegiescience.edu/peng/work/galfit/galfit.html</footnote>
morphological parameters. The $r_{\text{eff}}$ and $n$ are more robust to lens modeling uncertainties than $q$.

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 \approx 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\sigma$ 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_{\text{eff}} \approx 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 $\approx 4.3 \times$ the scatter of the relation. This places e1341 above the late-type galaxy mass–size relation.

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, $\Delta \log R_{\text{eff}} / \Delta \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 $\lesssim 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.

## 5. Discussion

In this section, we will discuss what our results imply for the overall evolution sequence of quenched galaxies. How did they come to be, and what would they evolve into? What is responsible for quenching star formation?

### 5.1. Tracking Evolution by Number Density

To figure out what our sample of lensed quenched galaxies would evolve into by $z \approx 0$, we estimate their halo mass ($M_{\text{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<fnref n="20" /> (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 \approx 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).

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_{\text{halo}}) \approx 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

<footnote n="20">²⁰ https://code.google.com/archive/p/nd-redshift/</footnote>
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.

## 5.2. Progenitors

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, parameteriza-
tion 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-consum-
ing, 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_{\odot}$ yr$^{-1}$. The most
massive ($\log(M_{\star}/M_{\odot}) > 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_{\odot}$ yr$^{-1}$ (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 an example of
 accelerated growth in dense environments. Compact star-
 forming galaxies at $z \gtrsim 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_{\odot}$ yr$^{-1}$). 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 \approx 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_{10}-t_{90}$, 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.

## 5.3. Descendants

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 \approx 1-2$
(Thomas et al. 2005; McDermid et al. 2015).<fnref n="21" /> 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
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

<footnote n="21">²¹ Although the stars in today’s massive galaxies formed early, galaxies could have assembled later through mergers long after stars formed.</footnote>
results from the MASSIVE and CALIFA surveys that provide
gradient measurements to the largest spatial extent, out to ≈2.5
r<sub>eff</sub>. 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.

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 rejuvena-
tion (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).

5.4. Implications for Quenching Mechanisms

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 examin-
ing 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.<fnref n="22" /> 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<sub>V</sub> that
declines with the SFR over time, starting from 2 to 0.4, and
assuming R<sub>V</sub> = 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<sub>10</sub> − t<sub>90</sub> = 0.07 − 0.47 Gyr. This conclusion
stands even for higher initial values of A<sub>V</sub> or adopting R<sub>V</sub>
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;

<footnote n="22">²² The conclusions of this comparison remain unchanged if we measure the
colors from the observed photometry instead. Template SEDs are commonly
used to interpolate between observed filters to obtain rest-frame colors, in any
case; see, e.g., Taylor et al. (2009).</footnote>
<table><caption>**Table 6**
Inferred Supermassive Black Hole Properties</caption><thead><tr><th></th><th>log(L<sub>[O III]</sub>)</th><th>log(L<sub>AGN</sub>)</th><th>log(M<sub>BH</sub>)</th><th>L<sub>AGN</sub>/L<sub>Edd</sub> (%)</th></tr></thead><tbody><tr><td>m1341</td><td>&lt;5.9</td><td>&lt;9.5</td><td>7.8 ± 0.1</td><td>&lt;0.1</td></tr><tr><td>m2129</td><td>7.9 ± 0.1</td><td>11.5 ± 0.1</td><td>9.5 ± 0.2</td><td>0.3 ± 0.2</td></tr><tr><td>m0451</td><td>8.1 ± 0.1</td><td>11.6 ± 0.1</td><td>7.9 ± 0.2</td><td>14.4 ± 6.2</td></tr><tr><td>m1423</td><td>8.0 ± 0.1</td><td>11.6 ± 0.1</td><td>9.1 ± 0.2</td><td>0.8 ± 0.3</td></tr></tbody></table>

**Notes.** The first three columns are quoted in logarithmic units of solar values. The uncertainties are propagated from measurement errors only. Upper limits are $1\sigma$ for the [O III] nondetection of e1341.

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 \approx 2$ quiescent galaxies (Kriek et al. 2016; Newman et al. 2018a; Zick et al. 2018; Belli et al. 2019).

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_\odot) = 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 \approx 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 \approx 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 \gtrsim 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 \sim 10^4$ 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 ($\lesssim 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$\alpha$ emission is present in m0451 (Jauzac et al. 2021, $z_{Ly\alpha} = 2.93$) across the three multiple images as detected with the VLT/MUSE spectrograph. The extended Ly$\alpha$ 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 $\sim 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 $\gtrsim 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 \gtrsim 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
<table><caption>**Table 7**
Morphological Constraints of Reconstructed Images</caption><thead><tr><th></th><th>Filter</th><th>λ<sub>rest</sub><br>(Å)</th><th>Angular Scale<br>(kpc/")</th><th>Sérsic Index (n)</th><th>Axis Ratio (b/a)</th><th>Position Angle<br>(°)</th><th>r<sub>eff</sub><br>(kpc)</th><th>r<sub>circ</sub><br>(kpc)</th></tr></thead><tbody><tr><td>e1341 (im2)</td><td>F140W</td><td>5366.1</td><td>8.471</td><td>1.68 ± 0.02</td><td>0.45 ± 0.09</td><td>35.0 ± 0.7</td><td>2.82 ± 0.07</td><td>1.89 ± 0.20</td></tr><tr><td>e1341 (im3)</td><td>—"—</td><td>—"—</td><td>—"—</td><td>1.47 ± 0.01</td><td>0.40 ± 0.06</td><td>14.8 ± 5.6</td><td>2.64 ± 0.11</td><td>1.67 ± 0.15</td></tr><tr><td>m2129</td><td>F160W</td><td>4881.7</td><td>8.294</td><td>1.17 ± 0.01</td><td>0.42 ± 0.01</td><td>-35.2 ± 0.5</td><td>2.12 ± 0.06</td><td>1.37 ± 0.04</td></tr><tr><td>m0451 (im1)</td><td>F160W</td><td>3919.7</td><td>7.762</td><td>1.09 ± 0.01</td><td>0.44 ± 0.03</td><td>-16.2 ± 2.3</td><td>0.99 ± 0.04</td><td>0.66 ± 0.03</td></tr><tr><td>m0451 (im2)</td><td>—"—</td><td>—"—</td><td>—"—</td><td>1.28 ± 0.90</td><td>0.68 ± 0.06</td><td>58.8 ± 22.8</td><td>1.05 ± 0.18</td><td>0.87 ± 0.15</td></tr><tr><td>m0451 (im3)</td><td>—"—</td><td>—"—</td><td>—"—</td><td>0.95 ± 0.03</td><td>0.63 ± 0.09</td><td>-37.3 ± 5.6</td><td>0.79 ± 0.04</td><td>0.63 ± 0.05</td></tr><tr><td>m1423</td><td>F160W</td><td>3651.9</td><td>7.542</td><td>1.22 ± 0.01</td><td>0.56 ± 0.03</td><td>-26.3 ± 1.6</td><td>1.09 ± 0.03</td><td>0.81 ± 0.03</td></tr></tbody></table>

Notes. Half-light effective radii ($r_{eff}$) and circularized radii ($r_{circ}$) are provided in kpc. The quoted values and errors refer to the mean and standard deviation of 100 Monte Carlo realizations perturbed by lensing model errors. Reconstruction of the im1 of e1341 is not possible, because the lensing arc is only a partial view of the source.

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 \gtrsim 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.

## 6. Conclusions

This paper presents the analysis of VLT/X-SHOOTER spectroscopic and HST imaging observations of four quenched galaxies at $z = 1.6-3.2$ that are gravitationally lensed by foreground clusters. Their magnification factors range from

$\approx 3$ to 30, affording us an exquisitely deep view of their stellar populations and morphologies. The photometry and spectra have been fitted with the Bagpipes stellar population synthesis fitting code. Our main findings are as follows:

1. The four quenched galaxies have intermediate stellar masses ($\log(M_*/M_\odot) = 10.2-11.2$), or $0.1-3 \times$ the characteristic stellar mass at their respective redshifts. The median ages of the stellar populations range between $t_{50} \approx 0.12$ and 1.2 Gyr, and they formed 80% of their stellar masses within 0.07-0.5 Gyr. Their specific SFR span a range of $\log(sSFR/yr^{-1})$ from -8.4 to -11.2. Three galaxies lie below the sequence of star-forming galaxies, where the youngest target is a rare example of a galaxy caught shortly after quenching, as its [O II] emission implies a SFR $\approx 2$ dex below that estimated from stellar populations fitting.

2. The targets of the sample have stellar metallicities of $\approx 0.5-1$ solar value, as inferred from Bagpipes fitting. However, systematic uncertainties inherent to stellar population synthesis modeling and metallicity calibration preclude a fair comparison with literature metallicity measurements across redshifts in a consistent manner at this stage. Further work is necessary to quantify such systematic effects and to self-consistently model the chemical evolution of galaxies.

3. All targets in our sample show Mg II $\lambda\lambda 2796,2804$ absorption. While stellar photospheres may partially account for the absorption, additional contribution from warm gas is required. Galactic outflows are present in the most recently quenched galaxies, as evidenced by their blueshifted absorption and/or redshifted emission.

4. Faint emission lines are detected in some of our targets. The [O II] $\lambda\lambda 3726,3729$ Å fluxes of our sample are consistent with having quiescent SF. The [O III] $\lambda 5007$ Å fluxes are at or below the faint end of $z \approx 2$ AGN surveys. We rule out the presence of type 1 AGN in our sample. The youngest two targets may host type 2 AGN, as inferred from their [O III]/[O II] ratio, and might be responsible for driving the Mg II outflow. If we attribute the [O III] emission to AGN only, this implies that the SMBH of the second youngest target ($t_{50} \approx 0.3$ Gyr) is accreting in radiative mode at a few percent of the Eddington limit, where the older two targets ($t_{50} \approx 0.6-1.2$ Gyr) accrete in jet-mode at less than 1% of the Eddington limit.

5. We use the lens models to reconstruct the rest-frame optical light profiles of our sample on their source planes.
<table><caption>**Table 8**
Photometry of the Targets</caption><thead><tr><th></th><th></th><th>e1341 (im1)</th><th>m2129</th><th>m0451 (im1)</th><th>m1423</th></tr></thead><tbody><tr><td rowspan="10">HST/ACS-WFC</td><td>F435W</td><td>...</td><td>25.809 ± 0.366</td><td>...</td><td>−99</td></tr><tr><td>F475W</td><td>...</td><td>25.494 ± 0.218</td><td>...</td><td>25.169 ± 0.086</td></tr><tr><td>F555W</td><td>...</td><td>25.313 ± 0.158</td><td>24.810 ± 0.069</td><td>24.426 ± 0.036</td></tr><tr><td>F606W</td><td>22.645 ± 0.037</td><td>24.836 ± 0.272</td><td>...</td><td>23.790 ± 0.022</td></tr><tr><td>F625W</td><td>...</td><td>24.567 ± 0.130</td><td>...</td><td>...</td></tr><tr><td>F775W</td><td>...</td><td>23.776 ± 1.095</td><td>23.289 ± 0.188</td><td>23.331 ± 0.038</td></tr><tr><td>F814W</td><td>20.791 ± 0.013</td><td>23.573 ± 0.028</td><td>23.125 ± 0.013</td><td>23.198 ± 0.013</td></tr><tr><td>F850LP</td><td>...</td><td>22.910 ± 0.046</td><td>23.060 ± 0.080</td><td>23.007 ± 0.023</td></tr><tr><td>F105W</td><td>19.300 ± 0.009</td><td>22.122 ± 0.017</td><td>...</td><td>22.875 ± 0.014</td></tr><tr><td>F110W</td><td>19.000 ± 0.009</td><td>21.339 ± 0.009</td><td>22.357 ± 0.014</td><td>22.751 ± 0.015</td></tr><tr><td rowspan="4">HST/WFC3-IR</td><td>F125W</td><td>...</td><td>20.944 ± 0.009</td><td>...</td><td>22.623 ± 0.013</td></tr><tr><td>F140W</td><td>18.493 ± 0.003</td><td>20.442 ± 0.005</td><td>...</td><td>22.465 ± 0.012</td></tr><tr><td>F160W</td><td>...</td><td>20.163 ± 0.003</td><td>20.893 ± 0.007</td><td>22.141 ± 0.006</td></tr><tr><td>CH1</td><td>...</td><td>19.062 ± 0.133</td><td>19.937 ± 0.133</td><td>20.994 ± 0.133</td></tr><tr><td>CH2</td><td>...</td><td>18.886 ± 0.133</td><td>19.716 ± 0.133</td><td>20.730 ± 0.133</td></tr><tr><td rowspan="2">IRAC</td><td>W1</td><td>17.652 ± 0.061</td><td>...</td><td>...</td><td>...</td></tr><tr><td>W2</td><td>18.414 ± 0.319</td><td>...</td><td>...</td><td>...</td></tr><tr><td>WISE</td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

**Notes.** Magnitudes are provided in the AB system. Correction for Galactic extinction has not been applied. Magnitude −99 indicates negative flux.

Their light profiles are best-fitted with low Sérsic indices of $n < 2$ and intermediate axis ratios of $q \approx 0.4–0.7$. The targets lie within the mass–size relation for early-type galaxies, except for e1341, which is as extended as late-type galaxies at its mass and redshift. Our sample is uniformly disky ($n < 2$), suggesting the need for additional morphological transformation if they are to evolve into metal-rich bulges.

6. Altogether, our results imply that star formation quenching at high redshift is a rapid process (<1 Gyr). Galactic-scale outflows of warm gas are detected in the most recently quenched targets, perhaps driven by faint, radiatively accreting type 2 AGN or a previous starburst. Contrary to some claims, quenching neither requires nor synchronizes with the formation of prominent stellar bulges.

This sample forms the backbone of the REsolving QUIEscent Magnified Galaxies Survey (REQUIEM-2D). HST grism observations will enable measurements of spatially resolved stellar populations, in order to quantify age and SFR gradients (Akhshik et al. 2020, 2021). ALMA surveys of the CO and dust continuum will provide molecular gas mass measurements (A. Man et al. 2021, in preparation; Whitaker et al. 2021). With the launch of the JWST, we anticipate huge leaps in the characterization of distant quiescent galaxies. NIR imaging surveys with the JWST/NIRCam will enable the identification of $z > 3$ quiescent galaxies. Their number density evolution is vital for understanding when the first galaxies become quenched. JWST/NIRCam can better resolve their light profiles, as they are currently barely resolved with HST. An unbiased census for the light distribution of $z \gtrsim 2$ quiescent galaxies will enable us to confirm whether morphological transformation and bulge formation are asynchronous with star formation quenching, as our work suggests. JWST/NIRSpec will provide resolved emission line maps of magnified quiescent galaxies, allowing us to identify the origin of the ionizing photons. All of these will provide valuable insights into the origins of star formation quenching.

<table><caption>**Table 9**
Rest-frame Colors and Covariance</caption><thead><tr><th></th><th>e1341</th><th>m2129</th><th>m0451</th><th>m1423</th></tr></thead><tbody><tr><td>(U-V)</td><td>1.70 ± 0.05</td><td>1.76 ± 0.04</td><td>1.31 ± 0.04</td><td>0.85 ± 0.05</td></tr><tr><td>(V - J)</td><td>0.95 ± 0.09</td><td>1.04 ± 0.06</td><td>0.63 ± 0.07</td><td>0.59 ± 0.11</td></tr><tr><td>cov(U - V, V - J)</td><td>0.0043</td><td>0.0020</td><td>0.0023</td><td>0.0029</td></tr></tbody></table>

We thank the referee, Andrew Newman, for timely and constructive reports that improved the quality of this manuscript. The authors are grateful to the ESO astronomers who assisted in planning and conducting the observations. We appreciate the CLASH and SURFSUP teams for making their images publicly available. A.M. is grateful to Andrew Zirm for assistance in the early stage of this work. A.M. thanks Sirio Belli and Anna Ferré-Mateu for providing data from their published figures for comparison. A.M. acknowledges helpful discussions with Alison Coil, Roberto Maiolino, and Tyrone Woods during the preparation of this work. We thank Adam Carnall for making Bagpipes publicly available and for clarifying the definitions of the solar abundance.

A.M. was supported by a Dunlap Fellowship at the Dunlap Institute for Astronomy & Astrophysics, funded through an endowment established by the David Dunlap family and the University of Toronto. The University of Toronto operates on the traditional land of the Huron-Wendat, the Seneca, and most recently, the Mississaugas of the Credit River; A.M. is grateful to have the opportunity to work on this land. J.Z. acknowledges support from ANR grant ANR-17-CE31-0017 (3DGasFlows). J.R. acknowledges support from the ERC Starting Grant 336736-CALENDS. G.B., S.T., and M.S. acknowledge support from the ERC Consolidator Grant funding scheme (project ConTExt, grant No. 648179). The Cosmic Dawn Center is funded by the Danish National Research Foundation under grant No. 140.

Based on observations collected at the European Southern Observatory under ESO programs 087.B-0812, 093.B-0815, 096.B-0994, 097.B-1064, and 099.B-0912. Based on observations made with the NASA/ESA Hubble Space Telescope,
obtained from the data archive at the Space Telescope Science Institute. STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555. This work is based in part on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, and NEOWISE, which is a project of the Jet Propulsion Laboratory/California Institute of Technology. WISE and NEOWISE are funded by the National Aeronautics and Space Administration.

*Facilities*: VLT:Kueyen(X-SHOOTER), HST, Spitzer, WISE.

*Software*: APLpy (Robitaille & Bressert 2012), AstroPy (Astropy Collaboration et al. 2013; Price-Whelan et al. 2018), Bagpipes (Carnall et al. 2018), Cloudy (Ferland et al. 2017), corner (Foreman-Mackey 2016), GALFIT (Peng et al. 2002, 2010, 2011), grizli (Brammer 2019), Lenstool (Jullo et al. 2007), Matplotlib

(Hunter 2007), Molecfit (Kausch et al. 2015; Smette et al. 2015), MultiNest (Feroz et al. 2019), ND-Redshift (Behroozi et al. 2013), Numpy (Harris et al. 2020), PyMultiNest (Buchner et al. 2014), SExtractor (Bertin & Arnouts 1996), TPHOT (Merlin et al. 2015, 2016).

# Appendix A: Photometry and Colors

The photometry of the targets obtained with the procedures described in Section 3.1 is listed in Table 8. The rest-frame U–V and V–J colors used in Section 5.4 are listed in Table 9.

# Appendix B: Stellar Population Fits

The covariance of the stellar population fitting parameters is shown in Figure 15. The model spectra and photometry are shown against the observed ones for each target in Figures 16–19.
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Covariance plot of the stellar population fitting. More details can be found in Section 3.4.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Photometry and spectrum of e1341. Photometric data are shown as large black circles. The binned VLT/X-SHOOTER spectrum is shown as small black squares. The best-fit stellar population model spectra and photometry are shown as colored lines and circles. The bottom panel shows the residual spectrum.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 17.** Photometry and spectrum of m2129. Refer to Figure 16 for explanations of the symbols.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure 18.** Photometry and spectrum of m0451. Refer to Figure 16 for explanations of the symbols.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 19.** Photometry and spectrum of m1423. Refer to Figure 16 for explanations of the symbols.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 20.** The stellar mass dependence of Sérsic indices and axis ratios. The symbols follow those of Figure 10.</figcaption></figure>
<table><caption>**Table 10**
Comparison of Derived Parameters of m2129</caption><thead><tr><th></th><th>This Work</th><th>Toft et al. (2017)</th><th>Newman et al. (2018a)</th></tr></thead><tbody><tr><td>z<sub>spec</sub></td><td>2.1487 ± 0.0002</td><td>2.1478 ± 0.0006</td><td>2.1487</td></tr><tr><td>μ</td><td>4.6 ± 0.2</td><td>4.6 ± 0.2</td><td>4.5</td></tr><tr><td>log(M<sub>*</sub>/M<sub>☉</sub>)</td><td>11.20 ± 0.05</td><td>11.10<sup>+0.26</sup><sub>-0.20</sub></td><td>10.96 ± 0.10</td></tr><tr><td>log(Z<sub>☉</sub>/Z<sub>☉</sub>)</td><td>-0.10 ± 0.11</td><td>-0.6 ± 0.5</td><td>0.16 ± 0.13</td></tr><tr><td>t<sub>age</sub> (Gyr)</td><td>0.61 ± 0.05 (t<sub>50</sub>)</td><td>1.05<sup>+0.81</sup><sub>-0.46</sub> (light-weighted)</td><td>0.80 ± 0.10</td></tr><tr><td>SFR (M<sub>☉</sub> yr<sup>-1</sup>)</td><td>0.03<sup>+0.61</sup><sub>-0.03</sub></td><td>0.0<sup>+0.2</sup><sub>-0.0</sub></td><td>0.6 ± 0.4</td></tr><tr><td>r<sub>eff</sub>''</td><td>0.26 ± 0.01</td><td>†0.27<sup>+0.05</sup><sub>-0.04</sub></td><td>0.29 ± 0.02</td></tr><tr><td>n</td><td>1.17 ± 0.01</td><td>†1.01<sup>+0.12</sup><sub>-0.06</sub></td><td>1</td></tr><tr><td>q</td><td>0.42 ± 0.01</td><td>†0.59<sup>+0.03</sup><sub>-0.09</sub></td><td>0.29 ± 0.03</td></tr></tbody></table>

**Note.** This work assumes the Kroupa (2001) IMF, while the Chabrier (2003) IMF was assumed in both Toft et al. (2017) and Newman et al. (2018a). Library B, Fit 1 (delayed-tau model and fitted over entire spectrum and photometry) is taken for Toft et al. (2017) values. Single Sérsic model value of Newman et al. (2018a) is shown here. The dagger symbols mark published values affected by an error as described in the text.

## Appendix C
Comparison of Parameters Derived for m2129

Table 10 shows a comparison of the stellar population and morphological parameters of m2129 derived in this work with published values from Toft et al. (2017) and Newman et al. (2018a). The stellar population parameters are in general agreement, except for the age and stellar metallicity. The values presented in this work lie between those in Toft et al. (2017) and Newman et al. (2018a). The substantial uncertainty in stellar metallicity across all three works implies that it is poorly constrained. The measurement presented in Toft et al. (2017) was derived using a different stellar population fitting procedure: the light-weighted stellar metallicity in the library varies with the star formation history as a function of stellar mass. The error bars in Toft et al. (2017) are larger because a large range of model parameters were allowed: for each parameter, the quoted uncertainties are given by the width of the probability distribution function marginalized over all the other parameters.

All three works uniformly point to a compact ($r_{eff} \approx 0''/3$) disk well-fitted by an exponential disk ($n \approx 1$). The GALFIT parameters reported in Toft et al. (2017) have been revised due to a mistake in the GALFIT input file: a PSF fine-sampling factor of 2 was erroneously used rather than 0.5. This work reports the rectified GALFIT results. The revised $r_{eff}$ and $n$ are similar to the values previously reported in Toft et al. (2017), and $q \approx 0.4$ is lower and thus closer to the $q \approx 0.3$ reported in Newman et al. (2018a). The axis ratio $q$ is not as well-constrained, and it shows a larger variation across the multiple images, available for two other targets, as discussed in Section 4.5. Despite the differences in modeling methods and priors, the resulting stellar population and morphological parameters are in good qualitative agreement. These differences do not affect the conclusions drawn in this work.

### ORCID iDs

Allison W. S. Man https://orcid.org/0000-0003-2475-124X
Johannes Zabl https://orcid.org/0000-0002-9842-6354

Gabriel B. Brammer https://orcid.org/0000-0003-2680-005X

Johan Richard https://orcid.org/0000-0001-5492-1049

Sune Toft https://orcid.org/0000-0003-3631-7176

Mikkel Stockmann https://orcid.org/0000-0001-5983-6273

Anna R. Gallazzi https://orcid.org/0000-0002-9656-1800

Stefano Zibetti https://orcid.org/0000-0003-1734-8356

Harald Ebeling https://orcid.org/0000-0001-8249-2739

## References

Abramson, L. E., Gladders, M. D., Dressler, A., et al. 2016, ApJ, 832, 7

Aird, J., Coil, A. L., & Georgakakis, A. 2019, MNRAS, 484, 4360

Akhshik, M., Whitaker, K. E., Brammer, G., et al. 2020, ApJ, 900, 184

Akhshik, M., Whitaker, K. E., Leja, J., et al. 2021, ApJL, 907, L8

Alexander, D. M., & Hickox, R. C. 2012, NewAR, 56, 93

Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Balogh, M. L., Gilbank, D. G., Muzzin, A., et al. 2017, MNRAS, 470, 4168

Barišić, I., van der Wel, A., Bezanson, R., et al. 2017, ApJ, 847, 72

Barnes, J. E., & Hernquist, L. E. 1991, ApJL, 370, L65

Baron, D., Netzer, H., Poznanski, D., Prochaska, J. X., & Förster Schreiber, N. M. 2017, MNRAS, 470, 1687

Baron, D., Netzer, H., Prochaska, J. X., et al. 2018, MNRAS, 480, 3993

Barro, G., Faber, S. M., Pérez-González, P. G., et al. 2014, ApJ, 791, 52

Barro, G., Faber, S. M., Koo, D. C., et al. 2017, ApJ, 840, 47

Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57

Belfiore, F., Maiolino, R., Maraston, C., et al. 2016, MNRAS, 461, 3111

Bell, S., Newman, A. B., & Ellis, R. S. 2017, ApJ, 834, 18

Bell, S., Newman, A. B., & Ellis, R. S. 2019, ApJ, 874, 17

Bell, S., Newman, A. B., Ellis, R. S., & Konidaris, N. P. 2014, ApJL, 788, L29

Bernardi, M., Domínguez Sánchez, H., Brownstein, J. R., Drory, N., & Sheth, R. K. 2019, MNRAS, 489, 5633

Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393

Best, P. N., & Heckman, T. M. 2012, MNRAS, 421, 1569

Bezanson, R., van Dokkum, P., van de Sande, J., Franx, M., & Kriek, M. 2013, ApJL, 764, L8

Bezanson, R., van der Wel, A., Pacifici, C., et al. 2018, ApJ, 858, 60

Birnboim, Y., & Dekel, A. 2003, MNRAS, 345, 349

Booth, C. M., & Schaye, J. 2009, MNRAS, 398, 53

Bothwell, M. S., Smail, I., Chapman, S. C., et al. 2013, MNRAS, 429, 3047

Bower, R. G., Schaye, J., Frenk, C. S., et al. 2017, MNRAS, 465, 32

Bradač, M., Ryan, R., Casertano, S., et al. 2014, ApJ, 785, 108

Brammer, G. 2019, Grizli: Grism Redshift and Line Analysis Software, Astrophysics Source Code Library, ascl:1905.001

Bruce, V. A., Dunlop, J. S., Cirasuolo, M., et al. 2012, MNRAS, 427, 1666

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125

Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682

Carnall, A. C., Leja, J., Johnson, B. D., et al. 2019a, ApJ, 873, 44

Carnall, A. C., McLure, R. J., Dunlop, J. S., & Davé, R. 2018, MNRAS, 480, 4379

Carnall, A. C., McLure, R. J., Dunlop, J. S., et al. 2019b, MNRAS, 490, 417

Casertano, S., de Mello, D., Dickinson, M., et al. 2000, AJ, 120, 2747

Chabrier, G. 2003, PASP, 115, 763

Chang, Y.-Y., van der Wel, A., Rix, H.-W., et al. 2013a, ApJ, 762, 83

Chang, Y.-Y., van der Wel, A., Rix, H.-W., et al. 2013b, ApJ, 773, 149

Chauke, P., van der Wel, A., Pacifici, C., et al. 2018, ApJ, 861, 13

Ciesla, L., Elbaz, D., & Fensch, J. 2017, A&A, 608, A41

Cimatti, A., Daddi, E., Renzini, A., et al. 2004, Natur, 430, 184

Cimatti, A., Cassata, P., Pozzetti, L., et al. 2008, A&A, 482, 21

Ciotti, L., D'Ercole, A., Pellegrini, S., & Renzini, A. 1991, ApJ, 376, 380

Circosta, C., Mainieri, V., Padovani, P., et al. 2018, A&A, 620, A82

Coccato, L., Gerhard, O., & Arnaboldi, M. 2010, MNRAS, 407, L26

Coil, A. L., Weiner, B. J., Holz, D. E., et al. 2011, ApJ, 743, 46

Conroy, C., & Gunn, J. E. 2010, ApJ, 712, 833

Conroy, C., Gunn, J. E., & White, M. 2009, ApJ, 699, 486

Conroy, C., & van Dokkum, P. 2012, ApJ, 747, 69

Conroy, C., van Dokkum, P. G., & Kravtsov, A. 2015, ApJ, 803, 77

Conroy, C., Villaume, A., van Dokkum, P. G., & Lind, K. 2018, ApJ, 854, 139

Cox, T. J., Dutta, S. N., Di Matteo, T., et al. 2006, ApJ, 650, 791
Crook, A. C., Huchra, J. P., Martimbeau, N., et al. 2007, ApJ, 655, 790

Croton, D. J., Stevens, A. R. H., Tonini, C., et al. 2016, ApJS, 222, 22

Cullen, F., McLure, R. J., Dunlop, J. S., et al. 2019, MNRAS, 487, 2038

Daddi, E., Renzini, A., Pirzkal, N., et al. 2005, ApJ, 626, 680

Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, A&A, 605, A70

D'Eugenio, C., Daddi, E., Gobat, R., et al. 2020, ApJL, 898, L2

Di Matteo, P., Pipino, A., Lehnert, M. D., Combes, F., & Semelin, B. 2009, A&A, 499, 427

Di Matteo, T., Springel, V., & Hernquist, L. 2005, Natur, 433, 604

Domínguez Sánchez, H., Pérez-González, P. G., Esquej, P., et al. 2016, MNRAS, 457, 3743

D'Souza, R., Kauffman, G., Wang, J., & Vegetti, S. 2014, MNRAS, 443, 1433

Dubois, Y., Devriendt, J., Slyz, A., & Teissier, R. 2012, MNRAS, 420, 2662

Dunlop, J., Peacock, J., Spinrad, H., et al. 1996, Natur, 381, 581

Ebeling, H., Barrett, E., Donovan, D., et al. 2007, ApJL, 661, L33

Ebeling, H., Edge, A. C., & Henry, J. P. 2001, ApJ, 553, 668

Ebeling, H., Edge, A. C., Mantz, A., et al. 2010, MNRAS, 407, 83

Ebeling, H., Stockmann, M., Richard, J., et al. 2018, ApJL, 852, L7

Egami, E., Rex, M., Rawle, T. D., et al. 2010, A&A, 518, L12

Elbaz, D., Leiton, R., Nagar, N., et al. 2018, A&A, 616, A110

Emsellem, E., Cappellari, M., Krajnović, D., et al. 2011, MNRAS, 414, 888

Estrada-Carpenter, V., Papovich, C., Momcheva, I., et al. 2019, ApJ, 870, 133

Fabian, A. C. 2012, ARA&A, 50, 455

Falkendal, T., De Breuck, C., Lehnert, M. D., et al. 2019, A&A, 621, A27

Fanelli, M. N., O'Connell, R. W., Burstein, D., & Wu, C.-C. 1992, ApJS, 82, 197

Feltre, A., Charlot, S., & Gutkin, J. 2016, MNRAS, 456, 3354

Ferland, G. J., Chatzikos, M., Guzmán, F., et al. 2017, RMXAA, 53, 385

Feroz, F., Hobson, M. P., Cameron, E., & Pettitt, A. N. 2019, OJAp, 2, 10

Ferré-Mateu, A., Trujillo, I., Martín-Navarro, I., et al. 2017, MNRAS, 467, 1929

Ferreras, I., Pasquali, A., Pirzkal, N., et al. 2019b, MNRAS, 486, 1358

Ferreras, I., Scott, N., La Barbera, F., et al. 2019a, MNRAS, 489, 608

Fitzpatrick, E. L. 1999, PASP, 111, 63

Foreman-Mackey, D. 2016, JOSS, 1, 24

Forrest, B., Annunziatella, M., Wilson, G., et al. 2020, ApJL, 890, L1

Förster Schreiber, N. M., Übler, H., Davies, R. L., et al. 2019, ApJ, 875, 21

Franx, M., Labbé, I., Rudnick, G., et al. 2003, ApJL, 587, L79

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1

Gallazzi, A., Bell, E. F., Zibetti, S., Brinchmann, J., & Kelson, D. D. 2014, ApJ, 788, 72

Gallazzi, A., Charlot, S., Brinchmann, J., & White, S. D. M. 2006, MNRAS, 370, 1106

Gallazzi, A., Charlot, S., Brinchmann, J., White, S. D. M., & Tremonti, C. A. 2005, MNRAS, 362, 41

Geier, S., Richard, J., Man, A. W. S., et al. 2013, ApJ, 777, 87

Gladders, M. D., Oemler, A., Dressler, A., et al. 2013, ApJ, 770, 64

Glazebrook, K., Schreiber, C., Labbé, I., et al. 2017, Natur, 544, 71

Gómez-Guijarro, C., Magdis, G. E., Valentino, F., et al. 2019, ApJ, 886, 88

González Delgado, R. M., García-Benito, R., Pérez, E., et al. 2015, A&A, 581, A103

Greene, J. E., Janish, R., Ma, C.-P., et al. 2015, ApJ, 807, 11

Greene, J. E., Murphy, J. D., Graves, G. J., et al. 2013, ApJ, 776, 64

Greene, J. E., Veale, M., Ma, C.-P., et al. 2019, ApJ, 874, 66

Halliday, C., Daddi, E., Cimatti, A., et al. 2008, A&A, 479, 417

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Heckman, T. M., & Best, P. N. 2014, ARA&A, 52, 589

Heckman, T. M., Kauffmann, G., Brinchmann, J., et al. 2004, ApJ, 613, 109

Hickox, R. C., Mullaney, J. R., Alexander, D. M., et al. 2014, ApJ, 782, 9

Hill, A. R., Muzzin, A., Franx, M., & van der Sande, J. 2016, ApJ, 819, 74

Hill, A. R., van der Wel, A., Franx, M., et al. 2019, ApJ, 871, 76

Hirschmann, M., Naab, T., Ostriker, J. P., et al. 2015, MNRAS, 449, 528

Hopkins, P. F., Cox, T. J., Kereš, D., & Hernquist, L. 2008, ApJS, 175, 390

Hopkins, P. F., Bundy, K., Croton, D., et al. 2010, ApJ, 715, 202

Huang, K.-H., Bradač, M., Lemaux, B. C., et al. 2016, ApJ, 817, 11

Hunter, J. D. 2007, CSE, 9, 90

Iyer, K. G., Gawiser, E., Faber, S. M., et al. 2019, ApJ, 879, 116

Jafariyazani, M., Newman, A. B., Mobasher, B., et al. 2020, ApJL, 897, L42

Jauzac, M., Klein, B., Kneib, J.-P., et al. 2021, MNRAS, in press

Johansson, J., Woods, T. E., Gilfanov, M., et al. 2016, MNRAS, 461, 4505

Jullo, E., Kneib, J.-P., Limousin, M., et al. 2007, NJPh, 9, 447

Kado-Fong, E., Marchesini, D., Marsan, Z. C., et al. 2017, ApJ, 838, 57

Kauffmann, G., Heckman, T. M., Tremonti, C., et al. 2003, MNRAS, 346, 1055

Kausch, W., Noll, S., Smette, A., et al. 2015, A&A, 576, A78

Kennicutt, R. C., Jr. 1998, ARA&A, 36, 189

Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531

Kewley, L. J., & Dopita, M. A. 2002, ApJS, 142, 35

Kewley, L. J., Groves, B., Kauffmann, G., & Heckman, T. 2006, MNRAS, 372, 961

Kewley, L. J., Nicholls, D. C., & Sutherland, R. S. 2019, ARA&A, 57, 511

Khochfar, S., & Ostriker, J. P. 2008, ApJ, 680, 54

Khochfar, S., Emsellem, E., Serra, P., et al. 2011, MNRAS, 417, 845

Khoperskov, S., Haywood, M., Di Matteo, P., Lehnert, M. D., & Combes, F. 2018, A&A, 609, A60

King, A. R., Zubovas, K., & Power, C. 2011, MNRAS, 415, L6

Krajnović, D., Alatalo, K., Blitz, L., et al. 2013, MNRAS, 432, 1768

Krajnović, D., Ural, U., Kuntschner, H., et al. 2020, A&A, 635, A129

Kriek, M., van Dokkum, P. G., Labbé, I., et al. 2009, ApJ, 700, 221

Kriek, M., Conroy, C., van Dokkum, P. G., et al. 2016, Natur, 540, 248

Kriek, M., Price, S. H., Conroy, C., et al. 2019, ApJL, 880, L31

Kroupa, P. 2001, MNRAS, 322, 231

Kuntschner, H., Emsellem, E., Bacon, R., et al. 2010, MNRAS, 408, 97

Lackner, C. N., Cen, R., Ostriker, J. P., & Joung, M. R. 2012, MNRAS, 425, 641

Lagos, C. d. P., Schaye, J., Bahé, Y., et al. 2018a, MNRAS, 476, 4327

Lagos, C. d. P., Stevens, A. R. H., Bower, R. G., et al. 2018b, MNRAS, 473, 4956

Leja, J., Carnall, A. C., Johnson, B. D., Conroy, C., & Speagle, J. S. 2019, ApJ, 876, 3

Leja, J., Johnson, B. D., Conroy, C., van Dokkum, P. G., & Byler, N. 2017, ApJ, 837, 170

Leung, G. C. K., Coil, A. L., Aird, J., et al. 2019, ApJ, 886, 11

Li, Y.-S., & White, S. D. M. 2008, MNRAS, 384, 1459

Limousin, M., Ebeling, H., Ma, C.-J., et al. 2010, MNRAS, 405, 777

Lonoce, I., Longhetti, M., Maraston, C., et al. 2015, MNRAS, 454, 3912

MacKenzie, T. P., Scott, D., Smail, I., et al. 2014, MNRAS, 445, 201

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Maíz Apellániz, J. 2006, AJ, 131, 1184

Maltby, D. T., Almaini, O., McLure, R. J., et al. 2019, MNRAS, 489, 1139

Man, A., & Belli, S. 2018, NatAs, 2, 695

Man, A. W. S., Lehnert, M. D., Vernet, J. D. R., De Breuck, C., & Falkendal, T. 2019, A&A, 624, A81

Man, A. W. S., Zirm, A. W., & Toft, S. 2016a, ApJ, 830, 89

Man, A. W. S., Greve, T. R., Toft, S., et al. 2016b, ApJ, 820, 102

Maraston, C., Nieves Colmenárez, L., Bender, R., & Thomas, D. 2009, A&A, 493, 425

Maraston, C., Pforr, J., Renzini, A., et al. 2010, MNRAS, 407, 830

Martig, M., Bournaud, F., Teissier, R., & Dekel, A. 2009, ApJ, 707, 250

Martín-Navarro, I., Vazdekis, A., Falcón-Barroso, J., et al. 2018, MNRAS, 475, 3700

Mathews, W. G. 1990, ApJ, 354, 468

McCarthy, P. J., Le Borgne, D., Crampton, D., et al. 2004, ApJL, 614, L9

McConnell, N. J., & Ma, C.-P. 2013, ApJ, 764, 184

McDermid, R. M., Alatalo, K., Blitz, L., et al. 2015, MNRAS, 448, 3484

Merlin, E., Bourne, N., Castellano, M., et al. 2016, A&A, 595, A97

Merlin, E., Fontana, A., Ferguson, H. C., et al. 2015, A&A, 582, A15

Merlin, E., Fontana, A., Castellano, M., et al. 2018, MNRAS, 473, 2098

Modigliani, A., Goldoni, P., Royer, F., et al. 2010, Proc. SPIE, 7737, 773728

Momcheva, I. G., Brammer, G. B., van Dokkum, P. G., et al. 2016, ApJS, 225, 27

Moster, B. P., Naab, T., & White, S. D. M. 2020, MNRAS, 499, 4748

Muzzin, A., Labbé, I., Franx, M., et al. 2012, ApJ, 761, 142

Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJ, 777, 18

Naab, T., Jesseit, R., & Burkert, A. 2006, MNRAS, 372, 839

Naab, T., Oser, L., Emsellem, E., et al. 2014, MNRAS, 444, 3357

Nesvadba, N. P. H., De Breuck, C., Lehnert, M. D., Best, P. N., & Collet, C. 2017, A&A, 599, A123

Netzer, H. 1990, in Active Galactic Nuclei, ed. R. D. Blandford (Saas-Fee: Swiss Society for Astrophysics and Astronomy), 57

Newman, A. B., Belli, S., & Ellis, R. S. 2015, ApJL, 813, L7

Newman, A. B., Belli, S., Ellis, R. S., & Patel, S. G. 2018a, ApJ, 862, 125

Newman, A. B., Belli, S., Ellis, R. S., & Patel, S. G. 2018b, ApJ, 862, 126

Oke, J. B., & Gunn, J. E. 1983, ApJ, 266, 713

Olsen, K. P., Rasmussen, J., Toft, S., & Zirm, A. W. 2013, ApJ, 764, 4

Onodera, M., Carollo, C. M., Renzini, A., et al. 2015, ApJ, 808, 161

Onodera, M., Daddi, E., Gobat, R., et al. 2010, ApJL, 715, L6

Onodera, M., Renzini, A., Carollo, M., et al. 2012, ApJ, 755, 26

Oser, L., Ostriker, J. P., Naab, T., Johansson, P. H., & Burkert, A. 2010, ApJ, 725, 2312

Oyarzún, G. A., Bundy, K., Westfall, K. B., et al. 2019, ApJ, 880, 111

Pacifici, C., Kassin, S. A., Weiner, B. J., et al. 2016, ApJ, 832, 79
Padovani, P., Alexander, D. M., Assef, R. J., et al. 2017, A&ARv, 25, 2

Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2002, AJ, 124, 266

Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2010, AJ, 139, 2097

Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2011, GALFIT: Detailed Structural Decomposition of Galaxy Images, Astrophysics Source Code Library, ascl:1104.010

Popping, G., Decarli, R., Man, A. W. S., et al. 2017, A&A, 602, A11

Postman, M., Coe, D., Benítez, N., et al. 2012, ApJS, 199, 25

Praderie, F., Talavera, A., & Lamers, H. J. G. L. M. 1980, A&A, 86, 271

Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M., et al. 2018, AJ, 156, 123

Robitaille, T., & Bressert, E. 2012, APLpy: Astronomical Plotting Library in Python, Astrophysics Source Code Library, ascl:1208.017

Rodríguez-Gomez, V., Pillepich, A., Sales, L. V., et al. 2016, MNRAS, 458, 2371

Rodríguez-Merino, L. H., Chavez, M., Bertone, E., & Buzzoni, A. 2005, ApJ, 626, 411

Rubin, K. H. R., Prochaska, J. X., Ménard, B., et al. 2011, ApJ, 728, 55

Rupke, D. S. N., Coil, A., Geach, J. E., et al. 2019, Natur, 574, 643

Saracco, P., La Barbera, F., Gargiulo, A., et al. 2019, MNRAS, 484, 2281

Schawinski, K., Virani, S., Simmons, B., et al. 2009, ApJL, 692, L19

Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103

Schönebeck, F., Puzia, T. H., Pasquali, A., et al. 2014, A&A, 572, A13

Schreiber, C., Glazebrook, K., Nanayakkara, T., et al. 2018, A&A, 618, A85

Sell, P. H., Tremonti, C. A., Hickox, R. C., et al. 2014, MNRAS, 441, 3417

Sérsic, J. L. 1963, BAAA, 6, 41

Sersic, J. L. 1968, Atlas de Galaxias Australes (Cordoba: Observatorio Astronomico)

Shen, J., Man, A. W. S., Zabl, J., et al. 2021, ApJ, 917, 79

Sijacki, D., Springel, V., Di Matteo, T., & Hernquist, L. 2007, MNRAS, 380, 877

Sijacki, D., Vogelsberger, M., Genel, S., et al. 2015, MNRAS, 452, 575

Silverman, J. D., Lamareille, F., Maier, C., et al. 2009, ApJ, 696, 396

Singh, R., van de Ven, G., Jahnke, K., et al. 2013, A&A, 558, A43

Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163

Smette, A., Sana, H., Noll, S., et al. 2015, A&A, 576, A77

Snow, T. P., Lamers, H. J. G. L. M., Lindholm, D. M., & Odell, A. P. 1994, ApJS, 95, 163

Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15

Spilker, J. S., Bezanson, R., Marrone, D. P., et al. 2016, ApJ, 832, 19

Spinrad, H., Dey, A., Stern, D., et al. 1997, ApJ, 484, 581

Springel, V., Di Matteo, T., & Hernquist, L. 2005, MNRAS, 361, 776

Stockmann, M., Toft, S., Gallazzi, A., et al. 2020, ApJ, 888, 4

Stockton, A., McGrath, E., Canalizo, G., Iye, M., & Maihara, T. 2008, ApJ, 672, 146

Storey, P. J., & Zeippen, C. J. 2000, MNRAS, 312, 813

Su, K.-Y., Hopkins, P. F., Hayward, C. C., et al. 2019, MNRAS, 487, 4393

Su, K.-Y., Hopkins, P. F., Hayward, C. C., et al. 2020, MNRAS, 491, 1190

Tacconi, L. J., Genzel, R., Saintonge, A., et al. 2018, ApJ, 853, 179

Takata, T., Kashikawa, N., Nakanishi, K., et al. 2003, PASJ, 55, 789

Taniguchi, Y., Kajisawa, M., Kobayashi, M. A. R., et al. 2015, ApJL, 809, L7

Taylor, E. N., Franx, M., van Dokkum, P. G., et al. 2009, ApJS, 183, 295

Taylor, P., & Kobayashi, C. 2017, MNRAS, 471, 3856

Thomas, D., Maraston, C., & Bender, R. 2003, MNRAS, 339, 897

Thomas, D., Maraston, C., Bender, R., & Mendes de Oliveira, C. 2005, ApJ, 621, 673

Thomas, D., Maraston, C., Schawinski, K., Sarzi, M., & Silk, J. 2010, MNRAS, 404, 1775

Toft, S., Gallazzi, A., Zirm, A., et al. 2012, ApJ, 754, 3

Toft, S., Smolčić, V., Magnelli, B., et al. 2014, ApJ, 782, 68

Toft, S., Zabl, J., Richard, J., et al. 2017, Natur, 546, 510

Tremonti, C. A., Moustakas, J., & Diamond-Stanic, A. M. 2007, ApJL, 663, L77

Trujillo, I., Ferré-Mateu, A., Balcells, M., Vazdekis, A., & Sánchez-Blázquez, P. 2014, ApJL, 780, L20

Trump, J. R., Hall, P. B., Reichard, T. A., et al. 2006, ApJS, 165, 1

Valentino, F., Tanaka, M., Davidzon, I., et al. 2020, ApJ, 889, 93

van de Sande, J., Kriek, M., Franx, M., et al. 2011, ApJL, 736, L9

van de Sande, J., Kriek, M., Franx, M., et al. 2013, ApJ, 771, 85

van der Wel, A., Franx, M., van Dokkum, P. G., et al. 2014, ApJ, 788, 28

van der Wel, A., Rix, H.-W., Wuyts, S., et al. 2011, ApJ, 730, 38

van Dokkum, P. G., Kriek, M., & Franx, M. 2009, Natur, 460, 717

van Dokkum, P. G., Nelson, E. J., Franx, M., et al. 2015, ApJ, 813, 23

Veale, M., Ma, C.-P., Greene, J. E., et al. 2017a, MNRAS, 471, 1428

Veale, M., Ma, C.-P., Greene, J. E., et al. 2018, MNRAS, 473, 5446

Veale, M., Ma, C.-P., Thomas, J., et al. 2017b, MNRAS, 464, 356

Vernet, J., Dekker, H., D'Odorico, S., et al. 2011, A&A, 536, A105

Volonteri, M., Capelo, P. R., Netzer, H., et al. 2015a, MNRAS, 449, 1470

Volonteri, M., Capelo, P. R., Netzer, H., et al. 2015b, MNRAS, 452, L6

Weiner, B. J., Coil, A. L., Prochaska, J. X., et al. 2009, ApJ, 692, 187

Wellons, S., Torrey, P., Ma, C.-P., et al. 2016, MNRAS, 456, 1030

Westfall, K. B., Cappellari, M., Bershady, M. A., et al. 2019, AJ, 158, 231

Whitaker, K. E., van Dokkum, P. G., Brammer, G., et al. 2013, ApJL, 770, L39

Whitaker, K. E., Williams, C. C., Mowla, L., et al. 2021, Natur, in press

Wild, V., Heckman, T., & Charlot, S. 2010, MNRAS, 405, 933

Williams, R. J., Quadri, R. F., Franx, M., van Dokkum, P., & Labbé, I. 2009, ApJ, 691, 1879

Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868

Wu, P.-F., van der Wel, A., Gallazzi, A., et al. 2018, ApJ, 855, 85

Wu, P.-F., van der Wel, A., Bezanson, R., et al. 2020, ApJ, 888, 77

Wuyts, S., Cox, T. J., Hayward, C. C., et al. 2010, ApJ, 722, 1666

Wuyts, S., Labbé, I., Franx, M., et al. 2007, ApJ, 655, 51

Yan, R., & Blanton, M. R. 2012, ApJ, 747, 61

Yan, R., Newman, J. A., Faber, S. M., et al. 2006, ApJ, 648, 281

Zabl, J., Nørgaard-Nielsen, H. U., Fynbo, J. P. U., et al. 2015, MNRAS, 451, 2050

Zakamska, N. L., Strauss, M. A., Krolik, J. H., et al. 2003, AJ, 126, 2125

Zhu, G. B., Comparat, J., Kneib, J.-P., et al. 2015, ApJ, 815, 48

Zibetti, S., Gallazzi, A. R., Hirschmann, M., et al. 2020, MNRAS, 491, 3562

Zick, T. O., Kriek, M., Shapley, A. E., et al. 2018, ApJL, 867, L16

Zitrin, A., Broadhurst, T., Barkana, R., Rephaeli, Y., & Benítez, N. 2011, MNRAS, 410, 1939

Zubovas, K., & King, A. R. 2014, MNRAS, 439, 400