Main

During the past 5 years, several studies have identified GS-9209 as a candidate high-redshift massive quiescent galaxy11,12, on the basis of its blue colours at wavelengths λ = 2–8 μm and non-detection at millimetre wavelengths13. GS-9209 is also not detected in X-rays14, at radio wavelengths15 or at λ = 24 μm (ref. 16). The faint, red nature of the source (with H and K-band apparent magnitudes of HAB = 24.7 and KAB = 23.6) means that near-infrared spectroscopy with ground-based instrumentation is prohibitively expensive. The James Webb Space Telescope (JWST) Near-Infrared Spectrograph (NIRSpec) data, shown in Fig. 1, reveal a full suite of extremely deep Balmer absorption features, with a Hδ equivalent width, as measured by the HδA Lick index, of 7.9 ± 0.3 Å, comparable to the most extreme values observed in the local Universe17. These spectral features strongly indicate that this galaxy has undergone a sharp decline in star-formation rate (SFR) during the preceding few hundred million years.

Fig. 1: JWST NIRSpec observations of GS-9209.
figure 1

The figure shows flux per unit wavelength (fλ) as a function of wavelength, λ. Data were taken on 16 November 2022, using the G235M and G395M gratings (R = 1,000) 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 galaxies40,41,42, clearly indicating that this galaxy experienced a substantial, rapid drop in SFR in 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 Fig. 2.

Full size image

The spectrum exhibits only the merest suspicion of [O ii] 3,727 Å and [O iii] 4,959 Å, 5,007 Å emission, and no apparent infilling of Hβ or any of the higher-order Balmer absorption lines. However, as can be seen in Fig. 2, both Hα and [Nii] 6,584 Å are clearly, albeit weakly, detected in emission, with Hα also exhibiting an obvious broad component. This broad component, along with the relative strength of [N ii] compared with the narrow Hα line, indicates the presence of an accreting supermassive black hole: an active galactic nucleus (AGN). However, the extreme equivalent widths of the observed Balmer absorption features indicate that the continuum emission must be strongly dominated by the stellar component.

Fig. 2: JWST NIRSpec observations of GS-9209 with a zoom in on Hβ and Hα.
figure 2

Data are shown in blue, with their associated (1σ) 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.

Full size image

To measure the stellar population properties of GS-9209, we performed full spectrophotometric fitting using the Bayesian Analysis of Galaxies for Physical Inference and Parameter EStimation (Bagpipes) code18 (Methods). Briefly, we first masked the wavelengths corresponding to [O ii], [O iii], narrow Hα and [N ii], due to likely AGN contributions. We then fitted a 22-parameter model for the stellar, dust, nebular and AGN components, as well as spectrophotometric calibration, to the spectroscopic data in combination with multiwavelength 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 substantially larger19,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 obtained a stellar mass of log10(M*/M⊙) = 10.58 ± 0.02, under the assumption of a Kroupa initial mass function (IMF)22. We also recovered a very low level of dust attenuation, with a V-band attenuation in magnitudes of AV = 0.02 ± 0.02. The SFR we measured averaged over the past 100 Myr is consistent with zero, with a very stringent upper bound, although this is largely a result of our chosen SFH parameterization19. We provide a detailed discussion of the SFR of GS-9209 in Methods.

The SFH we recovered is shown in Fig. 3. We found that GS-9209 formed its stellar population largely during an approximately 200 Myr period, from around 600–800 Myr after the Big Bang (z ≈ 7–8). We recovered a mass-weighted mean formation time, tform = 0.76 ± 0.03 Gyr after the Big Bang, corresponding to a formation redshift of zform = 6.9 ± 0.2. This is the redshift at which GS-9209 would have had half its current stellar mass, approximately log10(M*/M☉) = 10.3. We find that GS-9209 quenched (which we define as the time at which its specific star-formation rate (sSFR) fell below 0.2 divided by the Hubble time23) at time \({t}_{{\rm{quench}}}={0\,.\,83}_{-0.06}^{+0.08}\) Gyr after the Big Bang, corresponding to a quenching redshift of \({z}_{{\rm{quench}}}=6\,.\,{5}_{-0.5}^{+0.2}\).

Fig. 3: The SFR and stellar mass of GS-9209 as a function of time.
figure 3

a, The SFR as a function of time (the SFH). b, 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, zform = 6.9 ± 0.2 and a quenching redshift, \({z}_{{\rm{quench}}}=6\,.\,{5}_{-0.5}^{+0.2}\). The sample of massive z ≈ 8 galaxy candidates from the JWST CEERS Survey reported by Labbe et al.6 is also shown in b, demonstrating that these candidates are plausible progenitors for GS-9209. The uncertainties shown on the red points are 1σ standard deviation values.

Full size image

Our model predicts that the peak historical SFR for GS-9209 (at approximately zform) was within the range SFRpeak \(=\,49{0}_{-300}^{+680}\) M☉ yr−1. This is similar to the SFRs of bright submillimetre galaxies (SMGs). The number density of SMGs with a SFR of more than 300 M☉ yr−1 at 5 < z < 6 has been estimated to be around 3 × 10−6 Mpc−3 (ref. 24). Extrapolation then suggests that the SMG number density at z ≈ 7 is approximately 1 × 10−6 Mpc−3, which equates to roughly 1 SMG at z ≈ 7 over the roughly 400-arcmin2 area from which GS-9209 and one other z > 4 quiescent galaxy were selected12. This broadly consistent number density suggests that it is entirely plausible that GS-9209 went through a SMG phase at z ≈ 7, shortly before quenching.

In Fig. 3b, we show the positions of the massive, high-redshift galaxy candidates recently reported by Labbe et al.6 in the first imaging release from the JWST Cosmic Evolution Early Release Science (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 Survey imaging area12. It therefore seems unlikely that a large fraction of the candidates reported by Labbe et al.6 will evolve in a way similar to that of GS-9209 over the redshift interval z ≈ 5–8.

From our Bagpipes full spectral fit, we measured an observed broad Hα flux of fHα,broad = 1.26 ± 0.08 × 10−17 = erg s−1 cm−2 and full width at half maximum (FWHM) of 10,300 ± 700 km s−1 in the rest frame. This linewidth, while very broad, is consistent with rest-frame ultraviolet broad linewidths measured for some z ≈ 6 quasars25,26.

As visualized in Fig. 2, we fitted Gaussian components to the narrow Hα and [N ii] lines following subtraction of our best-fitting Bagpipes model (Methods). We obtained 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 log10([N ii]/Hα) = −0.01 ± 0.04. This line ratio is substantially 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 outflows27. Such outflows are commonly observed in post-starburst galaxies at z ≳ 1 (ref. 28). We discuss what can be learned about the SFR of GS-9209 from the observed Hα flux in Methods.

We estimated the black-hole mass for GS-9209, M•, from our combined Hα flux and broad-line width, using the relation presented in equation 6 of Greene and Ho29, obtaining log10(M•/M⊙) = 8.7 ± 0.1. From our Bagpipes full spectral fit, we inferred a stellar velocity dispersion of σ = 247 ± 16 km s−1 for GS-9209, after correcting for the intrinsic dispersion of our template set and instrumental dispersion. Given this measurement, the relationship between velocity dispersion and black-hole mass presented by Kormendy and Ho30 predicts log10(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 (ref. 31). The large amount of historical AGN accretion implied by this substantial black-hole mass suggests that AGN feedback may have been responsible for quenching this galaxy32.

GS-9209 is an extremely compact source, which is only marginally resolved in the highest-resolution available imaging data. We measured the size of GS-9209 using newly available JWST Near-Infrared Camera (NIRCam) F210M-band imaging, which has a FWHM of around 0.07″ (Methods). Accounting for the AGN point-source contribution, we measured an effective radius, re = 0.033 ± 0.003″ for the stellar component of GS-9209, along with a Sérsic index, n = 2.3 ± 0.3. At z = 4.658, this corresponds to re = 215 ± 20 parsecs. This is consistent with previous measurements by the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS)/3D-HST team33, and is about 0.7 dex below the mean relationship between re and stellar mass for quiescent galaxies at z ≈ 1 (refs. 33,34). This is interesting given that post-starburst galaxies at z ≈ 1 are known to be more compact than is typical for the wider quiescent population35. We calculate a stellar-mass surface density within re of log10(Σeff/M☉ kpc−2) = 11.15 ± 0.08, consistent with the densest stellar systems in the Universe36. We show the F210M data for GS-9209, along with our posterior median model, in Fig. 4.

Fig. 4: JWST NIRCam imaging of GS-9209.
figure 4

Each cut-out image shows an area of 1.5″ × 1.5″. a, RGB image, constructed with F430M as red, F210M as green and F182M as blue. b,c, The F210M image (b), with our posterior median PetroFit model shown in c. d, The residuals between model and data, on the same colour scale as b and c.

Full size image

We estimated the dynamical mass using our size and velocity dispersion measurements28, obtaining a value of log10(Mdyn/M☉) = 10.3 ± 0.1. This is about 0.3-dex lower than the stellar mass we measure. As GS-9209 is only marginally resolved, even in JWST imaging data, and owing to the presence of the AGN component, it is plausible that our measured re may be subject to systematic uncertainties. Furthermore, because 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 (ref. 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 around 800 Myr after the Big Bang. On the basis of the properties we measured, GS-9209 seems likely to be associated with the most extreme galaxy populations now known at z > 5, such as the highest-redshift SMGs and quasars26,38,39. GS-9209 and similar objects8 are also likely progenitors for the dense, ancient cores of the most massive galaxies in the local Universe.

Methods

Spectroscopic data and reduction

The spectroscopic data shown in Fig. 1 were taken on 16 November 2022. The target was acquired directly by means of 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 = 1,000. With each grating, a total of five integrations were taken at different dither positions along the slit. The read-out 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 reduced our NIRSpec data with the JWST Science Calibration Pipeline v.1.8.4, using v.1017 of the JWST calibration reference data. To improve the spectrophotometric calibration of our data, we also reduced observations of the A-type standard star 2MASS J18083474+6927286 (ref. 43), taken as part of JWST commissioning programme 1128 (principal investigator: Lützgendorf)44 using the same instrument modes. We compared the resulting stellar spectrum against a spectral model for this star from the CALSPEC library45 to construct a calibration function, which we then applied to our observations of GS-9209. The resulting spectrophotometry is well matched with the available near-infrared photometric data, and the calibration polynomial we fitted along with our Bagpipes model results only in further calibration changes at roughly the 10% level. We also find that the uncertainties output by the pipeline are only moderately underestimated, with the error bar expansion term in our Bagpipes model resulting in an increase of 50% to the pipeline-produced uncertainties, in agreement with other recent analyses (for example, https://github.com/spacetelescope/jwst/issues/7362).

Photometric data reduction

Most of our photometric data were taken directly from the CANDELS GOODS South catalogue46. We supplemented these data with new JWST NIRCam photometric data taken as part of the Ultra Deep Field Medium-Band Survey47 (Programme ID: 1963; PI: Williams). Data are available in the F182M, F210M, F430M, F460M and F480M bands. We reduced these data using the PRIMER Enhanced NIRCam Image-processing Library (PENCIL)7, a custom version of the JWST Science Calibration Pipeline (v.1.8.0), and using v.1011 of the JWST calibration reference data. We measured photometric fluxes for GS-9209 in large, 1″-diameter apertures to ensure we measured the total flux in each band (the object is isolated, with no other sources within this radius; Fig. 4). We measured uncertainties as the standard deviation of flux values in the nearest 100 blank-sky apertures, masking out nearby objects48.

Bagpipes full spectral fitting

We fitted the available photometry in parallel with our new spectroscopic data using the Bagpipes code18. 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 fitted our model to the data using the MultiNest nested sampling algorithm49,50,51. The full Bagpipes fitted to our combined dataset, along with residuals, is shown in Extended Data Fig. 1. Posterior percentiles for our fit to the data are given in Extended Data Table 2.

We used the 2016 revised version of the BC03 (refs. 52,53) stellar population models, using the MILES stellar spectral library54 and revised stellar evolutionary tracks55,56. We assumed a double power-law SFH model18,19. We allowed the logarithm of the stellar metallicity, Z*, to vary freely from log10(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) (ref. 57).

We masked out the narrow emission lines in our spectrum during our Bagpipes fitting because of likely AGN contributions, whereas Bagpipes is capable of modelling emission lines only from star-forming regions. We did however still include a nebular model in our Bagpipes fit to allow for the possibility of nebular continuum emission from star-forming regions. We assumed a stellar birth cloud lifetime of 10 Myr, and varied the logarithm of the ionization parameter, U, from log10(U) = −4 to −2. We also allowed the logarithm of the gas-phase metallicity, Zg, to vary freely from log10(Zg/Z☉) = −2.45 to 0.55. Because our eventual fitted model includes only an extremely small amount of star formation in the past 10 Myr for GS-9209, this nebular component makes a negligible contribution to the fitted model spectrum.

We modelled attenuation of the above components by dust using the model of Noll et al.58and Salim et al.59, which is parameterized as a power-law deviation from the Calzetti dust attenuation law60, and also includes a Drude profile to model the 2,175-Å bump. We allowed the V-band attenuation, AV to vary from 0 to 4 magnitudes. We further assumed that attenuation is multiplied by an extra factor for all stars with ages below 10 Myr and resulting nebular emission. This factor is commonly assumed to be 2; however, we allowed it to vary from 1 to 5.

We allowed redshift to vary, using a narrow Gaussian prior with a mean of 4.66 and standard deviation of 0.01. We also convolved 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 of 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 included a model for AGN continuum, broad Hα and Hβ emission. Following Vanden Berk et al.61, we modelled AGN continuum emission with a broken power law, with two spectral indices and a break at λrest = 5,000 Å in the rest frame. We varied the spectral index at λrest < 5,000 Å using a Gaussian prior with a mean value of αλ = −1.5 (αν = −0.5) and standard deviation of 0.5. We also varied the spectral index at λrest > 5,000 Å using a Gaussian prior with a mean value of αλ = 0.5 (αν = −2.5) and standard deviation of 0.5. We parameterized the normalization of the AGN continuum component using f5100, the flux at rest-frame 5,100 Å, which we allowed to vary with a linear prior from 0 to 10−19 erg s−1 cm−2 Å−1.

We modelled broad Hα with a Gaussian component, varying the normalization from 0 to 2.5 × 10−17 erg s−1 cm−2 using a linear prior, and the velocity dispersion from 1,000 to 5,000 km s−1 in the rest frame using a logarithmic prior. We also included a broad Hβ component in the model, which has the same parameters as the broad Hα line, but with normalization divided by the standard 2.86 ratio from Case B recombination theory. However, as shown in Fig. 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 included intergalactic medium absorption using the model of Inoue et al.62. To allow for imperfect spectrophotometric calibration of our spectroscopic data, we also included a second-order Chebyshev polynomial63,64,65, which the above components of our combined model were all divided by before being compared with our spectroscopic data. We finally fitted an extra 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 to 10.

Investigation of alternative SFH models

The functional forms used to model galaxy SFHs are well known to substantially affect physical parameter inferences19,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, first for the continuity non-parametric model20, and second 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, 1,000 and 1,260 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 Fig. 2. This figure shows two alternative versions of Fig. 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, while varying in their detailed shapes, are in good overall agreement with our original double power-law model. For the double power-law model, we recovered tform = 0.76 ± 0.03 Gyr and \({t}_{{\rm{quench}}}=0.8{3}_{-0.06}^{+0.08}\) 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 recovered \({t}_{{\rm{form}}}=0.7{4}_{-0.03}^{+0.02}\) Gyr and \({t}_{{\rm{quench}}}=0.8{6}_{-0.01}^{+0.19}\) Gyr. For the top-hat model we recovered tform = 0.74 ± 0.02 Gyr and \({t}_{{\rm{quench}}}=0.9{1}_{-0.06}^{+0.04}\) 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, although still consistent within 1σ. We conclude that our key results are not strongly dependent on our choice of SFH model.

AGN contribution and fitting of narrow emission lines

From our Bagpipes full spectral fit, we recovered an observed AGN continuum flux at rest-frame wavelength λrest = 5,100 Å of f5100 = 0.06 ± 0.01 × 10−19 erg s−1 cm−2 Å−1. This is approximately 7.5% of the total observed flux from GS-9209 at λ = 2.9 μm. We measured a power-law index for the AGN continuum emission of αλ = −0.5 ± 0.3 at λrest < 5,000 Å and αλ = 0.4 ± 0.3 at λrest > 5,000 Å. The AGN contribution to the continuum flux from GS-9209 rises to around 10% at the blue end of our spectrum (λ = 1.7 μm), and around 20% at the red end (λ = 5 μm). Just above the Lyman break at λ ≈ 7,000 Å, the AGN contribution is about 35% of the observed flux.

Given our measured fHα,broad, which is more direct than our AGN continuum measurement, the average relation for local AGN presented by Greene and Ho29 predicts f5100 to be roughly 0.2 dex brighter than we measure. However, given the intrinsic scatter of 0.2 dex that they report, our measured f5100 is only 1σ below the mean relation. The extreme equivalent widths of the observed Balmer absorption features firmly disfavour stronger AGN continuum emission.

We fitted the narrow Hα and [N ii] lines in our spectrum as follows. We first subtracted from our observed spectrum the posterior median Bagpipes model from our full spectral fitting. We then simultaneously fitted Gaussian components to both lines, assuming the same velocity width for both, which was allowed to vary. This process is visualized in Fig. 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 SFR 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 sSFRthreshold = 0.2/tH, where tH is the age of the Universe23. For GS-9209 at z = 4.658 and log10(M*/M☉) ≈ 10.6, this corresponds to log10(sSFRthreshold/yr−1) ≈ −9.8, or SFRthreshold ≈ 6 M☉ yr−1.

In Santini et al.66, the authors report that GS-9209 is undetected in the Atacama Large Millimeter/submillimeter Array (ALMA) band-6 data, with a flux of −0.05 ± 0.16 mJy per 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 that the average SFR for this sample is well below the individual-object detection limit. The extremely blue spectral shape of this object in the rest-frame red-optical to near-infrared (observed frame 2–8 μm; Extended Data Fig. 1) is also consistent with no substantial 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 that all the narrow Hα emission is driven by continuing star formation, neglecting dust attenuation, we obtain SFR = 1.9 ± 0.1 M☉ yr−1(ref. 67), corresponding to log10(sSFR/yr−1) = −10.3 ± 0.1. Measurements of the average dust attenuation on Hα emission, AHα, 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 attenuation68,69, with little evolution observed across this redshift interval. At z ≈ 2.3, the average AHα for galaxies with log10(M*/M☉) ≈ 10.6 is 1.25 magnitudes69, which would suggest that the SFR of GS-9209 is roughly 6 M☉ yr−1. However, given the multiple lines of evidence we uncovered for a substantial non-stellar component to the Hα line, combined with the fact that the extremely low stellar continuum AV implies that the gas-phase attenuation is also low70, it is probable 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/3D-HST team33 measured an effective radius of re = 0.029 ± 0.002″ for GS-9209 in the HST F125W filter by means of Sérsic fitting, along with a Sérsic index, n = 6.0 ± 0.8. At z = 4.658, this corresponds to re = 189 ± 13 parsecs. We revised this measurement using the newly acquired JWST NIRCam F210M imaging data discussed above. We modelled the light distribution of GS-9209 using PetroFit71. We fitted these PetroFit models to our data using the MultiNest nested sampling algorithm49,50,51. We used 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, thereby being more representative of the stellar mass present rather than any continuing star formation.

As our spectroscopic data contain strong evidence for an AGN, we fitted both Sérsic and delta-function components simultaneously, convolved by an empirically estimated point spread function (PSF), derived by stacking bright stars. In preliminary fitting, we found that the relative fluxes of these two components are entirely degenerate with the Sérsic parameters. We therefore predicted the AGN contribution to the flux in this band on the basis of our full spectral fitting result, obtaining a value of 8% ± 1%. We then imposed 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 normalization, which we fitted with a logarithmic prior; the effective radius, re; 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 PSF; and the fraction of light in the delta-function component, which we fitted with a Gaussian prior with a mean of 8% and standard deviation of 1%, on the basis of our full spectral fitting result.

Deeper imaging data in the F200W or F277W bands (for example, from the JWST Advanced Deep Extragalactic Survey) will provide a useful check on our size measurement for GS-9209, particularly given the lower AGN fraction in the F277W band.