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

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^2,3. 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^4,5. 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^6,7. We present integral-field spectroscopy (IFS) observations made by the near-infrared spectrograph (NIRSpec) onboard JWST during the Cosmic Assembly Near-infrared Deep Extragalactic 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¹⁰ M_⊙, half-light radius R_e = 0.5–2 kpc) thought to be the progenitors of compact quiescent galaxies at z = 2 (‘red nuggets’⁹). 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.

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α–[N ii] and Hβ–[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 π/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-spaxel spectrum centred on the Hβ–[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α–[N ii]. h, Same as f, but centred on Hδ and Hγ. i, Same as f, but centred on the Na I doublet. The horizontal bar in the top right corner of each panel is δv = 500 km s⁻¹. f and g show the complex, multi-component emission-line profile of [O iii] and Hα [N ii]. h and i show the same spaxel. The vertical dotted lines mark the position of the stellar Hδ, Hγ 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.

Full size image

Imaging with the near-infrared camera (NIRCam)^11,12,13 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-spaxel 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. ¹⁵) and comparable to the quiescence threshold of 19 M_⊙ yr⁻¹ (ref. ¹⁶). 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. ¹⁷) and with the SFR inferred from the emission-line analysis.

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 1σ 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, 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 1σ 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.

Full size image

The NIRSpec/IFS observation in the low spectral-resolution configuration reveals extended nebular emission in both the Hβ–[O iii] and Hα–[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α from [N ii] and Hβ (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α 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 annulus around GS-10578, [O iii] λ5007/Hβ increases from 4 to 10, whereas [N ii] λ6584/Hα 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] λ6584/Hα ratios (empty circles), although our values are like those of post-starburst (PSB) galaxies at z = 1–2.5 (ref. ¹⁹) and to the most extreme PSB galaxies in the local Universe (grey star). Assuming all Hα emission was due to star formation, the extinction-corrected SFR on timescales of 3–10 Myr is 9 M_⊙ yr⁻¹ (inside 2R_e, where R_e = 1.1 kpc from ref. ²⁰). On varying the aperture size, the SFR ranged between 4 M_⊙ yr⁻¹ (inside 1R_e) to 19 M_⊙ yr⁻¹ (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_⊙ yr⁻¹, respectively for GS-10578b, c and d. The SFRs based on Hα are, overall, like the result from the SED analysis. However, as we have seen, the Hα 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α SFR values are in agreement with the quiescent interpretation. The BPT diagram shows both high [O iii] λ5007/Hβ and [N ii] λ6584/Hα, 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 × 10⁴⁴ erg s⁻¹)²¹ 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 ± 0.9 × 10⁴⁵ erg s⁻¹ (Methods)²¹.

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

Full size image

The spatially resolved emission-line kinematics (Fig. 4) displays line broadening in excess of 500 km s⁻¹ 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⁻¹, whereas the extended emission is at velocities of ±100 km s⁻¹. Furthermore, the central regions have interpercentile line widths w₈₀ in excess 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^22,23. Overall, the fast outflow velocity, the low SFR and the AGN are evidence for SMBH-feedback-driven outflows^24,25,26; but see refs. ^27,28 for cases of fast starburst-driven outflows. We estimate a mass outflow rate \({\dot{M}}_{{\rm{out}}}^{[{\rm{O}}\,{\textsc{iii}}]}\) = 0.14–2.9 M_⊙ yr⁻¹ (Methods), meaning the ionized-gas outflow is much lower than the SFR (Fig. 2c).

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

Full size image

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δ_A = 6.6–8.1 Å and Hγ_A = 4.1–6.8 Å; 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 × 10²² cm⁻², ten times lower than the column density towards the X-ray AGN (6 × 10²³ cm⁻²; ref. ²¹). 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_{Na i} = −(250–350) km s⁻¹), 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. ²⁰, which gives a projected axis ratio q = 0.44 (we assumed an intrinsic axis ratio q₀ = 0.2). After correcting for inclination by a factor \(\cos i\), we obtained an outflow velocity of the order of −1,000 km s⁻¹, again too fast to be explained by star formation alone. Assuming a maximal escape velocity of v_esc = 1,200 km s⁻¹, we inferred a velocity ratio v_{Na i}/v_esc = 0.8 ± 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_out = 2.7 kpc, based on the largest projected distance between the Na i detection and the centre of the galaxy. With this R_out, the mass outflow rate is \({\dot{M}}_{\rm{out}}^{{\rm{Na}}\,{\rm{I}}}\approx 30{-}100\,{{\rm{M}}}_{\odot}\,{{\rm{yr}}}^{-1}\,\), depending on the assumed geometry (Methods). Even with the lowest estimate, \({\dot{M}}_{{\rm{out}}}^{{\rm{Na}}\,{\rm{I}}}\) 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. ²¹), lower than the value for coeval star-forming galaxies²⁹. Combining this upper limit to the gas fraction with the \({\dot{M}}_{{\rm{out}}}^{{\rm{Na}}\,{\rm{I}}}\) (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 σ_⋆ 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, σ_⋆ 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.

a–d, Voronoi-binned 2D maps, showing the empirical S/N (a), stellar velocity v_⋆ (d), instrument-deconvolved velocity dispersion σ_⋆ (b) and neutral-gas velocity v_{Na i} (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 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.

Full size image

By applying the virial theorem to the aperture velocity dispersion inside 1R_e, σ_⋆,e = 356 ± 36 km s⁻¹. Thus, we found a dynamical mass M_dyn = 2.0 ± 0.5 × 10¹¹ M_⊙ (using the calibration of ref. ³⁰), which is comparable to M_⋆. In a major merger, both σ_⋆,e and R_e would be severely overestimated, meaning the true M_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}}}_{{\rm{e}}}}\) spin parameter, a proxy for the angular momentum per unit mass³², defined by

where we calculate the radius r_i of each spatial element along the best-fitting ellipse³³. The aperture is the observed isophote of semimajor axis 1R_e, F_i is the continuum flux, and v_i and σ_i are the stellar velocity and velocity dispersion. We measured \({\lambda }_{{{{R}}}_{{\rm{e}}}}\,\text{(observed)}=0.46\), then applied a model-based point-spread function (PSF) correction (dependent on the galaxy size and light profile³⁴) for a Gaussian PSF with full-width at half-maximum FWHM = 0.09 arcsec to obtain \({\lambda }_{{{{R}}}_{{\rm{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}}}_{{\rm{e}}}}\) versus ϵ plane (where ϵ is the ellipticity of the isophote at 1R_e), GS-10578 is clearly a fast-rotator galaxy³², 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 high redshifts¹⁹. 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}}}_{{\rm{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¹⁴.

a, \({\lambda }_{{{{R}}}_{{\rm{e}}}}\) versus ϵ 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¹⁰⁰). It is like other quiescent galaxies at high redshift (MRG sample, diamonds)¹⁹ and unlike local quiescent galaxies (SAMI sample, circles)³³, including the most massive local quiescent galaxies (MASSIVE sample, squares)¹⁰¹. 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 1σ uncertainties were obtained by perturbing the measurements and bootstrapping 100 times.

Full size image

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 ineffective 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⁴¹, where the outflows have increasing mass loading going from the ionized phase to the neutral and molecular phases⁴². 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 quenching 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 suggests 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 stabilization against gravitational collapse or (AGN-induced) turbulence, as for local PSB galaxies⁴⁵ and AGN hosts⁴⁶. 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 turbulence 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 rotator 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.

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⁴⁷ (Principal Investigators (PIs) S. Arribas and R. Maiolino). Data for GS-10578 were obtained from programme ID 1216 (PI N. Lützgendorf) and consist of two disperser and filter combinations. Observations with the prism/clear set-up cover the full wavelength range of NIRSpec with spectral resolution R = 30–300 (ref. ⁴⁸), whereas the g235h/f170lp observations probe only 1.66–3.17 μm but with higher resolution R = 2,700. Both observations had a JWST position angle of −99°, resulting in a NIRSpec position angle of 39.5° (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/f170lp, 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. ⁴⁹. 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/f170lp) and context file v.1068 (for prism/clear). The pipeline was patched to include a number of corrections and improvements, as described in ref. ⁴⁹. Compared to the default reduction, we used an extra correction for the 1/f noise. For g235h/f170lp, 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. ⁵⁰. 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 ×3 the noise, whichever was highest), 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⁻¹, for both the prism/clear and the g235h/f170lp data cubes. For the PSF modelling (method III), we used a cube with 0.03 arcsec spaxel⁻¹.

We also use publicly available photometry for flux calibration (‘SED modelling and SFH’). At wavelengths 0.4–1.6 μ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 CANDELS^53,54. At wavelengths 0.9–4.8 μm, we used JWST/NIRCam data from the JWST Advanced Deep Extragalactic Survey (JADES; PID 1180, PI D. J. Eisenstein^11,12) and from the JWST Extragalactic Medium-band Survey (JEMS; PID 1963, PIs C. C. Williams, S. Tacchella and M. Maseda¹³).

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₃ = 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. ⁵⁵). 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 χ². The formal flux uncertainties for the wavelength-collapsed data cube were very small and probably dominated by systematic uncertainties. Our χ² 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:

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⁻¹) 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: σ, 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.⁵⁶). We used flat priors for the kernel parameters (0.003 < σ < 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. ⁵⁷). 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

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

(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 dithers meant that we could reconstruct the image only with a scale of 0.05 arcsec spaxel⁻¹, 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 }_{{{\rm{R}}}_{{\rm{e}}}}\), we applied an empirical correction to obtain the intrinsic value from the observed measurement³⁴. We assumed a PSF FWHM of 0.09 arcsec (at 2 μm). A larger FWHM value would increase the intrinsic λ_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⁻¹, 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_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_sys = 3.06404 ± 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 Hϵ).

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 ± 0.0001, z = 3.0618 ± 0.0001 and z = 3.0641 ± 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^11,12. 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 J 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. ⁶⁰. 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. ⁶², assuming a SFH described by a delayed exponential characterized by a timescale τ (ranging from 1 Myr to 1 Gyr in 0.1 dex steps) and age t₀ (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. ⁶² 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_V 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^64,66. 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 τ, t₀, Z and A_V.

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₀ 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. ⁶⁷ as well as the NIRCam photometry released by the JADES and JEMS collaborations^11,12,13. The MIR to far-infrared photometry was taken from ref. ⁶⁸. 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-τ (exponentially declining) SFH, the stellar-population models of ref. ⁶², 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. ⁷¹, where we assumed an AGN contribution equal to 0 to treat it separately. (3) The AGN emission component was reproduced using the models presented by ref. ⁷². We ran CIGALE v.2018.0.2, adopting the same input grid used by ref. ⁷³. 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 bolometric luminosity L_AGN = 2.6 ± 0.9 × 10⁴⁵ 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 × 10⁴⁴ erg s⁻¹ (ref. ²¹) and an X-ray bolometric correction from the literature (ref. ⁷⁴, their equation (3)), finding K_X(L_X) = 26. This gives an X-ray based L_AGN ≈ 1.7 × 10⁴⁶ 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 χ² 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⁷⁷ and C3K model atmospheres⁷⁸. 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⁷⁹, 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⁴⁸.

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δ and Hγ were present, σ_⋆ 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¹⁹, 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, σ_⋆,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δ–Hα (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] λ5200 and the [O i] λλ6300,6363 doublet. The second component consisted of the [O iii] λλ4949,5007, [N ii] λλ6548,6583 and [S ii] λλ6716,6731 doublets. The third component consisted of [O iii] and [N ii]. The fourth and last component was the resonant Na i D₁ and D₂ 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. ⁸⁰). The D₂/D₁ ratio for Na i was constrained between 0.89 and 1.94 (ref. ⁸¹). For nebular emission only, we used the Calzetti attenuation law⁶⁵.

From the elliptical aperture at 1R_e we found σ_⋆,e = 356 ± 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. ³⁰ and the light profile measurements from CANDELS⁵³ F160W imaging²⁰, we found a dynamical mass M_dyn = k(n)k(q)σ_⋆,e²R_e/G = (2.0 ± 0.1) × 10¹¹ 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δ_A and Hγ_A, using the Lick definition⁸². For Hδ_A, the Lick bands lay partially outside the wavelength range of the g235h/f170lp data cube. We, therefore, measured the Hδ absorption strength using the best-fitting model. The resulting values were Hδ_A = 6.63 ± 0.01 and 8.06 ± 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δ_A > 4–5 Å; for example, ref. ⁸³). For Hγ_A, we measured the index from the data and found, respectively, Hγ_A = 4.07 ± 0.01 and 6.83 ± 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 ≳ 10–15 Å⁻¹, unless specific constraints are applied⁷⁹. 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⁸⁴, but for data with low spatial resolution, they risk averaging over regions with different stellar populations and kinematics⁷⁹. 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 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⁸⁴. 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δ–Hα series, the [O iii] λλ4949,5007 doublet, [N i] λ5200, and the three doublets [O i] λλ6300,6363, [N ii] λλ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_e, 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 σ_⋆. 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 σ_⋆ values, weighted by the fractional overlap of the spaxel with the annuli. We fitted the spaxel spectrum using this restricted template selection and with σ_⋆ fixed to the average σ_⋆. 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 σ_⋆ free between ±30% from the input σ_⋆. 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. σ_⋆ 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 \({\lambda }_{{{{R}}}_{{\rm{e}}}}=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. ⁵⁶). 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₁₀, v₅₀, v₉₀ and w₈₀ as the 10th, 50th and 90th velocity percentiles of the emission line and the velocity width of the line 80% of the flux (v₉₀ to v₁₀). The systemic velocity of the galaxy (v₀) 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):

where L_[O iii] is the line luminosity and n_e is the mean electron density. We assumed that the electron temperature was 10,000 K, that all the oxygen was ionized to O²⁺ and that the metallicity [O/H] of the outflowing material was equal to the solar value [O/H]_⊙. We further assumed a uniformly filled biconical outflow and used the equation from ref. ⁸⁹ as

where v_out is the velocity of the outflows (defined as v_out = v₁₀ = 1,200 km s⁻¹), M_out is the outflowing mass and R_out is the radius of the outflow. The latter was estimated in the w₈₀ map as the radius of the region where w₈₀ > 600 km s⁻¹ (refs, ^88,90). This value corresponded to a radius R_out = 3.57 kpc. The largest source of uncertainty in this calculation was n_e, 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_⊙ yr⁻¹. This range is shown as the cyan error bars in Fig. 2c.

For the neutral-gas outflows, we used the methods of ref. ⁹³ based on the Na i absorption. These authors extend the methods of ref. ⁹⁴ to spatially resolved spectroscopy. We assumed a thin-shell, conical geometry orthogonal to the galaxy disk. Reference ⁹⁵ 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

where R_out is the extent of the outflow, Ω_k is the subtended solid angle, C_f = 0.37 is the cloud covering fraction (estimated from ref. ⁹⁴), N_k is the gas column density and v_out,k is the outflow physical velocity. The extent of the outflow R_out = 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. ⁹³, this size measurement relies on the presence of a sufficiently bright background illumination to detect foreground Na i absorption. Therefore, we can consider R_out 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_k was estimated from the EW of Na i using equation (1) from ref. ⁹³. 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 \({\dot{M}}_{{\rm{out}}}^{{\rm{Na}}\,{\rm{I}}}=100\) M_⊙ yr⁻¹.

To estimate the escape velocity v_esc, 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}_{\max }\); (refs. ^57,96,97). Using equation (4) from ref. ⁵⁷ with \({r}_{\max }=100\) kpc, we obtained an escape velocity v_esc = 1,200 ± 200 km s⁻¹. 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 σ_⋆,e was drawn from a normal distribution \({\mathcal{N}}(\mu ,\sigma )\) with mean μ = 356 and standard deviation σ = 36 km s⁻¹ (10%). R_e was drawn from \({\mathcal{N}}(0.14,0.02)\) (units of kiloparsec). The Sérsic index n was drawn from a truncated normal distribution \({\mathcal{T}}(\mu ,\sigma ;{\mu }_{0},{\mu }_{1})\) with mean μ = 1.96 and σ = 0.4 (20%) and truncated between μ₀ = 1 and μ₁ = 4. r_out was drawn from \({\mathcal{T}}(4.5,1,2.5,6)\) (units of spaxels). The galaxy intrinsic axis ratio q₀ was drawn from \({\mathcal{T}}(0.2,0.1,0.1,0.5)\), and the projected axis ratio q was drawn from \({\mathcal{T}}(0.44,0.1,0.3,0.5)\), subject to q ≥ q₀. q and q₀ were used to infer the inclination i.

Most of the uncertainty on v_esc was due to the outflow radius R_out (60%), followed by σ_⋆,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_{Na i}∣/v_esc = 0.8 ± 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).