RUBIES: Evolved Stellar Populations with Extended Formation Histories at z ∼ 7–8 in Candidate Massive Galaxies Identified with JWST/NIRSpec

The RUBIES targets in the EGS were selected from public JWST/NIRCam data from the Cosmic Evolution Early Release Science Survey (CEERS, JWST-GO-1345; PI Finkelstein; Bagley et al. 2023; Finkelstein et al. 2023), which provides imaging in seven bands (F115W, F150W, F200W, F277W, F356W, F410M, and F444W). Additionally, we use archival imaging in seven different filters (F435W, F606W, F814W, F105W, F125W, F140W, and F160W) from the Hubble Space Telescope (HST) from the CANDELS survey (Grogin et al. 2011; Koekemoer et al. 2011).

The image mosaics, with a pixel scale of 004 pixel⁻¹, are publicly available in the DAWN JWST Archive (DJA; version 7.2) and were reduced using grizli (Brammer 2023a). Fluxes are measured from point-spread-function-matched images in circular apertures with a radius of 016 and then Kron-corrected to the total flux, as described in Weibel et al. (2024).

Each target was observed for 48 minutes in the Prism/CLEAR mode and the G395M/F290LP mode, using a standard three-shutter slitlet and three-point nodding pattern. The spectra are reduced, combined, and extracted using the JWST calibration pipeline version 1.12.5 (Backhaus et al. 2024) and msaexp (Brammer 2023b), corresponding to the version 2 reduction on DJA. To account for wavelength-dependent flux calibration that is not yet captured well by the pipeline, we renormalize the Prism spectrum to match the NIRCam photometry using a dynamic high-order polynomial as described in Section 5. The G395M spectrum is subsequently rescaled by this polynomial. We find that there is a small systematic offset (≈10%–20%) in the flux calibration between the Prism and G395M spectra, which has recently also been reported by D'Eugenio et al. (2024). To determine the offset, we fit the [O iii] doublet in both the Prism and G395M spectra (Section 4) and hence rescale the full G395M spectrum by the ratio of the two to match the flux calibration of the Prism spectrum.

This Letter focuses on the three targets in this sample that exhibit clear Balmer breaks at z_spec = 6.6–8.4, found via visual inspection of the 2D and 1D spectra. These objects are extremely red in F277W–F444W (Figure 1), the characterization of red objects being one of the core targeting criteria of RUBIES. In what follows, we establish the detections of Balmer breaks and the broad emission lines, which are the two key characteristics of this sample.

The emission-line widths of the Balmer Hβ and [O iii] lines are modeled using the G395M spectra. Prior to fitting, we rescale the error spectrum, as we find that the noise estimated by the msaexp pipeline is underestimated compared to the observed pixel-to-pixel variation. To estimate the rescaling factor, we select the region outside of the Hβ and [O iii] emission lines and calculate the ratio between the pixel-to-pixel variance and the median of the error spectrum; this results in a rescaling factor in the range 1.3–2.0.

Our fitting methodology broadly follows that described in Wang et al. (2024a): given the compact morphology of the sources (Baggen et al. 2023), we use a point-source line-spread function (LSF) from de Graaff et al. (2024a) and the Markov Chain Monte Carlo ensemble sampler emcee (Foreman-Mackey et al. 2013) to estimate the posteriors. Crucially, our method explicitly includes a systematic uncertainty on the model LSF and therefore provides realistic measurement uncertainties for marginally resolved lines.

We select the region around the emission-line complex (±0.2 μm) and begin by masking the Hβ line (a region of 7000 km s⁻¹) to fit only the [O iii] doublet to estimate the narrow-line width and redshift of the source. We use a single Gaussian component, the dispersion of which is the same for both lines of the [O iii] doublet and constrained to be in the range σ_narrow ∈ [0, 500] km s⁻¹ using a uniform prior. The continuum is fit with a first-order polynomial. We then explore whether the [O iii] lines show evidence for a second broader component by fitting a two-component Gaussian model, where the width of the second Gaussian is constrained to be larger than that of the narrow component. The Bayesian information criterion (BIC) is computed to compare the two models (Jeffreys 1961; Liddle 2004). We find evidence (ΔBIC = 13) for a broader component in the [O iii] line for RUBIES-EGS-55604 with a dispersion ${\sigma }_{\mathrm{broad},\ \mathrm{OIII}}={248}_{-41}^{+61}\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and no evidence for a second component in the other two sources.

Next, we include the Hβ line in the fitting. We first fit a single Gaussian component for the Hβ line. Given the limited signal-to-noise ratio (S/N) of the data, the width of this component is tied to that of the [O iii] lines. We then include a second broad component for the Hβ line, where σ_broad ∈ [500, 2500] km s⁻¹. For RUBIES-EGS-55604 we find a blueshifted absorption feature in the Hβ line, which is modeled as an additional Gaussian component with a velocity offset Δv ∈ [0, 1000] km s⁻¹ and dispersion σ_outflow ∈ [0, 1000] km s⁻¹. Such an absorption feature in Balmer and/or He lines has been seen in recent studies of AGNs as well (Wang et al. 2024a; Kocevski et al. 2024; Matthee et al. 2024). The multicomponent fits and the residuals are shown in Figure 2.

Again, using the BIC to compare the models, we find that RUBIES-EGS-49140 and 55604 have unambiguous broad components in the Hβ line (ΔBIC > 100). The third source, RUBIES-EGS-966323, is fainter and at higher redshift, and it has only a marginal detection of a broad component (ΔBIC = 1.94; broad-line flux detected at 3.5σ, F_Hβ,broad = (2.1 ± 0.6) × 10⁻¹⁸ erg s⁻¹ cm⁻²). For this source we also perform an independent fit of the Prism spectrum, finding similarly weak evidence for a two-component model and posteriors that are consistent with the fit to the G395M spectrum.

Finally, we fit the Balmer Hα emission line, which is available in the G395M spectrum for RUBIES-EGS-49140 and only in the Prism spectrum for RUBIES-EGS-55604 (for RUBIES-EGS-966323 Hα is redshifted out of the wavelength range accessible with NIRSpec). Hα is modeled with a narrow and broad component and also fit simultaneously with the [N ii]_{λ λ6549,6585} doublet. We assume that [N ii] is narrow and set the width to be equal to that of the narrow Hα line, and we fix the flux ratio of the doublet to 1:2.94. Moreover, informed by the fits to Hβ and [O iii], we constrain the dispersion of the narrow lines to <150 km s⁻¹. For RUBIES-EGS-49140 we also find a blueshifted absorption feature in the Hα line, which is fit in the same manner as described previously for the Hβ line of RUBIES-EGS-55604. We note that these fits are only used to obtain an estimate of the broad-line flux in order to calculate the broad Hα equivalent width (EW) in Section 6.1.

Reverberation mapping is a method employed to determine the size of the broad-line region by analyzing the time delay between variations in the AGN continuum and corresponding changes in the broad permitted lines (e.g., Blandford & McKee 1982). It has enabled the establishment of empirical correlations between the size, line luminosities, and widths in the nearby Universe (e.g., Kaspi et al. 2000; Landt et al. 2013). These correlations facilitate the estimation of black hole masses from single-epoch observations.

Following Assef et al. (2011), we estimate the black hole mass based on Hβ and luminosity at rest 5100 Å as

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{log}({M}_{\mathrm{BH}}/{M}_{\odot })=0.895\\ +0.529\mathrm{log}\left(\displaystyle \frac{{L}_{5100}}{{10}^{44}\,\mathrm{erg}\,{{\rm{s}}}^{-1}}\right)+2.0\mathrm{log}\left(\displaystyle \frac{{\mathrm{FWHM}}_{H\beta }}{\mathrm{km}\,{{\rm{s}}}^{-1}}\right).\end{array}\end{eqnarray} \tag{ 1 }$$

Hβ is used since Hα is not available for all sources and suffers from uncertainty due to the blending of the [N ii] doublet with the broad component of Hα. L₅₁₀₀ is calculated from the intrinsic AGN spectrum inferred from the SED models.

We adopt this relation to better illustrate the uncertain in the black hole masses, as the different SED models predict a wide range of AGN luminosities. All line luminosities here are dereddened using the dust attenuation inferred from SED fitting. While subject to systematic uncertainties, this represents our best estimate of the dust content in the absence of a Balmer decrement. We note, however, that additional systematic uncertainties are likely introduced by the application of these methods at higher redshifts and in different physical conditions than where they are calibrated (e.g., Yue et al. 2024).

The available JWST and HST photometric data are jointly fitted with the full NIRSpec/Prism spectrum within the Prospector inference framework (Johnson et al. 2021), following Wang et al. (2024a). The setup for the stellar populations is detailed in Wang et al. (2024b). In brief, we adopt the MIST stellar isochrones (Choi et al. 2016; Dotter 2016) along with the MILES stellar spectral library (Sánchez-Blázquez et al. 2006) in FSPS (Conroy & Gunn 2010). Star formation history (SFH) follows the nonparametric Prospector-α description, characterizing mass formation in seven logarithmically spaced time bins (Leja et al. 2017), with a weakly informative prior assumption of a rising SFH from Wang et al. (2023). We also impose a joint prior on stellar mass and stellar metallicity, as introduced in Leja et al. (2019). This prior is a Gaussian approximation of the relationship measured from the Sloan Digital Sky Survey (SDSS; Gallazzi et al. 2005), but the confidence intervals are widened by a factor of 2 to account for potential systematics or redshift evolution. Dust attenuation is described by a Calzetti et al. (2000) curve with a flexible power-law slope (Noll et al. 2009). The fraction of starlight permitted outside the dust screen is allowed to vary, a nonzero fraction of which suggests the presence of blue stars possibly existing outside the dust or having created holes within it. Dust emission is incorporated based on the model by Draine & Li (2007).

Model spectra are convolved with the NIRSpec/Prism instrumental resolution curve, tailored for a point-source morphology using msafit (de Graaff et al. 2024a). We account for wavelength-dependent slit losses by scaling the normalization of the spectrum to match the photometry through a seventh-order polynomial calibration vector applied during the fitting process.

Emission lines are fit using a one-component Gaussian model and, where applicable, a two-component Gaussian model to account for the narrow and broad components. This approach means that the emission lines are not interpreted physically, merely described and included in the modeling of the photometry. To prevent the likelihood from being influenced by residuals from non-Gaussian line kinematics, we enforce a 10% error floor in the spectroscopic data, which is higher than the conventional 5% error floor applied to photometry. However, we introduce a multiplicative noise inflation term, with a prior range from 0.5 to 5. Typically, this value is found to be around 1.5, suggesting that an additional 50% inflation of the random noise produces a good fit. The posteriors are sampled via the dynamic nested sampler dynesty (Speagle 2020).

The composite galaxy and AGN model is presented in Wang et al. (2024a), in which we approximate the direct UV–optical emission from an AGN accretion disk as piecewise power laws, with varying normalization. The slopes are fixed to the best-fit values in Temple et al. (2021), which are calibrated to the median colors of quasars in SDSS, the UKIRT Infrared Deep Sky Survey, and unWISE. The light from the accretion disk experiences the same dust attenuation as the stars but is additionally reddened by a separate dust attenuation curve modeled as a power law with varying normalization and shape. In all, the AGN continuum is controlled by five free parameters (the AGN-to-galaxy flux ratio and four dust attenuation parameters, two of which are shared with the galaxy light).

The key uncertainty in the interpretation of these objects is the source of the red continuum: is it powered by AGNs, stars, or a mixture of both? The implied total stellar masses and formation histories are a strong function of this decomposition. A very red, luminous rest-frame optical stellar solution has a high M/L and an extended formation history, whereas a blue and/or less luminous stellar optical continuum can instead be fit with a flat or rising SFH and relatively little stellar mass. The Balmer break is a key constraint here, as this, alongside the resolved sizes in the blue (Baggen et al. 2023), suggests that the continuum blueward of rest 4000 Å is dominated by stars.

The central question is therefore the origin of the continuum redward of the Balmer break. We consider three models that bracket the possible range of inferred stellar and AGN contribution to the observed rest-optical continua.

5.3.1. Maximal Stellar Contribution

5.3.2. A Mixture of AGN and Stellar Contributions

5.3.3. Minimal Stellar Contribution

The models with minimal and medium stellar contribution infer that the AGN accretion disk dominates the rest-optical, while starlight dominates the rest-UV and also the Balmer break (by necessity). The AGN emission is heavily attenuated by dust, resulting in the observed red color. The transition from stellar- to AGN-dominated light can in principle be corroborated by the wavelength-dependent morphology—an unresolved morphology is more indicative of AGN-dominated light. The presence of broad emission lines is a suggestion of AGN activity, and if these broad lines are AGN powered as opposed to, e.g., a stellar-driven outflow, they can put a soft lower limit on the AGN continuum contribution. This is because the AGN-powered broad Hα EWs are typically observed to be ≲1000 Å (more discussion of this point in Section 7); decreasing the AGN contribution to the continuum underneath these broad lines requires correspondingly higher AGN EWs, potentially outside of the previously observed range, which would require different and new AGN physics.

We find unambiguous broad Hβ lines with FWHM > 2500 km s⁻¹ in RUBIES-EGS-49140 and 55604, which would be unusual for stellar-driven outflows (see Veilleux et al. 2005, for an overview; or see, e.g., Heckman et al. 2015; Wang et al. 2020, for local starbursts). Stellar-driven outflows would typically also be seen in the forbidden lines, which are narrow in these objects—though it is possible for the outflowing gas to be sufficiently dense that it is not luminous in the forbidden lines. A broad component is marginally detected in RUBIES-EGS-966323 (Section 4). Additionally, we compare the observed Hα EW to previous samples of LRDs. RUBIES-EGS-49140 and 55604 show large observed EWs of >600 Å. We reestimate the EW using the AGN continuum from the SED fits. Assuming that the broad Hα originates from an AGN, the intrinsic EW must be larger in the models where considerable galaxy light contributes to the continuum near Hα. As seen from Figure 4, the increased EWs at the AGN continuum inferred from the medium stellar model are still roughly consistent with the distributions measured from a grism-selected LRD sample at lower redshifts (Matthee et al. 2024), although at the higher end.

Zoom In Zoom Out Reset image size

We can also predict Hα fluxes from the galaxy, given the inferred stellar population parameters, which provide additional leverage on interpreting the continuum emission. For all three SED models, the predicted Hα fluxes are at least 1 order of magnitude lower than the observed flux. This means that the SFRs in all models are much lower than the observed total Hα flux, i.e., none of the stellar components in these models can be ruled out by requiring star formation that produces more Hα flux than is observed. Conversely, the observed large Hα EW suggests an abundance of ionizing photons, the evidence of which is not obvious from the spectra in the UV–optical.

Furthermore, we observe an emission line at ∼3869 Å in all three objects. One possibility is He I 3889, given that a strong He I emission line has been reported in spectra of LRDs (Wang et al. 2024a) and AGNs more generally (Riffel et al. 2006; Landt et al. 2008). More likely it corresponds to [Ne iii], in which case it would suggest a high ionization state, particularly considering the weak (nondetected) [O ii] λ3727 emission. A strong [Ne iii]/[O ii] ratio can also be indicative of AGN activity (Backhaus et al. 2024).

Finally, we cross-check the morphologies of RUBIES-EGS-55604 and 966323 as studied in Baggen et al. (2023). RUBIES-EGS-55604 (L23-38094) is among the most compact in L23, which is unresolved in the long-wavelength filters (<0.22 kpc) but clearly exhibits two components in F115W and F150W. RUBIES-EGS-966323 is also unresolved in F444W, and likely resolved in F200W, although the F200W data are too faint to have a robust size estimate. Preliminary analysis on RUBIES-EGS-49140 reveals the same trend. A thorough study on the morphologies will be presented in J. Baggen et al. (in preparation). These findings again indicate the presence of a transition from stellar to AGN light near the spectral break, thus favoring the minimum/medium stellar models.

The inferred SFHs based on the three SED models are shown in Figures 3(a)–(c). All models have significant SFHs extending over hundreds of millions of years, necessary to produce the evolved stellar populations responsible for the Balmer breaks. However, there is some variation in the timescales of star formation and significant variation in the amplitudes. The medium and maximal stellar models require significant star formation at z = 10 and a recent decline in the SFR. The decline may be representative of a very early termination in star formation necessary to produce the very old galaxies observed at lower redshifts (Carnall et al. 2023; Glazebrook et al. 2024) or a more temporary mini-quenching event (Looser et al. 2024); however, we note that at least two of these Balmer breaks are substantially stronger than the Looser et al. (2024) break, implying a correspondingly longer period of quiescence (≳100 Myr) than a mini-quenching event.

Conversely, the minimal stellar models decrease the strength of the stellar Balmer breaks by assuming a very red AGN continuum underneath, and in this way they can infer purely rising SFHs for RUBIES-EGS-55604 and 966323. The signature of a recent decline in the SFH is still preserved in the brightest source, RUBIES-EGS-49140, which also has the strongest Balmer break among the three. This is an expected behavior—the less prominent Balmer breaks may be fit with a weaker break and steeper AGN continuum emission as the prior drives to minimize the stellar mass. However, it is worth emphasizing that this scenario hinges on a peculiar AGN continuum shape, which coincidentally aligns perfectly with that of the host galaxy to produce the observed spectral break. Alternatively, if not all the UV light is from stars (e.g., scattered light from the AGN; Greene et al. 2024), this would remove an important constraint pushing the inferred stellar ages younger and perhaps relax the stringent requirement on the galaxy and the AGN conspiring to create a spectral break by allowing an older stellar population.

One potential concern may arise from the age–metallicity degeneracy. We note that the data do not constrain the metallicity, as the posteriors roughly follow the prior distributions. This effectively means that we marginalize over metallicity, i.e., the uncertainty in metallicity is propagated into the rest of the inferred stellar population parameters.

Interestingly, the maximal stellar interpretation can be more easily connected with several massive, very old quiescent galaxies recently discovered with JWST at z = 3–5 (Carnall et al. 2023; de Graaff et al. 2024b; Glazebrook et al. 2024), as illustrated in Figure 5. Such a link may serve as an argument in favor of the maximal stellar model.

Zoom In Zoom Out Reset image size

We note, though, that the above connection is still hard to establish given the low number densities of these massive objects (a lower limit of ∼10⁻⁵ Mpc⁻³ for the three objects in this work and, e.g., ∼5 × 10⁻⁶ Mpc⁻³ for RUBIES-EGS-QG-1 estimated in de Graaff et al. 2024b). Nevertheless, the most straightforward prediction stemming from the discoveries of these quiescent galaxies at z = 3–5 is the presence of Balmer breaks at z = 7–8, suggesting very early and/or rapid formation. While the number density may be low, they are indicative of a class of objects with early and rapid formation.

One of the key results from L23 is that the high stellar masses in these objects suggest that they dominate not just the high-mass end of the pre-JWST UV-selected galaxy stellar mass function (Stefanon et al. 2021) but indeed the entire stellar mass budget at this cosmic epoch. We revisit this point by estimating the cumulative stellar mass density using the stellar masses inferred in this work. Specifically, we bin the three objects into two redshift ranges (5.5 < z ≤ 7.0 and 7.0 < z ≤ 8.5) and sum the mass in each bin. The cosmic volume is estimated by integrating between the redshift limits over the areal coverage of CEERS (88.1 arcmin²; Finkelstein et al. 2023). As in L23, we only consider the Poisson uncertainty and cosmic variance, and we neglect corrections for incompleteness, given that any correction for incompleteness would increase the inferred stellar mass densities.

We note that there is some evidence indicating that EGS exhibits an overdensity of red objects. So far, nothing similar to the sample of this Letter, in terms of high-z Balmer breaks in the available spectra and, more broadly, of objects as bright and red in our parent photometric catalogs, has been identified in other RUBIES fields (e.g., UDS, which has an area of 224 arcmin²; Weibel et al. 2024). A naive estimate of the number density might include this larger empty volume. However, as will become evident subsequently, this factor of ∼3 decrease does not affect any of the main conclusions. A subtle point, though, is that for objects displaying very red colors in broadband photometry, their spectra may not necessarily be similar. Defining volumes using photometry might therefore be suboptimal, but it perhaps can serve as a reasonable approximation here given the strong correlation between stellar masses and luminosity/flux.

Figure 6 shows the cumulative stellar mass density in these objects compared to rest-UV-selected objects. Importantly, rest-UV selection would fail to select any of the objects in this Letter, which are extremely red: the total mass density is thus inferred to be the sum of the two (independent) curves. Although stellar mass functions incorporating the JWST observations have been estimated (e.g., Harvey et al. 2024; Weibel et al. 2024), the unclear physical nature of the reddest objects makes it especially difficult to both estimate and interpret the high-mass end of the high-redshift stellar mass functions. We thus compare only to the pre-JWST UV-selected mass functions, as a test of what a rest-UV selection function may have missed.

Zoom In Zoom Out Reset image size

Assuming the medium stellar model, our analysis suggests that at z ∼ 7–8 this small population contains 20%–50% of the stellar mass density of the entire rest-UV-selected galaxy population, implying that ∼20%–50% of star formation at z > 7 occurs in the progenitors of these optically red objects and that they would dominate the high-mass end of the mass function. Adapting the minimal stellar model, these objects constitute a much smaller ∼1% of the cosmic mass budget and live in much more typical galaxies. Given that the Balmer break objects may have slowing or declining SFHs, whereas UV-selected galaxies typically have rising formation histories by nature, this fraction is likely to be even higher at higher redshifts.

For reference, the number density of our sample is more than an order of magnitude lower than the total number density of the rest-UV-selected sample, obtained by integrating the Schechter fits down to a stellar mass limit of 10⁸ M_⊙ (Stefanon et al. 2021). This reinforces that the Balmer break objects may host disproportionately high levels of past star formation compared to the general galaxy population.

We now put the above stellar mass density into context with theoretical predictions from the standard ΛCDM cosmology. Given the suite of cosmological parameters, the matter power spectrum describing the density contrast of the Universe on large scales can be specified. As gravitational collapse becomes nonlinear on smaller scales that are relevant for dark matter halos that host galaxies, higher-order statistics becomes necessary. However, under the assumption that the initial density fluctuations are Gaussian and small, reasonable analytic predictions of the halo distribution can be made (Press & Schechter 1974).

Following Boylan-Kolchin (2023), we adopt the Sheth & Tormen (1999) halo mass function, which is an extension of the Press & Schechter (1974) formalism by assuming ellipsoidal collapse instead of the more simplified spherical collapse. The comoving number density of halos above a given halo mass threshold, M_halo, is estimated as

$$\begin{eqnarray}&&n(\gt {M}_{\mathrm{halo}},z)={\int }_{{M}_{\mathrm{halo}}}^{\infty }dM\displaystyle \frac{dn(M,z)}{dM},\end{eqnarray} \tag{ 2 }$$

where d n(M, z)/d M is the number of dark matter halos of mass M per unit mass per unit comoving volume at redshift z (i.e., the halo mass function). Then, the comoving mass density in halos more massive than M_halo can be computed straightforwardly as

$$\begin{eqnarray}&&{\rho }_{m}(\gt {M}_{\mathrm{halo}},z)={\int }_{{M}_{\mathrm{halo}}}^{\infty }d{MM}\displaystyle \frac{dn(M,z)}{dM}.\end{eqnarray} \tag{ 3 }$$

From the above, we can obtain the corresponding statistics of galaxies by assuming

$$\begin{eqnarray}&&{M}_{\star }=\epsilon {f}_{b}{M}_{\mathrm{halo}},\end{eqnarray} \tag{ 4 }$$

where is the efficiency of conversion of baryons into stars and f_b is the cosmic baryon fraction. Here = 1 sets the stringent upper limit on the stellar content that a halo can have. We also consider two more likely values of = 0.1 and 0.32 (e.g., Behroozi et al. 2013; Moster et al. 2013).

As illustrated in Figure 7, only in the limiting case of = 1.0 does the expected number density of galaxies with the maximal stellar mass approximately align with the ΛCDM prediction, corresponding to halos with cumulative comoving number densities ≲10^−5.4 Mpc⁻³). For = 0.32 (0.10), the implied number density is ≲10^−7.1 (10^−9.5) Mpc⁻³, suggesting that the maximal stellar mass model results in objects that are unexpectedly massive. Furthermore, similar to findings in Boylan-Kolchin (2023), the implied star formation efficiency from our maximal stellar mass estimate lies at the extreme end of ΛCDM expectations. These tensions are alleviated if assuming instead the minimal and medium stellar mass models. For completeness, we also include the second set of lower number density, based on the CEERS+UDS volume, as open gray symbols in Figure 7(b). Additional systematic uncertainties in SED fitting, e.g., the initial mass function, which can likewise decrease the tension with the standard model, are presented in Wang et al. (2024c).

Zoom In Zoom Out Reset image size

While the Balmer lines are very broad, the forbidden lines are relatively narrow: the [O iii] line widths for the three galaxies are FWHM ∼ 160, 113, and 162 km s⁻¹ (Table 1). As gas is viscous, the forbidden lines likely originate in disks and—assuming that those disks trace the gravitational potential of the galaxies—can be used to derive dynamical masses. As discussed in, e.g., Förster Schreiber et al. (2014), van Dokkum et al. (2015), and Price et al. (2016), the conversion from line width to dynamical mass depends on the spatial distribution of the gas, the orientation of the disks, and the orientation of the slit with respect to the rotation axis. While the gas is certainly compact (based on 2D spectra, where [O iii] lines have about the same pointlike morphology as the Balmer lines), the latter two dependencies are both unknown.

Nevertheless, it is instructive to estimate the dynamical mass. Following van Dokkum et al. (2015), we have

$$\begin{eqnarray}&&{M}_{\mathrm{dyn}}=2\displaystyle \frac{{V}_{\mathrm{rot}}^{2}{r}_{e}}{G},\end{eqnarray} \tag{ 5 }$$

where r_e is the spatial scale, taken to be ∼100 pc (Baggen et al. 2023), and V_rot is the rotation velocity derived via

$$\begin{eqnarray}&&{V}_{\mathrm{rot}}=\displaystyle \frac{{\sigma }_{\mathrm{gas}}}{\alpha {\sin }^{-1}(i)}.\end{eqnarray} \tag{ 6 }$$

Assuming α = 0.8 (Rix et al. 1997; Weiner et al. 2006) and an inclination angle, i, of 45°, the dynamical masses for RUBIES-EGS-49149, 55604, and 966323 are ∼10^8.6, 10^8.3, and 10^8.6 M_⊙, respectively, all at the low end of the stellar mass estimates, and actually in some tension with the black hole mass estimates alone (discussed in the next section).

Higher masses are possible if the gas has a larger spatial extent than the stars (van Dokkum et al. 2015), though the [O iii] emission lines appear to have a highly compact morphology in the 2D spectra. Importantly, the LSF of NIRSpec is strongly dependent on morphology (up to a factor of ≈2 difference in resolution between a point source and uniformly illuminated slit; de Graaff et al. 2024a). Therefore, if the source were to have r_e ∼ 1 kpc, then σ_gas would decrease to <30 km s⁻¹; that is, by making the source larger, the dynamical mass remains similarly small. Low disk inclinations can also increase the dynamical mass. This is unlikely to be the case by chance, but it could result from selection effects; the red color selection may preferentially select obscured AGNs at particular down-the-barrel (i.e., face-on) orientations. However, given the significant inferred dust attenuation, a down-the-barrel orientation would be in some tension with the standard AGN model.

Figure 8 compares the stellar mass–black hole mass relationship inferred for these objects to the local relationships. In the optically red regime where our sample resides, the M/L–color relationship used for the local relationships (Zibetti et al. 2009) agrees well with the Prospector M/L–color relationship (Li & Leja 2022).

Zoom In Zoom Out Reset image size

Based on our minimal/medium stellar masses, these objects host massive black holes ∼10 × above the z = 0 M_⋆−M_BH relationship (Reines & Volonteri 2015). The minimal stellar mass model is well above the typical scatter in this relationship, while the medium model is within ∼2σ. These results align with previous findings on single objects (Kokorev et al. 2023; Furtak et al. 2024), on quasars (Stone et al. 2024), as well as on a z ∼ 5 AGN sample curated in Pacucci et al. (2023, which is found to be ∼10−100 × above the local relationship). At face value, this implies that these black holes are overmassive relative to their stellar components and that the stellar components must continue to grow over the next 13 Gyr in order to produce the relationship observed today. However, we caution that the selection bias, where the luminous and massive black holes with emission sufficiently dominating their host galaxy are more readily observed in a flux-limited survey (Lauer et al. 2007), is likely exacerbated by the high-redshift frontier.

Meanwhile, the maximal stellar model predicts an AGN continuum orders of magnitude lower than the galaxy continuum and hence implies much lower black hole masses using the L₅₁₀₀ relation, although, as noted above, it is difficult to know how to calculate the black hole mass when the SEDs are so abnormal. Alternatively, since this model effectively predicts starlight dominating over the entire observed wavelength range, it may be possible that these galaxies do not host AGNs at all, which is reasonable given the nondetection in X-rays (Ananna et al. 2024; Yue et al. 2024). This, then, would perhaps require some yet-to-be-understood mechanism to produce the broad emission lines—for example, broad emission lines from supernovae have been previously mistaken for quasar activity (Filippenko 1997; Aretxaga et al. 1999; Baldassare et al. 2016), though it remains unclear why supernovae would be consistently associated with the other spectral characteristics of these objects, and such an interpretation may be in some tension with the lack of variability seen in LRDs (Maiolino et al. 2024).

While it is instructive to put our objects in the context of the M_⋆–M_BH relationship, despite substantial uncertainty in both the black hole and galaxy mass, we note that the z ≳ 7 objects do not necessarily need to lie on the z = 0 relationship. Peng (2007) proposed that a large number of mergers lead to a statistical convergence process, and thus the slope of the M_⋆–M_BH relationship always becomes ∼1, regardless of the initial condition. In any case, the above scenarios all point to very different pictures of the early M_⋆–M_BH relationship and the subsequent evolution required to match the local scaling laws.

The maximal stellar model, if correct, would have a remarkable impact on the first billion years of galaxy evolution. It proposes that the cumulative stellar mass densities in these objects are comparable to those of all the UV-selected objects at these redshifts (Figure 6). Such early and efficient formation would be in tension with the standard assumption on the baryon-to-stellar conversion efficiency (Figure 7; see also Boylan-Kolchin 2023). Yet such an interpretation is consistent with the very old stellar populations observed in high-mass galaxies at z = 3–5 (Figure 5)—these require rapid, early assembly of ≳10^10.5 M_⊙ in stellar mass, followed by an early cessation in star formation. The SEDs of their progenitors would look much like these objects—though these objects are ∼10 × smaller in physical size than these later quiescent galaxies.

In contrast, the models with minimal and medium stellar contributions (which, in turn, imply larger fractional AGN contributions) are supported by the observations of broad emission lines and compact morphologies at redder wavelengths. They require the black hole continuum to be steeply rising in the rest-optical, such that a rapid transition from a stellar-dominated continuum to a black-hole-dominated continuum occurs. Furthermore, barring line-of-sight or other arguments, the small dynamical mass is in distinct tension with all models except the minimal stellar contribution. Meanwhile, the medium stellar model leads to an M_⋆–M_BH relation in less tension with (or requiring fewer mergers to become consistent with) the local relations (Figure 8).

It has long been suspected that the cores of the most massive galaxies in the local Universe formed in a spectacular early burst of star formation at z ≳ 5, followed by an immediate and permanent quenching. This hypothesis is underpinned by their high α-element abundances and their old, but unresolved, ages (e.g., Thomas et al. 2005; Conroy et al. 2014). Further evidence comes from the number density analysis of massive, dense galaxy cores at redshifts z = 0–3 (van Dokkum et al. 2015).

Further evidence has risen to support this interpretation in the JWST era, in the form of "maximally old" quiescent galaxies identified at z = 3–5 with formation redshifts z > 10 (Glazebrook et al. 2024), z ∼ 10 (de Graaff et al. 2024b), and z ∼ 8 (Carnall et al. 2023). The higher redshift that these objects are observed at allows for better resolution of the formation timescales, due to the younger age of the Universe. This puts the formation of at least some of these massive galaxies in the first 500–600 Myr after the big bang. These objects have high stellar masses of (0.4–1) × 10¹¹ M_⊙ and must have quenched shortly after the formation of the bulk of their stellar mass, but likely progenitors have yet to be conclusively identified at this epoch.

The above gap in the observed evolution of massive quiescent galaxies could be bridged by these objects, if the formation history inferred from the maximal stellar model is correct (Figure 5). A key question is whether the number densities of these objects match the number densities of massive quiescent galaxies at later times. Answering this question requires a larger area and/or a more well-defined selection function; for now we simply note that the number densities of these objects are comfortably below the number density of z = 3 compact quiescent cores (van Dokkum et al. 2015). The presence of Balmer breaks at z = 7–8, indicating very early and/or rapid formation, is probably the most straightforward prediction from finding "maximally old" z = 2–5 quiescent galaxies. While the number density of these objects may be low, they are indicative of a class of objects with early and rapid formation. In addition, the strong spectroscopic selection effect, which we do not attempt to account for in this work, suggests that the value of ∼10⁻⁵ Mpc⁻³ (Figure 6) is a lower limit of this class of objects. This is a key and perhaps convincing argument that, in at least some of these objects, the stellar mass could be high (i.e., the AGN emission contribution to the continuum in the minimal/medium stellar models may be overestimated).

Another distinguishing property of the class of low-z massive quiescent galaxies is the remarkable compactness and high stellar densities, with effective radii of 1 kpc or less at z ∼ 2 (van Dokkum et al. 2008). Baggen et al. (2023) showed that a fully stellar interpretation as presented in L23 (or, similarly, the maximal stellar model of this Letter) yields similar stellar densities to the cores observed at later times, albeit with significantly smaller sizes, by a factor of 10.

In addition, the enhanced SFR inferred from the maximal stellar model would corroborate the surprising overabundance of luminous galaxies at z ≳ 10 (e.g., Atek et al. 2023; Finkelstein et al. 2023; Robertson et al. 2023; Casey et al. 2024). Their dense stellar structures would presumably be formed from very dense gas associated with highly efficient star formation (Schmidt 1959; Kennicutt 1998), perhaps in a similar manner to the super star clusters with unusually high cloud surface densities in the local Universe (Smith et al. 2006; Turner et al. 2015).

Finally, the greater stellar mass case may be preferred on the grounds that it does not require a conveniently located AGN continuum to account for the observed spectral break. As the stellar mass decreases, so does the Balmer break strength in the galaxy spectrum, meaning that the AGN continuum must precisely match the shape and intensity needed to replicate the spectral break.

Therefore, we conclude that one possible interpretation is that the high-redshift Balmer break sample presented in this work consists of the progenitors of the first massive galaxies, observed directly after their rapid co-formation with their supermassive black holes. Certainly, this interpretation leaves difficult problems. We discuss the outstanding questions below, in connection with an alternative interpretation of these objects being low-mass galaxies hosting AGNs.

To start, it is not intuitive to explain the broad emission lines with a non-AGN origin, since stellar-driven scenarios are seemingly inconsistent with the lack of velocity offset between the narrow and broad emission line components and the symmetry of the lines. The broad emission lines also suggest an abundance of ionizing photons, which is likewise difficult to make sense of under the stellar interpretation. However, all the spectroscopically confirmed LRDs, which share similar but nonidentical SED shapes to our sample, are underluminous in X-ray (Kocevski et al. 2023; Maiolino et al. 2023; Wang et al. 2024a; Furtak et al. 2024; Greene et al. 2024; Matthee et al. 2024). More recently, most LRDs are found to be underdetected in X-ray in Chandra observations (Ananna et al. 2024; Yue et al. 2024). It perhaps is then fair to speculate alternative, non-AGN-driven physical causes.

Second, the implied formation is uncomfortably early and efficient, compared to the conventional assumption on the baryon-to-stellar conversion efficiency (Figure 7; Boylan-Kolchin 2023). Paradoxically, the cores of local massive galaxies appear to have bottom-heavy stellar initial mass functions (e.g., Conroy & van Dokkum 2012), with some evidence that this persists or even strengthens at z ∼ 2 (van Dokkum et al. 2024; though see Mercier et al. 2023, for an alternate take). This would increase the inferred stellar masses (without changing the observed SED in any way; Wang et al. 2024c) by a factor of a few, further increasing the tension with the cosmic baryon fraction.

Third, the objects of this work are remarkably compact (r_e ≲ 0.1 kpc; Baggen et al. 2023). This is even more striking when comparing to RUBIES-EGS-QG-1, which has r_e ∼ 0.6 kpc (de Graaff et al. 2024b), and the z ∼ 2–3 compact quiescent galaxies, which are likewise much larger (by a factor of ∼10; van Dokkum et al. 2015). At face value, the stellar bodies of these objects would have to rapidly expand over the subsequent few hundred million years. Additionally, it is difficult to imagine a massive galaxy being less than 100 pc in size. This would indicate an increased importance for dynamical evolution effects normally reserved for dense globular clusters (Spitzer 1987; Vesperini & Heggie 1997), e.g., segregation by mass, where the massive stars and binaries tend to sink toward the cores and the low-mass stars move outward into the halo. Interestingly, the mass segregation, if it happened, would change the initial mass distribution and thus relieve the tension in inferring uncomfortably large stellar mass in the current models, without evoking a change in the form of the initial mass function.

While thus far we have opted to center our discussion on the intriguing implications of the maximal stellar model—a choice motivated by the newly discovered high-z Balmer breaks in this Letter—we additionally acknowledge another possibility here. An alternative interpretation, suggested by the minimal/medium stellar mass models, posits that these objects are low-mass galaxies hosting AGNs. Possible evolutionary tracks for this scenario have been discussed extensively in the literature (e.g., Maiolino et al. 2023; Greene et al. 2024). It is worthy emphasizing, however, that these objects need not necessarily all belong to the same category—the continuum composition may vary on an individual basis. Future deeper and redder observations are critical to distinguishing between galaxy and AGN contributions to the continuum in this population. We elaborate on the lingering questions affecting the interpretation of our sample, along with possible ways forward, in the subsection below.

One remaining puzzle is the contrast between the broad Balmer lines and the narrow forbidden emission lines. While the line widths have been shown to correlate well with mass in statistical samples (Wuyts et al. 2016), narrow line widths in Hα and CO have also been observed in some massive galaxies at z ∼ 2.3 when there is considerably less ambiguity about the stellar masses (van Dokkum et al. 2015; Mowla et al. 2019). Possibly, these narrow line widths can be explained by some combination of face-on disks and peculiar gas geometry (i.e., the emitting gas is not associated with the deep stellar potential well). The discrepancy between the significantly lower dynamical mass and the inferred stellar mass from the medium/maximal stellar SED models of all three objects, coupled with the absence of available evidence indicating the existence of disks within them, makes such an argument less satisfying. However, it remains possible that the color-based selection function selects for a particular face-on orientation angle.

As for the M_⋆–M_BH relation, we cautioned about the selection bias (Lauer et al. 2007) in Section 6.6. Recently it has been argued that overmassive black holes are not inconsistent with the local relation after taking into account the selection bias (Li et al. 2024). A well-defined selection function was one of the key pieces in designing the observation strategy for the RUBIES program, and we will perform a more complete population-level analysis in a future Letter.

A key missing piece of evidence is deep MIRI imaging—the stellar-dominated and AGN-dominated models can diverge at these wavelengths by a factor of ≳2. The MIRI filters probing the 1.6 μm stellar bump (Laurent et al. 2000; Sawicki 2002) and at the reddest end generally provide the most discriminating power. At longer wavelengths, stellar light is expected to decline at longer wavelengths, whereas AGN emission is expected to show a flat spectrum (if deficient in hot dust; Williams et al. 2023; Wang et al. 2024a; Pérez-González et al. 2024) or a rising spectrum (from a dusty torus). For instance, MIRI/F1280W happens to capture the peak of the 1.6 μm stellar bump in the spectra of RUBIES-EGS-49140 and 55604, where the fluxes from the maximal and minimal stellar models differ by ∼3. X-ray detections or firmer upper limits would also help to establish the AGN nature of these objects.

Finally, a word of caution: despite our exploration of several models, ranging from the minimal to maximal extremes of stellar mass, none of these models fully align with our current understanding of galaxy formation and evolution. This challenge arises from a combination of the difficulty in separating AGN and host galaxy light and the ongoing debate surrounding evolutionary scenarios sparked by JWST's discoveries of (potentially) overmassive black holes and massive quiescent systems. In light of all the complexities, the statements in this Letter should all be considered contingent. There certainly remains ample room for reinterpretation in the future.

JWST is revolutionizing our knowledge of the formation of early galaxies and their black holes. In this Letter, we report a remarkable discovery of prominent Balmer breaks as early as z_spec = 8.35 and, intriguingly, their coexistence with broad emission lines. The high-redshift Balmer breaks reveal unambiguously the presence of evolved stellar populations with extended formation histories within the first 600–800 Myr after the big bang. However, all of the examined explanations on the potential origins and evolutionary tracks for these objects leave key lingering questions. Deeper spectroscopic data revealing stellar absorption features and JWST/MIRI data sampling the red continuum would further elucidate the nature of these intriguing objects. We conclude by emphasizing that observed Balmer breaks establishing the existence of evolved stellar populations with extended formation histories, presented herein, mark an important development in understanding the origins and evolution of galaxies and their central supermassive black holes.