A massive quiescent galaxy at redshift 4.658

The extremely rapid assembly of the earliest galaxies during the first billion years of cosmic history is a major challenge for our understanding of galaxy formation physics1–5. The advent of the James Webb Space Telescope (JWST) has exacerbated this issue by confirming the existence of galaxies in substantial numbers as early as the first few hundred million years6–8. Perhaps even more surprisingly, in some galaxies, this initial highly efficient star formation rapidly shuts down, or quenches, giving rise to massive quiescent galaxies as little as 1.5 billion years after the Big Bang9,10. However, due to their faintness and red colour, it has proven extremely challenging to learn about these extreme quiescent galaxies, or to confirm whether any existed at earlier times. Here we report the spectroscopic confirmation of a massive quiescent galaxy, GS-9209, at redshift, z = 4.658, just 1.25 billion years after the Big Bang, using the JWST Near-Infrared Spectrograph (NIRSpec). From these data we infer a stellar mass of M* = 3.8 ± 0.2 × 1010 M⊙, which formed over a roughly 200 Myr period before this galaxy quenched its star-formation activity at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z={6.5}_{-0.5}^{+0.2}$$\end{document}z=6.5−0.5+0.2, when the Universe was approximately 800 Myr old. This galaxy is both a likely descendent of the highest-redshift submillimetre galaxies and quasars, and a likely progenitor for the dense, ancient cores of the most massive local galaxies.

The extremely rapid assembly of the earliest galaxies during the first billion years of cosmic history is a major challenge for our understanding of galaxy formation physics [1][2][3][4][5].The advent of JWST has exacerbated this issue by confirming the existence of galaxies in significant numbers as early as the first few hundred million years [6][7][8].Perhaps even more surprisingly, in some galaxies, this initial highly efficient star formation rapidly shuts down, or quenches, giving rise to massive quiescent galaxies as little as 1.5 billion years after the Big Bang [9,10], however, due to their faintness and red colour, it has proven extremely challenging to learn about these extreme quiescent galaxies, or to confirm whether any exist at earlier times.Here we report the spectroscopic confirmation of a massive quiescent galaxy, GS-9209, at redshift, z = 4.658, just 1.25 billion years after the Big Bang, using JWST NIRSpec.From these data we infer a stellar mass of M * = 3.8 ± 0.2 × 10 10 M ⊙ , which formed over a ≃ 200 Myr period before this galaxy quenched its star formation activity at z = 6.5 +0.2 −0.5 , when the Universe was ≃ 800 million years old.This galaxy is both a likely descendent of the highest-redshift submillimetre galaxies and quasars, and a likely progenitor for the dense, ancient cores of the most massive local galaxies.Rest-frame Wavelength / m The spectrum strongly resembles that of an A-type star, and is reminiscent of lower-redshift post-starburst galaxies [40][41][42], clearly indicating this galaxy experienced a significant, rapid drop in star-formation rate (SFR) within the past few hundred million years.The spectral region from λ = 2.6 − 4.0µm, containing Hβ and Hα, is shown at a larger scale in Figure 2.
During the past 5 years, several studies have identified GS-9209 as a candidate high-redshift massive quiescent galaxy [11,12], based on its blue colours at wavelengths, λ = 2 − 8 µm and non-detection at millimetre wavelengths [13].GS-9209 is also not detected in X-rays [14], at radio wavelengths [15], or at λ = 24 µm [16].The faint, red nature of the source (with magnitudes H AB = 24.7 and K AB = 23.6)means that near-infrared spectroscopy with ground-based instrumentation is prohibitively expensive.The JWST NIRSpec data, shown in Figure 1, reveal a full suite of extremely deep Balmer absorption features, with a Hδ equivalent width (EW), as measured by the Hδ A Lick index, of 7.9 ± 0.3 Å, comparable to the most extreme values observed in the local Universe [17].These spectral features strongly indicate this galaxy has undergone a sharp decline in star-formation rate (SFR) during the preceding few hundred Myr.
The spectrum exhibits only the merest suspicion of [O ii] 3727 Å and [O iii] 4959 Å, 5007 Å emission, and no apparent infilling of Hβ or any of the higherorder Balmer absorption lines.However, as can be seen in Figure 2, both Hα and [Nii] 6584 Å are clearly albeit weakly detected in emission, with Hα also exhibiting an obvious broad component.This broad component, along with  To measure the stellar population properties of GS-9209, we perform full spectrophotometric fitting using the Bagpipes code [18] (see Methods).Briefly, we first mask the wavelengths corresponding to [O ii], [O iii], narrow Hα and [N ii], due to likely AGN contributions.We then fit a 22-parameter model for the stellar, dust, nebular and AGN components, as well as spectrophotometric calibration, to the spectroscopic data in combination with multi-wavelength photometry.Throughout the paper we report only statistical uncertainties on fitted parameters.It should be noted however that systematic uncertainties in galaxy spectral energy distribution analyses can be significantly larger [19][20][21].We investigate the effect of our choice of star-formation history (SFH) model in the Methods section.
The resulting posterior median model is shown in black in Figs 1 and 2. We obtain a stellar mass of log 10 (M * /M ⊙ ) = 10.58±0.02,under the assumption of a Kroupa initial mass function (IMF) [22].We additionally recover a very low level of dust attenuation, with A V = 0.02±0.02.The SFR we measure averaged over the past 100 Myr is consistent with zero, with a very stringent upper bound, though this is largely a result of our chosen SFH parameterisation [19].We provide a detailed discussion of the SFR of GS-9209 in Methods.We find a formation redshift, z form = 6.9 ± 0.2 and a quenching redshift, z quench = 6.5 +0.2 −0.5 .The sample of massive z ≃ 8 galaxy candidates from JWST CEERS reported by [6] is also shown in the right panel, demonstrating that these candidates are plausible progenitors for GS-9209.The uncertainties shown on the red points are 1σ standard deviation values.
The SFH we recover is shown in Figure 3.We find that GS-9209 formed its stellar population largely during a ≃ 200 Myr period, from around 600 − 800 Myr after the Big Bang (z ≃ 7 − 8).We recover a mass-weighted mean formation time, t form = 0.76 ± 0.03 Gyr after the Big Bang, corresponding to a formation redshift, z form = 6.9 ± 0.2.This is the redshift at which GS-9209 would have had half its current stellar mass, approximately log 10 (M * /M ⊙ ) = 10.3.We find that GS-9209 quenched (which we define as the time at which its sSFR fell below 0.2 divided by the Hubble time, e.g., [23]) at time t quench = 0.83 +0.08 −0.06 Gyr after the Big Bang, corresponding to a quenching redshift, z quench = 6.5 +0.2 −0.5 .Our model predicts that the peak historical SFR for GS-9209 (at approximately z form ) was within the range SFR peak = 490 +680 −300 M ⊙ yr −1 .This is similar to the SFRs of bright submillimetre galaxies (SMGs).The number density of SMGs with SFR > 300 M ⊙ yr −1 at 5 < z < 6 has been estimated to be ≃ 3 × 10 −6 Mpc −3 [24].Extrapolation then suggests that the SMG number density at z ≃ 7 is ≃ 1 × 10 −6 Mpc −3 , which equates to ≃ 1 SMG at z ≃ 7 over the ≃ 400 square arcmin area from which GS-9209 and one other z > 4 quiescent galaxy were selected [12].This broadly consistent number density suggests it is entirely plausible that GS-9209 went through a SMG phase at z ≃ 7, shortly before quenching.
In panel b of Figure 3, we show the positions of the massive, high-redshift galaxy candidates recently reported by [6] in the first imaging release from the JWST CEERS survey.The positions of these galaxies are broadly consistent with the SFH of GS-9209 at z ≃ 7 − 8.It should however be noted that, as previously discussed, GS-9209 was selected as one of only two robustly identified z > 4 massive quiescent galaxies in an area roughly 10 times the size of the initial CEERS imaging area [12].It therefore seems unlikely that a large fraction of the candidates reported by [6] will evolve in a similar way to GS-9209 over the redshift interval from z ≃ 5 − 8.
From our Bagpipes full spectral fit, we measure an observed broad Hα flux of f Hα, broad = 1.26 ± 0.08 × 10 −17 = erg s −1 cm −2 and full width at half maximum (FWHM) of 10300 ± 700 km s −1 in the rest frame.This line width, whilst very broad, is consistent with rest-frame UV broad line widths measured for some z ≃ 6 quasars (e.g., [25,26]).
As visualised in Figure 2, we fit Gaussian components to the narrow Hα and [N ii] lines following subtraction of our best-fitting Bagpipes model (see Methods).We obtain a Hα narrow-line flux of 1.58 ± 0.10 × 10 −18 erg s −1 cm −2 and a [N ii] flux of 1.56 ± 0.10 × 10 −18 erg s −1 cm −2 , giving a line ratio of log 10 ([N ii]/Hα) = −0.01 ± 0.04.This line ratio is significantly higher than would be expected as a result of ongoing star formation, and is consistent with excitation due to an AGN or shocks resulting from galactic outflows [27].Such outflows are commonly observed in post-starburst galaxies at z ≳ 1 [28].We discuss what can be learned about the star-formation rate of GS-9209 from the observed Hα flux in Methods.
We estimate the black-hole mass for GS-9209, M • , from our combined Hα flux and broad-line width, using the relation presented in Equation 6of [29], obtaining log 10 (M • /M ⊙ ) = 8.7 ± 0.1.From our Bagpipes full spectral fit, we infer a stellar velocity dispersion, σ = 247 ± 16 km s −1 for GS-9209, after correcting for the intrinsic dispersion of our template set, as well as instrumental dispersion.Given this measurement, the relationship between velocity dispersion and black-hole mass presented by [30] predicts log 10 (M • /M ⊙ ) = 8.9 ± 0.1.
Given the broad agreement between these estimators, it seems reasonable to conclude that GS-9209 contains a supermassive black hole with a mass of approximately half a billion to a billion Solar masses.It is interesting to note that this is ≃ 4−5 times the black-hole mass that would be expected given the stellar mass of the galaxy, assuming this is equivalent to the bulge mass.This is consistent with the observed increase in the average black-hole to bulge mass ratio for massive galaxies from 0 < z < 2 [31].The large amount of historical AGN accretion implied by this significant black hole mass suggests that AGN feedback may have been responsible for quenching this galaxy [32].
GS-9209 is an extremely compact source, which is only marginally resolved in the highest-resolution available imaging data.We measure the size of GS-9209 using newly available JWST NIRCam F210M-band imaging, which has a FWHM of ≃ 0.07 ′′ (see Methods).Accounting for the AGN point-source .This is consistent with previous measurements by the CANDELS/3DHST team [33], and is ≃ 0.7 dex below the mean relationship between r e and stellar mass for quiescent galaxies at z ≃ 1 [33,34].This is interesting given that post-starburst galaxies z ≃ 1 are known to be more compact than is typical for the wider quiescent population [35].We calculate a stellar-mass surface density within r e of log 10 (Σ eff /M ⊙ kpc −2 ) = 11.15 ± 0.08, consistent with the densest stellar systems in the Universe [36].
We show the F210M data for GS-9209, along with our posterior-median model in Figure 4.
We estimate the dynamical mass using our size and velocity dispersion measurements (e.g., [28]), obtaining a value of log 10 (M dyn /M ⊙ ) = 10.3 ± 0.1.This is ≃ 0.3 dex lower than the stellar mass we measure.As GS-9209 is only marginally resolved, even in JWST imaging data, and due to the presence of the AGN component, it is plausible that our measured r e may be subject to systematic uncertainties.Furthermore, since the pixel scale of NIRSpec is 0.1 ′′ , our velocity dispersion measurement may not accurately represent the central velocity dispersion, leading to an underestimated dynamical mass.It should also be noted that the stellar mass we measure is strongly dependent on our assumed IMF.A final, intriguing possibility would be a high level of rotational support in GS-9209, as has been observed for quiescent galaxies at 2 < z < 3 [37].Unfortunately, the extremely compact nature of the source makes any attempt at resolved studies extremely challenging, even with the JWST NIRSpec integral field unit.Resolved kinematics for this galaxy would be a clear use case for the High Angular Resolution Monolithic Optical and Near-infrared Integral field spectrograph (HARMONI) planned for the Extremely Large Telescope (ELT).
GS-9209 demonstrates unambiguously that massive galaxy formation was already well underway within the first billion years of cosmic history and that the earliest onset of galaxy quenching was no later than ≃ 800 Myr after the Big Bang.Based on the properties we measure, GS-9209 seems likely to be associated with the most extreme galaxy populations currently known at z > 5, such as the highest-redshift submillimetre galaxies and quasars (e.g., [26,38,39]).GS-9209 and similar objects (e.g., [8]) are also likely progenitors for the dense, ancient cores of the most massive galaxies in the local Universe.

Spectroscopic data and reduction
The spectroscopic data shown in Figure 1 were taken on 16 th November 2022.The target was acquired directly via Wide Aperture Target Acquisition (WATA), meaning the object is extremely well centred.Spectroscopic observations were taken through the NIRSpec fixed slit (S200A1), which has a width of 0.2 ′′ .Data were taken using the G235M and G395M gratings, providing an average spectral resolution of R = 1000.With each grating, a total of 5 integrations were taken at different dither positions along the slit.The readout pattern used was NRSIRS2, with 30 and 20 groups per integration for the two gratings respectively, providing total integration times of 3 hours and 2 hours respectively.
We reduce our NIRSpec data with the JWST Science Calibration Pipeline v1.8.4,using version 1017 of the JWST calibration reference data.To improve the spectrophotometric calibration of our data, we also reduce observations of the A-type standard star 2MASS J18083474+6927286 [43], taken as part of JWST commissioning programme 1128 (PI: Lützgendorf) [44] using the same instrument modes.We compare the resulting stellar spectrum against a spectral model for this star from the CALSPEC library [45] to construct a calibration function, which we then apply to our observations of GS-9209.The resulting spectrophotometry is well matched with the available near-infrared photometric data, and the calibration polynomial we fit along with our Bagpipes model only results in further calibration changes at the ≃ 10 per cent level.We additionally find that the uncertainties output by the pipeline are only moderately underestimated, with the errorbar expansion term in our Bagpipes model resulting in an increase of 50 per cent to the pipeline-produced uncertainties, in agreement with other recent analyses 1 .

Photometric data reduction
The majority of our photometric data are taken directly from the CANDELS GOODS South catalogue [46].We supplement this with new JWST NIRCam photometric data taken as part of the Ultra Deep Field Medium-Band Survey [47] (Programme ID: 1963; PI: Williams).Data are available in the F182M, F210M, F430M, F460M and F480M bands.We reduce these data using the PRIMER Enhanced NIRCam Image-processing Library (PENCIL, e.g., [7]), a custom version of the JWST Science Calibration Pipeline (v1.8.0), and using version 1011 of the JWST calibration reference data.We measure photometric fluxes for GS-9209 in large, 1 ′′ -diameter apertures to ensure we measure the total flux in each band (the object is isolated, with no other sources within this radius, see Figure 4).We measure uncertainties as the standard deviation of flux values in the nearest 100 blank-sky apertures, masking out nearby objects (e.g., [48]).

Bagpipes full spectral fitting
We fit the available photometry in parallel with our new spectroscopic data using the Bagpipes code [18].Our model has a total of 22 free parameters, describing the stellar, dust, nebular and AGN components of the spectrum.A full list of these parameters, along with their associated priors, is given in Extended Data Table 1.We fit our model to the data using the Multi-Nest nested sampling algorithm [49][50][51].The full Bagpipes fit to our combined dataset, along with residuals, is shown in Extended Data Figure 1.Posterior percentiles for our fit to the data are given in Extended Data Table 2.
We use the 2016 updated version of the BC03 [52,53] stellar population models, using the MILES stellar spectral library [54] and updated stellar evolutionary tracks [55,56].We assume a double-power-law star-formation-history model (e.g., [18,19]).We allow the logarithm of the stellar metallicity, Z * to vary freely from log 10 (Z * /Z ⊙ ) = −2.45 to 0.55.These are the limits of the range spanned by the BC03 model grid relative to our adopted Solar metallicity value (Z ⊙ = 0.0142 [57]).
We mask out the narrow emission lines in our spectrum during our Bagpipes fitting due to likely AGN contributions, whereas Bagpipes is only capable of modelling emission lines from star-forming regions.We do however still include a nebular model in our Bagpipes fit to allow for the possibility of nebular continuum emission from star-forming regions.We assume a stellarbirth-cloud lifetime of 10 Myr, and vary the logarithm of the ionization parameter, U, from log 10 (U ) = −4 to −2.We also allow the logarithm of the gas-phase metallicity, Z g , to vary freely from log 10 (Z g /Z ⊙ ) = −2.45 to 0.55.Because our eventual fitted model only includes an extremely small amount of star formation within the last 10 Myr for GS-9209, this nebular component makes a negligible contribution to the fitted model spectrum.
We model attenuation of the above components by dust using the model of [58,59], which is parameterised as a power-law deviation from the Calzetti dust attenuation law [60], and also includes a Drude profile to model the 2175 Å bump.We allow the V −band attenuation, A V to vary from 0 − 4 magnitudes.We further assume that attenuation is multiplied by an additional factor for all stars with ages below 10 Myr, and resulting nebular emission.This factor is commonly assumed to be 2, however we allow this to vary from 1 to 5.
We allow redshift to vary, using a narrow Gaussian prior with a mean of 4.66 and standard deviation of 0.01.We additionally convolve the spectral model with a Gaussian kernel in velocity space, to account for velocity dispersion in our target galaxy.The width of this kernel is allowed to vary with a logarithmic prior across a range from 50 − 500 km s −1 .The resolution of our spectroscopic data is high enough that the total dispersion is dominated by stellar velocity dispersion within the target galaxy, which has a standard deviation of σ ≃ 250 km s −1 , compared with the average instrumental dispersion of σ ≃ 128 km s −1 .
Separately from the above components, we also include a model for AGN continuum, broad Hα and Hβ emission.Following [61], we model AGN continuum emission with a broken power law, with two spectral indices and a break at λ rest = 5000 Å in the rest frame.We vary the spectral index at λ rest < 5000 Å using a Gaussian prior with a mean value of α λ = −1.5 (α ν = −0.5)and standard deviation of 0.5.We also vary the spectral index at λ rest > 5000 Å using a Gaussian prior with a mean value of α λ = 0.5 (α ν = −2.5)and standard deviation of 0.5.We parameterise the normalisation of the AGN continuum component using f 5100 , the flux at rest-frame 5100 Å, which we allow to vary with a linear prior from 0 to 10 −19 erg s −1 cm −2 Å−1 .
We model broad Hα with a Gaussian component, varying the normalisation from 0 to 2.5 × 10 −17 erg s −1 cm −2 using a linear prior, and the velocity dispersion from 1000 − 5000 km s −1 in the rest frame using a logarithmic prior.We also include a broad Hβ component in the model, which has the same parameters as the broad Hα line, but with normalisation divided by the standard 2.86 ratio from Case B recombination theory.However, as shown in Figure 2, this Hβ model peaks at around the noise level in our spectrum, and the line is therefore plausible in not being obviously detected in the observed spectrum.
We include intergalactic medium (IGM) absorption using the model of [62].To allow for imperfect spectrophotometric calibration of our spectroscopic data, we also include a second-order Chebyshev polynomial (e.g., [63][64][65]), which the above components of our combined model are all divided by before being compared with our spectroscopic data.We finally fit an additional white noise term, which multiplies the spectroscopic uncertainties from the JWST pipeline by a factor, a, which we vary with a logarithmic prior from 1 − 10.

Investigation of alternative star-formation history models
The functional forms used to model galaxy star-formation histories are well known to significantly affect physical parameter inferences [19][20][21], with the degree of systematic uncertainty highly dependent on the physical parameter of interest, the type of data, and the galaxy being studied.In this section, we test the dependence of our inferred formation and quenching times for GS-9209 on the SFH model used.We re-run our Bagpipes full-spectral-fitting analysis, substituting the double-power-law SFH model described above, firstly for the continuity non-parametric model [20], and secondly for a simple top-hat (constant) SFH model.For the continuity model, we use 8 time bins, with bin edges at 0, 10, 100, 200, 400, 600, 800, 1000 and 1260 Myr before observation.For the top-hat model, we vary the time at which star-formation turned on with a uniform prior between the Big Bang and time of observation.We vary the time at which star formation then stopped with a uniform prior from the time at which star formation turned on to the time of observation.
The results of these alternative fitting runs are shown in Extended Data Figure 2.This figure shows two alternative versions of Figure 3, with the continuity non-parametric model results shown in panels a and b, and the top-hat model results shown in panels c and d.The SFH posteriors shown, whilst varying in their detailed shapes, are in good overall agreement with our original double-power-law model.For the double-power-law model, we recover t form = 0.76 ± 0.03 Gyr and t quench = 0.83 +0.08 −0.06 Gyr after the Big Bang.The values returned under the assumption of these other two models are consistent to within 1σ.For the continuity non-parametric model we recover t form = 0.74 +0.02 −0.03 Gyr and t quench = 0.86 +0.19 −0.01 Gyr.For the top-hat model we recover t form = 0.74 ± 0.02 Gyr and t quench = 0.91 +0.04 −0.06 Gyr.Both of these models also produce stronger constraints on the peak historical SFR of GS-9209 at a lower level than the double-power-law model, though still consistent within 1σ.We conclude that our key results are not strongly dependent on our choice of star-formation history model.

AGN contribution and fitting of narrow emission lines
From our Bagpipes full spectral fit, we recover an observed AGN continuum flux at rest-frame wavelength, λ rest = 5100 Å of f 5100 = 0.06 ± 0.01 × 10 −19 erg s −1 cm −2 Å−1 .This is approximately 7.5 per cent of the total observed flux from GS-9209 at λ = 2.9µm.We measure a power-law index for the AGN continuum emission of α λ = −0.5 ± 0.3 at λ rest < 5000 Å, and α λ = 0.4 ± 0.3 at λ rest > 5000 Å.The AGN contribution to the continuum flux from GS-9209 rises to ≃ 10 per cent at the blue end of our spectrum (λ = 1.7µm), and ≃ 20 per cent at the red end (λ = 5µm).Just above the Lyman break at λ ≃ 7000 Å, the AGN contribution is ≃ 35 per cent of the observed flux.
Given our measured f Hα, broad , which is more direct than our AGN continuum measurement, the average relation for local AGN presented by [29] predicts f 5100 to be ≃ 0.2 dex brighter than we measure.However, given the intrinsic scatter of 0.2 dex they report, our measured f 5100 is only 1σ below the mean relation.The extreme equivalent widths of the observed Balmer absorption features firmly disfavour stronger AGN continuum emission.
We fit the narrow Hα and [N ii] lines in our spectrum as follows.We first subtract from our observed spectrum the posterior median Bagpipes model from our full spectral fitting.We then simultaneously fit Gaussian components to both lines, assuming the same velocity width for both, which is allowed to vary.This process is visualised in Figure 2. We also show the broad Hβ line in our AGN model, for which we assume the same width as broad Hα, as well as Case B recombination.It can be seen that the broad Hβ line peaks at around the noise level in our spectrum, and is hence too weak to be clearly observed in our data.

The star-formation rate of GS-9209
In this section we discuss the available observational indicators for the SFR of GS-9209.The commonly applied sSFR threshold for defining quiescent galaxies is sSFR threshold = 0.2/t H , where t H is the age of the Universe [23].
In [66], the authors report that GS-9209 is undetected in ALMA band 6 data, with a flux of −0.05±0.16mJy/beam, from which they derive a 1σ upper limit on SFR of 41 M ⊙ yr −1 .They also perform a stacking experiment, with stacked ALMA band 6 data for a sample of 20 objects selected as 3 < z < 5 quiescent galaxies (including GS-9209) still yielding no detection, implying the average SFR for this sample is significantly below the individual-object detection limit.The extremely blue spectral shape of this object in the restframe red-optical to near-infrared (observed-frame 2−8 µm, see Extended Data Figure 1) is also consistent with no significant obscured star-forming or AGN component.Deeper ALMA data for this object would be of value for setting a more stringent direct upper bound on obscured star formation.
As discussed in the main text, the high [N ii]/Hα ratio in our observed spectrum strongly suggests that this line emission is powered by AGN activity or shocks.However, if we assume all the narrow Hα emission is driven by ongoing star formation, neglecting dust attenuation we obtain SFR = 1.9 ± 0.1 M ⊙ yr −1 [67], corresponding to log 10 (sSFR/yr −1 ) = −10.3± 0.1.Measurements of the average dust attenuation on Hα emission, A Hα , are not yet available at z ≃ 5, however, from 0 < z < 2, stellar mass is found to be the most important factor in predicting the level of dust attenuation (e.g., [68,69]), with little evolution observed across this redshift interval.At z ≃ 2.3, the average A Hα for galaxies with log 10 (M * /M ⊙ ) ≃ 10.6 is 1.25 magnitudes [69], which would suggest that the SFR of GS-9209 is roughly 6 M ⊙ yr −1 .However, given the multiple lines of evidence we uncover for a significant non-stellar component to the Hα line, combined with the fact that the extremely low stellar continuum A V implies the gas-phase attenuation is also low [70], it is likely that the sSFR of GS-9209 is considerably lower than the threshold normally applied for selecting quiescent galaxies.

Size measurement from F210M-band imaging
The CANDELS/3DHST team [33] measured an effective radius, r e = 0.029 ± 0.002 ′′ for GS-9209 in the HST F125W filter via Sérsic fitting, along with a Sérsic index, n = 6.0 ± 0.8.At z = 4.658, this corresponds to r e = 189 ± 13 parsecs.We update this measurement using the newly acquired JWST NIR-Cam F210M imaging data discussed above.We model the light distribution of GS-9209 using PetroFit [71].We fit these PetroFit models to our data using the MultiNest nested sampling algorithm [49][50][51].We use F210M in preference to the F182M band due to the smaller AGN contribution in F210M and the fact that it sits above the Balmer break, therefore being more representative of the stellar mass present rather than any ongoing star formation.
As our spectroscopic data contains strong evidence for an AGN, we fit both Sérsic and delta-function components simultaneously, convolved by an empirically estimated PSF, derived by stacking bright stars.In preliminary fitting, we find that the relative fluxes of these two components are entirely degenerate with the Sérsic parameters.We therefore predict the AGN contribution to the flux in this band based on our full-spectral-fitting result, obtaining a value of 8 ± 1 per cent.We then impose this as a Gaussian prior on the relative contributions from the Sérsic and delta function components.The 11 free parameters of our model are the overall flux normalisation, which we fit with a logarithmic prior, the effective radius, r e , Sérsic index, n, ellipticity and position angle of the Sérsic component, the x and y centroids of both components, the position angle of the point spread function, and the fraction of light in the delta-function component, which we fit with a Gaussian prior with a mean of 8 per cent and standard deviation of 1 per cent, based on our full spectral fitting result.
Deeper imaging data in the F200W or F277W bands (e.g., from the JWST Advanced Deep Extragalactic Survey; JADES) will provide a useful check on our size measurement for GS-9209, particularly given the lower AGN fraction in the F277W band.

Data Availability
The datasets analysed during the current study are available from the Mikulski Archive for Space Telescopes (MAST) repository at https://mast.stsci.edu.The spectrum for GS-9209 was observed under JWST Programme ID 2285 (PI: Carnall).This programme has a 12-month proprietary period, and data will automatically become publicly available via MAST on 16 th November 2023.Reduced data products are available from the corresponding author upon request.
Extended Data Table 2: Posterior percentiles for the 22 free parameters of the Bagpipes model we fit to our spectroscopic and photometric data.Full definitions of these parameters, along with their associated prior distributions, are given in Extended Data Table 1

Redshift
Extended Data Figure 2: The star-formation rate and stellar mass of GS-9209 as a function of time under the assumption of different models.Panels a and b show the star-formation rate and stellar mass respectively as a function of time under the assumption of the continuity non-parametric SFH model [20] (as an alternative to the double-power-law model used in the main analysis and shown in Figure 3).Panels c and d again show the star-formation rate and stellar mass respectively as a function of time, this time assuming a top-hat SFH model.Consistent formation times are recovered using all three models.

Figure 1 :
Figure 1: JWST NIRSpec observations of GS-9209.Data were taken on 16 th November 2022, using the G235M and G395M gratings (R = 1000) with integration times of 3 hours and 2 hours respectively, providing wavelength coverage from λ = 1.7 − 5.1µm.The galaxy is at a redshift of z = 4.6582 ± 0.0002, and exhibits extremely deep Balmer absorption lines.The spectrum strongly resembles that of an A-type star, and is reminiscent of lower-redshift post-starburst galaxies[40][41][42], clearly indicating this galaxy experienced a significant, rapid drop in star-formation rate (SFR) within the past few hundred million years.The spectral region from λ = 2.6 − 4.0µm, containing Hβ and Hα, is shown at a larger scale in Figure2.

Figure 2 :
Figure 2: JWST NIRSpec observations of GS-9209: zoom in on Hβ and Hα.Data are shown in blue, with their associated (1σ standard deviation) uncertainties visible at the bottom in purple.The full Bagpipes fitted model is shown in black, with the AGN component shown in red.The narrow Hα and [N ii] lines were masked during the Bagpipes fitting process, and subsequently fitted with Gaussian functions, shown in green.Key emission and absorption features are also marked.

Figure 3 :
Figure3: The star-formation rate and stellar mass of GS-9209 as a function of time.Panel a shows the star-formation rate (SFR) as a function of time (the star-formation history).Panel b shows the stellar mass as a function of time.The blue lines show the posterior medians, with the darker and lighter shaded regions showing the 1σ and 2σ confidence intervals respectively.We find a formation redshift, z form = 6.9 ± 0.2 and a quenching redshift, z quench = 6.5 +0.2 −0.5 .The sample of massive z ≃ 8 galaxy candidates from JWST CEERS reported by[6] is also shown in the right panel, demonstrating that these candidates are plausible progenitors for GS-9209.The uncertainties shown on the red points are 1σ standard deviation values.

Figure 4 :
Figure 4: JWST NIRCam imaging of GS-9209.Each cutout image shows an area of 1.5 ′′ × 1.5 ′′ .Panel a is a RGB image, constructed with F430M as red, F210M as green and F182M as blue.Panel b shows the F210M image, with our posterior median PetroFit model shown in Panel c. Panel d shows the residuals between model and data, on the same colour scale as panels b and c.

Extended Data Figure 1 :
Model fit to our full combined dataset.Panel a shows our spectroscopic dataset in blue, with the full Bagpipes fitted model shown in black, and residuals between model and data shown below in green.Regions of our spectroscopic dataset masked during our Bagpipes fitting are shaded gray.Panel b shows our photometric data, again in blue, with the full model again shown in black, and residuals below marked with green crosses.