Age-dating early quiescent galaxies: high star formation efficiency, but consistent with direct, higher-redshift observations

ABSTRACT

We present a detailed analysis of JWST/NIRSpec and NIRCam observations of ZF-UDS-7329, a massive, quiescent galaxy at redshift |$z=3.2$|⁠, which has been put forward to challenge cosmology and galaxy formation physics. We study on the impact of different star formation history (SFH) priors, stellar libraries, metallicity, and initial mass function assumptions. Our results show that ZF-UDS-7329, with a formed stellar mass of |$M_{\star } \approx 10^{11.4}~{\rm M}_\odot $| (surviving mass |$M_{\star \mathrm{,surv}} \approx 10^{11.2}~{\rm M}_\odot $|⁠) and a specific star formation rate of |$\mathrm{sSFR} \approx 0.03~{\rm Gyr} ^{-1}$|⁠, formed efficiently in the first billion years of the Universe. In agreement with previous work, we find that the spectrum is consistent with mass-weighted stellar ages of |$1.3{\!-\!}1.8$| Gyr, depending on the SFH prior used. A physically motivated rising SFH prior makes the formation history of ZF-UDS-7329 compatible with stellar mass and star formation rate estimates of |$z\gt 6$| galaxies. Using NIRCam imaging, we identify a colour gradient indicative of an old, quiescent bulge and a younger disc component, as expected from a complex formation history. The inferred SFH is consistent a high stellar fraction of |$f_{\star }=M_{\star }/(f_b \cdot M_{\rm h}) \approx 100{{\ \rm per\ cent}}$| at |$z=7{\!-\!}12$|⁠, implying an extremely high integrated star formation efficiency. However, when considering cosmic variance and possible mergers as expected in overdense environments – as traced by ZF-UDS-7329 – the stellar fractions could be reduced to |$f_{\star } \approx 50{{\ \rm per\ cent}}$|⁠, which is more consistent with galaxy formation models and the stellar-to-halo mass relation at lower redshifts. We conclude that ZF-UDS-7329 forms extremely efficient in the early universe, but does not necessitate unseen galaxies at higher redshifts since the inferred SFR of ancestors are consistent with those seen in |$z\gt 6$| galaxies.

1 INTRODUCTION

The sensitivity and near-infrared wavelength coverage of JWST enables us to trace the galaxy population in the first 500 Myr (redshift |$z\gt 10$|⁠). After two years of observation, numerous high-redshift dropouts have been discovered (Castellano et al. 2022; Naidu et al. 2022; Finkelstein et al. 2023; Adams et al. 2024; Hainline et al. 2024), with seven spectroscopically confirmed at |$z\gtrsim 12$| (Curtis-Lake et al. 2023; Robertson et al. 2023; D’Eugenio et al. 2024c; Fujimoto et al. 2023; Wang et al. 2023b; Castellano et al. 2024; Carniani et al. 2024a). Results show that the bright-end of the ultraviolet (UV) luminosity function only weakly evolves with redshift (Harikane et al. 2023; Donnan et al. 2024; Finkelstein et al. 2024; Robertson et al. 2024), thereby challenging theoretical models – such as numerical simulations (Davé et al. 2019; Shen et al. 2020; Vogelsberger et al. 2020; Kannan et al. 2023; Wilkins et al. 2023), semi-analytical galaxy formation (Dayal et al. 2014; Cowley et al. 2018; Yung et al. 2019), or empirical models (Tacchella, Trenti & Carollo 2013; Mason, Trenti & Treu 2015; Sun & Furlanetto 2016; Tacchella et al. 2018b; Behroozi et al. 2020).

The surprisingly high number density of galaxies at |$z\gt 10$|⁠, including the extreme UV brightness of individual systems, has sparked a range of theoretical ideas and developments. One way to overcome this tension between observations and theory is to modify cosmology, i.e. not assume Lambda cold dark matter (⁠|$\Lambda$|CDM), to boost the number density of dark matter haloes. This can be achieved by enhancing the matter power spectrum (Sabti, Muñoz & Kamionkowski 2024) or increasing the expansion of the Universe at early times with ‘early dark energy’ (Shen et al. 2024). In contrast, other solutions focus on modifying the baryonic physics of galaxy formation, which has mostly been calibrated with lower-redshift observations, and hence may need some refinement for the low-metallicity and dense environment of cosmic dawn. Such ideas include a higher star formation efficiency (Dekel et al. 2023; Li et al. 2024), a top-heavy stellar initial mass function (IMF; Inayoshi et al. 2022; Cueto et al. 2024; Trinca et al. 2024; Ventura et al. 2024; Yung et al. 2024), and a UV contribution from accreting black holes (Inayoshi et al. 2022; Trinca et al. 2024; Hegde, Wyatt & Furlanetto 2024). In addition to boosting the UV luminosity of the whole galaxy population, scatter of the UV luminosity at fixed halo mass can have a drastic effect on the bright end of the UV luminosity function because of the steepness of the dark matter halo mass function at these early times (Mason, Trenti & Treu 2023; Mirocha & Furlanetto 2023; Shen et al. 2023; Sun et al. 2023; Gelli, Mason & Hayward 2024; Kravtsov & Belokurov 2024).

These different ideas and scenarios can be tested. Current JWST observations of |$z\gt 10$| galaxies reveal a great diversity of stellar populations and morphologies. Galaxies such as GN-z11 probably host an active galactic nucleus (AGN; Maiolino et al. 2024) and a significant fraction (⁠|$\gt 50{{\ \rm per\ cent}}$|⁠) of the UV output stems from the central point source (Tacchella et al. 2023b). GN-z11, after removing the point source, has a stellar mass of |$M_{\star }\approx 10^{9}~{\rm M}_\odot $| and a star formation rate (SFR) of |$\mathrm{SFR}\approx 10~{\rm M}_\odot ~\mathrm{yr}^{-1}$|⁠. However, not all bright, |$z\gt 10$| galaxies are AGN: JADES-GS-z14-0, the most distant spectroscopically confirmed galaxy at |$z=14.3$|⁠, is unexpectedly luminous, is spatially resolved with a half-light radius of 260 pc, and consistent with being dominated by stellar continuum emission (Carniani et al. 2024a; Helton et al. 2024a). This object, with a stellar mass of |$M_{\star }\approx 10^{8.7}~\mathrm{M}_{\odot }$|⁠, has a |$\mathrm{SFR}\approx 22~\mathrm{M}_{\odot }~\mathrm{yr}^{-1}$|⁠, which is an order of magnitude higher than the dark matter accretion rate of typical haloes at this redshift (Tacchella et al. 2018b). The less luminous galaxies at |$z\gt 10$| are typically doubling their stellar mass every few tens of Myr and are morphologically compact, with sizes of the order of 100 pc, implying high stellar mass and SFR surface densities (Robertson et al. 2023; Wang et al. 2023b; Morishita et al. 2024). While these observations directly probe the earliest galaxies known, they are probably biased towards the brightest systems and cover only the rest-frame UV emission, giving us a limited view on the stellar populations of those galaxies (but see Helton et al. 2024a; Zavala et al. 2025).

A complementary approach to probe the first generation of galaxies is via archaeological lookback studies. JWST enables us to probe the crucial rest-frame optical wavelength range out to |$z\approx 9$|⁠. While star-forming galaxies are abundant, extracting lookback information such as the stellar mass, the stellar age and the star formation history (SFH) is extremely challenging, because the bursty SFHs of galaxies outshine older stellar populations (Papovich et al. 2023; Tacchella et al. 2023a; Endsley et al. 2024; Witten et al. 2024). In the case where young stars dominate the emission and older stellar populations do not leave any spectral signatures, the priors related to the SFHs play a crucial role (Tacchella et al. 2022b; Whitler et al. 2023). Furthermore, AGN introduce significant uncertainties in both stellar mass and SFH estimations because of unknown contribution of the AGN to continuum emission (Kocevski et al. 2023, 2024; Wang et al. 2024; Williams et al. 2024).

Contrarily, quiescent galaxies are less affected by outshining and offer a detailed view of the early Universe. Observing quiescent galaxies at |$z\gt 3$| is challenging, as older stellar populations are both rarer and fainter, and a signal-to-noise ratio (S/N) of |$\approx 20$| per Å in the spectral continuum is required to provide meaningful constraints on early star formation, metallicity and chemical abundance (Conroy 2013; Nanayakkara et al. 2022). JWST provided new constraints on quiescent galaxies at cosmic noon, |$z\approx 1-3$|⁠, highlighting the importance of supermassive black holes in rapidly suppressing star formation (‘quenching’) in massive galaxies by efficiently ejecting gas (Belli et al. 2024; D’Eugenio et al. 2024b; Beverage et al. 2024; Park et al. 2024b; Scholtz et al. 2024; Davies et al. 2024). Quiescent galaxies have also been probed at earlier times, both at low stellar masses (Looser et al. 2024; Strait et al. 2023) and high stellar masses (Carnall et al. 2023; de Graaff et al. 2024a; Glazebrook et al. 2024; Nanayakkara et al. 2024; Setton et al. 2024).

A particularly notable example of a massive quiescent galaxy at |$z=3.2$| is ZF-UDS-7329. This galaxy was initially selected by Schreiber et al. (2018) as a quiescent candidate using the UVJ colour–colour diagram (e.g. Williams et al. 2009). Using JWST spectroscopy, Glazebrook et al. (2024) confirmed this galaxy to be quiescent and inferred a stellar mass |$M_\star \gt 10^{11}~{\rm M}_\odot $| forming at |$z \gtrsim 11$|⁠. Such an extreme SFH challenges our understanding of early stellar populations, galaxy formation and/or the nature of dark matter, and may point to the presence of undetected populations of early galaxies. Carnall et al. (2024) obtained three and analysed two medium resolution grating/filter combinations of ZF-UDS-7329, emphasizing that this galaxy aligns with the most massive galaxies expected under the assumption of 100 per cent conversion of baryons to stars. Their results suggest extreme galaxy formation physics during the first billion years, but no conflict with |$\Lambda$|CDM cosmology.

In this paper, we re-reduce and analyse the JWST data of ZF-UDS-7329 (Section 2). Specifically, we extend the analysis of Glazebrook et al. (2024) by studying the impact of different SFH priors, stellar libraries, metallicity, and IMF assumptions (Section 3). Using NIRCam imaging data, we show clear indication for a colour gradient in ZF-UDS-7329, which is consistent with an old, quiescent bulge component and a younger disc component (Section 4). Our results highlight efficient star formation in the first billion years, though uncertainties remain regarding self-consistent modelling of |$\alpha$|-enhanced stellar population models, mergers diluting the archaeological signal and cosmic variance (Section 5). We conclude in Section 6.

Throughout this work, we assume a Planck18 flat |$\Lambda$|CDM cosmology, with |$\Omega _m = 0.315$|⁠, |$H_0 = 67.4$| |$\mathrm{km\, s^{-1}}$|⁠, and a baryon fraction of |$f_b=\Omega _b/\Omega _m=0.156$| (Planck Collaboration VI 2020). We assume a Kroupa IMF (Kroupa 2001) and a solar metallicity |$\mathrm{Z}_\odot = 0.014$| (Asplund et al. 2009).

2 OBSERVATIONS AND DATA

2.1 NIRCam and NIRSpec observations

ZF-UDS-7329 was observed with the JWST NIRSpec instrument (Jakobsen et al. 2022) in multiobject mode using the MSA (Micro-Shutter Assembly; Ferruit et al. 2022), as part of programme #2565 as described in Glazebrook et al. (2024). The MSA target ID is 12629 (which is also the 3D-HST ID; Brammer et al. 2012; Skelton et al. 2014; Momcheva et al. 2016). The MSA was configured with five-shutter slitlets and three nodded exposures. Each exposure consisted of a single integration of 45 groups in NRSIRS2RAPID mode (Rauscher et al. 2007), giving 667 s per exposure and 2001 s total observing time. The disperser used was the prism, which produced spectral data covering wavelengths |$0.6\!-\!5.2$||$\mu \mathrm{m}$| at a resolving power of |$\sim 50\!-\!350$| (Jakobsen et al. 2022).

In addition to spectroscopy, we use JWST NIRCam data from the public PRIMER survey¹ (programme #1837; PI J. Dunlop). These data consist of imaging in 8 NIRCam bands across |$1.0\!-\!5.0$||$\mu \mathrm{m}$| (F090W, F115W, F150W, F200W, F277W, F356W, F410M and F444W). A composite image using three of these is displayed in Fig. 1. The flux measurements in the NIRCam bands are from Glazebrook et al. (2024). In addition, we use cutouts obtained from the DAWN JWST Archive² (DJA; Valentino et al. 2023).

2.2 NIRSpec reduction

To assess the impact of possible systematics in the data reduction on the measurements, we compared the data reduction from two different pipelines. The spectrum from Glazebrook et al. (2024, |$f^{\mathrm{K}}_{\lambda }$| in Fig. A1) was reduced with the publicly available jwst pipeline v1.12.5 and the Calibration Reference Data System (CRDS) context 1149 (see Nanayakkara et al. 2024). As an alternative, we use the GTO pipeline (Carniani in preparation; |$f^{\mathrm{C}}_{\lambda }$| in Fig. A1), described in Curti et al. (2024). The reduced spectra are compared in Fig. A1 (see Appendix  A), alongside the NIRCam data highlighted with red circles. Nanayakkara et al. (2024) calibrated the |$f^{\mathrm{K}}_{\lambda }$| spectrum to match the shape and normalization of the photometry by fitting a second-order polynomial. In the case of |$f^{\mathrm{C}}_{\lambda }$|⁠, we applied a uniform scaling to match the photometry. The difference in the spectral shape of the rescaling is chiefly responsible for |$f^{\mathrm{K}}_{\lambda }$| appearing slightly bluer of the order of 0.05 mag in the F200W|$-$|F356W colour, which is discernible from the plot of |$\chi _{\mathrm{eff}} := |f^{\mathrm{K}}_{\lambda } - f^{\mathrm{C}}_{\lambda }| / \sqrt{\sigma _{\mathrm{C}}^2 + \sigma _{\mathrm{K}}^2}$|⁠. It should be noted that since these are not actually independent measurements, |$\chi _{\mathrm{eff}}$| merely quantifies the discrepancy between the reductions relative to the estimates of statistical uncertainty, |$\sigma _{\mathrm{C}}$| and |$\sigma _{\mathrm{K}}$|⁠.

Even without any spectral calibration, the scaled |$f^{\mathrm{C}}_{\lambda }$| is reasonably consistent with the photometry. Both spectra, however, contain artefacts larger than the estimated noise, |$\sigma$|⁠, such as the prominent bumps at |$\sim 6800$| Å rest frame, that do not correspond to any known emission features. In all subsequent fitting, we used the reduction from Curti et al. (2024), with the spectrum and photometry S/N capped to 20. Since ZF-UDS-7329 is neither a point source, nor covers the slit angle, the line-spread function (LSF) of Glazebrook et al. (2024) was chosen. This was derived using the NIRCam images and NIRSpec/MSA instrument model to estimate the wavelength-dependent LSF, parametrized by its Gaussian full-width half maximum (FWHM; de Graaff et al. 2024c). The wavelength-dependent spectral resolution ranges from 4000 km s|$^{-1}$| (around 1 |$\mu \mathrm{m}$|⁠) to 400 km s|$^{-1}$| (5 |$\mu \mathrm{m}$|⁠). We note that the overall normalization and large-scale wavelength dependence of the spectrum is not important in our analysis because we use a polynomial to match the spectrum to the photometry (see Section 3.1).

2.3 Preliminary comments

The spectrum (see Figs 2 and A1) displays a prominent 4000 Å break that is typical for old, |$\gt 1$| Gyr stellar populations, for which ionized metal absorption lines and an absence of hot, blue stars suppress the strong Balmer break seen in A-type stellar spectra. This is in contrast to so-called post-starburst quiescent galaxies, which have recently quenched their star formation (within |$\sim 500 \, {\rm Myr}$| prior to observation) following a strong burst (e.g. D’Eugenio et al. 2023; Looser et al. 2023). An old age is also consistent with the strength of the Mg ib 5175 Å triplet and TiO absorption bands at |$\sim$| 7150 Å and 9400 Å, which may indicate some degree of |$\alpha$|-element enhancement. Particularly notable is the strong Na d 5890 Å doublet absorption line, hypothesized to be either a result of cold interstellar medium (ISM) or an enhanced Na/Fe abundance (Jeong et al. 2013).³ No prominent emission lines are visible, as expected for a galaxy with little ongoing star formation.

3 STELLAR POPULATION ANALYSIS

3.1 Prospector set-up

The spectrum and photometry are simultaneously fit with the Bayesian code prospector (Johnson et al. 2019, 2021). prospector uses a python implementation (Speagle 2020; Koposov et al. 2022) of the Dynamic Nested Sampling algorithm (Higson et al. 2018) in which the posterior probability distribution is estimated by Monte Carlo sampling over nested shells bounded by iso-likelihood contours. The stellar population synthesis model is based on fsps code (Felixible Stellar Population Synthesis; Conroy, Gunn & White 2009; Conroy & Gunn 2010), which allows for a great deal of flexibility in the inputted spectral templates, isochrones, dust laws, and emission models. We explore both the empirical MILES (Sánchez-Blázquez et al. 2006) or the synthetic C3K (Conroy et al. 2019) stellar libraries, in conjunction with the MIST isochrones (Choi et al. 2016; Dotter 2016). In the former case, the limited resolution of the templates means that only the |$0.8\!-\!3.0$||$\mu \mathrm{m}$| observed range (⁠|$0.2\!-\!0.7$||$\mu \mathrm{m}$| in rest-frame) of the SED is fitted.

Within the prospector framework, we fit a model with 20 free parameters that describe the stellar populations of ZF-UDS-7329 in conjunction with dust and nebular effects. Table 1 summarizes all free parameters of the fiducial model, alongside their priors and resulting posterior probabilities.

We use a flexible (sometimes referred to as non-parametric) ‘continuity’ SFH prior, which penalizes drastic changes in SFR by fitting for the ratios of the SFRs between adjacent time bins (Leja et al. 2019a). Specifically, while no parametric form for the SFH as a function of time is assumed, time bins and a prior still need to be specified. As done in previous works (Leja et al. 2019a; Tacchella et al. 2022a; Wan et al. 2024), we assume the continuity prior, which assumes a Student-t distribution for the log ratio of the SFR in adjacent time bins, with |$\nu =2$| degrees of freedom and a scale of |$\sigma =0.3$|⁠. This choice was originally motivated by the SFH in the Illustris simulations (Leja et al. 2019a). The scale |$\sigma$|⁠, which we do not vary in this work, controls the amount of variable star formation, and has been increased for high-z galaxies in some works to account for possible bursty SFHs (Tacchella et al. 2022b). Importantly, this prior enforces a constant expectation value for the SFH, something we will address in Section 3.3. Furthermore, in this fiducial model, no star formation is allowed before |$z=20$|⁠, leaving nine bins that are approximately logarithmically spaced except between lookback times |$1\!-\!1.75$| Gyr, to properly characterize the formation period while capturing any recent star formation.

Dust attenuation is implemented using the two-component model proposed in Charlot & Fall (2000) and Noll et al. (2009), in which the attenuation due to birth-clouds and ambient ISM are approximated as two separate screens. The diffuse component follows an attenuation curve with a variable slope index, n, and a UV dust absorption bump,

where |$\tau _{\mathrm{2,V}}$| is the effective V-band optical depth, and |$k(\lambda)$| is the Calzetti curve (Calzetti et al. 2000), which has total-to-selective ratio |$A_{\mathrm{V}} / E(B-V) = 4.05$|⁠. |$D(\lambda)$| is a Lorentzian-like Drude profile approximating the UV bump, for which we tie the strength to the best-fitting value of n according to the results of Kriek & Conroy (2013). We adopt uniform, and clipped-Gaussian priors on n, and |$\tau _{\mathrm{2,V}}$|⁠, respectively. The birth-cloud component, |$\tau _{\mathrm{1}}(\lambda)$|⁠, only attenuates nebular and stellar emission from stars that are younger than |$10~{\rm Myr}$|⁠, and follows a simple power law:

We fit for the ratio |$\tau _{\mathrm{1,V}}/\tau _{\mathrm{2,V}}$| with a clipped Gaussian prior. Dust emission is typically negligible at rest-frame |$\lambda \lesssim 10 \, {\mu \mathrm{m}}$|⁠; hence, it is omitted from all the models.

Nebular emission is included using the default grids of fsps (Byler et al. 2019), which are based on the simulation code cloudy (Ferland et al. 2013). This is controlled by three parameters, the gas ionization parameter, U; the gas metallicity, |$Z_{\mathrm{gas}}$|⁠; and the gas velocity dispersion |$\sigma _{\mathrm{gas}}$| (which for us is a nuisance parameter). Additional flexibility, by the means of directly fitting for, and marginalizing over, emission line fluxes, was not considered given the lack of any observable emission in the data. As our templates cannot reproduce the strong Na d absorption, this line is masked from the fitting procedure.

An outlier model and noise jitter term are also included to mitigate bad pixels, possible data-model mismatches, and underestimated noise as follows. We fit for a fraction |$f_{\mathrm{out}}$| of the data points treated as outliers, which have errors inflated by a factor of 50. The outlier locations are marginalized over as described in Appendix D of Johnson et al. (2021). A multiplicative inflation of all the errors by factor |$j_\mathrm{spec}$| is fitted for by including kernels for uncorrelated noise in the likelihood calculations.

Since the spectrum, |$f^{\mathrm{C}}_{\lambda }$|⁠, is not flux-calibrated, we fit a Chebyshev polynomial to match the model continuum to the data at each likelihood call. Specifically, at each likelihood call, we match the model spectrum to the normalization of the spectroscopic data by fitting a polynomial in wavelength to their ratio. The net effect of our approach is that the large-scale continuum shape and normalization of the model are set by the photometry. By fitting a moderate-order polynomial to the ratio between the observed and physical model spectrum at each likelihood call, the spectroscopic calibration model basically removes all information content from the continuum shape of the spectroscopic data. This means that the continuum shape of the observed spectrum does not inform any of the galaxy’s physical parameters. Instead, information about physical parameters that affect the continuum shape is derived from photometry, which does not include any multiplicative calibration model. Therefore, there is no degeneracy between the spectroscopic calibration model and the galaxy’s parameters, such as the dust content or the SFH. Overall, the fitting of this polynomial is found to have only a small effect on the inferred population parameters (typically less than |$\sim 5{{\ \rm per\ cent}}$|⁠). The order of the polynomial was chosen as five with MILES, or seven with C3K templates, such that it is only sensitive to spectral variations across |$\gtrsim 500 \, \mathrm{\mathring{\rm A}}$|⁠.

3.2 Fiducial results

The maximum a posteriori (MAP) spectra of the fiducial C3K and MILES models are presented in Fig. 2 and the posterior distribution of the prospector model parameters are listed in Table 1. Setting |$j_\mathrm{spec}=1$|⁠, we find reduced |$\chi ^2$| values of 0.37 and 0.44 for C3K and MILES, respectively, which show good agreement given S/N of |$\sim 20$| across |$\sim 200\!-\!500$| spectral pixels with 20 free parameters. We find |$j_{\rm spec}=2.5-3.5$|⁠, which is possibly caused by data-model mismatches (in addition to what |$f_{\rm outlier}$| covers), such as |$\alpha$|-enhancement or TP-AGB phase contributions (Sections 5.4 and 5.5). The extrapolation, without convolution with the LSF, of the MILES best-fitting past 3 |$\mu \mathrm{m}$|⁠, is also displayed with a dotted red line (it should be noted that the NIR photometry is still included in the fitting). A corner plot of the posterior distributions of some of the key free parameters is shown in Fig. 3, alongside the inferred SFHs. There are markedly stronger degeneracies between the dust attenuation parameters and stellar metallicity in the MILES posterior, as well as a larger uncertainty in these values overall. We speculate that this is primarily due to the exclusion of the rest-frame NIR spectral data from the fit. The derived dust optical depths, |$\tau _{\mathrm{2,V}} \sim 0.2$|⁠, are relatively low, with a dust law steeper than Calzetti, |$n \approx -0.2$|⁠, which is consistent with the galaxy having little recent star formation.

We find a stellar age of |$t_{50} \approx 1.6 \, {\rm Gyr}$| (corresponding to |$z_{50}\approx 11$|⁠), defined as the look-back time at which half of the stellar mass has formed. The degeneracy of the other parameters with the age is relatively weak. As expected, all the nebular emission parameters are broadly unconstrained. The models preferentially add some weak emission for a small reduction in the overall |$\chi ^2$|⁠; however, the average SFR over the last 100 Myr remains |$\lesssim 10 \, {\rm M}_\odot ~\mathrm{yr}^{-1}$|⁠. This translates to a specific SFR (sSFR) of |$\sim 0.03 \, {\rm Gyr} ^{-1}$| and mass doubling time-scale of |$\sim 30 \, {\rm Gyr}$|⁠, which is 15 times larger than the age of the Universe at that epoch (⁠|$t_{\rm H} \approx 2 \, {\rm Gyr}$|⁠), implying that this galaxy is well quiescent. We note that the SFH in Fig. 3 shows a downward trend to |$\sim 2~\mathrm{M}_{\odot }/\mathrm{yr}$| in the past 10 Myr, which is caused by the absence of emission line in the spectrum. The SFR in the time period of |$10~\mathrm{Myr}\lt t - t_{\rm obs} \lt 1~\mathrm{Gyr}$| is generally low (⁠|$\lt 10~\mathrm{M}_{\odot }/\mathrm{yr}$|⁠) and driven by the low UV flux, but less constrained than the SFR in the past 10 Myr due to the degeneracy with the dust attenuation law and stellar metallicity.

The fits also give high formed stellar masses of |$\log _{10}{(M_\star /{\rm M}_\odot)} = 11.34^{+0.01}_{-0.02}$|⁠, and |$11.37^{+0.02}_{-0.02}$|⁠, or surviving stellar masses of |$\log _{10}{(M_{\star \mathrm{,surv}}/{\rm M}_\odot)} = 11.17^{+0.01}_{-0.02}$|⁠, and |$11.20^{+0.02}_{-0.02}$|⁠, for the MILES and C3K fits, respectively. This assembles in a strong burst with |$\sim 70 {{\ \rm per\ cent}}$| of the mass forming between |$z=10\!-\!20$| with a |$\mathrm{SFR}\approx 500-700~{\rm M}_\odot ~\mathrm{yr}^{-1}$|⁠. We note that the strong constraint on the SFR in the oldest time bin is likely influenced by our decision to disallow any star formation before |$z = 20$|⁠, which we explore in detail in Section 3.3.

Both stellar libraries lead to solutions with super-solar metallicities, [Z/H] |$= 0.03 \pm 0.10$| (MILES), and |$0.11 \pm 0.05$| (C3K), which is surprising for a galaxy formed at |$z \approx 10-20$|⁠. The MILES results are all consistent with the findings of Glazebrook et al. (2024). However, the fitting of spectral data beyond 3 |$\mu \mathrm{m}$| seems to push the best fit to adopt higher Z due to the inclusion of metallicity-sensitive absorption bands. This effect is explored in the Section 3.4, and strongly implies a degree of |$\mathrm{\alpha }$|-element enhancement, which we discuss in Section 5.4. Carnall et al. (2024) find an even higher value of [Z/H] |$= 0.35^{+0.07}_{-0.08}$|⁠, which is beyond the range of our templates, while all other population parameters are relatively consistent with our findings. We also note that the inconsistency in the derived redshifts between MILES and C3K likely arises from a small wavelength calibration issue across the full spectral range (e.g. D’Eugenio et al. 2024a, their Section 9.3).

Lastly, we corroborate these results with those from the Penalized Pixel-Fitting (ppxf) method (Cappellari 2023), which proceeds via a |$\chi ^2$|-minimization process, in lieu of prospector’s Bayesian framework. In this model, following the approach of Looser et al. (2023), we use simple stellar population (SSP) templates that are a combination of the synthetic C3K model atmospheres (Conroy et al. 2019) and Mesa Isochrones and Stellar Tracks (MIST) isochrones (Choi et al. 2016). Similarly to the Bayesian approach, we find a super-solar metallicity, high stellar-mass and age > 1.5 |${\rm Gyr}$|⁠. The full methodology and results are presented in Appendix  B, and the MAP spectrum in Fig. 2, alongside the prospector best fits.

3.3 Rising SFH prior and stellar age constraints

The flat continuity prior utilized in the fiducial model typically favours a median stellar age half that of the universe (Leja et al. 2019b). In our case, however, we argue it biases towards the oldest possible solution for the data by penalizing substantial deviations from a constant SFR. We explore this effect in Appendix  D by fitting a mock SSP with age |$= 1.5~{\rm Gyr}$|⁠. Our results demonstrate that the fiducial model prior overestimates the age by |$\sim 250~{\rm Myr}$| due to the strong degeneracy between solutions older than |$\sim 1.3~{\rm Gyr}$|⁠.

In contrast, both observations and theory point toward increasing SFHs with cosmic time at |$z\gt 3$|⁠. Observations show that SFHs are increasing on average, with a high fraction of bursty galaxies with high-EW emission lines (Looser et al. 2023; Tacchella et al. 2023a; Endsley et al. 2024; Simmonds et al. 2024a).

To account for the expectation for rising SFHs of high-z galaxies, we now motivate a rising SFH prior, which is still flexible enough to include quenching SFHs. We motivate this prior based on dark matter accretion, which goes beyond recent considerations by de Graaff et al. (2024a, considered a fixed ratio of increase between adjacent time bins) and Wang et al. (2023a, considered an increasing prior based on the cosmic SFR density). We consider the |$\Lambda$|CDM cosmology in the Einstein-deSitter regime, which is a useful approximation at |$z\gt 1$| and becomes more and more accurate at higher redshifts. In such a regime, following Dekel et al. (2013), the average accretion rate of mass into haloes of mass |$M_{\rm h}$| at z can be approximated by an expression of the form

where the power-law index of |$\mu =5/2$| and |$\alpha =4/5$| can be understood from the Press-Schechter (PS) and Extended-PS (EPS) approximations of gravitational structure formation in cosmology (Press & Schechter 1974; Bond et al. 1991). This functional form was indeed found to be a good fit to the halo growth in cosmological N-body simulations (Wechsler et al. 2002; Neistein & Dekel 2008; Dekel et al. 2013).

To implement this rising SFH prior within prospector, we shift the mean of the t-distribution for each SFR ratio such that the star formation in each bin obeys

where z refers to the mean redshift of each time bin, motivated directly by equation (3). For example, for |$z_{\rm obs}=3.2$|⁠, the expectation is that the SFR at |$z=10$| is reduced by a factor of |$\approx 20$|⁠. The 1|$\sigma$| regions of both priors are shown in the left panel (flat continuity prior) and right panel (rising continuity prior) of Fig. 4. Both priors have been normalized to a total stellar mass of |$M_{\star }=10^{11.37}~{\rm M}_\odot $|⁠. Repeating the 1.5 |${\rm Gyr}$| mock SSP test with the rising prior, we now find the best fit underestimates the age by |$\sim 100~{\rm Myr}$| (see Appendix  D). The test indicates the importance of the choice of prior when measuring stellar ages on an accuracy of |$\sim 10{{\ \rm per\ cent}}$| of quiescent galaxies.

In addition to the 1|$\sigma$| prior region, Fig. 4 shows the resulting best-fitting SFHs for the flat continuity SFH prior and the rising SFH prior in the left and right columns, respectively. We assume for both fits the C3K stellar library. Both fits reproduce the photometry and the spectrum equally well, with the flat continuity SFH prior and the rising SFH prior giving a reduced |$\chi ^2 = 0.37$| and 0.42 (calculated with |$j_\mathrm{spec}=1$|⁠), respectively. Counter-intuitively, the rising prior has a lower noise inflation factor |$j_\mathrm{spec}\sim 2$| compared to the fiducial |$j_\mathrm{spec}\sim 3$|⁠, but the values are consistent within |$1\sigma$|⁠. The posterior distributions for the key population parameters are shown in Fig. C1 in Appendix  C.

The rising prior instigates a best fit that has the vast majority of the stellar mass in the second-oldest bin. This leads to a (very constrained) median age |$\sim 200$||${\rm Myr}$| younger than the fiducial result. Repeating the fit with more time bins did not substantially widen the posterior on the age, as the model always favours a short burst of star formation at |$z \sim 8$|⁠. It is also worth noting that our use of a non-parametric SFH artificially gives a very constrained posterior for the specific SFR in the first time bin (see Fig. 4). While the other parameters remain generally consistent (approximately within |$1\sigma$|⁠), the rising prior favours overall a younger, more dusty solution. Importantly, the SFH obtained with the rising prior is more consistent with high-redshift observations of star-forming galaxies, which are overlaid in Fig. 4. It also significantly reduces the tension with the inferred star formation efficiency at |$z\gt 8$| and theoretical expectations, which we discuss in Section 5.1.

In addition to changing the base of the SFH prior, we now modify when the SFH starts. In our prospector models described above, we assume a negligible SFR at |$z\gt 20$|⁠. We run two additional models, one assuming the SFH starts at |$z=10$|⁠, the other starting at |$z=\infty$|⁠, corresponding to the big bang. We run these two additional models for both the C3K and the MILES stellar libraries. We find a stellar age of |$t_{50}=1.4$| Gyr (⁠|$t_{50}=1.3$| Gyr) and |$t_{50}=1.7$| Gyr (⁠|$t_{50}=1.8$| Gyr) for the |$z=10$| and |$z=\infty$| model, assuming C3K (MILES).

The MAP spectra for these tests are displayed in Appendix Fig.  E. These demonstrate that forcing a |$\sim 300 \, {\rm Myr}$| younger fit does not significantly disagree with the data, with most values remaining within the |$5 {{\ \rm per\ cent}}$| error margin. From oldest to youngest age, we find reduced |$\chi ^2$| values (setting |$j_\mathrm{spec}=1$| for comparison) of 0.43, 0.44, 0.49 for the MILES models and 0.37, 0.37, 0.38 for the C3K models, which cover a larger spectral range. We also note that the inferred noise inflation factor is relatively unchanged in these tests, with all values remaining consistent with |$j_{\mathrm{spec}}\sim 2.5$|⁠.

3.4 Metallicity investigation

To quantify the poorness of the subsolar models in fitting the data, the C3K fits were re-run with fixed metallicities of [Z/H] |$= -1.0, -0.3, 0.0, 0.18$|⁠, giving reduced |$\chi ^2$| values of |$0.58, 0.40, 0.37, 0.37$|⁠, respectively (calculated assuming |$j_{\mathrm{spec}}=1$|⁠). The resulting MAP spectra are presented with lines coloured blue to red in Fig. 5, which is zoomed into the rest-frame optical-to-NIR region for clarity. The figure demonstrates that the subsolar models fail to reproduce many of the absorption features of the SED, particularly the Mg ib triplet or TiO band at |$\sim$| 7150 Å.

The best-fitting SFHs are also plotted for each case. Unsurprisingly, there is a rough trend in which the lower-metallicity models compensate by having a greater proportion of the stellar mass form in the oldest time bin. Conversely, the subsolar fits also have a higher very recent SFR |$\sim 3 \, {\rm M}_\odot \mathrm{yr}^{-1}$| in the last 10 Myr, partially due to the model erroneously fitting the continuum with broadened emission. Overall, the median ages remain largely unchanged from the fiducial fits. The possibility of |$\mathrm{\alpha }$|-element enhancement is discussed in Section 5.4.

3.5 IMF investigation

Finally, the effect of IMF variation is explored by re-running the fits after modifying the piece-wise power-law slope of the Kroupa (2001) IMF,

where

which was used in our fiducial prospector models. Here, the C3K templates are used to ensure that any features sensitive to IMF variation in the mid-IR region are captured by the fit.

Several studies have found evidence that the IMF becomes more bottom-heavy in present-epoch, |$z \sim 0$|⁠, ETGs than the Milky Way IMF (Conroy & van Dokkum 2012; Cappellari et al. 2013; La Barbera et al. 2013, 2017; Spiniello, Trager & Koopmans 2015). An analysis of optical-to-NIR absorption lines by Conroy & van Dokkum (2012) found the deviation is stronger at larger |$\alpha$|-enhancements and velocity dispersions. Given that we have derived a low level of dust attenuation, a bottom-heavy IMF (as opposed to cold ISM) would help to reproduce the strength of the Na d line if combined with some degree of Na-enhancement (Jeong et al. 2013). While this would increase our already unexpectedly-high |$M_\star$|⁠, an increased prevalence of M-dwarf features in the NIR region might be slightly degenerate with the population’s age at our limited resolution (Conroy & van Dokkum 2012). We test these possibilities by raising |$\alpha _1$| to an extreme value of 3.0, above even the Salpeter slope of 2.35, to inflate the proportion of low-mass stars. This leads to a best-fit surviving stellar mass enlarged by |$\sim 0.3$| dex, |$\log _{10}{(M_{\star \mathrm{,surv}}/{\rm M}_\odot)} = 11.60^{+0.02}_{-0.02}$|⁠, with all other parameters largely unchanged. This is because faint, low-mass stars only contribute on the order of |$1{{\ \rm per\ cent}}$| to the bolometric luminosity.

In contrast, observations at |$z \sim 2$|⁠, such as an analysis of |$^{13}\mathrm{C}/^{18}\mathrm{O}$| transitions by Zhang et al. (2018), have suggested starburst galaxies emerge with top-heavy IMFs (Narayanan & Davé 2013). A top-heavy IMF is an attractive possibility for our galaxy because an overabundance of more luminous, high-mass stars could lower the |$M_\star /L$| ratio. However, in our case lowering the high-mass slope to |$\alpha _2 = \alpha _3 = 2.1$| (to increase the relative abundance of stars of |$M \gtrsim 4 {\rm M}_\odot $|⁠) gives only a |$\sim 0.1$| dex reduction in the surviving stellar mass, |$\log _{10}{(M_{\star \mathrm{,surv}}/{\rm M}_\odot)} = 11.13^{+0.02}_{-0.02}$|⁠, with |$M_\star $| largely unchanged. This is because the light from older, |$\gt 1 \, {\rm Gyr}$|⁠, stellar populations is dominated by stars closer to |$\sim 1 \, {\rm M}_\odot $| (Greggio & Renzini 2011; Esdaile et al. 2021). In fact, a much top-heavier IMF would only give a higher dynamical mass-to-light ratio due to the surplus of stellar remnants (Dabringhausen, Kroupa & Baumgardt 2009; Narayanan & Davé 2013). Overall, the effects of IMF variation are unlikely to significantly change the inferred SFH of this galaxy. In Appendix  F, we show the different IMF models give posteriors that are all in agreement, with the exception of the stellar mass.

4 MORPHOLOGY

Fig. 1 shows that ZF-UDS-7329 has an overall regular and smooth morphology, with a central concentration (similar to a bulge) and an extended, edge-on disc-like structure. The NIRSpec slit traces the central region of the galaxy. In this section, we quantify the morphology of ZF-UDS-7329 in different NIRCam bands in order to further shed light onto the physical nature of this galaxy and put the stellar population analysis into context.

Both Glazebrook et al. (2024) and Carnall et al. (2024) have fit a single Sérsic profile to ZF-UDS-7329. Glazebrook et al. (2024) used Galight (Ding et al. 2022) to fit the F444W image, finding a half-light semimajor axis |$R_{\rm e}=1.15_{-0.08}^{+0.08}$| kpc, |$n=2.41_{-0.29}^{+0.24}$| and an axis ratio of |$(b/a)=0.33_{-0.01}^{+0.01}$|⁠. Carnall et al. (2024) fitted the F277W band with the PetroFit code (Geda et al. 2022), finding |$R_{\rm e}=0.91\pm 0.01$| kpc and |$n=2.5\pm 0.1$|⁠. Both works commented on the consistency with the visual appearance of an edge-on disc. We expand here upon these studies by analysing the dependence of the morphology on wavelength in order to constrain radial colour gradients and thereby shedding light onto the formation of ZF-UDS-7329.

4.1 Single Sérsic fit

We perform our morphological analysis with pysersic (Pasha & Miller 2023). Pysersic is a code for fitting Sérsic profiles to galaxy images using Bayesian inference while accounting for the effect of the point spread function (PSF). It is written in python using jax with inference performed using numpyro (Phan, Pradhan & Jankowiak 2019; Bingham et al. 2019). In our analysis, we use model PSFs (mPSFs) which we derive by mosaicing WebbPSF models (Perrin et al. 2014), see Ji et al. (2024) for a description. We fit the F150W, F200W, F277W, F356W, F410M, and F444W with both a single Sérsic model and a double, bulge + disc Sérsic model. In the latter model, we fix one of the Sérsic models to have |$n=1$| (disc component) and assume that both Sérsic components have the same centre. We do not fit the F090W and F115W bands because of their low S/N.

For both the one-component and two-component models we perform a fit to the individual bands independently of each other as well as jointly. For the joint single Sérsic fit, we assume that the Sérsic index n, the ellipticity and half-light radius are linked across wavelength by a second-order polynomial in |$1/\lambda$|⁠. The position angle and centre are the same for all bands. We find that the joint and individual fits are overall consistent, with Sérsic index n and half-light size |$R_{\rm eff}$| decreasing with increasing wavelength (Fig. 6). The largest, most significant difference is found in the bluest band (F150W), where the joint fit gives |$n=3.53\pm 0.28$| and |$R_{\rm eff}=1.63\pm 0.07$| kpc, while the fit to F150W individually gives |$n=4.33\pm 0.17$| and |$R_{\rm eff}=2.18\pm 0.14$| kpc. Beyond 2 |$\mu$|m, both Sérsic index and size do not vary significantly, and we find |$n=1.55\pm 0.03$| and |$R_{\rm eff}=1.27\pm 0.02$| kpc for F444W (rest-frame |$\lambda _{\rm rest}\approx 1~\rm{\mu m}$|⁠). For the ellipticity, which shows only a marginal wavelength dependence, we find |$\varepsilon =0.71\pm 0.01$| for F444W, which translates to a semi-to-major axis ratio of |$0.29\pm 0.01$|⁠. In summary, the single Sérsic fit is consistent with an edge-on disc-like galaxy at |$\lambda _{\rm rest}\approx 1~\rm{\mu m}$|⁠, but with evidence for a larger size and Sérsic index for F150W (⁠|$\lambda _{\rm rest}\approx 0.356~\rm{\mu m}$|⁠). We discuss the implication of the galaxy size for the dark matter halo mass in Section 5.2.

4.2 Bulge-disc decomposition

This clear indication for a wavelength-dependent morphology, implying a colour gradient, motivates us to perform a bulge-disc decomposition. For the bulge-disc decomposition, we assume that the fluxes of both bulge and disc can vary as a function of wavelength, while the structural parameters are constant with wavelength. Furthermore, we assume that the bulge has a smaller half-light radius than the disc, and the disc’s Sérsic index is fixed to 1.

We find that the bulge-to-total ratio (B/T) increases with wavelength, |$0.65\pm 0.02$| to |$0.91\pm 0.01$| in F150W (⁠|$\lambda _{\rm rest}\approx 0.356~\rm{\mu m}$|⁠) to F444W (⁠|$\lambda _{\rm rest}\approx 1.043~\rm{\mu m}$|⁠). This is consistent with the results from above, where the bluer bands are better described by a more extended profile. Furthermore, the Sérsic index of the bulge is |$n=2.00\pm 0.01$|⁠, which indicates that we are running up against the lower-bound of the prior. This means that this central structure is not a ‘classical’, |$n=4$| bulge. We find the size of the bulge and disc to be |$1.04\pm 0.01$| kpc and |$3.45\pm 0.10$| kpc, respectively.

4.3 Colour gradient

The single Sérsic fit and the bulge-disc decomposition both imply a red centre and blue outskirt. While the NIRSpec spectrum mostly probes the central region of ZF-UDS-7329 (Fig. 1) and is consistent with an old, quiescent stellar population, we now address how much bluer the outskirt is relative to this central region.

We base our analysis on the rest-frame U|$-$| V versus V|$-$|J colour–colour (UVJ) diagram (Williams et al. 2009), which is a powerful diagram to differentiate between reddening due to older stellar populations and due to dust attenuation (though this can be complicated by the fact of varying attenuation law; Leja, Tacchella & Conroy 2019c). We take the NIRCam filters F150W (⁠|$\lambda _{\rm rest}\approx 0.356~\rm{\mu m}$|⁠), F200W (⁠|$\lambda _{\rm rest}\approx 0.471~\rm{\mu m}$|⁠), and F444W (⁠|$\lambda _{\rm rest}\approx 1.043~\rm{\mu m}$|⁠) to trace the rest-frame U, V, and J band, respectively. Since these filter curves differ from each other, we derive colour correction terms from the best-fit prospector model. These corrections terms are 0.35 mag and |$-0.17$| mag for the (U|$-$| V) and (V|$-$| J) colour, respectively, assuming the best-fitting prospector model.

Fig. 7 shows the UVJ diagram. The black line marks the UVJ-quiescent box, i.e. colours in this region are consistent with being quiescent (i.e. old and low sSFR). The total, integrated photometry using apertures (black cross) and the bulge-disc decomposition give a red colour, within the quiescent region. Interestingly, the bulge-disc decomposition clearly shows that the bulge is significantly redder than the disc, with the bulge being consistent with a quiescent, old stellar population (SSP age of |$1-3$| Gyr, with possibly some dust attenuation), while the disc lies in the star-forming region of the UVJ diagram with an SSP age of 0.1–0.5 Gyr. This is overall confirmed by the colour–colour profile that is implied from the single Sérsic fit, which is shown by the connected circles in Fig. 7.

In summary, we find that the central region of ZF-UDS-7329, which is probed by the NIRSpec spectrum, is consistent with being quiescent, while the outer disc component is consistent with being significantly younger, consistent with a rapidly and recently quenched galaxy (e.g. Belli, Newman & Ellis 2019) or with star formation. This is overall in agreement with the SFH shown in Fig. 4, which shows an extended SFH with a non-negligible SFR over the past 1 Gyr. As we discuss below, such a colour gradient is consistent with inside-out quenching (Tacchella et al. 2015, 2018a). The high central stellar density can be formed through a gas compaction event (Zolotov et al. 2015; Tacchella et al. 2016a), which is triggered by mergers or counter-rotating streams (Dekel & Burkert 2014; Lapiner et al. 2023). This nuclear starburst will consume and expel most of the gas, leading to a quenching episode, soon after which the re-accreted and newly accreted gas with higher angular momentum will form mostly stars in an extended disc configuration (Tacchella et al. 2016b).

5 DISCUSSION

Our analysis of ZF-UDS-7329 shows that this is a massive (⁠|$M_{\star }\approx 10^{11.4}~{\rm M}_\odot $|⁠), quiescent (⁠|$\mathrm{sSFR}\approx 0.03~{\rm Gyr} ^{-1}$|⁠) galaxy at |$z=3.2$|⁠, consistent with previous studies of this object (Glazebrook et al. 2024; Carnall et al. 2024). While the spectrum is consistent with a very old stellar population (mass-weighted stellar age of 1.6 Gyr at an epoch when the Universe was only 2 Gyr old), we show that the spectrum is similarly well fit by a range of SFHs that have stellar ages of |$1.3-1.8$| Gyr. Importantly, by introducing a new SFH prior that follows dark matter halo accretion and therefore is increasing with cosmic time, the best-fitting SFH is consistent with direct, high-redshift observations of star-forming galaxies.

In this section, we discuss the implications of ZF-UDS-7329 for the star formation efficiency of early galaxies (Sections 5.1 and 5.2), the environment of ZF-UDS-7329 (Section 5.3), and highlight systematics related to the possible |$\alpha$|-element enhancement (Section 5.4).

5.1 Early forming galaxies

ZF-UDS-7329–given its early formation–has been put forward as serious challenge to our current understanding of cosmology (Glazebrook et al. 2024) and galaxy formation (Carnall et al. 2024). Our analysis shows that the JWST/NIRSpec + NIRCam data also support solutions that are consistent with the current paradigms of cosmology and galaxy formation.

While the SFH of ZF-UDS-7329 indeed shows a formation redshift of |$z\gt 7$|⁠, the exact formation epoch is not well determined: the flat SFH prior gives a formation epoch of |$z\approx 11$|⁠, consistent with Glazebrook et al. (2024), while the rising SFH prior leads to |$z\approx 8$| (Fig. 4). While this redshift difference is significant and the Universe roughly doubles its age over this redshift range, the stellar age only changes from 1.6 to 1.4 Gyr, a difference of only 200 Myr, which leads to negligible changes in the NIRSpec PRISM spectrum (see Fig. E2). While we largely break the age-metallicity degeneracy (Figs 3 and 5), there is still the fundamental limitation that we cannot self-consistently model |$\alpha$|-element enhancement, which can lead to systematics in our stellar age estimates (see Section 5.4).

Importantly, by using the rising SFH prior, which is physically motivated, we obtain a SFH for ZF-UDS-7329 that is roughly consistent with the latest high-redshift observations of star-forming galaxies. Specifically, Fig. 4 shows SFRs from several spectroscopically confirmed galaxies, including GS-z14-0 and GS-z14-1 (⁠|$M_\star \approx 10^8~{\rm M}_\odot $|⁠, Carniani et al. 2024a), GN-z11 (⁠|$M_\star \approx 10^9~{\rm M}_\odot $|⁠, Tacchella et al. 2023b), GLASS-z12 (⁠|$M_\star \approx 10^9~{\rm M}_\odot $|⁠, Castellano et al. 2024), REBELS-04 (⁠|$M_\star \approx 10^9~{\rm M}_\odot $|⁠, Bouwens et al. 2022), REBELS-25 (⁠|$M_\star \approx 10^{10}~{\rm M}_\odot $|⁠, Bouwens et al. 2022), JADES-318593 (⁠|$M_\star \approx 10^{11}~{\rm M}_\odot $| Simmonds et al. 2024b), RUBIES-UDS-z7 (⁠|$M_\star \approx 10^{10}~{\rm M}_\odot $|⁠, Weibel et al. 2024), and protocluster SPT0311-58 (⁠|$M_\star \approx 10^{11}~{\rm M}_\odot $|⁠, Arribas et al. 2024). In contrast, the flat SFH prior gives an SFR of |$\sim 600~{\rm M}_\odot ~{\rm yr} ^{-1}$| at |$z\approx 10{\!-\!}20$|⁠, which is nearly two orders of magnitude above those direct SFR measurements. On the other hand, the rising SFH prior leads to a SFR of |$\sim 40~{\rm M}_\odot ~{\rm yr} ^{-1}$| at |$z\approx 10{\!-\!}20$|⁠, which is much closer to these direct observations. Another important consideration is the metallicity of higher-redshift galaxies. While ZF-UDS-7329 has solar metallicity and high |$\alpha$|-element abundance (see also Carnall et al. 2024), galaxies at |$z\gt 10$| have invariably subsolar gas metallicity. Given that gas metallicity probes the |$\alpha$| element oxygen, the mismatch between the stellar metallicity in ZF-UDS-7329 and in the gas of these higher-redshift systems may be problematic. Scenarios where gas metallicity is diluted by inflows are possible, but would imply that we preferentially observe this class of systems, given that there are essentially no current observations of solar–metallicity systems at sufficiently high redshift (with the closest known examples being at too low mass and redshift, i.e. |$M_\star \sim 10^{10}~{\rm M}_\odot $| at |$zBlhyp4.6$|⁠, D’Eugenio et al. 2025; Bsemi and |$zBlhyp6.7$|⁠; Shapley et al. 2024). We can envision a period of rapid metal enrichment between |$z=10$| and the peak of the SFH, where most of the high-SFR sources are detected by ALMA; among these, only SPT0311-58 has a measurement of the gas metallicity, which is still only half solar (Arribas et al. 2024). Characterizing the gas metallicity of these systems would be a powerful consistency test for the scenario presented here.

Still, even disregarding the uncertainties due to metallicity evolution, our inferred SFHs lies roughly a factor of 2 above those direct estimates. A possible way to alleviate this tension is by considering mergers: we are measuring an SFH of a galaxy that assembled through both in situ star formation and ex situ mass accretion (mergers of galaxies). It is therefore not surprising that the SFR we estimate at early times in ZF-UDS-7329 is higher than what observed in any single high-redshift galaxy, since the SFH of ZF-UDS-7329 includes both the SFR of the main progenitor and all of its accreted satellites. While we cannot measure the ex situ mass fraction for this specific galaxy, we find that the extended morphology (relatively large size of |$\gt 1$| kpc) with its colour gradient is consistent with having undergone (multiple) mergers. Specifically, the colour gradient is consistent with inside-out quenching (Tacchella et al. 2015, 2018a). The high central stellar density could be formed through a gas compaction event (Springel & Hernquist 2005; Robertson et al. 2006; Hopkins et al. 2013; Zolotov et al. 2015; Tacchella et al. 2016a), which is triggered by mergers or counter-rotating streams (Dekel & Burkert 2014; Lapiner et al. 2023). This nuclear starburst will consume and expel most of the gas, leading to a quenching episode, soon after which the re-accreted and newly accreted gas with higher angular momentum will form mostly stars in an extended disc configuration (Tacchella et al. 2016b). In addition, this galaxy lies probably in an overdensity (see Section 5.3), where mergers are prevalent, in particular at higher redshifts (Hopkins et al. 2006).

In order to quantitatively assess the implication of the SFH on the efficiency of star formation in the early Universe, we connect ZF-UDS-7329 to the underlying dark matter halo mass function to estimate the stellar fraction, which is defined as |$f_{\star }=M_{\star }/(f_{\rm b}\times M_{\rm h})$|⁠, where |$f_{\rm b}=0.16$| is the cosmic baryon fraction and |$M_{\rm h}$| is the halo mass. The stellar fraction is also sometimes referred to as the integrated star formation efficiency. For example, |$f_{\star }=1$| corresponds to the case in which all baryons available to the galaxy have been converted into stars. It is challenging to estimate the halo mass |$M_{\rm h}$| for a single galaxy because of cosmic variance, stochastic sampling of the halo mass and galaxy stellar mass functions, and the limited volumes of observational surveys and simulations. To overcome the problem of characterizing a precise selection function for population studies, we use the extreme value statistics (EVS; e.g. Harrison & Coles 2011) approach introduced by Lovell et al. (2023) and discussed in Carnall et al. (2024), which provides a full probability distribution of halo masses given a survey volume and stellar mass of a galaxy. We adopt the Evstats package⁴ provided by Lovell et al. (2023). We assume a survey area of PRIMER from which ZF-UDS-7329 has been selected, which is a total of 378.2 arcmin|$^2$| (effective comoving volume |$\approx 1.27 \times 10^{-3}~\text{Gpc} ^3$| for |$3 \lt z \lt 4$|⁠).

The left panel of Fig. 8 presents the stellar mass growth histories obtained from our SFHs for the flat SFH prior (blue line) and the rising SFH prior (orange line). We compare these growth histories with EVS-based estimates for the most-massive galaxy expected in our survey volume as a function of redshift. Specifically, the dotted line shows the fiducial model for |$f_{\star }$|⁠, where |$f_{\star }$| is a truncated log-normal following |$f_{\star }=\ln (\mathcal {N}(\mu ,\, \sigma ^{2}))$|⁠, with |$\mu =0.135$| and |$\sigma =1$| across the interval |$0 \le f_{\star } \le 1$|⁠, motivated by a variety of theoretical and observational constraints (Lovell et al. 2023). The shaded regions indicate the 1|$\sigma$|⁠, 2|$\sigma$|⁠, and 3|$\sigma$| contours of this fiducial EVS model. The dashed line is the 3|$\sigma$| upper limit of the maximal case, in which |$f_{\star }=1$| and all available baryons are converted into stars.

Comparing the mass growth history of ZF-UDS-7329 with the EVS model estimates, we find that ZF-UDS-7329 is indeed consistent with being one of the most massive systems at the redshift of observations (⁠|$z=3.2$|⁠), but well within the fiducial model. However, at earlier times, because of the early formation and the associated high stellar mass, the growth history of both priors is in a |$\gt 3\sigma$| tension with the fiducial EVS model at |$z=8-9$|⁠. Even more challenging, the flat SFH prior leads to a growth history that is in significant tension with the |$f_{\star }=1$| model at |$z\gt 10$|⁠. However, since the growth history with the rising SFH prior shows most of the growth at |$z=7-10$|⁠, this growth history lies always below the |$3\sigma$| line of the |$f_{\star }=1$| model (shown as a dashed line), implying consistency with a more reasonable stellar fraction at early cosmic times.

The right panel of Fig. 8 shows the implied stellar fraction assuming that ZF-UDS-7329 is the most massive system in the survey volume at each redshift. We find a stellar fraction of |$f_{\star }=0.08_{-0.03}^{+0.05}$| at |$z=3.2$|⁠. For comparison, at |$z=0$|⁠, the stellar fraction peaks at a halo mass of |$\approx 10^{12}~\mathrm{M}_{\odot }$| with a value of |$f_{\star }=0.2$|⁠, and there seems to be only little evolution in the stellar fraction towards |$z=3$|⁠, with some studies finding values of |$f_{\star }=0.1$| (Moster et al. 2010; Behroozi, Wechsler & Conroy 2013; Rodríguez-Puebla et al. 2017; Behroozi & Silk 2018). We conclude that the stellar fraction at |$z=3.2$| of ZF-UDS-7329 is consistent with lower-redshift estimates from empirical models and abundance matching approaches.

At |$z=7{\!-\!}12$| (the right panel of Fig. 8), we find a stellar fraction of |$f_{\star }\approx 1{\!-\!}2$| and |$2{\!-\!}10$| for the rising and flat SFH prior, respectively. This shows that the flat SFH prior implies an unphysically high stellar fraction, while the rising prior is at maximal efficiency of star formation. Consistently, as discussed in Appendix  E, starting the SFH of the galaxy at |$z=10$| (instead of |$z=20$|⁠) is able to reproduce the observational data similarly well and significantly reduces the star formation activity at early (⁠|$z\gt 10$|⁠) times, thereby reducing the need for extreme star formation efficiencies at high redshifts.

Importantly, mergers can significantly reduce the high star formation efficiency. In the right panel of Fig. 8, we also show the stellar fraction evolution by assuming that ZF-UDS-7329 underwent one, two or three major mergers (1:1, 1:1:1, and 1:1:1:1) at |$z=6$|⁠, in which the progenitors have the same SFHs. We find that such a single major merger, double major merger and triple major merger can reduce the stellar fraction to |$f_{\star }=1.0$|⁠, |$f_{\star }=0.7$|⁠, and |$f_{\star }=0.5$| at |$z=8-9$|⁠. This highlights that mergers can further reduce the inferred stellar fraction at early times, further lowering the necessity for extreme galaxy formation physics.

But are such mergers expected at |$z\gt 3$|? Considering |$\Lambda$|CDM, Fakhouri, Ma & Boylan-Kolchin (2010, see their fig. 3) obtained a mean merger rate per halo per unit redshift that is nearly independent of z and only weakly dependent on |$M_h$| (see also Genel et al. 2010; Rodriguez-Gomez et al. 2015). Considering major mergers (mass ratio of |$\gt 1/3$|⁠), the mean merger rate is |$dN_m/dz\approx 0.6$| and 0.3 for a halo mass of |$M_h=10^{13}~\mathrm{M}_{\odot }$| and |$M_h=10^{11}~\mathrm{M}_{\odot }$|⁠, respectively. Considering the halo mass evolution obtained by EVS model, we estimate a total of |$4\pm 1$| major mergers to take place from |$z=12$| to |$z=3.2$|⁠. This indicates that our scenarios outlined above of one, two or three major mergers lies within the expectation of |$\Lambda$|CDM.

The EVS model provides a simple description of how to relate the stellar masses of galaxies to their host dark matter halo masses. Specifically, within the EVS framework, we assume that ZF-UDS-7329 traces the most massive halo in the probed cosmic volume at each redshift. A caveat of this approach is that each redshift is treated independently, i.e. there is a lack of continuity in the model. Expanding on this, we compare our growth histories to the cosmological hydrodynamic simulation flamingo (Schaye et al. 2023; Kugel et al. 2023). Specifically, we use the box with a volume of (1 Gpc)|$^3$| and a baryonic particle mass of |$1\times 10^8~{\rm M}_\odot $|⁠. We select the 100 oldest galaxies with |$M_{\star }\gt 10^{11}~{\rm M}_\odot $| and |$\mathrm{SFR}\lt 20~{\rm M}_\odot ~{\rm yr} ^{-1}$| at the redshift of observation. We plot the median SFH of those galaxies in the left panel of Fig. 8. We find that the flamingo galaxies’ SFHs have a significantly lower formation redshift (⁠|$z \sim 5-6$|⁠) with respect to ZF-UDS-7329, even when adopting the rising prior. Similar results are find in other theoretical and numerical models (e.g. Lagos et al. 2024). While part of this difference could be due to the numerical resolution, it might also point towards a higher star formation efficiency at early cosmic times than what current numerical models are assuming. In particular, as put forward by Dekel et al. (2023, see also Li et al. 2024; Renzini 2023), the high densities and low metallicities at these early epochs could lead to a high star formation efficiency in the most massive dark matter haloes, because the free-fall time is shorter than |$\sim 1$| Myr, below the time for low-metallicity massive stars to develop winds and supernovae. In such a regime, the integrated star formation efficiency lies in the range of |$0.2{\!-\!}1.0$|⁠, consistent with our results presented here. A more detailed comparison, including a larger set of galaxies, is needed to further constrain the numerical models and the star formation efficiency.

Finally, cosmic variance could be an important factor as well. In the above EVS-based discussion, we assume that ZF-UDS-7329 is the most massive system in the PRIMER survey. However, we do not know how representative this survey area is of the universe as a whole, which can lead to biased estimates of the abundance of rare objects such as ZF-UDS-7329. Cosmic variance may reach |$\gtrsim 100{{\ \rm per\ cent}}$| of the stellar mass density at the high-mass end (see e.g. Chen et al. 2023; Xiao et al. 2024). This would roughly double our 1|$\sigma$| uncertainties for the implied |$f_{\star }$| (indicated by the error bar in Fig. 8), further reducing the tension.

In summary, while ZF-UDS-7329 is an early forming galaxy, we estimate a stellar fraction (i.e. integrated star formation efficiency) of |$f_{\star } \approx 0.5-1.0$| at |$z\approx 7-12$|⁠, in particular when considering the impact of mergers and cosmic variance. While this stellar fraction at early times is higher than lower-redshift estimates, we find that ZF-UDS-7329 potentially evolves to a stellar fraction of 0.1 by |$z=3.2$|⁠, which is consistent with the peak of the stellar-to-halo mass ratio at |$z\approx 0-3$|⁠. We therefore conclude that ZF-UDS-7329 formed very efficiently at high redshifts, but is consistent with direct, high-redshift observations and realistic galaxy formation physics.

5.2 Dark matter halo mass from the galaxy size

In addition to assess the integrated star formation efficiency based on the most massive halo in the survey volume as done in the previous section, we now perform an independent crosscheck of the halo mass estimation by using the theoretically motivated and observationally confirmed relation between galaxy sizes and their virial radius.

To first order, the physical origin of the galaxy size-mass relation and its evolution can be understood by considering that both dark matter and diffuse gas acquire angular momentum via tidal torques and mergers (Peebles 1969; White 1984; Porciani, Dekel & Hoffman 2002; Vitvitska et al. 2002). The specific angular momentum is often written using the dimensionless spin parameter:

where J is the total angular momentum, E is the total energy, G is Newton’s gravitational constant and M is the total mass. In the classical picture, diffuse gas acquires about the same amount of specific angular momentum as the dark matter, and conserves most of this angular momentum as it cools, collapses, and forms stars. The Mo, Mao & White (1998) model predicts

where |$j_d$| and |$m_d$| are fractions of baryon angular momentum and mass budget within halo in the central disc, |$f_c$| is a function of halo concentration, and |$f_R$| is a function that takes into account baryonic contraction of halo in response to halo formation. At |$z=3$|⁠, the expectation is that |$R_{\rm eff} \approx 0.01 R_{200}$|⁠, which is supported by observations based on abundance matching (Kravtsov 2013; Somerville et al. 2018; Shibuya et al. 2019).

Our measured half-light size in F444W of |$R_{\rm eff} = 1.27 \pm 0.02$| kpc translates to |$R_{200}\approx 127 \pm 25$| kpc (assuming an uncertainty of 20 per cent for the spin parameter). At |$z=3$|⁠, this translate to a dark matter halo mass of |$M_{\rm h} = (4/3) \pi R_{200}^3 200\rho _{\rm crit}(z=3.2) = (8.6_{-4.2}^{+6.2})\times 10^{12}~{\rm M}_\odot $|⁠, i.e. |$\log (M_h/{\rm M}_\odot) = 12.9\pm 0.3$|⁠. This implies a stellar fraction of |$f_{\star }=M_{\star }/(f_b \cdot M_h)=0.16_{-0.07}^{+0.15}$|⁠, which is consistent within the uncertainties with our estimate discussed in the previous section.

5.3 Environment of ZF-UDS-7329

The expectation in |$\Lambda$|CDM structure formation is that the first collapsed structures are highly clustered (Springel, Frenk & White 2006). Observationally, several examples of overdensities have been put forward. For example, JADES-GS-z14-0–the most distant galaxy known at |$z = 14.32$| – has indeed a neighbour (JADES-GS-z14-1 at |$z = 13.9$|⁠) and it seems likely that these galaxies are at least mildly associated in an extended large-scale structure (Carniani et al. 2024a). GN-z11 has nine galaxies within |$\sim 5$| comoving Mpc transverse with photometric redshifts consistent with |$z = 10.6$|⁠, which is consistent with GN-z11 being hosted by a massive dark matter halo of |$\approx 8\times 10^{10}~{\rm M}_\odot $| (Tacchella et al. 2023b). Furthermore, several galaxy overdensities and protoclusters have been identified in the Epoch of Reionzation (Higuchi et al. 2019; Helton et al. 2024b; Sun et al. 2024). Theoretically, semi-analytical models and numerical simulations show that that the cosmic SFR density at |$z\gt 5$| is dominated by galaxies in protoclusters (Chiang et al. 2017; Lim et al. 2024). Recently, Rennehan (2024) shows that galaxies such as ZF-UDS-7329 indeed have SFHs and formation epochs that are consistent with the most massive galaxies within protocluster cores of galaxy clusters at |$z\sim 2$|⁠.

As early forming galaxy, we expect ZF-UDS-7329 to be part of an overdensity, and possibly protocluster. We find that ZF-UDS-7329 has two companions at redshift |$z=3.2$|⁠. These galaxies are ZF-UDS-7542 (MSA ID 13079 from programme #2565; Schreiber et al. 2018; Nanayakkara et al. 2024) and RUBIES-UDS-46261 from programme #4233⁵ (de Graaff et al. 2024b), and lie within a projected distances of 12 and 6.5 arcsec, respectively, corresponding to 90 and 50 kpc. While RUBIES-UDS-46261 is star forming, ZF-UDS-7542 is a massive, post-starburst galaxy (Nanayakkara et al. 2024), which suggests a connection with the massive quiescent nature of ZF-UDS-7329 itself. The presence and proximity of another massive quiescent galaxy in the vicinity of ZF-UDS-7329 suggests a link between the quiescent nature of these galaxies and their surrounding environment. Furthermore, it supports the idea that ZF-UDS-7329 might have undergone a merger-rich history. Other massive, quiescent galaxies have been reported to occupy overdensities (de Graaff et al. 2024a); however, the exact nature of the environment surrounding ZF-UDS-7329 requires further investigation (de Graaff et al., in preparation).

5.4 Metallicity: |$\alpha$|-element enhancement

The best-fitting metallicity values for ZF-UDS-7329 can only be possible for its stellar age if there is some deviation from a solar composition, because some of the solar chemical abundances are due to enrichment processes that act on long time-scales (e.g. Maiolino & Mannucci 2019), whereas the SFH of ZF-UDS-7329 is necessarily short. We may derive loose bounds on [|$\mathrm{\alpha /Fe}$|] using the simplified relation |$[Z/\mathrm{H}] = [\mathrm{Fe/H}] + A[\mathrm{\alpha /Fe}]$| introduced in Trager et al. (2000), where A is a constant dependent on the detailed abundance pattern (= 0.93 if all |$\mathrm{\alpha }$|-elements are varied identically). Taking the C3K best-fitting value |$[Z/\mathrm{H}] = 0.11$| to be genuine, we find |$[\mathrm{Fe/H}] \approx -0.1$| to |$-0.3$|⁠, for values of |$[\mathrm{\alpha /Fe}] = 0.2$| to 0.5. Such results would be concordant with studies of massive quiescent galaxies at |$z \sim 0.7$| from the LEGA-C survey (van der Wel et al. 2016; Beverage et al. 2021). Furthermore, both chemical evolution models (Thomas et al. 2005), and analyses of early-type galaxies from the Sloan Digital Sky Survey (de La Rosa et al. 2011), have provided strong evidence for a correlation between short star formation time-scales and |$[\mathrm{\alpha /Fe}]$|⁠. We adopt the relation of Thomas et al. (2005), |$[\mathrm{\alpha /Fe}] \approx 1/5 - 1/6\log _{10}(\Delta t)$|⁠, which assumes a Gaussian SFH of FWHM |$\Delta t$|⁠. Taking |$\Delta t$| to be the time between |$16 {{\ \rm per\ cent}}$| and |$84 {{\ \rm per\ cent}}$| of the mass forming, we find |$[\mathrm{\alpha /Fe}] \approx$| 0.3 from our best-fitting SFHs. Finally, Carnall et al. (2024) fit a higher-resolution version of the spectrum using the absorption line fitter (alf) code (Conroy & van Dokkum 2012; Conroy et al. 2018) to obtain |$[\mathrm{Mg/Fe}] = 0.42^{+0.19}_{-0.17}$|⁠, consistent with our estimations. Given our treatment has assumed a solar abundance pattern throughout, how do these results curtail our confidence in the inferred population parameters?

In Fig. 5, we display the spectra of synthetic |$\alpha$|-enhanced SSPs based on the MILES library from Vazdekis et al. (2015) at two separate ages, convolved to the same resolution and binning as our data and offset for clarity. The SEDs demonstrate that a 1.0 Gyr-old population with |$[\mathrm{\alpha /Fe}] = 0.4$| can appear similar to a 1.5 Gyr-old one with a solar composition. This effect is most prominent at precisely the ages between |$1 \!-\! 1.5 \, {\rm Gyr}$|⁠, with the TiO band at |$7150 \, \mathrm{\mathring{\rm A}}$| becoming particularly sensitive to the degree of |$\alpha$|-enhancement. This highlights that considering |$\alpha$|-element enhancement is important to derive more accurate stellar ages for quiescent galaxies. It should be noted, however, that these |$\alpha$|-enhanced SSP templates do not self-consistently account for the effect of [|$\alpha$|/Fe] on the stellar isochrones, which may change the lifetimes on the stellar main sequence, and hence the colour of the spectra. Thanks to recent advancements (Park et al. 2024a), this will however become possible to self-consistently model.

5.5 Possible TP-AGB contribution

The importance and contribution of thermally pulsing asymptotic giant branch (TP-AGB) stars has been debated for decades (Maraston 2005; Conroy et al. 2009; Kriek et al. 2010). The contribution peaks in a stellar age range of |$0.2{\!-\!}2$| Gyr, but large uncertainties in the modelling of the TP-AGB phase existed due to its double-shell burning regime, leading to instabilities and ‘thermal pulses’ on short time-scales alongside a strong mass loss. Since we are fitting the spectrum until |$\sim 1~\mu$|m, where molecular features are present, TP-ABG could contribute to or even dominate the emission. Lu et al. (2024) suggest the presence of a strong TP-AGB component in a galaxy at |$z=1.08$|⁠, which seems to have a similar stellar age and metallicity as ZF-UDS-7329, though possible a longer formation period. Depending on the TP-AGB models used, the authors find that the inferred ages of the delayed-|$\tau$| model SFH can change by several 100 Myr. While we cannot exclude that a similar effect on our inferred ages, we note that we do not fit the MILES model beyond |$\sim 0.7~\mu$|m due to resolution and find consistent age constraints with the C3K models (Table 1). This is probably the case because the |$D_n4000$| strength together with the age of the Universe derive our stellar age constraint.

6 CONCLUSIONS

Our analysis of the JWST data for ZF-UDS-7329 provides new insights into the formation and evolution of massive quiescent galaxies in the early Universe. ZF-UDS-7329, with a stellar mass of |$M_{\star }\approx 10^{11.4}~{\rm M}_\odot $| and a sSFR of |$\approx 0.03~{\rm Gyr} ^{-1}$| at |$z=3.2$|⁠, has a strong 4000 Å break, implying an extremely high formation redshift. There are only a few spectroscopically confirmed massive quiescent galaxies at |$z \gtrsim 3$| with such an extreme SFH, with most appearing as post-starbursts (D’Eugenio et al. 2020), or having quenched at |$z \sim 4$| (Nanayakkara et al. 2024). Therefore, ZF-UDS-7329 presents an interesting object for understanding the complexities of early galaxy formation, and has been used for challenging cosmology (Glazebrook et al. 2024) and claiming extreme galaxy formation physics (Carnall et al. 2024).

We extended the analysis by Glazebrook et al. (2024) to include different SFH priors, stellar libraries (MILES versus C3K), metallicity, and IMF assumptions. We re-reduce the NIRSpec/PRISM spectrum with the NIRSpec GTO pipeline (e.g. Curtis-Lake et al. 2023; Curti et al. 2023, 2024; Bunker et al. 2024). Our findings show that the spectrum of ZF-UDS-7329 is consistent with a range of SFHs, suggesting stellar ages between 1.3 and 1.8 Gyr. This range reflects uncertainties but also indicates the robustness of our methodology in constraining the age and formation history of this galaxy. We also discuss an enhanced [|$\alpha$|/Fe] abundance, which we argue could reduce the real age by a |$\sim 200$| Myr. Options for fitting SPS models with |$\alpha$|-enhanced libraries are currently limited, but upcoming models will hopefully shed light on this hypothesis.

Importantly, by employing a physically motivated rising SFH prior, which tracks dark matter accretion histories but allows for significant deviation from this, we derived a formation history for ZF-UDS-7329 that aligns well with direct high-redshift observations, indicating SFRs and sSFRs more consistent with known galaxies at |$z\gt 10$| (see Fig. 4). We find that ZF-UDS-7329 is an extremely early forming galaxy, consistent with a stellar fraction (i.e. integrated star formation efficiency) of |$f_{\star } \approx 0.5-1.0$| at |$z\approx 7-12$|⁠, in particular when considering the impact of mergers (Fig. 8). We stress that we infer the stellar fraction from the EVS model, where we assume that ZF-UDS-7329 traces the most massive halo in the probed cosmic volume at each redshift. A caveat of this approach is that each redshift is treated independently; that is, there is a lack of continuity in the model. Nevertheless, this stellar fraction at early times is higher than lower-redshift estimates, but we find that ZF-UDS-7329 potentially evolves to a stellar fraction of 0.1 by |$z=3.2$|⁠, which is consistent with the peak of the stellar-to-halo mass ratio at |$z\approx 0-3$|⁠. We therefore conclude that ZF-UDS-7329 formed very efficiently at high redshifts, but does not necessitate unseen galaxies at redshift |$z\gt 10$|⁠.

The presence of a colour gradient, as observed in NIRCam imaging data, suggests a complex morphology with an older, quiescent bulge component and a younger disc component (Fig. 7). This supports the notion that ZF-UDS-7329 may have undergone mergers, contributing to its mass assembly and extended SFH.

Our analysis also highlights the significance of the environment in shaping galaxy evolution. ZF-UDS-7329’s possible association with an overdensity or protocluster environment at |$z=3.2$| is consistent with theoretical predictions of early structure formation in the |$\Lambda$|CDM cosmology: massive, early forming galaxies are high-|$\sigma$| peaks and therefore highly clustered.

In summary, ZF-UDS-7329 represents an early forming galaxy. Its efficient star formation at high redshifts, coupled with its significant stellar mass and complex morphology, provides a valuable benchmark for future studies of galaxy evolution in the early Universe. Specifically, upcoming high spectral resolution JWST/NIRSpec IFU observations of ZF-UDS-7329 (JWST Proposal Cycle 3, ID #5069) will constrain its dynamical mass, chemical abundances, gas properties, and spatially resolved SFH, helping us to test IMF assumptions, to constrain star formation time-scales, and to look for supermassive black hole feedback induced quenching.

ACKNOWLEDGEMENTS

We thank the referee for their constructive comments that helped to improve the paper. We are grateful to Andrea Ferrara, Avishai Dekel, and Andrey Kravtsov for insightful discussions. We thank Joop Schaye, Matthieu Schaller, Roi Kugel, and John Helly for creating the FLAMINGO project and providing us with the data. This work has been conducted as part of a Part III Project at the Cavendish Laboratory. CT thanks for insightful feedback from Didier Queloz. ST acknowledges support by the Royal Society Research Grant G125142. FDE acknowledges support by the Science and Technology Facilities Council (STFC), by the ERC through Advanced Grant 695671 ‘QUENCH’, and by the UKRI Frontier Research grant RISEandFALL. KG and TN acknowledge support from Australian Research Council Laureate Fellowship FL180100060. This research was supported in part by grant NSF PHY-2309135 to the Kavli Institute for Theoretical Physics (KITP). The NIRCam cutouts presented herein were retrieved from the Dawn JWST Archive (DJA). DJA is an initiative of the Cosmic Dawn Center (DAWN), which is funded by the Danish National Research Foundation under grant DNRF140. These data were reduced using a combination of the pipeline jwst and grizli.⁶ NIRSpec spectra from the RUBIES survey were obtained from the DJA, and reduced using the pipeline jwst and msaexp.⁷ This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K00042X/1, ST/P002293/1, ST/R002371/1, and ST/S002502/1, Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

REFERENCES

APPENDIX A: NIRSpec REDUCTION COMPARISON

We show the reduced spectrum adopted in this paper compared to that of Glazebrook et al. (2024) in Fig. A1.

APPENDIX B: ppxf SPECTRAL FITTING

B1 Methodology

The full spectrum is fitted using the ppxf method (Cappellari 2023) with a library of SSP templates that are a combination of the synthetic C3K model atmospheres (Conroy et al. 2019) and MIST isochrones (Choi et al. 2016) with solar abundance patterns and a Salpeter IMF. In particular, the scheme minimizes the function |$\chi ^2 + R \, \beta$|⁠, where |$\chi ^2$| are the weighted model residuals and |$\beta$| is a function quantifying the smoothness of the weights (i.e. their second derivative) with adjustable coefficient R. SSPs older than the age of the Universe at |$z=3.2$| were excluded, leaving a logarithmic grid of 86 ages from |$10^{5.0}$| to |$10^{9.2}$| yr, across 12 metallicities, |$\log _{10}(Z_\mathrm{SSP}/{\rm Z}_\odot) = -2.5$| to 0.5. While noise could bias the fit to younger ages, it is reasoned that allowing older SSPs would not be suitable given the extreme SFH of this galaxy.

Although the spectrum has not been directly calibrated, it agrees reasonably well with the photometric uncertainties, so a multiplicative polynomial is not applied to the fit. This is to ensure the dust attenuation, which is included using the Calzetti et al. (2000) curve, could be properly constrained, as our current pipeline does not fit for photometry. The velocity dispersions of the kinematic components are allowed to vary from 1 to 1000 km s⁻¹, but it should be noted that the PRISM disperser resolution is far too low to constrain this quantity. The SSP weights are additionally assessed using bootstrapping of the initial ppxf best fit. Here, the fitting is repeated 100 times, with artificial noise from random samples of the original best-fitting residuals added to the spectrum at each iteration. In all cases, second-order regularization is used with the ppxf regularization parameter set to a (modest) value of |$R=5$|⁠.

B2 Results

The inferred SSP light-weights and mass-weights after the bootstrapping procedure are shown in Fig. B1. We find the vast majority of stellar mass is formed over 1 |${\rm Gyr}$| prior to observation, with both the median mass-weighted and light-weighted age |$\simeq 1.6$| Gyr. The effective dust extinction is found as |$\tau _\mathrm{V} = 0.35$|⁠. While most of the inferred population is supersolar, there is a marked subsolar solution, accounting for |$\sim 10 {{\ \rm per\ cent}}$| of the total flux, that vanishes at greater regularization (⁠|$R \sim 20$|⁠). The mean mass-weighted metallicity of [Z/H]  = 0.07 is consistent with that obtained from full-spectrum fitting with prospector (see the next section) and implies some degree of template mismatch to be consistent with such an old age. We find a high stellar mass, |$\log _{10}{(M_\star /{\rm M}_\odot)} = 11.38$|⁠, |$\sim$| 0.2 dex greater than with prospector, which is expected given the usage of the bottom-heavier Salpeter IMF. The best-fitting spectrum and residuals are plotted in Fig. 2.

APPENDIX C: SFH PRIOR COMPARISON

We show that the results from the prospector models with a flat versus ‘rising’ continuity prior are largely consistent. The posteriors for six key parameters are presented in Fig. C1; the tight constraint on the median age arises from the rising model placing almost all the stellar mass in a single time bin.

APPENDIX D: MOCK SSP TESTS

We quantify the effect of changing the SFH prior on the derived stellar ages by fitting mock SSP spectra. Two SSPs are generated, one with age |$= 1.5~{\rm Gyr}$|⁠, and one with age |$= 1.9~{\rm Gyr}$|⁠. The other parameters are identically chosen as |$\log _{10}{(M_\star /{\rm M}_\odot)} = 11.3$|⁠, [Z/H]|$=0.01$|⁠, |$\tau _{1,V}=\tau _{2,V}=0.3$|⁠, |$n=-0.1$|⁠, |$z=3.2$|⁠, and no nebular emission, to roughly mimic the results from our best fits of ZF-UDS-7329. The spectra are binned to the same resolution as our data, and convolved with our estimated LSF. In addition, we add Gaussian noise according to the same (fractional) noise vector.

Fig. D1 shows the best-fitting SFHs produced by our models. For the younger SSP, the flat prior overestimates the age by |$\sim 250~{\rm Myr}$|⁠, while the rising prior underestimates it by |$\sim 100~{\rm Myr}$|⁠. Both posteriors on the median age are |$\gtrsim 3\sigma$| inconsistent with the true value. This is compounded by the use of time binning, for which the placement of bin edges will also influence the results. Hence, the posterior distribution on the age should not be taken at face value for these non-parametric models. For the older SSP, the flat prior now correctly recovers the age, but the rising one gives a MAP SFH virtually unchanged from the previous test, hence underestimating the age by |$\sim 500~{\rm Myr}$|⁠. In all cases, the other stellar parameters are recovered within |$1\!-\!2 \sigma$|⁠, showing the uniqueness of this problem.

APPENDIX E: SFH CUT-OFFS

We present the best-fitting spectra for all our prospector models with no star formation allowed past |$z_0 = 10, 20, \infty$| in Fig. E2, with the corresponding SFHs in Fig. E1. These demonstrate that forcing a |$\sim 300 \, {\rm Myr}$| younger poor fit does not significantly disagree with the data, with most values remaining within the |$5 {{\ \rm per\ cent}}$| error margin. From oldest to youngest age, we find reduced |$\chi ^2$| values (calculated with |$j_{\mathrm{spec}}=1$|⁠) of 0.43, 0.44, 0.49 for the MILES models and 0.37, 0.37, 0.38 for the C3K models (which cover a larger spectral range). We also note that the inferred noise inflation factor, |$j_{\mathrm{spec}}$|⁠, is relatively unchanged in these tests.

APPENDIX F: IMF TESTS

Top-heavy and bottom-heavy IMFs were implemented in the C3K prospector model as described in Section 3.5. The posteriors for six key parameters are shown in Fig. F1 alongside the fiducial values, which used a Kroupa (2001) IMF. The results demonstrate that while a top-heavy IMF gives a |$\sim 0.1$| dex reduction in |$M_{\star \mathrm{,surv}}$|⁠, it is unable to reduce the total formed mass, |$M_\star $|⁠. All other population parameters are largely unchanged by the IMF variation.