# A fast-rotator post-starburst galaxy quenched by supermassive black-hole feedback at z = 3

Francesco d'Eugenio, Pablo G. Pérez-González, Roberto Maiolino, Jan Scholtz, Michele Perna, Chiara Circosta, Hannah Übler, Santiago Arribas, Torsten Böker, Andrew J. Bunker, et al.

► To cite this version:

Francesco d'Eugenio, Pablo G. Pérez-González, Roberto Maiolino, Jan Scholtz, Michele Perna, et al. A fast-rotator post-starburst galaxy quenched by supermassive black-hole feedback at z = 3. Nature Astronomy, 2024, 8, pp.1443-1456. 10.1038/s41550-024-02345-1. insu-04851487

HAL Id: insu-04851487

https://insu.hal.science/insu-04851487v1

Submitted on 21 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# A fast-rotator post-starburst galaxy quenched by supermassive black-hole feedback at z = 3

Received: 11 August 2023

Accepted: 22 July 2024

Published online: 16 September 2024

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

Check for updates

A list of authors and their affiliations appears at the end of the paper

The most massive galaxies in the Universe stopped forming stars due to the time-integrated feedback from central supermassive black holes (SMBHs). However, the exact quenching mechanism is not yet understood, because local massive galaxies were quenched billions of years ago. Here we present JWST/NIRSpec integral-field spectroscopy observations of GS-10578, a massive, quiescent galaxy at redshift z = 3.064 ± 0.002. From its spectrum, we measure a stellar mass M* = 1.6 ± 0.2 × 10^11 M⊙ and a dynamical mass Mdyn = 2.0 ± 0.5 × 10^11 M⊙. Half of its stellar mass formed at z = 3.7–4.6, and the system is now quiescent, with a current star-formation rate of less than 19 M⊙ yr⁻¹. We detect ionized- and neutral-gas outflows traced by [O III] emission and Na I absorption, with mass outflow rates 0.14–2.9 and 30–100 M⊙ yr⁻¹, respectively. Outflow velocities reach vout ≈ 1,000 km s⁻¹, comparable to the galaxy escape velocity. GS-10578 hosts an active galactic nucleus, evidence that these outflows are due to SMBH feedback. The neutral outflow rate is higher than the star-formation rate. Hence, this is direct evidence for ejective SMBH feedback, with a mass loading capable of interrupting star formation by rapidly removing its fuel. Stellar kinematics show ordered rotation, with spin parameter λR_e = 0.62 ± 0.07, meaning GS-10578 is rotation-supported. This study presents direct evidence for ejective active galactic nucleus feedback in a massive, recently quenched galaxy, thus helping to clarify how SMBHs quench their hosts. The high value of λR_e implies that quenching can occur without destroying the stellar disk.

The Universe today is not what it used to be¹. Local, massive, quiescent galaxies stand like colossal wrecks of glorious but remote star-formation histories (SFHs) and mighty and rapid quenching, the likes of which have no present-day equals²,³. The James Webb Space Telescope (JWST) has enabled us for the first time to witness these monumental galaxies during the long-gone epoch when they arose and fell. By redshift z = 1.5–2, 3–4 Gyr after the Big Bang, massive quiescent galaxies have little to no cold gas, the fuel for star formation⁴,⁵. However, the question whether the missing fuel was consumed by starbursts or whether it was removed by ‘ejective’ feedback from supermassive black

holes (SMBHs) remains open⁶,⁷. We present integral-field spectroscopy (IFS) observations made by the near-infrared spectrograph (NIRSpec) onboard JWST during the Cosmic Assembly Near-infrared Deep Extra-galactic Legacy Survey (CANDELS) of GS-10578, a massive, quiescent galaxy at redshift z = 3.064 (look-back time 11.7 Gyr). Observed as part of the Galaxy Assembly with NIRSpec/IFS (GA-NIFS) guaranteed time observations (GTOs), this galaxy (Fig. 1a) was selected as a ‘blue nugget’⁸, a class of massive, extremely compact galaxies (stellar mass M* = 3–30 × 10^10 M⊙, half-light radius R_e = 0.5–2 kpc) thought to be the progenitors of compact quiescent galaxies at z = 2 (‘red nuggets’⁹).

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

e-mail: francesco.deugenio@gmail.com
<figure><img src="image_5.png" /></figure>

**Fig. 1** | Summary of GS-10578 observations. **a**, NIRCam false-colour image. The grey dashed square is the field of view of NIRSpec/IFS (with dithering). The solid white contour is a stellar isophote from NIRSpec/IFS prism spectroscopy. The red dashed ellipse is the best-fitting deconvolved model, with semimajor axis $1R_e$. The cyan contours are rest-frame UV from NIRCam F090W, which display the complex morphology, probably due to foreground dust. The orange and blue contours are isosignificance lines for the $H\alpha$-[N II] and $H\beta$-[O III] blended lines, from the prism data. Nebular emission is significantly more extended than the stellar light, indicating the ability of mergers and feedback to enrich the circumgalactic medium. **b**, Detail of the central region of **a**; the dashed black square outlines a single NIRSpec spaxel, whose spectrum is reproduced in panels **f–i**, **c**, Comparing the model isophote (red dashed ellipse) to NIRCam isophote morphology, there is a clear wavelength-dependent offset to the north-west. **d**, Brightest [O III] emission from g235h/f170lp data, showing a slight

offset from the centre of the stellar emission (red dot, distance <1 spaxel) and a notably different position angle, almost $\pi/2$ compared to the stars. **e**, Neutral gas Na I absorption, which is clearly asymmetric and located where the rest-UV photometry is faintest. **f**, Single-spxel spectrum centred on the $H\beta$-[O III] region, showing the data (grey), best-fitting continuum and nebular-line model (red and blue, respectively), and their sum (orange). **g**, Same as **f**, but centred on $H\alpha$-[N II]. **h**, Same as **f**, but centred on $H\delta$ and Hy. **i**, Same as **f**, but centred on the Na I doublet. The horizontal bar in the top right corner of each panel is $\delta v = 500 \text{ km s}^{-1}$. **f** and **g** show the complex, multi-component emission-line profile of [O III] and $H\alpha$ [N II]. **h** and **i** show the same spaxel. The vertical dotted lines mark the position of the stellar $H\delta$, Hy and Na I-doublet absorption. The deep Na I absorption, blueshifted with respect to the stars, betrays the presence of a neutral-gas outflow. Dec., declination; FOV, field of view; RA, right ascension.

Blue nuggets are believed to be undergoing 'gas-rich compaction', which means that a central starburst driven by disk instability or gas-rich major mergers is followed by rapid quenching¹⁰ to leave a compact, quiescent, red-nugget galaxy. As we will show, GS-10578 is, instead, already a red nugget in an advanced stage of quenching. The system is merging with several low-mass satellites and is undergoing powerful, ejective feedback from its SMBH.

Imaging with the near-infrared camera (NIRCam)¹¹⁻¹³ onboard JWST revealed a system with several smaller companions, all within 8 kpc in projection (labelled GS-10578b–f in Fig. 1a, of which GS-10578b–d have been spectroscopically confirmed to be satellites). It is common for massive quiescent galaxies at high redshift to have several low-mass satellites¹⁴. GS-10578 has a regular, elliptical shape at 2.77–4.44 µm, but the light distribution is increasingly lopsided blueward of 2 µm (Fig. 1b), with the brightest region at 0.9 µm clearly offset to the north-west (by 0.5 kpc) and displaying two peaks.

We used the NIRSpec/IFS data with the low-resolution prism disperser to perform spaxel-by-spxel full spectral energy distribution (SED) modelling, which allowed us to measure the surface stellar mass density, spatially resolved stellar-population properties and SFH (Methods). The resulting mass distribution is symmetric, with no strong signs of asymmetry (Fig. 2a). The asymmetric light distribution at rest-frame ultraviolet (UV) wavelengths arises from younger stellar populations (luminous but with low stellar mass) or from an asymmetric dust distribution. By integrating the two-dimensional (2D) properties inside a circular aperture of 1.15 arcsec radius and correcting for aperture losses, we measured the total stellar mass $M_*$ = 1.6 ± 0.2 × 10¹¹ M☉ (surviving stellar mass). The total SFH from all spaxels shows that the main star-formation episode happened at $z = 3.7–4.6$, 0.4–0.8 Gyr before

the observation. Afterwards, the star-formation rate (SFR) rapidly declined. There is some evidence of a recent upturn, but our models do not include an active galactic nucleus (AGN), therefore this upturn could be due to the nebular emission being interpreted as due to star formation. The SFR averaged over the last 100 Myr is 40 ± 20 M☉ yr⁻¹. This is five times below the star-forming main sequence at $z = 3$ (ref. 15) and comparable to the quiescence threshold of 19 M☉ yr⁻¹ (ref. 16). However, the nebular emission is mostly powered by the AGN, and part of the SFR is due to the satellites (as explained later), meaning that the true SFR is probably much lower. This implies that GS-10578 is presently quenched and on its way to quiescence, which is consistent with the empirical UVJ colour diagnostic ($U - V = 1.45 ± 0.05$ and $V - J = 0.75 ± 0.08$ mag; ref. 17) and with the SFR inferred from the emission-line analysis.

The NIRSpec/IFS observation in the low spectral-resolution configuration reveals extended nebular emission in both the $H\beta$-[O III] and $H\alpha$-[N II] spectrally blended complexes (Fig. 1a). This nebular emission is elongated in a different direction compared to the stellar major axis but in the same direction as GS-10578b, which suggests that the ionized gas comes from interstellar medium (ISM) material stripped by the interaction with the satellite. An outflow origin is also possible, but in this case, the alignment with the satellite would be a coincidence.

Unlike the low-resolution observations, NIRSpec high-resolution data (Fig. 1c,d,f,g) are able to separate $H\alpha$ from [N II] and $H\beta$ (which in GS-10578 we observe mostly in absorption) from [O III]. The emission lines show a broad, multi-component profile, indicating complex kinematics. The BPT diagnostic diagram (Fig. 3) shows that GS-10578 has everywhere a high [N II]/$H\alpha$ ratio (emission-line ratios are measured inside concentric elliptical annuli, following the shape and position angle of the stellar isophotes). From the outermost to the innermost
<figure><img src="image_6.png" /><figcaption>**Fig. 2 | Summary of SED fitting.** a,b, Stellar-population surface mass-density map (a) and mass-weighted age map (b) from full spectral fitting of the low-resolution spectroscopy data, showing that most of GS-10578 consists of 0.5-Gyr-old stars. Merging satellites have lower mass-to-light ratios than GS-10578 and do not contribute to the mass-weighted properties. c, The SFH from the same models, integrated over the 2D map, shows the main episode of star formation occurred 0.5 Gyr ago. The solid line and shaded region are the median and $\sigma$ uncertainty from the 658 SFHs of each individual spaxel. The cyan diamond and black circle are the fiducial mass outflow rates from [O III] and Na I,</figcaption></figure>

respectively, with the error bars representing systematic uncertainties. Unlike the ionized-gas outflow, the neutral-gas outflow has a mass loading high enough to stop star formation. d, SED of the individual spaxels summed inside a circular aperture of radius 1.15 arcsec, for both the NIRSpec/IFS data (black) and the best-fitting model (red). We also show the aperture photometry as error bars (representing the $\sigma$ measurement uncertainty). The model does not include an AGN. Therefore, it does not capture the emission lines or the flux redder than 6 µm, which is dominated by host-dust emission from the AGN torus. For a model including an AGN torus, see Methods.

annulus around GS-10578, [O III] $\lambda$5007/H$\beta$ increases from 4 to 10, whereas [N II] $\lambda$6584/H$\alpha$ remains approximately constant at 5–8. All annuli lie outside the star-forming region of the BPT diagram, well beyond the demarcation curve between starburst-driven and AGN-driven photoionization¹⁸. Moreover, all annuli lie far outside the 99th percentile of the distribution of local galaxies and AGNs (grey and black contours). Star-forming and AGN galaxies at z = 2 also have much smaller [N II] $\lambda$6584/H$\alpha$ ratios (empty circles), although our values are like those of post-starburst (PSB) galaxies at z = 1–2.5 (ref. 19) and to the most extreme PSB galaxies in the local Universe (grey star). Assuming all H$\alpha$ emission was due to star formation, the extinction-corrected SFR on timescales of 3–10 Myr is 9 M$_{\odot}$ yr$^{-1}$ (inside $2R_e$, where $R_e = 1.1$ kpc from ref. 20). On varying the aperture size, the SFR ranged between 4 M$_{\odot}$ yr$^{-1}$ (inside $1R_e$) to 19 M$_{\odot}$ yr$^{-1}$ (summing over the entire map, including the satellites). However, the satellites contribute substantially to the SFR. We measure $7.4_{-2.1}^{+3.4}$, $4.6_{-2.6}^{+5.9}$ and $0.5_{-0.2}^{+0.9}$ M$_{\odot}$ yr$^{-1}$, respectively for GS-10578b, c and d. The SFRs based on H$\alpha$ are, overall, like the result from the SED analysis. However, as we have seen, the H$\alpha$ emission is not powered by star formation. Therefore, the true SFR is probably much lower, and we treat all these values as upper limits within their apertures. Overall, the low H$\alpha$ SFR values are in agreement with the

quiescent interpretation. The BPT diagram shows both high [O III] $\lambda$5007/H$\beta$ and [N II] $\lambda$6584/H$\alpha$, which may be due to AGN, shocks or a combination of the two. In addition to circumstantial evidence from the BPT diagram, direct evidence for an AGN comes from the X-ray detection ($L_X = 6.4 \times 10^{44}$ erg s$^{-1}$)²¹ and from the SED excess at medium-infrared (MIR) wavelengths, which can be explained by the presence of a dusty torus around the accreting SMBH (Fig. 2d, compare with the model and photometry at 24 µm). Modelling the photometry using both stellar and AGN emission models gives an AGN bolometric luminosity $L_{AGN} = 2.6 \pm 0.9 \times 10^{45}$ erg s$^{-1}$ (Methods)²¹.

The spatially resolved emission-line kinematics (Fig. 4) displays line broadening in excess of 500 km s$^{-1}$ and centroids that are blueshifted relative to the stellar absorption (the stellar distribution is traced by the red dashed ellipse in Fig. 4). Most of the nebular flux comes from the central regions, but is offset from the centre of the stellar emission and is elongated almost at a right angle compared to the stellar isophote (Fig. 1d). There is a clear kinematic transition between this bright central region and the more diffuse emission at larger radii. The central region has a velocity offset of –400 km s$^{-1}$, whereas the extended emission is at velocities of ±100 km s$^{-1}$. Furthermore, the central regions have interpercentile line widths $w_{80}$ in excess
<figure><img src="image_7.png" /><figcaption>**Fig. 3** | Photoionization diagnostic diagram for GS-10578. The emission-line ratios of our target place its ionization source firmly beyond the range that can be explained by star formation (dashed line)¹⁸. For the main target, we show the emission-line ratios from spectra integrated inside elliptical apertures (colour-coded by the semimajor axis of the aperture). The yellow circle is the eastern satellite. We show a number of samples for comparison. The contours and the small empty circles trace the distribution of galaxies at z = 0.1 (from the Sloan Digital Sky Survey (SDSS)) and z = 2 (from KLEVER)⁹⁸, respectively. The sand-coloured symbols and upper limits are a sample of PSB galaxies at z = 1–2 (ref. 99). The grey star is one of the most extreme PSB galaxies in SDSS. The extreme values we find for GS-10578 are like what are found for other PSB galaxies.</figcaption></figure>

of 1,000 km s⁻¹. The extent and elongation of the diffuse nebular emission suggest a merger origin, but most of the flux is probably associated with the blueshifted component, which we interpret as an outflow. The large velocity offset and velocity width of the nebular lines—and that the contours of highest flux are directed along the galaxy minor axis (Fig. 1d and Fig. 4a–c)—are evidence for an outflow normal to the disk²²,²³. Overall, the fast outflow velocity, the low SFR and the AGN are evidence for SMBH-feedback-driven outflows²⁴⁻²⁶; but see refs. 27,28 for cases of fast starburst-driven outflows. We estimate a mass outflow rate $M_{\text{out}}^{[\text{III}]} = 0.14–2.9 M_{\odot} \text{ yr}^{-1}$ (Methods), meaning the ionized-gas outflow is much lower than the SFR (Fig. 2c).

What sets GS-10578 apart from other high-redshift sources, however, is the clear detection of spatially resolved, deep absorption features. These include the prominent hydrogen absorption lines characteristic of 0.5-Gyr-old stellar populations (equivalent widths (EWs) $H\delta_A = 6.6–8.1 \AA$ and $H\gamma_A = 4.1–6.8 \AA$; Figs. 1h and 5e–f) and Na I resonant absorption (Fig. 1i). The Balmer absorption confirms the results of the SED analysis (which is based on spectroscopy with much lower spectral resolution). The Na I absorption is offset both spatially and kinematically from the stars. Therefore, we ruled out a stellar origin. The highest Na I EW translates into a column density of $6 \times 10^{22} \text{ cm}^{-2}$, ten times lower than the column density towards the X-ray AGN ($6 \times 10^{23} \text{ cm}^{-2}$; ref. 21). The Na I EW maps (Fig. 1e) show that the absorption extends from the centre of the galaxy to the south-east, avoiding completely the north-western half, where we detected the brightest rest-UV emission. This fact, combined with the blueshifted line-of-sight velocities ($v_{\text{Na I}} = -(250–350) \text{ km s}^{-1}$), indicates that the absorption is due to a neutral-gas outflow with a cone-like geometry possibly tilted with respect to the stellar disk, whose existence is inferred based on stellar kinematics as discussed below. We estimated the disk inclination from the best-fitting Sérsic model of ref. 20, which gives a projected axis ratio $q = 0.44$ (we assumed an intrinsic axis ratio $q_0 = 0.2$). After correcting for inclination by a factor $\cos i$, we obtained an outflow velocity of the order of $-1,000 \text{ km s}^{-1}$, again too fast to be explained by star formation alone. Assuming a maximal escape

velocity of $v_{\text{esc}} = 1,200 \text{ km s}^{-1}$, we inferred a velocity ratio $v_{\text{Na I}}/v_{\text{esc}} = 0.8 \pm 0.2$, so that the outflows are certainly capable of removing material from the star-forming disk, with a fraction possibly escaping the galaxy altogether.

We assumed an outflow extent $R_{\text{out}} = 2.7 \text{ kpc}$, based on the largest projected distance between the Na I detection and the centre of the galaxy. With this $R_{\text{out}}$, the mass outflow rate is $\dot{M}_{\text{out}}^{\text{NaI}} \approx 30–100 M_{\odot} \text{ yr}^{-1}$, depending on the assumed geometry (Methods). Even with the lowest estimate, $\dot{M}_{\text{out}}^{\text{NaI}}$ is still higher than the Hα-based SFR, particularly if we recall that the latter is overestimated. This indicates that the neutral-gas outflow is capable of halting star formation by removing the necessary fuel. Indeed, non-detection from existing submillimetre observations of CO emission give a molecular-gas mass fraction $f < 10\%$ (ref. 21), lower than the value for coeval star-forming galaxies²⁹. Combining this upper limit to the gas fraction with the $\dot{M}_{\text{out}}^{\text{NaI}}$ (and assuming negligible mass contribution from the ionized phase), we estimated an upper limit for the depletion time of 200 Myr, which is very short compared to galaxy evolution timescales.

Figure 5 shows the stellar kinematics of GS-10578. We measured a clear velocity gradient between the east and west sides (Fig. 5e,f). Could this gradient be due to a major merger? In favour of this hypothesis, the two sides of GS-10578 have different continuum shapes, with the south-east side (Fig. 5e) displaying a flatter continuum slope and broader, less deep absorption. The rest-frame UV photometry shows several peaks, and the stellar velocity dispersion $\sigma_*$ is highest in the south-east rather than at the centre of the galaxy. However, contrary to the merger expectation, both UV peaks are in the same (receding) region of the velocity map (cyan contours in Fig. 5d). Therefore, these two peaks are unlikely to trace two distinct galaxies merging. If there is a merger, $\sigma_*$ would be highest between the two galaxies, where absorption lines of different redshift blend. Moreover, the mass maps and the red-wavelength NIRCam imaging show an elliptical morphology centred at—or very near to—the kinematic centre of the maps, which is also in contrast with the merger hypothesis (Fig. 2a and orange contours in Fig. 5d; see also Fig. 1c). Finally, the putative rotation axis is aligned with the minor axis of the isophotes (and isodensity contours), as expected from regular rotation. Under the rotation hypothesis, the different spectral shapes between the receding and approaching halves of the velocity map could be due to a combination of more dust in the south-east (where we also detected the neutral-gas outflow, as shown in Figs. 1e and 5c), a slightly younger stellar population in the north-west (Fig. 2b) or a stronger absorption-line infill due to nebular emission in the south-east. Most probably, all three factors contribute to the observed spectral differences.

By applying the virial theorem to the aperture velocity dispersion inside $1R_e$, $\sigma_{*,e} = 356 \pm 36 \text{ km s}^{-1}$. Thus, we found a dynamical mass $M_{\text{dyn}} = 2.0 \pm 0.5 \times 10^{11} M_{\odot}$ (using the calibration of ref. 30), which is comparable to $M_*$. In a major merger, both $\sigma_{*,e}$ and $R_e$ would be severely overestimated, meaning the true $M_{\text{dyn}}$ would be much lower than what is inferred from the virial theorem, in tension with $M_*$. Based on all this evidence, we classify GS-10578 as a regular (fast) rotator, the most distant yet found from stellar kinematics. An alternative possibility is that GS-10578 is a system in an advanced stage of merger, with the orbital angular momentum of the two progenitors locked in the stars³¹. To measure the degree of rotation support, we use the $\lambda_{R_e}$ spin parameter, a proxy for the angular momentum per unit mass³², defined by

$$ \lambda_{R_e} \equiv \frac{\sum_{r<R_e} F_i r_i |v_i|}{\sum_{r<R_e} F_i r_i \sqrt{v_i^2 + \sigma_i^2}} \qquad (1) $$

where we calculate the radius $r_i$ of each spatial element along the best-fitting ellipse³³. The aperture is the observed isophote of semi-major axis $1R_e$, $F_i$ is the continuum flux, and $v_i$ and $\sigma_i$ are the stellar
<figure><img src="image_8.png" /><figcaption>**Fig. 4** | Flux and kinematic maps of the [O III] λλ4949,5007 (Hα and [N II] λλ6548,6583 nebular lines). **a–c.** Flux maps obtained by fitting one or two Gaussians to the emission lines (Methods) for [O III] λλ4949,5007 (a), Hα (b) and [N II] λλ6548,6583 (c). **d–f.** Non-parametric median velocities for [O III] λλ4949,5007 (d), Hα (e) and [N II] λλ6548,6583 (f). **g–i.** Interpercentile line widths for [O III] λλ4949,5007 (g), Hα (h) and [N II] λλ6548,6583 (i). **j–l.** 10th</figcaption></figure>

percentile velocity for [O III] λλ4949,5007 (j), Hα (k) and [N II] λλ6548,6583 (l). **m–o.** 90th percentile velocity for [O III] λλ4949,5007 (m), Hα (n) and [N II] λλ6548,6583 (o). For GS-10578 (whose continuum emission is highlighted by the red dashed ellipse), most of the emission is from a broad, blueshifted component. Together with information from the BPT diagram (Fig. 3), this is evidence for an AGN-powered outflow.

<figure><img src="image_9.png" /><figcaption>**Fig. 5** | Absorption-line kinematics of GS-10578. **a–d.** Voronoi-binned 2D maps, showing the empirical S/N (a), stellar velocity $v_*$ (d), instrument-deconvolved velocity dispersion $\sigma_*$ (b) and neutral-gas velocity $v_{\text{NaI}}$ (c). In these four panels the dashed ellipses trace the shape and $R_e$ of the best-fitting Sérsic model to the stellar continuum. **a** also shows the contours of the Voronoi bins (white dotted lines). **d** shows an enlarged region of the RGB image in Fig. 1a, for comparison. **d** shows</figcaption></figure>

clear evidence for ordered rotation. The rotation axis is along the minor axis of the stellar distribution. **e.** Spectrum of a Voronoi bin, illustrating the Hδ and Hγ stellar absorption features; the spatial location of this bin is highlighted by the dotted white square labelled (e) in **d.** **f.** Same as **e**, but for the Voronoi bin labelled (f) in **d**. The two Voronoi bins (e) and (f) were chosen to highlight the velocity offset that we attribute to stellar rotation. The line colours are the same as Fig. 1f–i.

velocity and velocity dispersion. We measured $\lambda_{R_e}$ (observed) = 0.46, then applied a model-based point-spread function (PSF) correction (dependent on the galaxy size and light profile$^{34}$) for a Gaussian PSF with full-width at half-maximum FWHM = 0.09 arcsec to obtain $\lambda_{R_e} = 0.62 \pm 0.07$ (the uncertainties include the uncertainties due to the correction).

Figure 6a compares GS-10578 to other quiescent galaxies. Based on its location on the $\lambda_{R_e}$ versus $\epsilon$ plane (where $\epsilon$ is the ellipticity of the isophote at $1R_e$), GS-10578 is clearly a fast-rotator galaxy$^{32}$, which is consistent with it being an oblate spheroid with mild anisotropy, dominated by rotation support. Fast rotators are rare among galaxies of similar mass in the local Universe$^{35,36}$ but are much more common at
<figure><img src="image_10.png" /></figure>

**Fig. 6** | Kinematics and structure of GS-10578. a, $\lambda_{R_e}$ versus $\epsilon$ diagram showing that, by interpreting the velocity field as stellar rotation, GS-10578 is supported by rotation so that it is a so-called fast-rotator galaxy (slow rotators fall inside or near the dashed box$^{100}$). It is like other quiescent galaxies at high redshift (MRG sample, diamonds)$^{19}$ and unlike local quiescent galaxies (SAMI sample, circles)$^{33}$, including the most massive local quiescent galaxies (MASSIVE sample, squares)$^{101}$. b, Mass versus size diagram, illustrating that GS-10578 is extremely compact. To join the $z=0$ distribution, our galaxy must grow in size by a factor of 10. Our $10$ uncertainties were obtained by perturbing the measurements and bootstrapping 100 times.

high redshifts$^{19}$. Figure 6b shows GS-10578 on the mass–size plane. Its compact size again makes it a strong outlier at $z=0$ and even possibly at $z=2$. Taken together, Fig. 6a,b implies that galaxies like GS-10578 must have undergone a significant dynamical and morphological evolution between $z=3$ and the present. This evolution has to decrease the angular momentum and increase the size while not exceeding the stellar mass of local massive galaxies. It has been postulated that this kind of growth can be achieved only through gas-poor minor mergers$^{37,38}$, because gas-rich mergers trigger starbursts and disk (re-) formation, which would increase the light-weighted $\lambda_{R_e}$ parameter$^{39,40}$. Our target is surrounded by from three to five faint, low-mass satellites, which may represent the fuel for the postulated size growth$^{14}$.

GS-10578 represents a unique opportunity to study how the most massive galaxies in the Universe became—and stayed—quiescent. Even though we cannot draw general conclusions from a single target, we show that AGN feedback is capable of powering neutral-gas outflows with high velocity and high mass loading, sufficient to interrupt star formation by removing its cold-gas fuel. By contrast, the ionized-gas outflows appear to have much lower mass loading and may be inef- fective at halting star formation, in agreement with studies of other AGNs at high redshift. The stark difference in mass loading between the neutral and ionized phases had previously been observed only in star-forming galaxies and quasars$^{41}$, where the outflows have increasing mass loading going from the ionized phase to the neutral and molecular phases$^{42}$. More recently, other quiescent and quenching systems have been found with similar properties to GS-10578 (refs. 43,44), suggesting that neutral-phase outflows are a key mechanism for efficiently quench- ing massive galaxies. How exactly these outflows are coupled with the AGN is not yet clear. A possibility is that the AGN is transitioning from a dust-obscured phase by blowing out its surrounding material. It is unclear how much cold gas an AGN like the one in GS-10578 may leave behind. The upper limit on the molecular-gas mass fraction $f<10\%$ is not sufficiently constraining. Current evidence from redshifts $z=1.5–2$ sug- gests that quiescent galaxies have low (3%) molecular-gas fractions$^{4,5}$, an order of magnitude lower than coeval star-forming galaxies. Any leftover gas may have low star-forming efficiency, perhaps due to sta- bilization against gravitational collapse or (AGN-induced) turbulence, as for local PSB galaxies$^{45}$ and AGN hosts$^{46}$. The ISM and circumgalactic medium around GS-10578 suggest that a combination of AGN feedback

and mergers may be providing the energy required to sustain this tur- bulence against dissipation. The presence of a stellar disk suggests that the feedback was ‘gentle’ enough to preserve the dynamical structure of GS-10578 throughout its quenching phase. Therefore, a dynamical transformation from a disk-like fast rotator to a spheroid-like slow rota- tor can happen after quenching. A census of the emission-line properties of high-redshift, massive, quiescent galaxies is essential for clarifying whether strong, ejective SMBH feedback is episodic or whether it is a key mechanism for driving and maintaining quiescence.

## Methods

### NIRSpec observations and data reduction

Our observations are part of the GA-NIFS survey, a large NIRSpec GTO survey targeting extended galaxies and AGNs at $z>3$ with the NIRSpec/ IFS mode$^{47}$ (Principal Investigators (PIs) S. Arribas and R. Maiolino). Data for GS-10578 were obtained from programme ID 1216 (PIN. Lützgendorf) and consist of two disperser and filter combinations. Observations with the prism/clear set-up cover the full wavelength range of NIR- Spec with spectral resolution $R=30–300$ (ref. 48), whereas the g235h/ fl170lp observations probe only $1.66–3.17 \ \mu m$ but with higher resolution $R=2,700$. Both observations had a JWST position angle of $-99^\circ$, result- ing in a NIRSpec position angle of $39.5^\circ$ (measured north to east). For prism/clear, we used the IRS2RAPID read-out pattern with 33 groups per integration and eight exposures, for a total of 1.1 h on source. For g235h/ fl170lp, we used the IRS2 read-out pattern with 25 groups per integration and eight exposures, for a total of 4.1 h on source. The exposures were dithered following the medium cycling pattern to improve the PSF sampling and to marginalize against detector defects or leakage and open shutters from the NIRSpec micro-shutter assembly.

For the data reduction, we followed the procedure described in ref. 49. Here we present a brief summary and highlight any differences. We reduced the data with the publicly available JWST pipeline, v.1.8.2, using static calibrations from the context file v.1014 (for g235h/fl170lp) and context file v.1068 (for prism/clear). The pipeline was patched to include a number of corrections and improvements, as described in ref. 49. Compared to the default reduction, we used an extra correc- tion for the 1/f noise. For g235h/fl170lp, flux calibration was done in postprocessing using data from the standard star TYC 4433-1800-1 (PID 1128, observation 9). To detect and flag outliers, we used a custom algorithm, like that in ref. 50. That algorithm was developed specifically for (spatially) undersampled Hubble Space Telescope (HST) images. However, NIRSpec/IFS data were undersampled only in the spatial direction. Therefore, we calculated the derivative of the count-rate maps only along the detector x axis, which is the approximate direction of the dispersion axis (the trace curvature can be neglected on scales of two detector pixels). The absolute value of the derivative was then normalized by the local flux (or by $\times3$ the noise, whichever was high- est), and we rejected pixels for which the normalized derivative was higher than the 95th percentile of the resulting distribution. For the cube-building stage, we used the drizzle algorithm to produce cubes with a scale of 0.05 arcsec spaxel$^{-1}$, for both the prism/clear and the g235h/fl170lp data cubes. For the PSF modelling (method III), we used a cube with 0.03 arcsec spaxel$^{-1}$.

We also use publicly available photometry for flux calibration ('SED modelling and SFH'). At wavelengths 0.4–1.6 $\mu m$, we used data from the Great Observatories Origins Deep Survey captured by the Advanced Camera for Surveys onboard HST$^{51,52}$ and HST/WFC3 IR data from CAN- DELS$^{53,54}$. At wavelengths 0.9–4.8 $\mu m$, we used JWST/NIRCam data from the JWST Advanced Deep Extragalactic Survey (JADES; PID 1180, PID J. Eisenstein$^{11,12}$) and from the JWST Extragalactic Medium-band Survey (JEMS; PID 1963, PIs C. C. Williams, S. Tacchella and M. Maseda$^{13}$).

### PSF determination

We present three independent measurements of the shape and size of the NIRSpec/IFS PSF. The first measurement used observations of
a standard star with the g235h/f170lp grating and filter combination. The second method used a serendipitous star inside the NIRSpec/IFS field of view of GS-10578. The third method used NIRCam imaging from JADES and JEMS (for which the PSF is well understood) to infer the NIRSpec/IFS PSF. Note that all these methods use reconstructed data cubes. Therefore, they depend on both the data as well as the cube-build algorithm.

**Method I: standard star observations.** We used the standard star TYC 4433-1800-1 (V$_{3}$ = 188° and position angle of 146.5°), reduced using the same pipeline and context file as the galaxy data, and modelled the PSF using the multi-Gaussian expansion algorithm MGEFIT (ref. 55). We divided the data cube into 100 wavelength intervals between 1.66 and 3.17 µm, and then we found the PSF position angle from the second moment of the light distribution and the best-fitting multi-Gaussian expansion model by minimizing $\chi^2$. The formal flux uncertainties for the wavelength-collapsed data cube were very small and probably dominated by systematic uncertainties. Our $\chi^2$ minimization, therefore, used flux weighting instead of error weighting. The resulting FWHM is shown in Supplementary Fig. 1a. Filled and empty stars are the model FWHM along the major and minor axes of the PSF, respectively. The solid lines are the running median (note the detector gap at 2.5 µm). In Supplementary Fig. 1b, the solid line is the position angle of the PSF with respect to the NIRSpec position angle. A value of 0° indicates that the PSF major axis has the same direction as the long axes of the slices (hereafter, along slices whereas 90° is across slices). We found the NIRSpec/IFS PSF to be preferentially aligned close to along slices. At short wavelengths, the PSF was larger than the NIRCam PSF (coloured dots), and there seemed to be a drop in the FWHM around 3 µm.

**Method II: serendipitous star.** For this method, we used the relatively bright star inside the NIRSpec field of view (Fig. 1a; F775W = 25 mag), modelled as a Gaussian. The prism/clear data cube for GS-10578 was summed over six wavelength elements, and each synthetic image was fitted independently. We then modelled the resulting FWHM as a function of wavelength:

$$ \text{FWHM}(\lambda) = (0.129 + 1.61\lambda [\mu\text{m}] \exp(-21.43/\lambda [\mu\text{m}])) \text{arcsec.} \quad (2) $$

This model is shown as a dotted-dashed grey line in Supplementary Fig. 1a. Note that this determination may suffer from edge effects, particularly at the longest wavelengths.

**Method III: matching NIRCam imaging.** For this method, we selected a NIRCam wideband or medium-band filter, integrated the NIRSpec/IFS prism/clear cube (0.03 arcsec spaxel$^{-1}$) under the filter transmission curve to create a synthetic image in that band and then proceeded to find the kernel convolution that best matched the available NIRCam image to the synthetic NIRSpec image. The NIRSpec/IFS PSF was then estimated as the convolution of the NIRCam PSF with the best-fitting kernel. The model consists of six free parameters. The kernel is assumed to be a pixel-integrated, bivariate Gaussian, with three free parameters: $\sigma$, axis ratio $q$ and position angle. In addition, we allowed a small pixel adjustment to align the NIRCam and NIRSpec images and a flux scaling factor. The model was optimized using the Markov chain Monte Carlo algorithm emcee (ref. 56). We used flat priors for the kernel parameters (0.003 < $\sigma$ < 0.5 arcsec, -90° < position angle < 90° and 0 < $q$ < 1) and informative, Gaussian priors for matching the NIRCam image position and flux (the pixel shifts were Gaussians with mean 0 and standard deviation three pixels, and the flux scaling factor was a Gaussian with a mean equal to the mean flux ratio and a standard deviation of 10%). The example posterior distribution in Supplementary Fig. 2 shows that, as for method I above, the PSF is not circular. Supplementary Fig. 1a shows the results from method III as coloured ellipses, with the filled and empty symbols denoting the PSF FWHM along the major and minor

axes of the PSF, respectively. The position angle is once again close to 0, so that the PSF tends to be aligned with the IFS slices (Supplementary Fig. 1b). This method also found a drop in the PSF FWHM between the NIRCam F210M and F277W filters. The origin of this drop is unclear, and we did not model it in the subsequent analysis. We tested that our results are consistent when modelling the convolution kernel as a superposition of two Gaussians with the same axis ratio and position angle. In addition to GS-10578, we also modelled GS-3073, a AGN-host galaxy at z = 5.55 (ref. 57). The results are shown with GS-3073-shaped markers, with filled and empty symbols having the same meaning as for GS-10578. Using both GS-10578 and GS-3073, we modelled the PSF FWHM along the major axis (along slices) as

$$ \text{FWHM}(\lambda) = (0.125 + 1.88\lambda [\mu\text{m}] \exp(-24.35/\lambda [\mu\text{m}])) \text{arcsec.} \quad (3) $$

(dotted grey line in Supplementary Fig. 1a) and along the minor axis (across slices) as

$$ \text{FWHM}(\lambda) = (0.094 + 0.20\lambda [\mu\text{m}] \exp(-12.48/\lambda [\mu\text{m}])) \text{arcsec} \quad (4) $$

(dashed grey line in Supplementary Fig. 1a).

None of the three methods adopted is ideal. For method I, we used observations with a different number of dither positions than GS-10578, which may affect the image and PSF reconstruction. The limited number of dither meant that we could reconstruct the image only with a scale of 0.05 arcsec spaxel$^{-1}$, which may have limited the accuracy of the PSF modelling. As the serendipitous star is close to the edge of the NIRSpec/IFS field of view, method II may, therefore, suffer from edge effects. For method III, the starting NIRCam image has a different PSF than the NIRSpec/IFS image. Therefore, there is no guarantee that a Gaussian (or two-Gaussian) kernel can model the resulting image adequately. However, that the three methods agree on the elongated shape of the PSF leads credence to our conclusion. Furthermore, in the meantime, the elongation of the PSF along the instrument slicers has been reported also by other authors$^{58,59}$. This feature was expected from the instrument design and had already been noted during ground-based evaluations by the engineering team (Böker, T., private communication).

When measuring the spin parameter $\lambda_{R_e}$, we applied an empirical correction to obtain the intrinsic value from the observed measurement$^{34}$. We assumed a PSF FWHM of 0.09 arcsec (at 2 µm). A larger FWHM value would increase the intrinsic $\lambda_R$ and make our results more significant.

## Redshift determination

One would think that spectroscopy with a spectral resolution of R = 2,700 ought to yield redshift uncertainties better than 10 km s$^{-1}$, that is, a fraction of a pixel for well-sampled spectra (as is the case for NIRSpec). However, our data consist of two kinds of spectral features, neither of which are conducive to a precise redshift determination. Nebular emission consists of broad lines and several narrow lines. Even when a precise redshift could be determined, it was unclear which narrow lines could be used to trace the redshift of the galaxy. For the stellar features, the physical association with the systemic redshift was unambiguous, but all the absorption features belonged to the Balmer series of hydrogen and were, therefore, subject to nebular-line infill. The line infill may consist of several components too.

To take into account the different ways the absorption lines can be filled, we used a Monte Carlo approach. We considered the central spectrum as the unweighted sum of the spectra within the elliptical aperture of semimajor axis equal to $R_e/2$, axis ratio $q = 0.75$ and position angle 111.5°. Spaxels partially inside the aperture were weighted according to the fractional overlap between the spaxel and the ellipse. We then created 1,000 random-noise realizations of this spectrum by adding white noise from the uncertainty spectrum.
This procedure yielded an ensemble of spectra that had a 40% lower signal-to-noise ratio (S/N) than the actual data but crucially preserved any systematics compared to, for example, bootstrapping based on a best-fitting model spectrum. We fitted each random realization with the same procedure used in 'Stellar kinematics', but we also used a random starting guess for the redshift, chosen from 3.04 < z < 3.09. The resulting redshift distribution was bimodal, with the strongest peak at z = 3.06404 and a secondary peak at z = 3.06568 (the two distributions contain 63% and 37% of the total probability). We, therefore, adopted a systemic redshift $z_{\text{sys}} = 3.06404$. The random uncertainties (from a Gaussian fit to the posterior distribution) were of the order of 0.0004, but this value is overestimated because we added random noise to the already noisy spectrum. However, the systematic uncertainty (defined as the distance between the strongest peak and the secondary peak) was much larger at 0.0016. The random and systematic redshift uncertainties correspond to velocity uncertainties of 31 and 120 km s⁻¹, respectively. Our final redshift determination is $z_{\text{sys}} = 3.06404 \pm 0.0004$ (rand) ± 0.0016 (systematic).

Spectroscopy in band 1 (g140m and g140h gratings) may remove the systematic uncertainty entirely by targeting the higher-order Balmer lines where the stellar absorption is much stronger than the nebular emission. However, even this region may present complex features (for example, ISM absorption with a high EW in the Ca II doublet, which would overlap other stellar features, such as stellar Ca II and He).

In addition to the main target, we identified five smaller galaxies in the NIRCam image (Fig. 1a; labelled GS-10578b-f). Of these, targets b-d were clearly detected in the NIRSpec high-resolution data through the bright [O III] λλ4949,5007 doublet. We fitted a simple Gaussian to the emission-line complex Hβ-[O III] and found redshifts $z = 3.0639 \pm 0.0001$, $z = 3.0618 \pm 0.0001$ and $z = 3.0641 \pm 0.0001$ for GS-10578b, c and d, respectively. All these values are consistent with their being gravitationally bound to GS-10578. These satellites have also been reported in other works⁵⁸, with the names AGN-B, Clump1 and Clump2 for GS-10578b, c and d, respectively.

## UVJ colour-colour diagram

A UVJ diagram is a powerful diagnostic tool of star-forming activity¹⁷. To measure rest-frame *U*–*V* and *V*–*J* colours, we considered the prism/clear spectrum extracted from an aperture with 1.15 arcsec radius. The spectrum was recalibrated in flux (with correction factors ranging from 0.86 at 1 µm to 1.29 at 2.0–2.5 µm) to match the integrated magnitude of the NIRCam wideband and medium-band data obtained by JADES¹¹,¹². We considered a 5% constant error added in quadrature to the formal observational uncertainties to account for the absolute flux-calibration uncertainties.

We calculated the *U*–*V* and *V*–*J* rest-frame AB colours from the NIRSpec prism spectrum (Supplementary Fig. 3a). Our data do not fully cover the rest-frame/band, missing roughly 25% of the filter band-pass. We completed the missing data with two alternative methods, either by fitting a line to the reddest end of the NIRSpec spectrum (Extrap. 1 in Supplementary Fig. 3a) or by fitting a line to the photometry (Extrap. 2 in Supplementary Fig. 3a). These two methods give different *V*–*J* colours of 0.67 and 0.82, respectively. The observational uncertainties were negligible (<0.001 mag), but we added a 0.05 mag error in quadrature to reflect the uncertain calibration. Systematics in *V*–*J* are the largest source of uncertainty (0.15 mag). Regardless of which extrapolation was used, both methods placed GS-10578 in the quiescent region of the UVJ diagram (circles in Supplementary Fig. 3b). Our galaxy lies fairly close to the demarcation between quiescent and rapidly quenched galaxies.

## SED modelling and SFH

To study the physical properties of the stellar populations in GS-10578, the prism/clear NIRSpec/IFS observations were analysed in a spaxel-by-spaxel basis (0.05 arcsec on a side) following a similar

method to that described in ref. 60. Before comparing with models, the NIRSpec prism/clear data were PSF-matched using a Gaussian kernel to match the FWHM of the PSF at 4.5 µm, 0.16 arcsec FWHM (dotted line in Supplementary Fig. 1a), and were recalibrated in flux as outlined in 'UVJ colour-colour diagram'. As for the aperture spectrum, we added in quadrature a 5% constant error to the formal uncertainties and we analysed only spaxels inside a 1.15 arcsec aperture (9 kpc) centred on GS-10578. This aperture was from CANDELS⁶¹, and it captured the total flux of the source, including some of the nearest satellites.

The spectrum for each spaxel counting with more than ten wavelengths with S/N > 5 was compared to a grid of stellar-population models constructed with the library of ref. 62, assuming a SFH described by a delayed exponential characterized by a timescale $\tau$ (ranging from 1 Myr to 1 Gyr in 0.1 dex steps) and age $t_0$ (ranging from 1 Myr to the age of the Universe at the redshift of the galaxy, 2 Gyr). The stellar metallicity $Z$ was left as a free parameter, allowing us to take all the discrete values in the library of ref. 62 from 2% solar to 2.5 times solar. Nebular emission (continuum and line) was taken into account as described in refs. 63,64. The attenuation of the stellar and nebular emission was modelled with Calzetti's law⁶⁵, with $A_\nu$ ranging from 0 to 5 mag in 0.1 mag steps. The (surviving, that is, without including remnants and gas re-injected by winds and supernova explosions into the ISM) stellar mass was obtained by scaling the (mass-normalized) stellar model to the spectrum (the stellar mass is not a directly fitted quantity).

The stellar-population synthesis method includes a Monte Carlo algorithm to estimate uncertainties and degeneracies of the modelling⁶⁴,⁶⁶. This consists in varying the data points according to their uncertainties, refitting the spectrum and analysing the clusters of different solutions identified in the multi-dimensional parameter space defined by $\tau$, $t_0$, $Z$ and $A_\nu$.

Based on this method, we constructed the stellar mass and mass-weighted age maps presented in Fig. 2a,b. The mass-weighted age was calculated by integrating the SFH multiplied by the time from 0 to the age $t_0$ of the stellar population.

The integrated SFH and stellar mass evolution as well as the global SED, all presented in Fig. 2, were constructed by adding all the spaxels included in the stellar-population synthesis analysis. To account for low surface brightness pixels that were not included in the modelling, we applied an aperture correction obtained by comparing the flux of the analysed spaxels to the integrated photometric aperture.

### Modelling the panchromatic SED with CIGALE.
The panchromatic SED of the target was built by considering the multi-wavelength UV-to-NIR photometry available in ref. 67 as well as the NIRCam photometry released by the JADES and JEMS collaborations¹¹⁻¹³. The MIR to far-infrared photometry was taken from ref. 68. The counterpart in the different catalogues was found by using a positional matching radius of 1 arcsec. The SED was then analysed using the galaxy SED-fitting tool CIGALE (Code Investigating Galaxy Emission⁶⁹). This code was to used perform a multi-component fitting to separate the AGN contribution from the host galaxy emission. We used the following models to reproduce the different emission components: (1) The stellar emission was modelled by adopting a delayed-$\tau$ (exponentially declining) SFH, the stellar-population models of ref. 62, a Chabrier initial mass function⁷⁰ and solar metallicity. Dust attenuation was taken into account through the modified version of the Calzetti curve⁶⁵ available within CIGALE. (2) The emission from dust heated by star formation was reproduced using the model library in ref. 71, where we assumed an AGN contribution equal to 0 to treat it separately. (3) The AGN emission component was reproduced using the models presented in ref. 72. We ran CIGALE v.2018.0.2, adopting the same input grid used by ref. 73. The resulting best-fitting model is shown in Supplementary Fig. 4. The presence of a dust-obscured AGN (already suggested by the NIR data; Fig. 2d) is clear at 24 µm (~6 µm in the rest frame). From this model, we found a bolo-metric luminosity $L_{\text{AGN}} = 2.6 \pm 0.9 \times 10^{45}$ erg s⁻¹, a factor of four fainter
than previously reported values²¹. This discrepancy can be explained by the combination of different models, as well as new and updated photometry. We compared our value of $L_{AGN}$ (which was driven by the MIR photometry) to the bolometric luminosity estimated from the X-ray detection. We used the absorption-corrected $L_X = 6.4 \times 10^{44}$ erg s⁻¹ (ref. 21) and an X-ray bolometric correction from the literature (ref. 74, their equation (3)), finding $K_X(L_X) = 26$. This gives an X-ray based $L_{AGN} \approx 1.7 \times 10^{46}$ erg s⁻¹, six times higher than the MIR-based $L_{AGN}$, which is acceptable given the large uncertainties in this conversion.

## Stellar kinematics

We next present spatially resolved stellar kinematics beyond redshift $z=2.5$, and 2D stellar kinematics beyond $z=0.8$. For these measurements, we used PPXF, which simultaneously models the stellar continuum and nebular lines with a nonlinear $\chi^2$ minimization$^{55,75,76}$, all while taking into account the instrument line-spread function. The stellar continuum was modelled as a (non-negative) linear superposition of simple stellar population (SSP) spectra. As input, we used a subset of the SSP library based on the Mesa Isochrones and Stellar Tracks (MIST) isochrones$^{77}$ and C3K model atmospheres$^{78}$. The SSP templates span a logarithmically spaced grid of ages and metallicities, covering the ranges 10–2,500 Myr with 0.2 dex sampling and [Z/H] = −2.00–0.25 dex with 0.25 dex sampling. The maximum age was set to 500 Myr older than the age of the Universe at the redshift of GS-10578, to avoid edge effects (changing this setting did not alter our conclusions). The nebular emission and absorption lines were modelled with several Gaussians. It is also common practice to use additive Legendre polynomials when fitting stellar kinematics$^{79}$, but we found that PPXF tended to fit a strongly negative and almost featureless polynomial, akin to an oversubtracted background as bright as the target, which—given our flux calibration—is unphysical. We, therefore, used a 15th-order multiplicative Legendre polynomial to model residual flux-calibration issues. For the line-spread function, we applied a uniform spectral resolution of $R=2,700$, using the nominal resolution curve for g235h$^{48}$.

We found that our results depended critically on the wavelength range adopted. This was due to the already mentioned degeneracy between Balmer emission and absorption ('Redshift determination'). For example, on restricting the fit to the range where only $H\delta$ and $H\gamma$ were present, $\sigma_*$ was 25% lower (in this case the degree of the multiplicative polynomial was only 3). It is common practice to assess this uncertainty by repeating the fit and masking the Balmer lines$^{19}$, but in our case, these lines were the only prominent stellar absorption features in the observed spectroscopic configuration. Therefore, as for the redshift, a conclusive test may require new observations at shorter wavelengths.

### Aperture-integrated kinematics.
We used two sets of aperture-integrated kinematics: inside an ellipse with semimajor axis equal to $1R_e$ (to calculate the second moment of the line-of-sight velocity distribution of the stars, $\sigma_{*,e}$) and inside elliptical annuli (for use in the BPT diagram, Fig. 3). In both cases, we used emission lines coupled in four sets called 'kinematic components', which share the same line-of-sight velocity and intrinsic velocity dispersion (before convolution with the NIRSpec g235h line-spread function). The first kinematic set consisted of the Balmer series of hydrogen $H\delta-H\alpha$ (using the intrinsic line ratios appropriate for case B recombination, electron density $n_e = 100$ cm⁻³ and temperature $T_e = 10,000$ K; note that the relative line fluxes were then modified with a dust attenuation model), [N I] $\lambda$5200 and the [O I] $\lambda\lambda$6300,6363 doublet. The second component consisted of the [O III] $\lambda\lambda$4949,5007, [N II] $\lambda\lambda$6548,6583 and [S II] $\lambda\lambda$6716,6731 doublets. The third component consisted of [O III] and [N II]. The fourth and last component was the resonant Na I $D_1$ and $D_2$ absorption features. The line ratios of the [O III], [O I] and [N II] doublets were fixed to the appropriate values from atomic physics. The line ratio of the [S II] doublet was constrained between 0.44 and 1.47 (ref. 80). The $D_2/D_1$ ratio for Na I was

constrained between 0.89 and 1.94 (ref. 81). For nebular emission only, we used the Calzetti attenuation law$^{65}$.

From the elliptical aperture at $1R_e$ we found $\sigma_{*,e} = 356 \pm 10$ km s⁻¹. The uncertainties were probably dominated by systematics, so we considered the fiducial uncertainties to be 10%. Using the virial theorem calibration of ref. 30 and the light profile measurements from CANDELS$^{53}$ F160W imaging$^{20}$, we found a dynamical mass $M_{dyn} = k(n) k(q)\sigma_{*,e}^2 R_e/G = (2.0 \pm 0.1) \times 10^{11}$ M☉. This was dominated by the scatter in the calibration, which was of order 25%. For a single-burst stellar population of age 1 Gyr, the surviving stellar mass fraction was 0.7 (assuming a Chabrier initial mass function and solar metallicity). Downscaling the stellar mass formed by this factor, we obtained a dynamical-to-stellar mass ratio of 70%.

We also used the $1R_e$ apertures to measure the Balmer EW indices $H\delta$ and $H\gamma$, using the Lick definition$^{82}$. For $H\delta$, the Lick bands lay partially outside the wavelength range of the g235h/f170lp data cube. We, therefore, measured the $H\delta$ absorption strength using the best-fitting model. The resulting values were $H\delta = 6.63 \pm 0.01$ and $8.06 \pm 0.08$ Å (rest frame). For the second value, we subtracted the nebular emission. The uncertainties were from bootstrapping the uncertainties on the data. Clearly, the true stellar EW must be higher than the nonsubtracted measurement, which itself lay well into the range that defines PSB galaxies ($H\delta > 4-5$ Å; for example, ref. 83). For $H\gamma$, we measured the index from the data and found, respectively, $H\gamma = 4.07 \pm 0.01$ and $6.83 \pm 0.09$ Å for the unsubtracted and subtracted spectra. Both values are again in the PSB regime.

### Spatially resolved stellar kinematics.
Measuring spatially resolved stellar kinematics with PPXF typically requires $S/N \ge 10-15$ Å⁻¹, unless specific constraints are applied$^{79}$. Our data have a continuum $S/N$ as low as $3-5$ Å⁻¹. Therefore, we needed to adopt a careful approach. We measured the spatially resolved stellar kinematics with two different set-ups: Voronoi-binned data (method I) and spaxel-by-spaxel measurements (method II). Briefly, Voronoi bins have the advantage of a higher $S/N$$^{84}$, but for data with low spatial resolution, they risk averaging over regions with different stellar populations and kinematics$^{79}$. Spaxel-by-spaxel fitting overcomes the low-$S/N$ limitations by restricting the set of input templates$^{79,85}$. The results in this article are based on method I. The other method is described only briefly to show that our findings do not depend on the procedure adopted.

### Method I: Voronoi-binned data.
We measured the $S/N$ per spaxel with a preliminary fit to each spaxel independently using the best-fitting template from the fit to the aperture spectrum as the only input template ('Aperture-integrated kinematics'). The $S/N$ was then measured from the standard deviation of the residuals between the best-fitting spaxel model and the data. Having determined the $S/N$ empirically, we created the Voronoi bins using VORBIN$^{84}$. We selected a target $S/N = 12$ Å⁻¹. The resulting bins are outlined in Fig. 5a and are colour-coded by their empirical $S/N$, which was estimated from the residuals of the PPXF fit (see below). Note that in the outer bins ($1<R<2R_e$), the empirical $S/N$ was lower than the target $S/N$. This outcome was, in fact, expected from the correlated noise, which becomes more relevant as the average size of the Voronoi bins increases outward from the galaxy centre.

To fit each spectrum, we used the same procedure as for the integrated spectrum, but the choice of SSP spectra was restricted only to the spectra that had a non-zero weight in the linear combination forming the best-fitting aperture spectrum.

Moreover, the selection of Gaussian templates representing nebular emission and absorption was also simplified compared to what we used for the aperture spectrum, reflecting the simpler kinematics of the spatially resolved spectra. The first kinematic component consisted of the $H\delta-H\alpha$ series, the [O III] $\lambda\lambda$4949,5007 doublet, [N I] $\lambda$5200, and the three doublets [O I] $\lambda\lambda$6300,6363, [N II] $\lambda\lambda$6548,6583
and [S II] λλ6716,6731. The second kinematic component consisted of [O III] λλ4949,5007 and [N II] λλ6548,6583, and the third kinematic component consisted of Na I absorption. All doublets were fixed or constrained, as described for the aperture kinematics. Example fits of Voronoi bins are shown in Fig. 5e,f.

**Method II: Spaxel-by-spaxel data.** To measure the spaxel-by-spaxel kinematics, we used a threshold S/N > 3 pixel⁻¹. In general, such a low value is considered suboptimal, but if the templates are carefully preselected, it is still possible to obtain accurate measurements⁷⁹. To optimize the selection of templates, we first created 11 annular-binned spectra between 0 and 2.5R<sub>e</sub>, spaced by 0.025 arcsec and with the same shape and position angle as the observed stellar continuum. We then de-rotated the spectrum in each spaxel according to the line-of-sight velocity measured in the relevant Voronoi bins, which avoided the spectral broadening due to rotation. The resulting bin spectra were fitted with the same procedure used for the Voronoi bins. From these fits, we took the weights of the SSP templates and the measured σ<sub>*</sub>. We then proceeded to fit the individual spaxels. Given a spaxel, we identified the annular bins intersected by the spaxel and then also considered two more bins (if they existed): the annulus adjacent to the innermost selected annulus and the annulus adjacent to the outermost selected annulus. The templates selected were simply the union of all non-zero templates from the best-fitting annular spectra (typically seven to nine templates were selected). The velocity dispersion was the average of the annular σ<sub>*</sub> values, weighted by the fractional overlap of the spaxel with the annuli. We fitted the spaxel spectrum using this restricted template selection and with σ<sub>*</sub> fixed to the average σ<sub>*</sub>. If the resulting fit had S/N < 5, the fit was concluded and we proceeded to another spaxel. If the resulting fit had S/N > 5, we repeated the fit leaving σ<sub>*</sub> free between ±30% from the input σ<sub>*</sub>. The uncertainties were measured again using bootstrapping. The resulting maps are shown in Fig. 5, where they are compared with the Voronoi-binned maps. The resulting velocity values were quantitatively different, but the results were qualitatively the same. σ<sub>*</sub> has higher values in the central regions, probably due to a combination of lower average S/N and smaller template set compared to the Voronoi bins. Using this method, we measured the spin parameter λ<sub>R<sub>e</sub></sub> = 0.40, which rose to 0.60 after the PSF correction. This latter value is consistent with the measurement from method I.

## Emission-line measurements

To measure the line-of-sight velocity and the velocity width of the nebular emission lines, we used the QUBESPEC algorithm. We simultaneously modelled [O III] λλ4949,5007, Hα and [N II] λλ6548,6583. Each emission line was modelled with a single- or a double-Gaussian profile, with the redshift and FWHM of each Gaussian component being the same between the individual emission lines. We fixed the flux ratio of the [O III] and [N II] doublets as described in 'Stellar kinematics'. We estimated the posterior probability distribution of both the single- and double-Gaussian models with the Markov chain Monte Carlo algorithm emcee (ref. 56). We selected the best model (either one or two Gaussians per emission line) based on the Bayesian information criterion⁸⁶. For each spaxel, we took a median spectrum within the PSF (±2 spaxels) to increase the S/N of each fitted spectrum, while maintaining the diffraction-limited spatial resolution.

For each fitted spaxel, we calculated v<sub>10</sub>, v<sub>50</sub>, v<sub>90</sub> and w<sub>80</sub> as the 10th, 50th and 90th velocity percentiles of the emission line and the velocity width of the line 80% of the flux (v<sub>90</sub> to v<sub>10</sub>). The systemic velocity of the galaxy (v<sub>0</sub>) was set to the systemic velocity of the stellar component ('Redshift determination'). Although each of the Gaussian components was fitted as having the same redshift and FWHM, the final kinematics could be different due to the different flux ratios of the Gaussian components in each emission line. The resulting flux and kinematic maps are presented in Fig. 4.

## Outflow measurements

To estimate the outflow rate of the ionized gas, we used the [O III] emission line. Because Hβ was seen only in absorption, we calculated the mass of the outflowing ionized gas from the [O III] line as follows (see, for example, refs. 87,88):

$$M_{\text{out}}^{[\text{O III}]} = 8 \times 10^7 \left( \frac{1}{10^{[\text{O/H}]} - [\text{O/H}]_{\odot}} \right) \left( \frac{L_{[\text{O III}]}}{10^{44} \text{ erg s}^{-1}} \right) \left( \frac{n_e}{500 \text{ cm}^{-3}} \right)^{-1} M_{\odot}, \quad (5)$$

where L<sub>[O III]</sub> is the line luminosity and n<sub>e</sub> is the mean electron density. We assumed that the electron temperature was 10,000 K, that all the oxygen was ionized to O<sup>2+</sup> and that the metallicity [O/H] of the outflowing material was equal to the solar value [O/H]<sub>☉</sub>. We further assumed a uniformly filled biconical outflow and used the equation from ref. 89 as

$$\dot{M}_{\text{out}} = 3v_{\text{out}} \frac{M_{\text{out}}}{R_{\text{out}}}, \quad (6)$$

where v<sub>out</sub> is the velocity of the outflows (defined as v<sub>out</sub> = v<sub>10</sub> = 1,200 km s⁻¹), M<sub>out</sub> is the outflowing mass and R<sub>out</sub> is the radius of the outflow. The latter was estimated in the w<sub>80</sub> map as the radius of the region where w<sub>80</sub> > 600 km s⁻¹ (refs. 88,90). This value corresponded to a radius R<sub>out</sub> = 3.57 kpc. The largest source of uncertainty in this calculation was n<sub>e</sub>, which can range from 500 to 10⁴ cm⁻³ (for example, refs. 91,92). Therefore, we estimated the outflow rate for the range of density values, resulting in an outflow rate in the ionized gas of 0.14–2.9 M<sub>☉</sub> yr⁻¹. This range is shown as the cyan error bars in Fig. 2c.

For the neutral-gas outflows, we used the methods of ref. 93 based on the Na I absorption. These authors extend the methods of ref. 94 to spatially resolved spectroscopy. We assumed a thin-shell, conical geometry orthogonal to the galaxy disk. Reference 95 remarks that assuming a conical geometry with uniform density would increase the outflow rate by a factor of 3, which is the main source of uncertainty. Our disk inclination is 66° from the axis ratio of the best-fitting Sérsic model²⁰. The outflow rate in the kth spaxel is

$$\dot{M}_{\text{out,k}}^{\text{NaI}} = 11.5 \left( \frac{R_{\text{out}}}{10 \text{ kpc}} \right) \left( C_f \frac{\Omega_k}{0.4} \right) \left( \frac{N_k}{10^{21} \text{ cm}^2} \right) \left( \frac{v_{\text{out,k}}}{200 \text{ km s}^{-1}} \right) M_{\odot} \text{yr}^{-1}, \quad (7)$$

where R<sub>out</sub> is the extent of the outflow, Ω<sub>k</sub> is the subtended solid angle, C<sub>f</sub> = 0.37 is the cloud covering fraction (estimated from ref. 94), N<sub>k</sub> is the gas column density and v<sub>out,k</sub> is the outflow physical velocity. The extent of the outflow R<sub>out</sub> = 2.7 kpc, which was estimated as the projected distance between the centre of the galaxy and the outermost detected Na I absorption, corrected for inclination. As noted by, for example, ref. 93, this size measurement relies on the presence of a sufficiently bright background illumination to detect foreground Na I absorption. Therefore, we can consider R<sub>out</sub> a lower limit to the true extent of the outflow. The solid angle subtended by each spaxel was calculated from the spherical shell geometry. N<sub>k</sub> was estimated from the EW of Na I using equation (1) from ref. 93. Our EW measurements were corrected for the contribution due to the stellar atmospheric absorption. The outflow velocity was taken from the PPXF fits and was the (absolute value) of the centroid Gaussian-line velocity, corrected for the outflow inclination by cos i.

From this formula, summing over all spaxels where we detect blueshifted Na I absorption, we obtained M<sub>out</sub><sup>NaI</sup> = 100 M<sub>☉</sub> yr⁻¹.

To estimate the escape velocity v<sub>esc</sub>, we need the gravitational potential of the galaxy, including its dark matter halo. Using the relation between stellar mass and halo mass to infer the mass of the dark matter halo is impractical. This is because, for the stellar mass of GS-10578, the slope of the relation is unfavourable for a precise inference. We adopted the alternative method of assuming that the gravitational potential of the galaxy is an isothermal sphere, truncated at a radius r<sub>max</sub>; (refs. 57,96,97). Using equation (4) from ref. 57 with r<sub>max</sub> = 100 kpc,
we obtained an escape velocity $v_{esc} = 1,200 \pm 200 \text{ km s}^{-1}$. The uncertainties were estimated with a simple Monte Carlo analysis, for which we randomly drew the various observables that go into the equation of $v_{esc}$ as follows. The aperture velocity dispersion $\sigma_{*,e}$ was drawn from a normal distribution $N(\mu, \sigma)$ with mean $\mu = 356$ and standard deviation $\sigma = 36 \text{ km s}^{-1}(10\%)$. $R_e$ was drawn from $N(0.14, 0.02)(\text{units of kiloparsec})$. The Sérsic index $n$ was drawn from a truncated normal distribution $T(\mu, \sigma; \mu_0, \mu_1)$ with mean $\mu = 1.96$ and $\sigma = 0.4(20\%)$ and truncated between $\mu_0 = 1$ and $\mu_1 = 4$. $r_{\text{out}}$ was drawn from $T(4.5, 1, 2.5, 6)(\text{units of spaxels})$. The galaxy intrinsic axis ratio $q_0$ was drawn from $T(0.2, 0.1, 0.1, 0.5)$, and the projected axis ratio $q$ was drawn from $T(0.44, 0.1, 0.3, 0.5)$, subject to $q \ge q_0$. $q$ and $q_0$ were used to infer the inclination $i$.

Most of the uncertainty on $v_{esc}$ was due to the outflow radius $R_{\text{out}}(60\%)$, followed by $\sigma_{*,e}(30\%)$ and $R_e(10\%)$. The other three variables account for a few per cent each. Comparing our $v_{esc}$ to the map of deprojected Na I velocities (Fig. 5c), we found an average ratio $|v_{\text{Na I}}|/v_{esc} = 0.8 \pm 0.2$. Given the uncertainties on the outflow geometry and escape velocity, we can conclude that at least part of the neutral-gas outflows is able to escape the galaxy’s gravitational potential. At the very least, the gas will be temporarily removed from the star-forming region.

Comparing the ionized and neutral outflow phases, we conclude that the ionized-gas outflow rate is negligible, being up to two orders of magnitude lower than the neutral-gas outflow rate. In contrast, the mass-loading factor of the neutral outflow is higher than the recent (Hα-based) SFR. Therefore, this outflow phase is sufficient to stop star formation. Similar events in the recent past of GS-10578 may have caused the galaxy to become quiescent (SFH in Fig. 2b).

## Data availability

The observational data that support the findings of this study are publicly available from https://archive.stsci.edu/.

## Code availability

This work made extensive use of the freely available operating system Debian (http://www.debian.org). We used the programming language Python (www.python.org), which is maintained and distributed by the Python Software Foundation¹⁰². We further acknowledge direct use of the software Astropy (https://pypi.org/project/astropy/)¹⁰³, DYNESTY (https://pypi.org/project/dynesty/)¹⁰⁴,¹⁰⁵, FSPS (https://github.com/cconroy20/fspsv)¹⁰⁶,¹⁰⁷, Matplotlib (https://pypi.org/project/matplotlib)¹⁰⁸, NumPy (https://pypi.org/project/numpy)¹⁰⁹, Prospector (https://github.com/bd-j/prospector)¹¹⁰, QFitsView (https://www.mpe.mpg.de/-ott/QFitsView/) and SciPy (https://pypi.org/project/scipy/)¹¹¹.

## References

1. Madau, P. & Dickinson, M. Cosmic star-formation history. Annu. Rev. Astron. Astrophys. 52, 415–486 (2014).

2. Thomas, D., Maraston, C., Schawinski, K., Sarzi, M. & Silk, J. Environment and self-regulation in galaxy formation. Mon. Not. R. Astron. Soc. 404, 1775–1789 (2010).

3. McDermid, R. M. et al. The ATLAS3D Project. XXX. Star formation histories and stellar population scaling relations of early-type galaxies. Mon. Not. R. Astron. Soc. 448, 3484–3513 (2015).

4. Whitaker, K. E. et al. Quenching of star formation from a lack of inflowing gas to galaxies. Nature 597, 485–488 (2021).

5. Williams, C. C. et al. ALMA measures rapidly depleted molecular gas reservoirs in massive quiescent galaxies at z~1.5. Astrophys. J. 908, 54 (2021).

6. Man, A. & Belli, S. Star formation quenching in massive galaxies. Nat. Astron. 2, 695–697 (2018).

7. Somerville, R. S. & Davé, R. Physical models of galaxy formation in a cosmological framework. Annu. Rev. Astron. Astrophys. 53, 51–113 (2015).

8. Barro, G. et al. CANDELS: the progenitors of compact quiescent galaxies at z~2. Astrophys. J. 765, 104 (2013).

9. Damjanov, I. et al. Red nuggets at z~1.5: compact passive galaxies and the formation of the Kormendy relation. Astrophys. J. 695, 101–115 (2009).

10. Dekel, A. & Burkert, A. Wet disc contraction to galactic blue nuggets and quenching to red nuggets. Mon. Not. R. Astron. Soc. 438, 1870–1879 (2014).

11. Eisenstein, D. J. et al. Overview of the JWST advanced deep extragalactic survey (JADES). Preprint at https://arxiv.org/abs/2306.02465 (2023).

12. Rieke, M. et al. JADES initial data release for the Hubble ultra deep field: revealing the faint infrared sky with deep JWST NIRCam imaging. Astrophys. J. Suppl. Ser. 269, 16 (2023).

13. Williams, C. C. et al. JEMS: a deep medium-band imaging survey in the Hubble ultra-deep field with JWST NIRCam and NIRISS. Astrophys. J. Suppl. Ser. 268, 64 (2023).

14. Suess, K. A. et al. Minor merger growth in action: JWST detects faint blue companions around massive quiescent galaxies at $0.5 < z < 3.0$. Astrophys. J. Lett. 956, L42 (2023).

15. Popesso, P. et al. The main sequence of star-forming galaxies across cosmic times. Mon. Not. R. Astron. Soc. 519, 1526–1544 (2023).

16. Pacifici, C. et al. The evolution of star formation histories of quiescent galaxies. Astrophys. J. 832, 79 (2016).

17. Williams, R. J., Quadri, R. F., Franx, M., van Dokkum, P. & Labbé, I. Detection of quiescent galaxies in a bicolor sequence from $z=0-2$. Astrophys. J. 691, 1879–1895 (2009).

18. Kauffmann, G. et al. The host galaxies of active galactic nuclei. Mon. Not. R. Astron. Soc. 346, 1055–1077 (2003).

19. Newman, A. B., Belli, S., Ellis, R. S. & Patel, S. G. Resolving quiescent galaxies at $z \ge 2$. II. Direct measures of rotational support. Astrophys. J. 862, 126 (2018).

20. van der Wel, A. et al. 3D-HST+CANDELS: the evolution of the galaxy size–mass distribution since $z=3$. Astrophys. J. 788, 28 (2014).

21. Circosta, C. et al. X-ray emission of $z>2.5$ active galactic nuclei can be obscured by their host galaxies. Astron. Astrophys. 623, A172 (2019).

22. Cresci, G. & Maiolino, R. Observing positive and negative AGN feedback. Nat. Astron. 2, 179–180 (2018).

23. Harrison, C. M. et al. AGN outflows and feedback twenty years on. Nat. Astron. 2, 198–205 (2018).

24. Förster Schreiber, N. M. et al. The KMOS3D survey: demographics and properties of galactic outflows at $z=0.6-2.7$. Astrophys. J. 875, 21 (2019).

25. Nelson, D. et al. First results from the TNG50 simulation: galactic outflows driven by supernovae and black hole feedback. Mon. Not. R. Astron. Soc. 490, 3234–3261 (2019).

26. Concas, A. et al. Being KLEVER at cosmic noon: ionized gas outflows are inconspicuous in low-mass star-forming galaxies but prominent in massive AGN hosts. Mon. Not. R. Astron. Soc. 513, 2535–2562 (2022).

27. Diamond-Stanic, A. M. et al. High-velocity outflows without AGN feedback: Eddington-limited star formation in compact massive galaxies. Astrophys. J. Lett. 755, L26 (2012).

28. Rupke, D. S. N. et al. A 100-kiloparsec wind feeding the circumgalactic medium of a massive compact galaxy. Nature 574, 643–646 (2019).

29. Tacconi, L. J., Genzel, R. & Sternberg, A. The evolution of the star-forming interstellar medium across cosmic time. Annu. Rev. Astron. Astrophys. 58, 157–203 (2020).

30. van der Wel, A. et al. The mass scale of high-redshift galaxies: virial mass estimates calibrated with stellar dynamical models from LEGA-C. Astrophys. J. 936, 9 (2022).
31. Bois, M. et al. Formation of slowly rotating early-type galaxies via major mergers: a resolution study. *Mon. Not. R. Astron. Soc.* **406**, 2405–2420 (2010).

32. Emsellem, E. et al. The SAURON project. IX. A kinematic classification for early-type galaxies. *Mon. Not. R. Astron. Soc.* **379**, 401–417 (2007).

33. Brough, S. et al. The SAMI galaxy survey: mass as the driver of the kinematic morphology–density relation in clusters. *Astrophys. J.* **844**, 59 (2017).

34. Harborne, K. E. et al. Recovering $\lambda_R$ and V/σ from seeing-dominated IFS data. *Mon. Not. R. Astron. Soc.* **497**, 2018–2038 (2020).

35. Trujillo, I. et al. Superdense massive galaxies in the nearby Universe. *Astrophys. J. Lett.* **692**, L118–L122 (2009).

36. Comerón, S. et al. The massive relic galaxy NGC 1277 is dark matter deficient. From dynamical models of integral-field stellar kinematics out to five effective radii. *Astron. Astrophys.* **675**, A143 (2023).

37. Naab, T., Johansson, P. H. & Ostriker, J. P. Minor mergers and the size evolution of elliptical galaxies. *Astrophys. J. Lett.* **699**, L178–L182 (2009).

38. Bezanson, R. et al. The relation between compact, quiescent high-redshift galaxies and massive nearby elliptical galaxies: evidence for hierarchical, inside-out growth. *Astrophys. J.* **697**, 1290–1298 (2009).

39. Lagos, Cd. P. et al. Quantifying the impact of mergers on the angular momentum of simulated galaxies. *Mon. Not. R. Astron. Soc.* **473**, 4956–4974 (2018).

40. Lagos, Cd. P. et al. The diverse nature and formation paths of slow rotator galaxies in the EAGLE simulations. *Mon. Not. R. Astron. Soc.* **509**, 4372–4391 (2022).

41. Veilleux, S., Maiolino, R., Bolatto, A. D. & Aalto, S. Cool outflows in galaxies and their implications. *Astron. Astrophys. Rev.* **28**, 2 (2020).

42. Fluetsch, A. et al. Properties of the multiphase outflows in local (ultra)luminous infrared galaxies. *Mon. Not. R. Astron. Soc.* **505**, 5753–5783 (2021).

43. Belli, S. et al. Star formation shut down by multiphase gas outflow in a galaxy at a redshift of 2.45. *Nature* **630**, 54–58 (2024).

44. Davies, R. L. et al. JWST reveals widespread AGN-driven neutral gas outflows in massive z~2 galaxies. *Mon. Not. R. Astron. Soc.* **528**, 4976–4992 (2023).

45. Smercina, A. et al. After the fall: resolving the molecular gas in post-starburst galaxies. *Astrophys. J.* **929**, 154 (2022).

46. Venturi, G. et al. MAGNUM survey: compact jets causing large turmoil in galaxies. Enhanced line widths perpendicular to radio jets as tracers of jet–ISM interaction. *Astron. Astrophys.* **648**, A17 (2021).

47. Böker, T. et al. The near-infrared spectrograph (NIRSpec) on the James Webb Space Telescope. III. Integral-field spectroscopy. *Astron. Astrophys.* **661**, A82 (2022).

48. Jakobsen, P. et al. The near-infrared spectrograph (NIRSpec) on the James Webb Space Telescope. I. Overview of the instrument and its capabilities. *Astron. Astrophys.* **661**, A80 (2022).

49. Perna, M. et al. GA-NIFS: The ultra-dense, interacting environment of a dual AGN at z~3.3 revealed by JWST/NIRSpec IFS. *Astron. Astrophys.* **679**, A89 (2023).

50. van Dokkum, P. G. Cosmic-ray rejection by Laplacian edge detection. *Publ. Astron. Soc. Pac.* **113**, 1420–1427 (2001).

51. Giavalisco, M. et al. The Great Observatories Origins Deep Survey: initial results from optical and near-infrared imaging. *Astrophys. J. Lett.* **600**, L93–L98 (2004).

52. Illingworth, G. et al. The Hubble legacy fields (HLF-GOODS-S) v1.5 data products: combining 2442 orbits of GOODS-S/CDF-S region ACS and WFC3/IR images. Preprint at https://arxiv.org/abs/1606.00841 (2016).

53. Grogin, N. A. et al. CANDELS: The Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey. *Astrophys. J. Suppl. Ser.* **197**, 35 (2011).

54. Whitaker, K. E. et al. The Hubble legacy field GOODS-S photometric catalog. *Astrophys. J. Suppl. Ser.* **244**, 16 (2019).

55. Cappellari, M. & Emsellem, E. Parametric recovery of line-of-sight velocity distributions from absorption-line spectra of galaxies via penalized likelihood. *Publ. Astron. Soc. Pac.* **116**, 138–147 (2004).

56. Foreman-Mackey, D., Hogg, D. W., Lang, D. & Goodman, J. emcee: the MCMC hammer. *Publ. Astron. Soc. Pac.* **125**, 306 (2013).

57. Übler, H. et al. GA-NIFS: a massive black hole in a low-metallicity AGN at z~5.55 revealed by JWST/NIRSpec IFS. *Astron. Astrophys.* **677**, A145 (2023).

58. Perna, M. et al. A surprisingly high number of dual active galactic nuclei in the early Universe. Preprint at https://arxiv.org/abs/2310.03067 (2023).

59. Parlanti, E. et al. GA-NIFS: early-stage feedback in a heavily obscured AGN at z=4.76. *Astron. Astrophys.* **684**, A24 (2023).

60. Pérez-González, P. G. et al. CEERS key paper. IV. A triality in the nature of HST-dark galaxies. *Astrophys. J. Lett.* **946**, L16 (2023).

61. Koekemoer, A. M. et al. CANDELS: The Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey—The Hubble Space Telescope observations, imaging data products, and mosaics. *Astrophys. J. Suppl. Ser.* **197**, 36 (2011).

62. Bruzual, G. & Charlot, S. Stellar population synthesis at the resolution of 2003. *Mon. Not. R. Astron. Soc.* **344**, 1000–1028 (2003).

63. Pérez-González, P. G. et al. Stellar populations in local star-forming galaxies. I. Data and modelling procedure. *Mon. Not. R. Astron. Soc.* **338**, 508–524 (2003).

64. Pérez-González, P. G. et al. The stellar mass assembly of galaxies from z=0 to z=4: analysis of a sample selected in the rest-frame near-infrared with Spitzer. *Astrophys. J.* **675**, 234–261 (2008).

65. Calzetti, D. et al. The dust content and opacity of actively star-forming galaxies. *Astrophys. J.* **533**, 682–695 (2000).

66. Domínguez Sánchez, H. et al. Pathways to quiescence: SHARDS view on the star formation histories of massive quiescent galaxies at 1.0<z<1.5. *Mon. Not. R. Astron. Soc.* **457**, 3743–3768 (2016).

67. Merlin, E. et al. The ASTRODEEP-GS43 catalogue: new photometry and redshifts for the CANDELS GOODS-South field. *Astron. Astrophys.* **649**, A22 (2021).

68. Shirley, R. et al. HELP: the Herschel extragalactic legacy project. *Mon. Not. R. Astron. Soc.* **507**, 129–155 (2021).

69. Boquien, M. et al. CIGALE: a Python code investigating galaxy emission. *Astron. Astrophys.* **622**, A103 (2019).

70. Chabrier, G. Galactic stellar and substellar initial mass function. *Publ. Astron. Soc. Pac.* **115**, 763–795 (2003).

71. Dale, D. A. et al. A two-parameter model for the infrared/ submillimeter/radio spectral energy distributions of galaxies and active galactic nuclei. *Astrophys. J.* **784**, 83 (2014).

72. Fritz, J., Franceschini, A. & Hatziminaoglou, E. Revisiting the infrared spectra of active galactic nuclei with a new torus emission model. *Mon. Not. R. Astron. Soc.* **366**, 767–786 (2006).

73. Circosta, C. et al. SUPER. I. Toward an unbiased study of ionized outflows in z~2 active galactic nuclei: survey overview and sample characterization. *Astron. Astrophys.* **620**, A82 (2018).

74. Duras, F. et al. Universal bolometric corrections for active galactic nuclei over seven luminosity decades. *Astron. Astrophys.* **636**, A73 (2020).

75. Cappellari, M. Improving the full spectrum fitting method: accurate convolution with Gauss–Hermite functions. *Mon. Not. R. Astron. Soc.* **466**, 798–811 (2017).

76. Cappellari, M. Full spectrum fitting with photometry in pXF: stellar population versus dynamical masses, non-parametric star formation history and metallicity for 3200 LEGA-C galaxies at redshift z~0.8. *Mon. Not. R. Astron. Soc.* **526**, 3273–3300 (2023).
77. Choi, J. et al. Mesa Isochrones and Stellar Tracks (MIST). I. Solar-scaled models. Astrophys. J. 823, 102 (2016).

78. Conroy, C. et al. Resolving the metallicity distribution of the stellar halo with the H3 survey. Astrophys. J. 887, 237 (2019).

79. van de Sande, J. et al. The SAMI galaxy survey: revisiting galaxy classification through high-order stellar kinematics. Astrophys. J. 835, 104 (2017).

80. Kewley, L. J. et al. Theoretical ISM pressure and electron density diagnostics for local and high-redshift galaxies. Astrophys. J. 880, 16 (2019).

81. Poznanski, D., Prochaska, J. X. & Bloom, J. S. An empirical relation between sodium absorption and dust extinction. Mon. Not. R. Astron. Soc. 426, 1465-1474 (2012).

82. Worthey, G. Comprehensive stellar population models and the disentanglement of age and metallicity effects. Astrophys. J. Suppl. Ser. 95, 107-149 (1994).

83. Wu, P.-F. et al. The colors and sizes of recently quenched galaxies: a result of compact starburst before quenching. Astrophys. J. 888, 77 (2020).

84. Cappellari, M. & Copin, Y. Adaptive spatial binning of integral-field spectroscopic data using Voronoi tessellations. Mon. Not. R. Astron. Soc. 342, 345-354 (2003).

85. Foster, C. et al. The MAGPI survey: science goals, design, observing strategy, early results and theoretical framework. Publ. Astron. Soc. Aust. 38, e031 (2021).

86. Schwarz, G. Estimating the dimension of a model. Ann. Stat. 6, 461-464 (1978).

87. Carniani, S. et al. Ionised outflows in z~2.4 quasar host galaxies. Astron. Astrophys. 580, A102 (2015).

88. Kakkad, D. et al. BASS. XXXI. Outflow scaling relations in low redshift X-ray AGN host galaxies with MUSE. Mon. Not. R. Astron. Soc. 511, 2105-2124 (2022).

89. Fiore, F. et al. AGN wind scaling relations and the co-evolution of black holes and galaxies. Astron. Astrophys. 601, A143 (2017).

90. Harrison, C. M. et al. The KMOS AGN Survey at High redshift (KASHz): the prevalence and drivers of ionized outflows in the host galaxies of X-ray AGN. Mon. Not. R. Astron. Soc. 456, 1195-1220 (2016).

91. Baron, D. & Netzer, H. Discovering AGN-driven winds through their infrared emission. II. Mass outflow rate and energetics. Mon. Not. R. Astron. Soc. 486, 4290-4303 (2019).

92. Davies, R. et al. Ionized outflows in local luminous AGN: what are the real densities and outflow rates? Mon. Not. R. Astron. Soc. 498, 4150-4177 (2020).

93. Cazzoli, S., Arribas, S., Maiolino, R. & Colina, L. Neutral gas outflows in nearby [U]LIRGs via optical NaD feature. Astron. Astrophys. 590, A125 (2016).

94. Rupke, D. S., Veilleux, S. & Sanders, D. B. Outflows in infrared-luminous starbursts at z<0.5. II. Analysis and discussion. Astrophys. J. Suppl. Ser. 160, 115-148 (2005).

95. Maiolino, R. et al. Evidence of strong quasar feedback in the early Universe. Mon. Not. R. Astron. Soc. 425, L66-L70 (2012).

96. Veilleux, S., Cecil, G. & Bland-Hawthorn, J. Galactic winds. Annu. Rev. Astron. Astrophys. 43, 769-826 (2005).

97. Arribas, S., Colina, L., Bellocchi, E., Maiolino, R. & Villar-Martín, M. Ionized gas outflows and global kinematics of low-z luminous star-forming galaxies. Astron. Astrophys. 568, A14 (2014).

98. Curti, M. et al. The KLEVER survey: spatially resolved metallicity maps and gradients in a sample of 1.2<z<2.5 lensed galaxies. Mon. Not. R. Astron. Soc. 492, 821-842 (2020).

99. Newman, A. B., Belli, S., Ellis, R. S. & Patel, S. G. Resolving quiescent galaxies at z≥2. I. Search for gravitationally lensed sources and characterization of their structure, stellar populations, and line emission. Astrophys. J. 862, 125 (2018).

100. Cappellari, M. Structure and kinematics of early-type galaxies from integral field spectroscopy. Annu. Rev. Astron. Astrophys. 54, 597-665 (2016).

101. Veale, M. et al. The MASSIVE survey. VII. The relationship of angular momentum, stellar mass and environment of early-type galaxies. Mon. Not. R. Astron. Soc. 471, 1428-1445 (2017).

102. van Rossum, G. Python Tutorial (Universe, 1995).

103. Astropy Collaboration. et al. Astropy: A community Python package for astronomy. Astron. Astrophys. 558, A33 (2013).

104. Speagle, J. S. DYNESTY: a dynamic nested sampling package for estimating Bayesian posteriors and evidences. Mon. Not. R. Astron. Soc. 493, 3132-3158 (2020).

105. Johnson, B. et al. dfm/python-fsps: v0.4.4. Zenodo https://doi.org/10.5281/zenodo.8230430 (2023).

106. Conroy, C., Gunn, J. E. & White, M. The propagation of uncertainties in stellar population synthesis modeling. I. The relevance of uncertain aspects of stellar evolution and the initial mass function to the derived physical properties of galaxies. Astrophys. J. 699, 486-506 (2009).

107. Conroy, C. & Gunn, J. E. FSPS: flexible stellar population synthesis. Astrophysics Source Code Library, record ascl:1010.043 (ASCL, 2010).

108. Hunter, J. D. Matplotlib: a 2D graphics environment. Comput. Sci. Eng. 9, 90-95 (2007).

109. Harris, C. R. et al. Array programming with NumPy. Nature 585, 357-362 (2020).

110. Johnson, B. D., Leja, J., Conroy, C. & Speagle, J. S. Stellar population inference with Prospector. Astrophys. J. Suppl. Ser. 254, 22 (2021).

111. Jones, E., Oliphant, T., Peterson, P. et al. SciPy: open source scientific tools for Python (2001); www.scipy.org/

Acknowledgements

F.D.E., T.J.L., R.M. and J.S. acknowledge support from the Science and Technology Facilities Council, the European Research Council (ERC; Advanced Grant No. 695671, QUENCH) and UK Research and Innovation (Frontier Research Grant RISEandFALL). R.M. is also supported by a research professorship from the Royal Society. F.D.E. is grateful to C. Conroy for sharing the high-resolution C3K/MIST library and to C. Ramos Almeida for improving the quality of the manuscript. P.G.P.-G. acknowledges support from the Spanish Ministry for Science and Innovation/State Research Agency (MICIN/AEI; Grant Nos PGC2018-093499-B-100 and PID2022-139567NB-100 under project MCIN/AEI/10.13039/501100011033) and the European Regional Development Fund. S.A., B.R.P. and M.P. acknowledge support from MICIN/AEI (Research Project PID2021-127718NB-100). M.P. is also supported by the Programa Atracción de Talento de la Comunidad de Madrid (Grant No. 2018-T2/TIC-11715). H.Ü. gratefully acknowledges support from the Isaac Newton Trust and the Kavli Foundation through a Newton-Kavli Junior Fellowship. A.J.B., J.C. and G.C.J. acknowledge funding from a FirstGalaxies Advanced Grant from the ERC under the European Union's Horizon 2020 research and innovation programme (Grant Agreement No. 789056). S. Carniani, E.P. and G.V. acknowledge support from Horizon Europe and the ERC (Starting Grant No. 101040227, WINGS). G.C. acknowledges support from the Italian National Institute for Astrophysics (Large Grant 2022, The metal circle: a new sharp view of the baryon cycle up to Cosmic Dawn with the latest generation IFU facilities). I.L. acknowledges support from AEI (Project No. PID2022-140483NB-C22 under 10.13039/501100011033), MCIN through the Recovery, Transformation and Resilience Plan from the Spanish State (Grant No. BDC 20221289) and NextGenerationEU from the European Union through the Recovery and Resilience Facility. B.E.R. acknowledges support from the NIRCam Science Team contract to the University of Arizona (Contract No. NAS5-02015).
We acknowledge use of the lux supercomputer at University of California, Santa Cruz, which is funded by the Major Research Instrumentation Program of the National Science Foundation (Grant No. AST 1828315).

## Author contributions

F.D.E. measured the stellar and gas kinematics, wrote the main text, created the figures, and measured the dynamics and outflow properties. P.G.P.-G. designed the observations and measured the spatially resolved stellar-population properties. R.M. designed the survey and observations and helped with the interpretation of the results. J.S. measured the spatially resolved gas kinematics. M.P. and H.Ü. reduced the data. C.C. measured the integrated SED properties including the AGN component. S.A. designed the survey and observations and contributed to the measurement and interpretation of the outflows. B.R. and S.T. provided the imaging data. C.J.W. helped with data reduction. In addition, all authors helped at various stages of the survey design and execution, read the manuscript, and provided comments and insight for the interpretation.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41550-024-02345-1.

**Correspondence and requests for materials** should be addressed to Francesco D'Eugenio.

**Peer review information** *Nature Astronomy* thanks the anonymous reviewers for their contribution to the peer review of this work.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© Crown 2024

Francesco D'Eugenio Ⓓ1,2,16 Ⓒ, Pablo G. Pérez-González Ⓓ3,16, Roberto Maiolino Ⓓ1,2,4, Jan Scholtz Ⓓ1,2, Michele Perna³, Chiara Circosta⁴,⁵, Hannah Übler Ⓓ1,2, Santiago Arribas Ⓓ3, Torsten Böker Ⓓ6, Andrew J. Bunker Ⓓ7, Stefano Carniani Ⓓ8, Stephane Charlot Ⓓ9, Jacopo Chevallard⁷, Giovanni Cresci Ⓓ10, Emma Curtis-Lake Ⓓ11, Gareth C. Jones Ⓓ7, Nimisha Kumari Ⓓ12, Isabella Lamperti³, Tobias J. Looser Ⓓ1,2, Eleonora Parlanti Ⓓ8, Hans-Walter Rix Ⓓ13, Brant Robertson Ⓓ14, Bruno Rodríguez Del Pino³, Sandro Tacchella Ⓓ1,2, Giacomo Venturi⁸ & Chris J. Willott¹⁵

<footnote n="1">¹Kavli Institute for Cosmology, University of Cambridge, Cambridge, UK. ²Cavendish Laboratory - Astrophysics Group, University of Cambridge, Cambridge, UK. ³Centro de Astrobiología (CAB), CSIC-INTA, Torrejón de Ardoz, Spain. ⁴Department of Physics and Astronomy, University College London, London, UK. ⁵European Space Agency (ESA), European Space Astronomy Centre (ESAC), Villanueva de la Cañada, Spain. ⁶European Space Agency, Baltimore, MD, USA. ⁷University of Oxford, Department of Physics, Oxford, UK. ⁸Scuola Normale Superiore, Pisa, Italy. ⁹Sorbonne Université, UPMC-CNRS, UMR7095, Institut d'Astrophysique de Paris, Paris, France. ¹⁰INAF - Osservatorio Astrofisico di Arcetri, Florence, Italy. ¹¹Centre for Astrophysics Research, Department of Physics, Astronomy and Mathematics, University of Hertfordshire, Hatfield, UK. ¹²AURA for European Space Agency, Space Telescope Science Institute, Baltimore, MD, USA. ¹³Max Planck Institute for Astronomy, Heidelberg, Germany. ¹⁴Department of Astronomy and Astrophysics University of California, Santa Cruz, CA, USA. ¹⁵NRC Herzberg, Victoria, British Columbia, Canada. ¹⁶These authors contributed equally: Francesco D'Eugenio, Pablo G. Pérez-González. e-mail: francesco.deugenio@gmail.com</footnote>