The identification of red apparently massive galaxies at z gt 7 in early James Webb Space Telescope JWST photometry suggests a strongly accelerated time line compared to standard models of galaxy growth. A major uncertainty in the interpretation is whether the red colors are caused by evolved stellar populations dust or other effects such as emission lines or active galactic nuclei AGNs. Here we show that three of the massive galaxy candidates at z  6.78.4 have prominent Balmer breaks in JWSTNIRSpec spectroscopy from the RUBIES program. The Balmer breaks demonstrate unambiguously that stellar emission dominates at  SUBrestSUB  0.4 m and require formation histories extending hundreds of millions of years into the past in galaxies only 600800 Myr after the big bang. Two of the three galaxies also show broad Balmer lines with H FWHM gt 2500 km sSUP1SUP suggesting that dustreddened AGNs contribute to or even dominate the spectral energy distributions of these galaxies at  SUBrestSUB  0.6 m. All three galaxies have relatively narrow O III lines seemingly ruling out a highmass interpretation if the lines arise in dynamically relaxed inclined disks. Yet the inferred masses also remain highly uncertain. We model the highquality spectra using Prospector to decompose the continuum into stellar and AGN components and explore limiting cases in stellarAGN contribution. This produces a wide range of possible stellar masses spanning M SUBSUB  10SUP9SUP10SUP11SUP M SUBSUB. Nevertheless all fits suggest a very early and rapid formation most of which follow with a truncation in star formation. Potential origins and evolutionary tracks for these objects are discussed from the cores of massive galaxies to lowmass galaxies with overmassive black holes. Intriguingly we find all of these explanations to be incomplete deeper and redder data are needed to understand the physics of these systems.
['Active galactic nuclei', 'AGN host galaxies', 'Galaxy evolution', 'Galaxy formation', 'High-redshift galaxies', 'Spectral energy distribution', '16', '2017', '594', '595', '734', '2129', 'Astrophysics - Astrophysics of Galaxies']
RUBIES Evolved Stellar Populations with Extended Formation Histories at z  78 in Candidate Massive Galaxies Identified with JWSTNIRSpec
{'No Header': '9', 'RUBIES: Evolved Stellar Populations with Extended Formation Histories at z ∼ 7 -8 in Candidate Massive Galaxies Identified with JWST/NIRSpec': "Bingjie Wang ( 王 冰 洁 ), 1, 2, 3 Joel Leja, 1, 2, 3 Anna de Graaff, 4 Gabriel B. Brammer, 5, 6 Andrea Weibel, 7 Pieter van Dokkum, 8 Josephine F.W. Baggen, 8 Katherine A. Suess, 9, ∗ Jenny E. Greene, 10 Rachel Bezanson, 11 Nikko J. Cleri, 12, 13 Michaela Hirschmann, 14 Ivo Labb'e, 15 Jorryt Matthee, 16 Ian McConachie, 17 Rohan P. Naidu, 18, ∗ Erica Nelson, 19 Pascal A. Oesch, 7, 5 David J. Setton, 10, † and Christina C. Williams 20, 21  \n- 1 Department of Astronomy & Astrophysics, The Pennsylvania State University, University Park, PA 16802, USA\n- 2 Institute for Computational & Data Sciences, The Pennsylvania State University, University Park, PA 16802, USA\n- 3 Institute for Gravitation and the Cosmos, The Pennsylvania State University, University Park, PA 16802, USA  \n4 Max-Planck-Institut fur Astronomie, Konigstuhl 17, D-69117, Heidelberg, Germany  \n5 Cosmic Dawn Center (DAWN), Copenhagen, Denmark  \n6 Niels Bohr Institute, University of Copenhagen, Jagtvej 128, Copenhagen, Denmark  \n7 Department of Astronomy, University of Geneva, Chemin Pegasi 51, 1290 Versoix, Switzerland  \n8 Department of Astronomy, Yale University, New Haven, CT 06511, USA  \nKavli Institute for Particle Astrophysics and Cosmology and Department of Physics, Stanford University, Stanford, CA 94305, USA  \n10 Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA  \n- 11 Department of Physics & Astronomy and PITT PACC, University of Pittsburgh, Pittsburgh, PA 15260, USA  \n12 Department of Physics and Astronomy, Texas A&M University, College Station, TX, 77843-4242 USA  \n- 13 George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, Texas A&M University, College Station, TX, 77843-4242 USA  \n14 Institute of Physics, Lab for Galaxy Evolution, EPFL, Observatoire de Sauverny, Chemin Pegasi 51, 1290 Versoix, Switzerland 15 Centre for Astrophysics and Supercomputing, Swinburne University of Technology, Melbourne, VIC 3122, Australia 16 Institute of Science and Technology Austria (ISTA), Am Campus 1, 3400 Klosterneuburg, Austria 17 Department of Physics and Astronomy, University of California, Riverside, Riverside, CA 92521, USA 18 MIT Kavli Institute for Astrophysics and Space Research, Cambridge, MA 02139, USA 19 Department for Astrophysical and Planetary Science, University of Colorado, Boulder, CO 80309, USA 20 NSF's National Optical-Infrared Astronomy Research Laboratory, Tucson, AZ 85719, USA 21 Steward Observatory, University of Arizona, Tucson, AZ 85721, USA", 'ABSTRACT': 'The identification of red, apparently massive galaxies at z > 7 in early James Webb Space Telescope (JWST) photometry suggests a strongly accelerated timeline compared to standard models of galaxy growth. A major uncertainty in the interpretation is whether the red colors are caused by evolved stellar populations, dust, or other effects such as emission lines or AGN. Here we show that three of the massive galaxy candidates at z = 6 . 7 -8 . 4 have prominent Balmer breaks in JWST/NIRSpec spectroscopy from the RUBIES program. The Balmer breaks demonstrate unambiguously that stellar emission dominates at λ rest = 0 . 4 µ m, and require formation histories extending hundreds of Myr into the past in galaxies only 600-800 Myr after the Big Bang. Two of the three galaxies also show broad Balmer lines, with H β FWHM > 2500 kms -1 , suggesting that dust-reddened AGN contribute to, or even dominate, the SEDs of these galaxies at λ rest ≳ 0 . 6 µ m. All three galaxies have relatively narrow [O iii ] lines, seemingly ruling out a high-mass interpretation if the lines arise in dynamically-relaxed, inclined disks. Yet, the inferred masses also remain highly uncertain. We model the high-quality spectra using Prospector to decompose the continuum into stellar and AGN components, and explore limiting cases in stellar/AGN contribution. This produces a wide range of possible stellar masses, spanning M ⋆ ∼ 10 9 -10 11 M ⊙ . Nevertheless, all fits suggest a very early and rapid formation, most of which follow with a truncation in star formation. Potential origins and evolutionary tracks for these  \nobjects are discussed, from the cores of massive galaxies to low-mass galaxies with over-massive black holes. Intriguingly, we find all of these explanations to be incomplete; deeper and redder data are needed to understand the physics of these systems.  \nKeywords: Active galactic nuclei (16) - AGN host galaxies (2017) - Galaxy evolution (594) - Galaxy formation (595) - High-redshift galaxies (734) - Spectral energy distribution (2129)', '1. INTRODUCTION': "In the cores of the most massive galaxies in the local universe, stars have inferred stellar age of ∼ 13 Gyr and high α -element abundance, suggesting their stellar components are formed at z ≳ 5 in a spectacular and short burst of star formation (e.g., Thomas et al. 2005). This hypothesis was bolstered by the discovery of their putative z ∼ 2 progenitors, compact galaxies with high stellar masses ∼ 10 11 M ⊙ and small effective radii of ∼ 1 kpc, which importantly have stellar densities similar to the cores of z ∼ 0 ellipticals (Bezanson et al. 2009). When exactly their stellar bodies formed, however, has not yet been clearly established. While some likely compact star-forming progenitors have been identified at z = 2 -3 (Nelson et al. 2014; Barro et al. 2014), simulations and number density arguments suggest that at least some of these massive cores must have formed earlier (Wellons et al. 2015). This has since been buttressed by the James Webb Space Telescope (JWST) discovering and characterizing massive quiescent galaxies at z = 2 -5, with stellar bodies of 10 10 . 5 -11 M ⊙ and inferred formation redshifts of z ≳ 7 (Carnall et al. 2023; Glazebrook et al. 2024; de Graaff et al. 2024a; Park et al. 2024). Taken at face value, this implies a very rapid formation and quenching of very massive galaxies in the first billion years of the universe-a spectacular event that should produce a huge amount of observable light. Yet, while candidates have been found (e.g., Hashimoto et al. 2018; Williams et al. 2023), the progenitors of these quenched galaxies existing at t univ ∼ 1 Gyr have yet to be conclusively identified.  \nContemporaneously, one of the most surprising early discoveries made with the JWST is the identification of a population of seemingly massive ( M ⋆ ≳ 10 10 M ⊙ ), compact (effective radii ≲ 1 kpc), and rest-optical red galaxies at redshift z > 6 via a double-break color selection (Labb'e et al. 2023a; hereafter L23). This selection targets distinct spectral energy distributions (SEDs) which include a drop-out at ∼ 1 µ m and a very red color at ∼ 3 µ m in observed frame. This is a highly efficient  \nselection for 7 ≲ z ≲ 9 objects, but the red rest-frame optical color may be driven by very different underlying physics: it could come from a Balmer break, from strong emission lines, from a very red continuum, or from some combination of these features.  \nIn their main analysis, L23 interpreted the red color as a combination of a Balmer break and a very red rest-optical stellar continuum. This interpretation yields very high stellar mass-to-light ratios (M/L) and implies extreme stellar masses, up to M ⋆ ∼ 10 11 M ⊙ . While massive galaxies must form early and quickly, this would imply that these objects both emerged earlier and hosted more mass than expected. Indeed, it soon became clear that such early massive galaxies are difficult to make in the standard model of cosmology, as the amount of baryons inferred in stars is comparable to the cosmic baryon abundance in these early halos (Boylan-Kolchin 2023). Star formation feedback processes typically limit the fraction of baryons locked up in stars to far below the cosmic baryon fraction.  \nThese uncomfortably high stellar masses prompted a wave of alternative explanations for the observed fluxes and colors, including extreme emission line galaxies (Endsley et al. 2023), a top-heavy stellar initial mass function (Boylan-Kolchin 2023; Steinhardt et al. 2023), or obscured active galactic nuclei (AGN) (Kocevski et al. 2023; Barro et al. 2024). Obscured AGNs significantly lower the inferred stellar masses by contributing to the red optical continuum flux. They can also lower the photometric redshift, and so decrease the inferred number densities of massive galaxies at high redshifts. In addition to their red color, the L23 sample exhibits compact morphology, supporting the idea that AGN could contribute to the fluxes.  \nIndeed, one of the bright red objects in L23 has been confirmed to host a broad-line AGN at a lower redshift of z = 5 . 62 (Kocevski et al. 2023), and the more broadly defined red compact sources (often dubbed as little red dots, or LRDs; e.g., Labb'e et al. 2023b; Matthee et al. 2024) have shown a prevalence of broad-line AGNs (Greene et al. 2024). However, as found in Baggen et al. (2023), most of the L23 sample are marginally resolved at rest-UV wavelengths-so their rest-UV flux at least is not dominated by a point source. Furthermore, the  \ninferred stellar densities, assuming that the L23 objects are galaxies, are consistent with the central regions of today's elliptical galaxies.  \nTherefore, these objects may still be massive galaxies with evolved populations at high redshifts-in fact, it would be surprising if star formation slowed sufficiently to show prominent Balmer breaks at high redshift without influence from an AGN. The two massive quiescent systems currently known at z ∼ 5 both show indications of AGN activity (Carnall et al. 2023; de Graaff et al. 2024a). It has long been suspected that the stellar cores of galaxies must have co-formed with their supermassive black holes; such a formation scenario is the simplest interpretation of the observed tight correlation between the mass of the bulge and the mass of the central supermassive black hole (Ferrarese & Merritt 2000). The key questions regarding these red objects are therefore threefold: (i) are the photometric redshifts accurate; or more generally, what is the source of the red colors-emission lines, continua, or spectral breaks? (ii) do these objects host evolved stellar populations? and (iii) how much of the continuum is powered by stellar versus AGN emission? These core questions are unanswerable without rest-frame optical spectra.  \nHere we conduct a follow-up study on double-break candidates selected from the RUBIES JWST/NIRSpec spectroscopic program (JWST-GO-4233; PIs de Graaff & Brammer). These targets partially overlap with the L23 sample, and were observed with high priority for their red colors (F150W-F444W > 2), and bright apparent magnitudes (F444W < 26 mag). In this paper, we present three objects with detected Balmer breaks, existing as early as z spec = 8 . 35, suggesting a formation history extending hundreds of Myr into the past in galaxies only 600-800 Myr after the Big Bang. Unambiguous broad emission lines are also observed in 2/3 objects, motivating a deeper look at the source of the red continua.  \nThe structure of this paper is as follows. Section 2 provides an overview of the data. Section 3 presents the key observational features of these objects, including observed Balmer breaks and broad lines. Section 4 focuses on the analysis of the emission lines, while Section 5 details the AGN and host galaxy composite SED modeling. Section 6 presents evidence for/against the different proposed physical models, including kinematics, formation histories, and population-level characteristics. We conclude in Section 7 with discussion of possible interpretations and evolutionary scenarios, and outstanding questions.  \nWhere applicable, we adopt the best-fit cosmological parameters from the WMAP 9 yr results: H 0 = 69 . 32  \nkms -1 Mpc -1 , Ω M = 0 . 2865, and Ω Λ = 0 . 7135 (Hinshaw et al. 2013). Unless otherwise mentioned, we assume the Chabrier (2003) initial mass function, and report the median of the posterior, with associated 1 σ error bars being the 16th and 84th percentiles.", '2. DATA': 'The spectroscopic survey RUBIES uses JWST/NIRSpec to observe approximately 4000-5000 sources selected from public NIRCam imaging in the EGS and UDS fields. It utilizes the NIRSpec microshutter array (MSA; Ferruit et al. 2022) with both the low-resolution Prism/CLEAR and medium-resolution G395M/F290LP disperser/filter combinations. The sample in this paper was targeted in 6 masks observed in March 2024. Figure 1 presents an overview of the sample in color space.  \nReduction of the imaging data is presented in Valentino et al. (2023) while the spectroscopic reductions are described in Heintz et al. (2024); Wang et al. (2024a). Full details of the RUBIES observing program and data reduction will be detailed in A. de Graaff et al. in prep. The subsequent sections provide a brief summary.', '2.1. Imaging': "The RUBIES targets in the EGS were selected from public JWST/NIRCam data from the Cosmic Evolution Early Release Science Survey (CEERS, JWST-GO1345; PI Finkelstein; Bagley et al. 2023; Finkelstein et al. 2023), which provides imaging in 7 bands (F115W, F150W, F200W, F277W, F356W, F410M, and F444W). Additionally, we use archival imaging in 7 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).  \nThe image mosaics, with a pixel scale of 0 . 04 '' pix -1 , are publicly available in the DAWN JWST Archive (DJA; version 7.2) and were reduced using grizli (Brammer 2023a). Fluxes are measured from PSFmatched images in circular apertures with a radius of 0 . 16 '' , and then Kron-corrected to the total flux, as described in Weibel et al. (2024).", '2.2. Spectroscopy and Sample Selection': "Each target was observed for 48 min in the Prism/CLEAR mode and the G395M/F290LP mode, using a standard 3-shutter slitlet and 3-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  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 1. Sample of this paper. (Left) This work focuses on continuum-detected sources with Balmer breaks, shown as solid diamonds on the color plane. The L23 sample is denoted as plus signs, and all sources in CEERS brighter than 27 AB mag are grey hexagons included for reference. (Right) We show the 2 × 2 '' color images of the three Balmer-break-detected sources, colors from JWST/NIRCam F115W, F277W, and F444W. They are remarkably compact at red wavelengths, with some evidence for spatial structure at blue wavelengths.  \n<!-- image -->  \n2023b), 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 re-normalize 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 ( § 4), and hence rescale the full G395M spectrum by the ratio of the two to match the flux calibration of the Prism spectrum.  \nThis paper focuses on the 3 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.", '3. THE COEXISTENCE OF BALMER BREAKS AND BROAD EMISSION LINES AT z = 6 . 6 -8 . 4': 'The most striking discovery is the seeming appearance of Balmer breaks observed between z = 6 . 7 -8 . 4. These  \nare produced in the atmospheres of older stars and typically only appear in evolved stellar populations after a significant reduction in star formation rate (SFR) lasting at least ∼ 100 Myr (Bruzual A. 1983; Hamilton 1985; Worthey et al. 1994; Balogh et al. 1999). It is surprising to see SEDs dominated by evolved stellar populations at these times-the age of the universe for the highest redshift object ( z = 8 . 35) is 610 Myr, giving very little time to form the stellar populations. This suggests an extremely rapid and early formation for these objects. The highest redshift candidate at z = 8 . 35 represents the highest-redshift Balmer break identified to date, with the next-highest redshift being a z = 7 . 3 object with a low mass of 10 8 . 5 M ⊙ (Looser et al. 2024).  \nTo verify that the observed spectral breaks are indeed Balmer breaks, as opposed to observational artifacts or a mis-identification of spectral breaks driven by other physics, we take the average of the three observed spectra after scaling each spectrum by the median of the continuum flux close to the break (between 4150 and 4250 ˚ A in rest frame). The same average is then calculated for the best-fit galaxy model spectra ( § 5). We show the normalized individual Prism spectra, over-plotted with the averaged model spectra in Figure 2. It is evident that all the breaks are located at the expected wavelength, and have similar shapes, buttressing the Balmer break interpretation. Critically, the stacked galaxy model spectra agree well with the data, suggesting that the model fits  \nFigure 2. Characteristics of the sample of this paper. (a) Balmer breaks detected at z spec = 6 . 7 -8 . 4. The scaled individual spectra which exhibit potential Balmer breaks are plotted in light colors, and the averaged spectrum is shown in black. Overplotted in red is the averaged best-fit model spectrum ( § 5). All breaks are located at the expected wavelength, and show similar shapes, supporting the Balmer break interpretation. The break strength for the average observed spectrum, quantified by the flux ratio in the wavelength windows illustrated in blue and red shades, is 2.0. (b) Emission line fits to the H β and [O iii ] complex ( § 4). Unambiguous broadened H β lines are found in the G395M spectrum for two sources (RUBIES-EGS-49140 and 55604), indicative of an AGN. A broad component of H β is only marginally detected in both the G395M and Prism fits for RUBIES-EGS-966323.  \n<!-- image -->  \nare also properly interpreting the light as evolved stellar populations. Section 5 elaborates on the SED modeling. A commonly used measure of the Balmer break strength is D 4000 , as originally defined in Bruzual A. (1983); Balogh et al. (1999). The wavelength windows in these definitions locate redward of the Balmer limit (3645 ˚ A), as they are intended for measuring the 4000 ˚ A break. However, we are more interested in the Balmer break. We thus define a spectral break strength, taking the median of the fluxes in two windows at [3620,3720] ˚ A and [4000,4100] ˚ A instead. These windows avoid contamination from strong nebular line emission and are close enough to be minimally impacted by dust attenuation, or the overall slope of the spectrum. Our spectral break strength differs primarily from D 4000 for galaxies with ages ≲ 1 . 5 Gyr, where the Balmer break is prominent but the 4000 ˚ A break is less visible. It is therefore a better age indicator for highredshift galaxies. A comparison between our new index with the Balogh et al. (1999) definition is supplemented in Appendix A. Meanwhile, the break strengths for RUBIES-EGS-49140, 55604, and 966323 are 2 . 44 ± 0 . 10, 2 . 18 ± 0 . 11, and 1 . 96 ± 0 . 14, respectively.  \nIn addition to hosting Balmer breaks, we note that the broad emission lines are evident in the spectra of RUBIES-EGS-49140 and 55604, as shown in Figure 2, while the existence of broad H β is more ambiguous for the final object RUBIES-EGS-966323. These observations motivate detailed emission line fitting to decompose them into broad and (potentially) narrow components, and also present an opportunity to estimate black hole masses. We do so in Section 4.', '4.1. Emission Line Decomposition': '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.  \nOur fitting methodology broadly follows that described in Wang et al. (2024a): given the compact mor-  \nphology of the sources (Baggen et al. 2023), we use a point source line-spread function (LSF) from de Graaff et al. (2024b), and the Markov Chain Monte Carlo (MCMC) 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.  \nWe select the region around the emission line complex ( ± 0 . 2 µ m), and begin by masking the H β line (a region of 7000 km s -1 ) 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] kms -1 using a uniform prior. The continuum is fit with a 1st order polynomial. We then explore whether the [O iii ] lines show evidence for a second broader component by fitting a 2-component Gaussian model, where the width of the second Gaussian is constrained to be larger than that of 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 σ broad , OIII = 248 +61 -41 kms -1 , and no evidence for a second component in the other two sources.  \nNext, we include the H β line in the fitting. We first fit a single Gaussian component for the H β line. Given the limited 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] kms -1 . For RUBIES-EGS-55604 we find a blue-shifted absorption feature in the H β line, which is modeled as an additional Gaussian component with a velocity offset ∆ v ∈ [0 , 1000] kms -1 and dispersion σ outflow ∈ [0 , 1000] km s -1 . Such absorption feature in Balmer and/or He lines has been seen in recent studies of AGNs as well (Kocevski et al. 2024; Matthee et al. 2024; Wang et al. 2024a). The multi-component fits and the residuals are shown in Figure 2.  \nAgain, 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 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 -18 erg s -1 cm -2 ). For this source we also perform an independent fit of the Prism spectrum, finding similarly weak evidence for a two-component model and  \nposteriors that are consistent with the fit to the G395M spectrum.  \nFinally, we fit the Balmer H α emission line, which is available in the G395M spectrum for RUBIES-EGS49140 and only in the Prism spectrum for RUBIESEGS-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 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 kms -1 . For RUBIES-EGS-49140 we also find a blue-shifted 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.', '4.2. Single-epoch Black Hole Mass Estimates': '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.  \nFollowing Assef et al. (2011), we estimate the black hole mass based on H β and luminosity at rest 5100 ˚ A as  \nlog( M BH /M ⊙ ) = 0 . 895+ 0 . 529 log ( L 5100 10 44 erg s -1 ) +2 . 0 log ( FWHM H β km s -1 ) . (1)  \nH β is used since H α is not available for all sources and suffers from uncertainty due the blending of the [N ii ] doublet with the broad component of H α . L 5100 is calculated from the intrinsic AGN spectrum inferred from the SED models.  \nWe 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  \nthese methods at higher redshifts and in different physical conditions than where they are calibrated (e.g., Yue et al. 2024).', '5. SPECTRAL ENERGY DISTRIBUTION MODELING': 'The clear detection of the Balmer breaks means that the rest-frame 3500 ˚ A wavelength range is dominated by starlight. However, Balmer breaks are indicative only that evolved stellar populations exist and are relatively bright-they can appear even with ongoing star formation. Therefore, unambiguous Balmer breaks leave some room for interpretation when inferring formation histories. The mixture of stellar (as suggested by the Balmer break) and AGN continuum emission (as suggested by the broad emission lines) motivate the need for detailed spectrophotometric modeling which includes contributions from both components.  \nGiven that decomposing galaxy/AGN light is a known challenge, we consider three models to capture the systematic uncertainties in the inferred properties and formation histories. Below we describe the free and fixed parameters of the SED model, and then introduce three different priors which lead to three different interpretations of the observed light.', '5.1. Core Model Components': "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'anchez-Bl'azquez et al. 2006) in FSPS (Conroy & Gunn 2010). Star formation history (SFH) follows the non-parametric Prospector -α description, characterizing mass formation in 7 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 non-zero fraction of which suggests the presence of blue stars possibly existing outside the dust or having created holes within  \nit. Dust emission is incorporated based on the model by Draine & Li (2007).  \nModel spectra are convolved with the NIRSpec/Prism instrumental resolution curve, tailored for a point source morphology using msafit (de Graaff et al. 2024b). We account for wavelength-dependent slit losses by scaling the normalization of the spectrum to match the photometry through a 7th-order polynomial calibration vector applied during the fitting process.  \nEmission 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 an additional 50% inflation of the random noise produces a good fit. The posteriors are sampled via the dynamic nested sampler dynesty (Speagle 2020).", '5.2. Continuum Contribution from A Dusty Active Galactic Nucleus': 'The composite galaxy and AGN model is presented in Wang et al. (2024a), in which we approximate the direct ultraviolet (UV)-optical emission from an AGN accretion disk as piece-wise 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, UKIDSS, 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 5 free parameters (the AGN-to-galaxy flux ratio, and four dust attenuation parameters, two of which are shared with the galaxy light).', '5.3. Modeling A Range of Possible Stellar and AGN Contribution to the Rest-optical Continua': 'The key uncertainty in the interpretation of these objects is the source of the red continuum: is it powered by AGN, 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  \nFigure 3-a. Spectrophotometric modeling for RUBIES-EGS-49140, with models including maximal / medium / minimal stellar contribution shown. The other two objects with detected Balmer breaks are shown in subsequent figures. (Left panels) The photometric and spectroscopic data are shown in gray. The best-fit model spectrum, which includes the marginalized emission lines as annotated, is plotted in light brown. The emission at ∼ 3869 ˚ A is likely [Ne iii ] 3869, although He I 3889 is also possible. The galaxy and the AGN continuum components are over-plotted in blue and red, respectively. The spectral break strength, SBS, predicted by the galaxy model spectrum is indicated to the upper left corner. The spectral regions that are masked due to low S/N or detector gap are shaded in gray. (Right panels) The inferred SFRs are plotted as a function of the age of the universe. The gray and light gray shading indicates 1 σ and 2 σ uncertainties, respectively. The post-starburst feature is primarily driven by the Balmer break.  \n<!-- image -->  \nless 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), suggest the continuum blueward of rest 4000 ˚ A is dominated by stars.  \nThe central question is therefore the origin of the continuum redward of the Balmer break. We consider three models which bracket the possible range of inferred stellar and AGN contribution to the observed rest-optical continua.', '5.3.1. Maximal Stellar Contribution': "We begin by considering a model which maximizes the stellar contributions. This is achieved by placing a lognormal prior on the pre-dust-attenuated AGN-to-galaxy flux ratio at rest 5500 ˚ A, f AGN , 5500A , with mean µ = -3,  \nand standard deviation σ = 1. A log-uniform (i.e., flat in log-space) prior is used on the galaxy mass.  \nThis model down-weights the AGN contribution at wavelengths where it is naturally brightest, and effectively leads to stellar masses similar to those from a 'galaxy-only' fit. While this serves as a useful benchmark, it also implies the broad line EWs, assuming they are driven by AGN, are extremely, perhaps unphysically, high.", '5.3.2. A Mixture of AGN and Stellar Contributions': "We also define a middle ground model, where loguniform priors are assumed on the total mass formed and f AGN , 5500A . While it is impossible to write down an agnostic prior, the intent here is to let the data inform the inference process to the maximal extent. This model typically falls between the 'maximum' and 'minimum'  \nFigure 3-b. Spectrophotometric modeling for RUBIES-EGS-55604 (identified as 38094 in L23), with format as described in Figure 3-a.  \n<!-- image -->  \nstellar contributions. However, this does not mean that this model is to be taken as the fiducial model or in some way better motivated than the other choices. The available spectroscopic data are not constraining enough to break the degeneracies among the stellar populations, black hole properties, and dust attenuation.", '5.3.3. Minimal Stellar Contribution': 'For the minimal stellar model, we impose a prior on the galaxy mass, with the probability P ( M ) ∝ M α . The slope, α = -1 . 7, is taken to be the low-mass slope of the theoretical stellar mass function at z = 7 (Tacchella et al. 2018). This prior is not intended to exactly replicate the mass function and serves its purpose of producing a low-mass solution even if the true mass function has a different low-mass slope.  \nA more intuitive prior down-weighting the galaxy contribution would be on the fractional light contribution instead of directly on the stellar mass, as is done with the maximal stellar model. However, the complex translation from light to mass means that, with a prior on the light only, the dust and stellar M/L can adjust to keep  \nthe stellar masses similar; a direct prior on the stellar mass avoids this degeneracy.', '6. RESULTS': 'Here we describe the properties resulting from the individual fits to the three objects with Balmer breaks detected at z spec = 6 . 7 -8 . 4. As alluded to in Section 5, even in the presence of an unambiguous Balmer break, there can be a wide range of inferred star formation histories. The Balmer break strengths, as measured by the spectral break strength defined in Section 3, vary from 2.2 to 1.7 (Figure 3-a), 2.1 to 1.2 (Figure 3-b), and 1.7 to 1.1 (Figure 3-c), by altering the AGN contribution from ∼ 0 to ∼ 100% for RUBIES-EGS-49140, 55604, and 966323, respectively.  \nThe range of spectral break strengths also suggests that the adopted sets of SED models bracket the uncertainty in continuum decomposition. As summarized in Table 1, the medium stellar model typically leads to M ⋆ ∼ 10 9 . 8 M ⊙ . If the red optical light is instead driven by an older dusty stellar population, M/L can rise dra-  \nFigure 3-c. Spectrophotometric modeling for RUBIES-EGS-966323 (identified as 14924 in L23), with format as described in Figure 3-a.  \n<!-- image -->  \nlly, producing masses up to M ⋆ ∼ 10 11 M ⊙ . In the other direction, a prior that minimizes the inferred stellar masses can still reproduce the light while lowering the inferred masses by up to 1 dex.  \nIn this section we lay out the physical properties and the implications of the formation histories inferred by the three physical models, including both observational and evolutionary considerations. The aim here is to provide an unbiased, comprehensive view into the very different possible physical interpretations of these objects-all of which result in statistically acceptable fits to the observed spectra.', '6.1. Apparent Wavelength-dependent Contribution of AGN and Stellar Light': '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 AGNdominated light can in principle be corroborated by the wavelength-dependent morphology-an unresolved mor-  \nphology is more indicative of AGN-dominated light. The presence of broad emission lines is both a suggestion of AGNactivity 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 ˚ A (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.  \nWe find unambiguous broad H β lines with FWHM > 2500 kms -1 in RUBIES-EGS-49140 and 55604, which would be unusual for stellar-driven outflows (see Veilleux et al. 2005 for an overview, or 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-  \nTable 1. Inferred AGN and Host Galaxy Properties, and Emission Line Kinematics  \n966323 ( § 4). Additionally, we compare the observed H α EW to previous samples of LRDs. RUBIES-EGS49140 and 55604 show large observed EWs of > 600 ˚ A. We re-estimate the EW using the AGN continuum from the SED fits. Assuming the broad H α is originated 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.  \nWe 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 which produces more H α flux than is observed. Conversely, the observed  \nlarge H α EWsuggests an abundance of ionizing photons, the evidence of which is not obvious from the spectra in the UV-optical.  \nFurthermore, we observe an emission line at ∼ 3869 ˚ A in all three objects. One possibility is He I 3889, given that 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 high ionization state, particularly considering the weak (non-detected) [O ii ] 3727 emission. A strong [Ne iii ]/[O ii ] ratio can also be indicative of AGN activity (Backhaus et al. 2024).  \nFinally, we cross-check the morphologies of RUBIESEGS-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 LW 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 is too faint to have a robust size estimate. Preliminary analysis on RUBIES-EGS-49140 reveals the same trend. A thorough study on the mor-  \nFigure 4. Rest-frame equivalent widths of H α . EWs of the broad component of H α measured for the sample of this paper at the total continuum and AGN continuum (assuming the minimal stellar model) are represented in black and red, respectively. For reference, the EWs of a grism-selected LRD sample are plotted in gray (Matthee et al. 2024). Assuming the broad lines trace the AGN, then the implied EW at the AGN continuum in the case of starlight dominating the rest optical would be orders of magnitudes higher than the observed EW. Conversely, an AGN-dominated interpretation puts the EWs of our sample in a similar range as the grism-selected LRD sample.  \n<!-- image -->  \nphologies will be presented in J. Baggen et al. in prep. 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.', '6.2. Inferred Formation Histories': 'The inferred SFHs based on the three SED models are shown in Figure 3-a-3-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  \nold galaxies observed at lower redshifts (Carnall et al. 2023; Glazebrook et al. 2024) or a more temporary miniquenching 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.  \nConversely, the minimal-stellar models decrease the strength of the stellar Balmer breaks by assuming a very red AGN continuum underneath, and in this way 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-EGS49140, 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.  \nOne potential concern may arise from the agemetallicity 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.  \nInterestingly, 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; Glazebrook et al. 2024; de Graaff et al. 2024a), as illustrated in Figure 5. Such a link may serve as an argument in favor of the maximalstellar model.  \nWe note, though, the above connection is still hard to establish given the low number densities of these massive objects (a lower limit of ∼ 10 -5 Mpc -3 for the three objects in this work, and e.g., ∼ 5 × 10 -6 Mpc -3 for RUBIES-EGS-QG-1 estimated in de Graaff et al. 2024a). 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  \nFigure 5. Formation history inferred for RUBIES-EGS-49140, compared to those inferred for maximally old quiescent galaxies at z ∼ 3 -5. From left to right, we show the SFHs of RUBIES-EGS-QG-1 (at z = 4 . 9 with 10 11 M ⊙ ; de Graaff et al. 2024a), GS-9209 (at z = 4 . 7 with 4 × 10 10 M ⊙ ; Carnall et al. 2023), and ZF-UDS-7329 (at z = 3 . 2 with 10 11 M ⊙ ; Glazebrook et al. 2024). For the the z = 4 . 9 galaxy, we include the fiducial low-metallicity solution in black, with 1 σ and 2 σ uncertainties shaded in grey, and also an alternative solution assuming a solar metallicity in blue. For the other two quiescent galaxies, we only show the 1 σ uncertainty. The top panels show the formation history from our medium-stellar model while the bottom panels show the formation history from our maximal-stellar model. A stellar interpretation of the light in RUBIES-EGS-49140 produces greater consistency with the massive quiescent galaxies, while the medium model fails to predict sufficient stellar mass to connect RUBIES-EGS-49140 to these massive old populations at lower redshift.  \n<!-- image -->  \nare indicative of a class of objects with early and rapid formation.', '6.3. Contribution to the Cosmic Stellar Mass Density': 'One of the key results from L23 is that the high stellar masses in these objects suggest 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 2 ; Finkelstein et al. 2023). As in L23, we only consider the Poisson uncertainty and cosmic variance, and neglect corrections for incompleteness, given that any correction for incom-  \neness would increase the inferred stellar mass densities.  \nWe note that there is some evidence indicating that EGS exhibits an over-density of red objects. So far, nothing similar to the sample of this paper, in terms of highz 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 (i.e., UDS, which has an area of 224 arcmin 2 ; 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.  \nFigure 6. Implied number density and mass density. (a) Number density of different samples. The red and orange curves are taken from van Dokkum et al. (2015). The blue curve is the total number density for UV-selected sample, obtained by integrating the Schechter fits down to 10 8 M ⊙ (Stefanon et al. 2021). Number density for the sample of this paper assuming the EGS area for the volume, in two redshift bins (5 . 5 < z ≤ 7 . 0, and 7 . 0 < z ≤ 8 . 5), are shown in black; uncertainties reflect Poisson statistics and cosmic variance (Gehrels 1986). This illustrates, preliminarily, that the number density of the Balmer break sample are comfortably below that of the z = 3 compact quiescent cores. (b) Cumulative cosmic stellar mass density, ρ ⋆ . The curves are derived from the same Schechter fits, with the extrapolated regions indicated in dash. ρ ⋆ from this work are reported for the three SED models. Uncertainties again reflect Poisson statistics and cosmic variance. To facilitate a direct comparison to the Stefanon et al. (2021) stellar mass functions, we adjust all our Chabrier stellar masses by +0 . 24 dex to Salpeter masses (Salpeter 1955). While the number density of the Balmer break sample is low, the mass density in these objects at least comprise ∼ 1% of the cosmic mass budget (albeit with large uncertainties; ∼ 0 . 1% within 1 σ ). With the medium-stellar model, their contribution is roughly equal to the mass density in all UV-selected objects combined. Note that we do not attempt to account for selection effects, so all the reported densities for our sample should be taken as lower limits (see § 6.3 for details).  \n<!-- image -->  \nFigure 7. Comparison to theoretical limits. (a) Limits on the abundance of galaxies as a function of redshift. Curves illustrate M ⋆ as functions of z at fixed cumulative halo abundance, assuming the physically maximal ϵ = 1 . 0 to the left, and the more likely case of ϵ = 0 . 32 to the right. These plots suggest that it would be rare to find galaxies as massive as the maximal stellar mass case in these redshifts. (b) Stellar mass density limits. Similar to Figure 6-b, but we now include the theoretical co-moving stellar mass density contained within galaxies more massive than M ⋆ in the two redshift bins for three values of ϵ . The maximal stellar mass case implies an an unrealistic limit that all available baryons in the haloes are converted into stars. For completeness, the naive estimates of the lower limits on the number density, based on the CEERS+UDS volume (instead of considering the CEERS field alone as in the fiducial case), are included here as unfilled gray markers.  \n<!-- image -->  \n/circledot  \nFigure 8. Stellar-black hole mass scaling relationship. The z ∼ 0 (solid line) and z ∼ 0 ellipticals (dashed line) relations correspond to the black hole mass and host galaxy total stellar mass scaling, and non-active ellipticals with dynamical black hole mass measurements, respectively (Reines & Volonteri 2015). The gray shading indicates the intrinsic scatter. Over-plotted are the black hole masses of our sample estimated from L 5100 and the broad H β lines, and their corresponding stellar mass inferred from three SED models. The semi-transparent points indicate RUBIES-EGS-966323 for which the presence of a broad line is ambiguous. The black hole masses are uncertain as the inferred AGN luminosity from the SED models vary.  \n<!-- image -->  \nFigure 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 paper, 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 UVselected mass functions, as a test of what a rest-UV selection function may have missed.  \nAssuming 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%  \nof 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 comprise 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 star formation histories, whereas UV-selected galaxies typically have rising formation histories by nature, this fraction is likely to be even higher at higher redshifts.  \nFor 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 8 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.', '6.4. Comparison to Theoretical Limits': '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 which are relevant for dark matter haloes 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).  \nFollowing 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 co-moving number density of haloes above a given halo mass threshold, M halo , is estimated as  \nn ( > M halo , z ) = ∫ ∞ M halo d M d n ( M,z ) d M , (2)  \nwhere d n ( M,z ) / d M is the number of dark matter haloes of mass M per unit mass per unit co-moving volume at redshift z (i.e., the halo mass function). Then the comoving mass density in haloes more massive than M halo can be computed straightforwardly as  \nρ m ( > M halo , z ) = ∫ ∞ M halo d MM d n ( M,z ) d M . (3)  \nFrom the above, we can obtain the corresponding statistics of galaxies by assuming  \nM ⋆ = ϵf b M halo , (4)  \nwhere ϵ is the efficiency of conversion of baryons into stars, and f b is the cosmic baryon fraction. ϵ = 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 , 0 . 32 (e.g., Behroozi et al. 2013; Moster et al. 2013).  \nAs 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 haloes with cumulative co-moving number densities ≲ 10 -5 . 4 Mpc -3 ). For ϵ = 0 . 32 (0 . 10), the implied number density is ≲ 10 -7 . 1 (10 -9 . 5 ) Mpc -3 , 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 unfilled gray markers 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).', '6.5. Dynamical Mass': '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 -1 (Table 1). As gas is viscous, the forbidden lines likely originate in disks, and-assuming those disks trace the gravitational potential of the galaxies-can be used to derive dynamical masses. As discussed in, e.g., Forster Schreiber et al. 2014; van Dokkum et al. 2015; 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 point-like morphology as the Balmer lines), the latter two dependencies are both unknown.  \nNevertheless, it is instructive to estimate the dynamical mass. Following van Dokkum et al. (2015), we have  \nM dyn = 2 V 2 rot r e G , (5)  \nwhere r e is the spatial scale, taken to be ∼ 100 pc (Baggen et al. 2023), and V rot is the rotation velocity derived via  \nV rot = σ gas α sin -1 ( i ) . (6)  \nAssuming α = 0 . 8 (Rix et al. 1997; Weiner et al. 2006), and an inclination angle, i , of 45 degrees, 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).  \nHigher 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. 2024b). Therefore, if the source were to have r e ∼ 1 kpc, then σ gas would decrease to < 30 kms -1 ; 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-thebarrel (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.', '6.6. Galaxy-Black Hole Scaling Relationship': '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) agree well with the Prospector M/L-color relationship (Li & Leja 2022).  \nBased 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), quasars (Stone et al. 2024), as well as that a z ∼ 5 AGN sample which is located ∼ 10 -100 × above the local relationship (Pacucci et al. 2023). At face value, this implies that these black holes are over-massive 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  \nin a flux-limited survey (Lauer et al. 2007), is likely exacerbated by the high-redshift frontier.  \nMeanwhile, the maximal stellar model predicts an AGN continuum orders of magnitudes lower than the galaxy continuum, and hence implies much lower black hole masses using the L 5100 relation, although as noted above, it is difficulty to know how to calculate the black hole mass when the SEDs are so abnormal. Alternatively, since this model effectively predicts starlight dominating over all the observed wavelength range, it may be possible that these galaxies do not host AGNs at all, which is reasonable given the non-detection in Xrays (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 SNe would be consistently associated with the other spectral characteristics of these objects, and such interpretation may be in some tension with the lack of variability seen in LRDs (Maiolino et al. 2024).  \nWhile 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.', '7. DISCUSSION AND CONCLUSIONS': 'Having discussed the implications of the physical properties and the formation histories resulting from each model, we now attempt to tie all the pieces together. We begin this section by briefly summarizing the modeling results, and then examine two potential interpretations of the nature of these objects.', '7.1. A Brief Summary of the Implications of Different Stellar Contributions': '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 is comparable to that of all the UV-selected objects at these redshifts (Figure 6). Such early and efficient formation would be in tension  \nwith the standard assumption on the baryon to stellar conversion efficiency (Figure 7; also see 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.  \nIn 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 a M ⋆ -M BH relation in less tension with (or requiring fewer mergers to become consistent with) the local relations (Figure 8).', '7.2. Progenitors of Massive Quiescent Galaxies, or Low-mass Galaxies with Bright Black Holes?': "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 is underpinned by their high α -element abundances and their old (but unresolved) ages (e.g., Thomas et al. 2005; Conroy et al. 2014), and from number density analysis of massive, dense galaxy cores at z = 0 -3 (van Dokkum et al. 2015).  \nFurther 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. 2024a) 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 11 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.  \nThe above gap in the observed evolution of massive quiescent galaxies could be bridged by these objects, if  \nthe 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 -5 Mpc -3 (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).  \nAnother distinguishing property of the class of lowz 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 paper) yields similar stellar densities to the cores observed at later times, albeit with significantly smaller sizes, by a factor of 10.  \nIn addition, the enhanced SFR inferred from the maximal stellar model would corroborate to 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 as the superstar clusters with unusually high cloud surface densities in the local universe (Smith et al. 2006; Turner et al. 2015).  \nFinally, 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.  \nTherefore, we conclude that one possible interpretation is that the high-redshift Balmer break sample presented in this work are the progenitors of the first massive galaxies, observed directly after their rapid coformation 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.  \nTo 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 as our sample, are under-luminous in X-ray (Kocevski et al. 2023; Maiolino et al. 2023; Furtak et al. 2024; Greene et al. 2024; Matthee et al. 2024; Wang et al. 2024a). More recently, most LRDs are found to be under-detected 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.  \nSecond, 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.  \nThird, 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. 2024a), 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.,  \nsegregation 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.  \nWhile 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 highz Balmer breaks in this paper-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 categorythe 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 section below.", '7.3. Key Remaining Questions About the Nature of the Balmer Break Sample': "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 compared to 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, make such an argument less satisfying. However, it remains possible that the color-based selection function selects for a particular face-on orientation angle.  \nAs for the M ⋆ -M BH relation, we cautioned about the selection bias (Lauer et al. 2007) in Section 6.6. Recently  \nit has been argued that over-massive black holes are not inconsistent with the local relation after taking in 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 paper.  \nA 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; P'erez-Gonz'alez et al. 2024; Wang et al. 2024a) 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 the minimal stellar models differ by ∼ 3. X-ray detections or firmer upper limits would also help to establish the AGN nature of these objects.  \nFinally, 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) over-massive black holes and massive quiescent systems. In light of all the complexities, the statements in this paper should all be considered contingent. There certainly remains ample room for reinterpretation in the future.", '7.4. Final Remarks': 'JWST is revolutionizing our knowledge of the formation of early galaxies and their black holes. In this paper, we report a remarkable discovery of prominent Balmer breaks as early as z spec = 8 . 35, and intriguingly, its 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  \nby 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.', 'ACKNOWLEDGMENTS': "We thank the anonymous referee for the helpful comments. B.W. and J.L. acknowledge support from JWSTGO-04233.009-A. The Cosmic Dawn Center is funded by the Danish National Research Foundation (DNRF) under grant #140. This research was supported by the International Space Science Institute (ISSI) in Bern, through ISSI International Team project #562 (First Light at Cosmic Dawn: Exploiting the James Webb Space Telescope Revolution). This work is based in part on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. Computations for this research were performed on the Pennsylvania State University's Institute for Computational and Data Sciences' Roar supercomputer. This publication made use of the NASA Astrophysical Data System for bibliographic information.  \nFacilities: HST (ACS, WFC3), JWST (NIRCam, NIRSpec).  \nSoftware: Astropy (Astropy Collaboration et al. 2013, 2018, 2022), dynesty (Speagle 2020), EAzY (Brammer et al. 2008), emcee (Foreman-Mackey et al. 2013), Matplotlib (Hunter 2007), msaexp (Brammer 2023b), msafit (de Graaff et al. 2024b), NumPy (Harris et al. 2020), Prospector (Johnson et al. 2021), Python-FSPS (Johnson et al. 2023).", 'A. AN AGE INDICATOR FOR HIGH-REDSHIFT GALAXIES': 'The commonly used age indicator, D 4000 (Bruzual A. 1983; Balogh et al. 1999), measures the 4000 ˚ A break, which results from the blanket absorption of high energy radiation from metals in stellar atmospheres and  \na deficiency of hot, blue stars. In this paper, we introduce a new definition, spectral break strength, to quantify the Balmer break strength using two windows at [3620,3720] ˚ A and [4000,4100] ˚ A. While Binggeli et al. (2019) proposed using the fluxes in [3400,3600] ˚ A and [4150,4250] ˚ A, the larger separation means that this definition is more sensitive to the overall slope of the spectrum. The wavelength windows used in this paper, however, avoid contamination from strong nebular line emission and are close enough to minimize the impact of dust attenuation and the overall spectrum slope.  \nTo gain more intuition on our new definition, we measure the D 4000 as defined in Balogh et al. (1999) on a set of model spectra, drawn from FSPS (Conroy & Gunn 2010). We assume solar metallicity, and exclude nebular emission lines. We also smooth the spectra to the JWST/Prism resolution. The comparison, illustrated in Figure A.1, shows that the difference is age-dependent, and differ mostly for galaxies with ages ≲ 1 . 5 Gyr. This is expected, as the Balmer break is more prominent in these younger galaxies than the 4000 ˚ A break. Our new age indicator is thus particularly suited for the highredshift universe.', 'B. INFERRED SPECTRAL BREAK STRENGTHS': 'As a test for the constraining power of the data on the formation history, we examine the residuals around the spectral break resulting from the maximal, medium, and minimal stellar mass models, and compare the corresponding spectral break strength measured on the total model spectra (as opposed to on the galaxy spectra as is done in Section 5). As seen from Figure B.1, the minimal stellar mass model predicts a weaker break strength than the data suggest in all cases, implying that the data are sufficiently constraining for the purpose of inferring an extended formation history.  \nThat being said, estimating a tight lower limit on the age of the stellar population is difficult at the current stage. First, the available data, as demonstrated throughout this paper, cannot break the degeneracies among the stellar populations, black hole properties, and dust attenuation. Second, while we have developed models that can describe the observed UV-optical SED, a full understanding of the underlying physical picture is still lacking. The most promising prospect may come from building a complete multi-wavelength view, leveraging the instrumental capability of e.g., Chandra, JWST/MIRI, and ALMA.', 'REFERENCES': 'Ananna, T. T., Bogd\'an, \' A., Kov\'acs, O. E., Natarajan, P., & Hickox, R. C. 2024, arXiv e-prints, arXiv:2404.19010. https://arxiv.org/abs/2404.19010  \nAretxaga, I., Joguet, B., Kunth, D., Melnick, J., & Terlevich, R. J. 1999, ApJL, 519, L123, doi: 10.1086/312114  \nFigure A.1. Comparison between D 4000 (Balogh et al. 1999) probing the 4000 ˚ A break and the new spectral break strength (SBS) definition proposed in this paper, which is designed for the Balmer break. The two panels to the left show normalized model spectra at t age = 300 Myr and 1000 Myr, respectively. The red color illustrates the wavelength window ([4000,4100] ˚ A) used by both Balogh et al. (1999) and this paper. Our definition differs in the bluer window. The blue color indicates the wavelength window ([3850,3950] ˚ A) used by Balogh et al. (1999), whereas the cyan color indicates the wavelength window ([3620,3720] ˚ A) used by this paper. The values of the two age indicators are annotated in the upper left corner. The right panel contrasts the two definitions, color-coded by the age of the model spectra. The age-dependent variation is driven by the 4000 ˚ A break being less visible than the Balmer break in galaxies with ages ≲ 1 . 5 Gyr.  \n<!-- image -->  \nAssef, R. J., Denney, K. D., Kochanek, C. S., et al. 2011,  \nApJ, 742, 93, doi: 10.1088/0004-637X/742/2/93 Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068 Astropy Collaboration, Price-Whelan, A. M., Sip"ocz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74 Atek, H., Chemerynska, I., Wang, B., et al. 2023, MNRAS, 524, 5486, doi: 10.1093/mnras/stad1998 Backhaus, B. E., Trump, J. R., Pirzkal, N., et al. 2024, ApJ, 962, 195, doi: 10.3847/1538-4357/ad1520 Baggen, J. F. W., van Dokkum, P., Labb\'e, I., et al. 2023, ApJL, 955, L12, doi: 10.3847/2041-8213/acf5ef Bagley, M. B., Finkelstein, S. L., Koekemoer, A. M., et al. 2023, ApJL, 946, L12, doi: 10.3847/2041-8213/acbb08 Baldassare, V. F., Reines, A. E., Gallo, E., et al. 2016, ApJ, 829, 57, doi: 10.3847/0004-637X/829/1/57 Balogh, M. L., Morris, S. L., Yee, H. K. C., Carlberg, R. G., & Ellingson, E. 1999, ApJ, 527, 54, doi: 10.1086/308056 Barro, G., Trump, J. R., Koo, D. C., et al. 2014, ApJ, 795, 145, doi: 10.1088/0004-637X/795/2/145 Barro, G., P\'erez-Gonz\'alez, P. G., Kocevski, D. D., et al. 2024, ApJ, 963, 128, doi: 10.3847/1538-4357/ad167e Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57, doi: 10.1088/0004-637X/770/1/57 Bezanson, R., van Dokkum, P. G., Tal, T., et al. 2009, ApJ, 697, 1290, doi: 10.1088/0004-637X/697/2/1290  \nBinggeli, C., Zackrisson, E., Ma, X., et al. 2019, MNRAS, 489, 3827, doi: 10.1093/mnras/stz2387 Blandford, R. D., & McKee, C. F. 1982, ApJ, 255, 419, doi: 10.1086/159843 Boylan-Kolchin, M. 2023, Nature Astronomy, 7, 731, doi: 10.1038/s41550-023-01937-7 Brammer, G. 2023a, grizli, 1.9.11, Zenodo, doi: 10.5281/zenodo.1146904 -. 2023b, msaexp: NIRSpec analyis tools, 0.6.17, Zenodo, doi: 10.5281/zenodo.7299500 Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503, doi: 10.1086/591786 Bruzual A., G. 1983, ApJ, 273, 105, doi: 10.1086/161352 Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682, doi: 10.1086/308692 Carnall, A. C., McLure, R. J., Dunlop, J. S., et al. 2023, Nature, 619, 716, doi: 10.1038/s41586-023-06158-6 Casey, C. M., Akins, H. B., Shuntov, M., et al. 2024, ApJ, 965, 98, doi: 10.3847/1538-4357/ad2075 Chabrier, G. 2003, PASP, 115, 763, doi: 10.1086/376392 Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102, doi: 10.3847/0004-637X/823/2/102 Conroy, C., Graves, G. J., & van Dokkum, P. G. 2014, ApJ, 780, 33, doi: 10.1088/0004-637X/780/1/33 Conroy, C., & Gunn, J. E. 2010, ApJ, 712, 833, doi: 10.1088/0004-637X/712/2/833 Conroy, C., & van Dokkum, P. G. 2012, ApJ, 760, 71,  \ndoi: 10.1088/0004-637X/760/1/71  \nFigure B.1. Zoom-in of the spectral break region. In each row, we show the regions near the spectral break assuming the maximal, medium, and minimal stellar mass model from left to right. The residuals show minor differences, suggesting that all models provide a statistically acceptable solution. The spectral break strength (SBS), calculated as the flux ratio in the two wavelength windows in the blue and red shading, is indicated in the upper left corner in each panel. The minimal stellar mass model always predicts a weaker break strength than the data suggest, implying that the inferred age of the stellar population can be approximately taken as a lower limit.  \n<!-- image -->  \nWeibel, A., Oesch, P. A., Barrufet, L., et al. 2024, arXiv e-prints, arXiv:2403.08872,  \ndoi: 10.48550/arXiv.2403.08872  \nWeiner, B. J., Willmer, C. N. A., Faber, S. M., et al. 2006, ApJ, 653, 1027, doi: 10.1086/508921  \nWellons, S., Torrey, P., Ma, C.-P., et al. 2015, MNRAS,  \n449, 361, doi: 10.1093/mnras/stv303  \nWilliams, C. C., Alberts, S., Ji, Z., et al. 2023, arXiv e-prints, arXiv:2311.07483,  \ndoi: 10.48550/arXiv.2311.07483  \nWorthey, G., Faber, S. M., Gonzalez, J. J., & Burstein, D.  \n1994, ApJS, 94, 687, doi: 10.1086/192087 Wuyts, S., Forster Schreiber, N. M., Wisnioski, E., et al. 2016, ApJ, 831, 149, doi: 10.3847/0004-637X/831/2/149 Yue, M., Eilers, A.-C., Annana, T. T., et al. 2024, arXiv e-prints, arXiv:2404.13290, doi: 10.48550/arXiv.2404.13290  \nZibetti, S., Charlot, S., & Rix, H.-W. 2009, MNRAS, 400, 1181, doi: 10.1111/j.1365-2966.2009.15528.x'}