We study from both a theoretical and observational perspective the physical origin and spectroscopic impact of extreme nebular emission in highredshift galaxies. The nebular continuum which can appear during extreme starbursts is of particular importance as it tends to redden UV slopes and has a significant contribution to the UV luminosities of galaxies. Furthermore its shape can be used to infer the gas density and temperature of the ISM. First we provide a theoretical background showing how different stellar populations SPS models IMFs and stellar temperatures and nebular conditions impact observed galaxy spectra. We demonstrate that for systems with strong nebular continuum emission 1 UV fluxes can increase by up to 0.7magnitudes or more in the case of hotmassive stars above the stellar continuum which may help reconcile the surprising abundance of bright highredshift galaxies and the elevated UV luminosity density at zgt10 2 at high gas densities UV slopes can redden from betalesssim2.5 to betasim1 3 observational measurements of xiion are grossly underestimated and 4 UV downturns from twophoton emission can masquerade as DLAs. Second we present a dataset of 58 galaxies observed with NIRSpec on JWST at 2.5ltzlt9.0 that are selected to have strong nebular continuum emission via the detection of the Balmer jump. Five of the 58 spectra are consistent with being dominated by nebular emission exhibiting both a Balmer jump and a UV downturn consistent with twophoton emission. For some galaxies this may imply the presence of hot massive stars and a topheavy IMF. We conclude by exploring the properties of spectroscopically confirmed zgt10 galaxies finding that UV slopes and UV downturns are in some cases redder or steeper than expected from SPS models which may hint at more exotic e.g. hottermore massive stars or AGN ionizing sources.
['Astrophysics - Astrophysics of Galaxies', 'Astrophysics - Cosmology and Nongalactic Astrophysics']
21 Balmer Jump Street The Nebular Continuum at High Redshift and Implications for the Bright Galaxy Problem UV Continuum Slopes and Early Stellar Populations
{'21 BALMER JUMP STREET: THE NEBULAR CONTINUUM AT HIGH REDSHIFT AND IMPLICATIONS FOR THE BRIGHT GALAXY PROBLEM, UV CONTINUUM SLOPES, AND EARLY STELLAR POPULATIONS': 'Harley Katz 1 , Alex J. Cameron 2 , Aayush Saxena 2 , Laia Barrufet 3 , Nicholas Choustikov 2 , Nikko J. Cleri 4 , 5 , 6 , Anna de Graaff 7 , Richard S. Ellis 8 , Robert A. E. Fosbury 9 , Kasper E. Heintz 10 , 11 , 12 , Michael Maseda 13 , Jorryt Matthee 14 , Ian McConachie 13 , 15 , Pascal A. Oesch 12 , 10 , 11  \n- 1 Department of Astronomy & Astrophysics, University of Chicago, 5640 S Ellis Avenue, Chicago, IL 60637, USA 2 Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford, OX1 3RH, UK 3 Institute for Astronomy, University of Edinburgh, Royal Observatory, Edinburgh, EH9 3HJ 4 Department of Astronomy and Astrophysics, The Pennsylvania State University, University Park, PA 16802, USA 5 Institute for Computational and Data Sciences, The Pennsylvania State University, University Park, PA 16802, USA 6 Institute for Gravitation and the Cosmos, The Pennsylvania State University, University Park, PA 16802, USA 7 Max-Planck-Institut fur Astronomie, Konigstuhl 17, D-69117, Heidelberg, Germany 8 Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK 9 European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching, Germany 10 Cosmic Dawn Center (DAWN), Denmark 11 Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200 Copenhagen N, Denmark. 12 Department of Astronomy, University of Geneva, Chemin Pegasi 51, 1290 Versoix, Switzerland. 13 Department of Astronomy, University of Wisconsin-Madison, 475 N. Charter St., Madison, WI 53706 USA 14 Institute of Science and Technology Austria (ISTA), Am Campus 1, 3400 Klosterneuburg, Austria and 15\n- Department of Physics and Astronomy, University of California, Riverside, 900 University Avenue, Riverside, CA 92521, USA Version August 8, 2024', 'ABSTRACT': 'We study, from both a theoretical and observational perspective, the physical origin and spectroscopic impact of extreme nebular emission in high-redshift galaxies. The nebular continuum, which can appear during an extreme starburst, is of particular importance as it tends to redden UV slopes and has a significant contribution to the UV luminosities of galaxies. Furthermore, its shape can be used to infer the gas density and temperature of the interstellar medium. First, we provide a theoretical background, showing how different stellar populations (SPS models, initial mass functions (IMFs), and stellar temperatures) and nebular conditions impact observed galaxy spectra. We demonstrate that, for systems with strong nebular continuum emission, 1) UV fluxes can increase by up to 0.7 magnitudes (or more in the case of hot/massive stars) above the stellar continuum, which may help reconcile the surprising abundance of bright high-redshift galaxies and the elevated UV luminosity density at z ≳ 10, 2) at high gas densities, UV slopes can redden from β ≲ -2 . 5 to β ∼ -1, 3) observational measurements of ξ ion are grossly underestimated, and 4) UV downturns from two-photon emission can masquerade as damped Ly α systems. Second, we present a dataset of 58 galaxies observed with NIRSpec on JWST at 2 . 5 < z < 9 . 0 that are selected to have strong nebular continuum emission via the detection of the Balmer jump. Five of the 58 spectra are consistent with being dominated by nebular emission, exhibiting both a Balmer jump and a UV downturn consistent with two-photon emission. For some galaxies, this may imply the presence of hot massive stars and a top-heavy IMF. We conclude by exploring the properties of spectroscopically confirmed z > 10 galaxies, finding that UV slopes and UV downturns are in some cases redder or steeper than expected from SPS models, which may hint at more exotic (e.g. hotter/more massive stars or AGN) ionizing sources.', '1. INTRODUCTION': 'The successful launch of JWST (Gardner et al. 2006) has opened a new frontier to study galaxy formation at the earliest cosmic epochs. Understanding the physical properties of these early galaxies, how they differ from those in the low-redshift Universe, and how they impacted subsequent galaxy evolution remain key open goals of modern extragalactic astronomy. One of the primary advantages of JWST compared to HST is its spectroscopic capabilities in the rest-frame UV and optical at high-redshift. This tool can be used not only to confirm the redshifts estimated from photometry, but also to obtain deep insights into the stellar populations, star formation histories, and the conditions of the interstellar medium (ISM) of galaxies in the early Universe (e.g. Sanders et al. 2023; Cameron et al. 2023b).  \nNumerous galaxies have now been spectroscopically  \nconfirmed at z > 10 (e.g. Curtis-Lake et al. 2023; Wang et al. 2023; Arrabal Haro et al. 2023; Carniani et al. 2024) and have been shown to exhibit a diverse set of properties. For example GNz11 (Bunker et al. 2023b) and GHZ2 (Castellano et al. 2024; Zavala et al. 2024) are characterized by extreme UV emission lines, including those rarely seen at lower redshifts (e.g. N IV ] λ 1486), while others exhibit much lower EW (and undetected) UV lines (e.g. Carniani et al. 2024; Curtis-Lake et al. 2023). Compilations of high-redshift galaxies (e.g. Heintz et al. 2024b; Roberts-Borsani et al. 2024) show that UV slopes can range from close to -3, which would be consistent with the intrinsic values of massive stars, to redder than -2, which is typically reconciled by assuming dust attenuation, clearly indicating a difference in ISM properties or stellar populations. With JWST, we are seeing unprecedented spectral features at high-redshift.  \nMore generally, the photometric detection and spectroscopic confirmation of numerous bright galaxies at z ≳ 10 as well as the elevated UV luminosity density at similar redshifts is in tension with most theoretical models of high-redshift galaxy formation (e.g. Finkelstein et al. 2023; Harikane et al. 2024; Leung et al. 2023; Chemerynska et al. 2023, although c.f. Willott et al. 2023). Various physical mechanisms have been proposed to explain this overabundance of bright (M UV ≳ -19) galaxies such as:  \n- 1. an increased scatter in the relation between dark matter halo mass and UV luminosity compared to the lower-redshift Universe, possibly driven by bursty star formation (e.g. Ren, Trenti & Mason 2019; Mason, Trenti & Treu 2023; Shen et al. 2023; Sun et al. 2023; Kravtsov & Belokurov 2024),\n- 2. an increased star formation efficiency at high redshift (e.g. McCaffrey et al. 2023), potentially caused by weak feedback (e.g. Dekel et al. 2023),\n- 3. a lack of dust in massive galaxies due to strong radiation pressure (Ferrara, Pallottini & Dayal 2023),\n- 4. or a change in the mass-to-light ratio of early stellar populations, possibly due to a top-heavy IMF (e.g. Yung et al. 2024, although c.f. Cueto et al. 2024).  \nThese solutions address the problem by either changing the observed mass-to-light ratio or forming stars more efficiently. The true origin of the bright galaxy problem and excess UV luminosity density is likely a combination of many physical effects. For example, nearly all high-resolution numerical simulations of high-redshift galaxies that model cold gas in the ISM predict that star-formation is bursty (e.g. Ma et al. 2018; Katz et al. 2023; Pallottini et al. 2022). However, as this scatter is mass-dependent (e.g. Gelli, Mason & Hayward 2024; Kravtsov & Belokurov 2024) it remains an open question whether the level of SFR burstiness generically expected in galaxy formation models is sufficient to produce the level of UV luminosity variability required to match observations (e.g. Pallottini & Ferrara 2023).  \nNumerous 1D (e.g. Omukai et al. 2005) and 3D models (e.g. Chon et al. 2022) predict that the stellar initial mass function (IMF) becomes more top-heavy at lower metallicities and in the early Universe when the cosmic microwave background (CMB) temperature is higher. A more top-heavy IMF is a generic prediction in the extreme limit of metal-free star formation at high redshift (e.g. Bromm, Coppi & Larson 2002; Abel, Bryan & Norman 2002; Hirano et al. 2014). The mass-to-light ratio is lower at young stellar ages for IMFs with flatter high-mass slopes, which would lead to an increase in UV luminosity at fixed star formation rate.  \nWhat has perhaps been considered less in the literature is the role of the nebular continuum, which is observed to become important in lower mass, low-metallicity environments (e.g. Izotov, Guseva & Thuan 2011). When ionizing photons are absorbed by gas they can be reprocessed into continuum emission through three different processes: free-free, free-bound, and two-photon 1 . For  \n1 Note that throughout this work we will only consider the nebular continuum contribution as we will analyze spectroscopic data. For photometric data, line emission can also contribute to the inferred luminosity integrated over a photometric passband.  \nFig. 1.Spectral shape of the three components of the nebular continuum as a function of wavelength. We show the shape of the nebular continuum for gas at 2 × 10 4 K and an electron density of 10 2 cm -3 assuming Case B recombination. Annotated are the locations of the Balmer jump caused by free-bound emission and the UV downturn from two-photon emission.  \n<!-- image -->  \na typical star-forming galaxy, nebular continuum emission is generally subdominant at rest-frame wavelengths of 1500 ˚ A, where UV luminosities are typically measured. However, strong nebular continuum emission has now been detected in both low (e.g. Guseva, Izotov & Thuan 2006; Guseva et al. 2007) and high-redshift galaxies (e.g. Cameron et al. 2023a; Roberts-Borsani et al. 2024; Welch et al. 2024) undergoing bursts of star formation.  \nFree-bound emission is arguably the most easily recognizable due to the Balmer jump, which is a sharp discontinuity at rest-frame 3645 ˚ A (see the pink line in Figure 1). Two-photon emission, particularly the 2 s → 1 s transition of neutral hydrogen peaks very close to 1500 ˚ A (in f λ ) and smoothly drops towards both shorter and longer wavelengths, going to zero at 1216 ˚ A (see the green line in Figure 1). In addition to being emitted by H II regions, two-photon emission can also arise from warm cooling flows that may be present in the high-redshift Universe (e.g. Dijkstra 2009). Two-photon emission is rarely directly detected because for Population II star formation, the power-law stellar continuum increases towards the blue end of the spectrum and typically outshines the nebula at short wavelengths. However, in low-metallicity environments where H II region temperatures are elevated, Case B departures can occur which increases both the two-photon and Ly α emission via collisional excitation (e.g. Raiter, Schaerer & Fosbury 2010; Mas-Ribas, Dijkstra & Forero-Romero 2016).  \nIn contrast to two-photon emission, free-bound emission increases towards longer wavelengths and can thus outshine the stellar emission. Hence it can be directly detected in the rest-frame optical. The key aspect of the nebular continuum is that it should only become important when galaxies are undergoing a burst of star formation. The natural UV magnitude fluctuations driven by deviations from the star-forming main-sequence can be further enhanced by continuum emission from the gas. Identifying a sample of high-redshift galaxies with strong  \nnebular continuum emission may provide an ideal testbed for elucidating the physical mechanisms that drive extreme bursts of star formation, which may explain the over-abundance of bright high-redshift galaxies and the excess UV luminosity density.  \nUnfortunately, a sample of high-redshift galaxies with detected nebular continuum has yet to be compiled. Only one individual galaxy with a Balmer jump has been reported at z > 5 (Cameron et al. 2023a) and this same feature has appeared in a stack of z ∼ 6 Ly α emitters (Roberts-Borsani et al. 2024).  \nThe lack of reported Balmer jumps is perhaps surprising because high-resolution simulations of early galaxies generically predict bursty star formation histories such that the SFRs can fluctuate by multiple orders of magnitude on tens of Myr timescales (e.g. Ma et al. 2018; Katz et al. 2023; Pallottini et al. 2022), which is exactly the scenario where the nebular continuum is expected to be strong. SED fitting of the photometry of highredshift galaxy candidates often predicts extreme nebular emission (e.g. Endsley et al. 2023; Bradaˇc et al. 2024). Moreover, the flux-limited nature of high-redshift surveys preferentially biases observations towards galaxies that are highly star-forming, which is when the nebular continuum should begin to appear. The spectral resolution of the NIRSpec prism and sensitivity of the grating may also impact our ability to detect the Balmer jump if present. As the prevalence and characteristics of the nebular continuum at high-redshift are poorly known, its importance for the luminosities of early galaxies remains to be quantified.  \nIn this work, we aim to better understand the physics of galaxies undergoing extreme bursts of star formation by focusing on the spectral impact of extreme nebular continuum emission and the underlying physics that may be driving its observability. We begin in Section 2 by providing a theoretical background for how nebular emission modifies the observational properties of starbursting galaxies, focusing on changes to the properties of the underlying stellar populations and ISM conditions. In Section 3, we test the theoretical results in the context of a newly compiled sample of 58 spectroscopically confirmed high-redshift galaxies that show Balmer jumps. Finally, in Sections 4 and 5, we discuss the implications our findings in the context of z > 10 galaxies and present our conclusions.', '2. THEORETICAL BACKGROUND': 'We begin by providing a theoretical overview of the expected nebular continuum contribution of different stellar populations to the observed spectra of star-bursting galaxies and how this emission impacts our ability to infer the underlying properties of these same galaxies.', '2.1. Spectral synthesis modeling': 'SEDs of galaxies including the nebular continuum are computed with CLOUDY version C23 (Chatzikos et al. 2023). Unless otherwise specified, models assume a slab configuration with a dimensionless ionization parameter of log 10 U = -2 . 0, a metallicity of 1% solar by scaling down abundance ratios from Asplund et al. (2009), and a hydrogen gas density of 10 3 cm -3 , consistent with that of galaxies in the early Universe (e.g. Isobe et al.  \n2023). Calculations are stopped when the electron fraction reaches 1% 2 . The exact details of the models have little impact on our conclusions. We consider five different stellar population synthesis (SPS) models, Starburst99 (Leitherer et al. 1999) 3 , BPASS v2.2.1 (Stanway & Eldridge 2018) 4 , FSPS (Conroy, Gunn & White 2009; Conroy & Gunn 2010) 5 , C&B 2019 models 6 , and individual, metal-poor, massive star models of Larkin, Gerasimov & Burgasser (2023) 7 . Within each subset of models, we choose the SEDs that have the closest metallicity to the fiducial value of 1% that we adopt for the gas 8 . All models assume an instantaneous burst of star formation and we consider stellar ages up to 20 Myr. Beyond this, the ionizing output is so low that the nebular emission is negligible. Note that the upper mass limits for all of our chosen SPS models differ. The fiducial upper mass cutoffs are 120 M ⊙ , 300 M ⊙ , and 100 M ⊙ for Starburst99 , BPASS , and FSPS , respectively. For the C&B19 models, we consider three different upper mass limits of 100 M ⊙ , 300 M ⊙ , and 600 M ⊙ .  \nThe output from these models are the transmitted stellar emission as well as the continuum and line emission from the gas. It should be noted that in real galaxies, the turbulent nature of H II regions can modify the gas emission from expectations of spherical models (e.g. Jin, Kewley & Sutherland 2022) and interpretation biases can arise if multiple H II regions with different properties combine to produce the spectrum of a galaxy (e.g. Cameron, Katz & Rey 2023). Nevertheless, the 1D models used here provide both quantitative and qualitative insights into how galaxy spectra may change under different conditions under a restrictive set of assumptions.', '2.2. How strong can the nebular continuum be?': 'We first address the question of how much the nebular continuum can contribute to the UV luminosity of a galaxy.', '2.2.1. Standard SPS models': 'In Figure 2 we show the fraction of the total 1500 ˚ A or 3640 ˚ A continuum emission from the galaxy comes from the nebula as a function of stellar age for four SPS models. For young stellar ages, ∼ 1 Myr, all models predict that the nebula contributes at least 30% to the 1500 ˚ Aluminosity, which corresponds to brightening the galaxy by ∼ 0 . 4 magnitudes. The nebular continuum contribution can be > 70% in the rest-frame optical. All SPS models predict that after ∼ 10 Myr, the nebular contribution at 1500 ˚ A is negligible, a consequence of the massive stars  \n2 We consider this to be the edge of an ionization bounded nebula. Adopting a density bounded model would reduce the contribution of nebular emission with respect to the stellar continuum as the escape fraction will be non-zero. Moreover, extending the calculation beyond the Stromgren sphere will not impact our results unless the gas is hot enough so that it can be collisionally excited. This possibility is further discussed below.  \n- 3 https://www.stsci.edu/science/starburst99/docs/default.htm\n- 4 https://bpass.auckland.ac.nz/9.html\n- 5 https://github.com/cconroy20/fsps\n- 6 http://www.bruzual.org\n- 7 https://atmos.ucsd.edu  \n8 Note that we do not assume any α enhancement in stellar SEDs which is expected both at low metallicities and at high redshifts (e.g. Kobayashi, Karakas & Lugaro 2020; Cullen et al. 2021).  \nFig. 2.Nebular continuum fraction of the total SED at either 1500 ˚ A (magenta), where the two-photon continuum becomes important, or at 3640 ˚ A (cyan), where the free-bound emission reaches its maximum. We show results for four different SPS models as a function of stellar age. For the C&B19 models, different line styles represent different upper mass limits for the IMF.  \n<!-- image -->  \nevolving off the main sequence and the lower ionizing photon production efficiencies of less massive stars. The timescale over which the nebular continuum remains important is slightly longer for the free-bound emission due to the fact that it is stronger with respect to the stellar continuum compared to two-photon emission. We note that burst models maximize the nebular contribution to each SSP and the nebular continuum contribution for realistic SFHs will be explored in Section 2.3.1.  \nWhile the general trends between the SPS curves are similar, in detail the shapes of the curves and the maximum nebular contributions are different. For example the 1500 ˚ A nebular contribution of the BPASS models peak at ∼ 40% at 1 Myr, much higher than the Starburst99 models, due to the increased ionizing photon production efficiencies caused by binary stellar evolution. The C&B19 models with maximum stellar masses of 300 M ⊙ and 600 M ⊙ can reach a 1500 ˚ A nebular contribution as high as 50% at similarly young stellar ages, which corresponds to a brightening the UV luminosity by 0.75 magnitudes. In contrast to the other three SPS models, the nebular contribution of the FSPS models peaks between 3 -5 Myr at nearly 45%, which results from different assumptions regarding Wolf-Rayet stars.', '2.2.2. A top-heavy IMF': "Numerous theoretical models predict that the stellar IMF becomes more top-heavy at high redshift due to lower gas metallicities and an elevated CMB temperature (e.g. Chon et al. 2022; Bate 2023). This has two important effects. First, as the upper-mass IMF slope becomes  \nflatter 9 , at young stellar ages, more UV luminosity is produced per unit star formation, which would help explain the bright galaxy problem by decreasing the mass-tolight ratio at young ages. Second, the ionizing photon production per unit star formation also increases as the IMF becomes more biased towards high-mass stars. To test the impact of a varying IMF, we adopt Starburst99 models 10 , where we have systematically varied the uppermass slope, α , from -2 . 6 (i.e. more bottom heavy than the canonical -2 . 3) to -0 . 5 (where nearly all of the mass is locked in massive stars). We define the IMF such that:  \nξ ( m ) ∝ { m -1 . 3 m ≤ 0 . 5 M ⊙ m α m> 0 . 5 M ⊙ (1)  \nThe lower and upper mass limits of the IMF are set to 0.08 M ⊙ and 120 M ⊙ , respectively.  \nIn the top panel of Figure 3, we show the fractional nebular contribution at 1500 ˚ A and 3640 ˚ A as a function of stellar age. The general trend of the nebular contribution is the same as above - the peak in nebular emission comes at young stellar ages and begins to subside after a few Myr, with a small secondary peak from Wolf-Rayet stars. Unsurprisingly, a more top-heavy IMF  \n9 Here we only consider modifications to the IMF where the power-law slope of the high-mass end is varied. However, the shape of the IMF may be fundamentally different and deviate from a power-law at high-redshift (e.g. Chon et al. 2022). For our work, the important parameter is the fraction of mass in high-mass stars with respect to that at lower masses rather than the exact shape.  \n10 Starburst99 models were adopted because it is trivial to vary the properties of the IMF, in contrast to some of the other SPS models. Our primary conclusions likely generalize to the other SPS models.  \nFig. 3.Contribution of the nebular continuum to the total SED at either 1500 ˚ A (solid), where the two-photon continuum becomes important, or at 3640 ˚ A (dashed), where the free-bound emission reaches its maximum. We show results for Starburst99 models assuming a gas density of 10 3 cm -3 with different stellar initial mass functions where the upper-mass slope is varied as listed in the legend (top) as well as for low-metallicity massive star models of different stellar temperatures (bottom) where we vary the gas density of the nebula. For models with varying stellar temperature, we only consider the zero-age main-sequence spectra.  \n<!-- image -->  \n1  \n.  \n0  \n0  \n.  \n8  \n0  \n.  \n6  \n0  \n.  \n4  \n0  \n.  \n2  \n0  \n.  \n0  \n40000  \n60000  \n80000  \n100000  \nStellar Temperature [K]  \nresults in an increased fractional nebular contribution to the 1500 ˚ A luminosity. While the fiducial model predicts that ∼ 30% of the 1500 ˚ A luminosity at an age of 1 Myr comes from nebular emission, by flattening the upper-mass IMF slope to -0 . 5, the nebular contribution increases to ∼ 38% at the same stellar age. This increase is rather modest, especially compared to the total increase in UV luminosity which is more than a factor of ten for an upper-mass IMF slope to -0 . 5 (see below).  \nIt is straightforward to understand why extreme IMF variations have only a minimal impact on the nebular continuum fraction at 1500 ˚ A. A massive O star has temperature of ∼ 40 , 000 K (Pecaut & Mamajek 2013), which according to Wien's law has a spectral peak at ∼ 725 ˚ A. For a standard IMF, these stars dominate the ionizing  \n.  \nneb  \nf  \n2  \n10  \n3  \n10  \n4  \n10  \n5  \n10  \ncm  \n-  \ncm  \n-  \ncm  \n-  \ncm  \n-  \n3  \n3  \n3  \n3  \noutput as going slightly lower in temperature pushes the peak redward of the hydrogen ionizing wavelength. Since hotter blackbodies are brighter at all wavelengths, the most massive stars also contribute significantly to the 1500 ˚ A UV luminosity, despite both their lower numbers and their spectral peak at shorter wavelengths. The absolute maximum nebular contribution is then that which comes from the most massive star in the IMF. Since it is these stars that provide both the ionizing photons and a significant amount of the 1500 ˚ A luminosity, flattening the upper-mass slope of the IMF (i.e. making it more 'top-heavy') does not substantially increase the fractional nebular contribution because the upper limit to this value is set by the most massive star in the IMF.", '2.2.3. Hotter and more massive stars': 'We have shown that for an IMF with a fixed upper mass, the fractional nebular contribution at 1500 ˚ A does not increase substantially if the upper-mass IMF slope flattens from -2 . 3 to -0 . 5. However, increasing the maximum mass sampled by the IMF represents a possible alternative. This was already hinted at in Section 2.2.1 where we showed for the C&B19 models that increasing the maximum mass from 100 M ⊙ to 600 M ⊙ increased the fractional nebular contribution at both 1500 ˚ A and 3640 ˚ A.  \nTo explore this effect further, we adopt the lowmetallicity, massive star models of Larkin, Gerasimov & Burgasser (2023) which are available for a range of stellar temperatures. We assume that these are generally representative of low-metallicity massive stars in the early Universe for which there is limited public data of stellar atmospheres. The effect on nebular contribution strength due to very massive stars was recently discussed in Schaerer et al. (2024). Our work differs from Schaerer et al. (2024) in that they chose a very specific model for massive stars at a metallicity that is higher than typically measured for high-redshift galaxies and hence the effective temperatures of those models are lower than considered here. Moreover, they apply a more simplistic model for the nebular continuum by assuming a constant temperature, which both ignores the temperature structure of the nebula and neglects potential Case B departures that may become important at low metallicity (e.g. Raiter, Schaerer & Fosbury 2010; Mas-Ribas, Dijkstra & Forero-Romero 2016).  \nIn the bottom panel of Figure 3 we show the fractional nebular contribution at 1500 ˚ A and 3640 ˚ A as a function of stellar temperature. While at a temperature of 40,000 K, the nebular fraction at 1500 ˚ A is only ∼ 14%, for stars with T = 100 , 000 K, nebular emission represents 80% of the total 1500 ˚ A UV luminosity. For these hot stars, the UV luminosity could then be brightened by 1.75 magnitudes, purely due to nebular emission.  \nThis prediction is not unique to our work. In fact, a generic characteristic of hot star models is that they have nebular dominated spectra across nearly the entire UV and optical (e.g. Panagia 2002; Schaerer 2003; Trussler et al. 2023), and Raiter, Schaerer & Fosbury (2010) have already shown that the nebular continuum can dominate the 1500 ˚ A UV luminosity in the context of Pop. III stars under various assumptions of IMF. Raiter, Schaerer & Fosbury (2010) also pointed out that at low metallici-  \nies, one expects departures from standard Case B assumptions due to the fact that the 2 s and 2 p states can be collisionally populated. Case B departures have the effect of enhancing the two-photon emission (and Ly α ), both due to the fact that the collisionally excited atom can fluoresce down to the ground state and that in the excited state, lower energy photons are needed to ionize the gas, although this latter effect is likely subdominant (Mas-Ribas, Dijkstra & Forero-Romero 2016). Case B departures due to collisional excitation will also increase the intrinsic ratio of H α /H β to values > 2 . 86.', '2.2.4. Other physics that impacts the nebular contribution': "Cooling Flows -Just as collisions impact the level populations within H II regions, gas accreting onto haloes can undergo the same process as long as the virial temperature of the halo is high enough ( ≳ 10 4 K) (e.g. Fardal et al. 2001; Faucher-Gigu'ere et al. 2010). This process then cools the gas and is particularly efficient in the temperature range 15 , 000 -20 , 000 K. Numerical simulations predict that cold-mode accretion, where cooling is strongest is much more prevalent at highredshift compared to the local Universe (Kereˇs et al. 2005). A nebular-only spectrum dominated by Ly α and two-photon emission is then expected when cooling is sufficiently strong (e.g. Dijkstra 2009). This process can boost the nebular contribution, particularly Ly α and two-photon emission, which has the effect of decreasing the overall EWs of UV emission lines, reducing the perceived ionizing photon production efficiency, and increasing the UV magnitude of a galaxy. Cooling flows have been used, for example, to explain the presence of Ly α blobs at intermediate redshift (e.g. Goerdt et al. 2010; Faucher-Gigu'ere et al. 2010; Rosdahl & Blaizot 2012) and the relative contribution of cooling radiation to the total Ly α (and also two-photon) emission is predicted to increase substantially at higher redshift (e.g. Dayal, Ferrara & Saro 2010; Yajima et al. 2012). Given the fact that catastrophic cooling is required to drive a starburst it is reasonable to expect that cooling radiation occurs simultaneously with star formation, although the fractional contribution of cooling radiation is likely maximized immediately preceding or just after the starburst (e.g. Smith et al. 2019). The exact details of how much cooling radiation can contribute to the total nebular continuum warrants further investigation under the conditions likely present in early galaxies and under conditions that may increase the electron fraction in the cooling gas (e.g. X-ray radiation).  \nHigh Density Gas -While we have primarily focused on the physical effects that can increase the nebular contribution to the UV luminosity, we also highlight the fact that certain nebular conditions can inhibit two-photon emission. For example, at high gas densities, l -changing collisions can shift electrons between the 2 s to 2 p states. However, because the Einstein coefficient is orders of magnitude larger for Ly α compared to two-photon emission, at high densities, Ly α is enhanced at the expense of two-photon emission. This is demonstrated in the bottom panel of Figure 3 where we show the fractional contribution of nebular emission at 1500 ˚ A at four densities from 10 2 -10 5 cm -3 . Between 10 3 cm -3 and 10 4 cm -3 the nebular contribution begins to decrease and at a den-  \ny of 10 5 cm -3 the decrease is much more significant. Nevertheless, at high stellar temperatures, even at extreme densities of 10 5 cm -3 , which have only been observed in a select few high-redshift galaxies (e.g. Topping et al. 2024c,b; Senchyna et al. 2024), the nebular continuum can still increase the UV luminosity of the galaxy by a factor of two. The electron densities of typical galaxies in the reionization epoch are such that collisional de-excitation is not expected to be very important (e.g. Isobe et al. 2023).", '2.3. Observational implications': 'It is clear that during a starburst with ages ≲ 5 Myr, the nebular continuum can represent a significant fraction of the observed SED of a galaxy. Under specific circumstances, e.g. the presence of very hot massive stars, the nebular continuum can even dominate over the stellar emission at all wavelengths observable by NIRSpec on JWST. Here we explore the observational implications for galaxies with strong nebular continuum emission.', '2.3.1. UV luminosities and implications for the bright galaxy problem': "As we have shown in Section 2.2.1, in some of the SPS models, e.g. the C&B19 models with maximum stellar masses of 300 M ⊙ and 600 M ⊙ , the nebular contribution at 1500 ˚ A can reach as high as 50%. This corresponds to a UV luminosity increase of 0.75 magnitudes. Even in the optimistic scenario of the C&B19 models, the excess UV luminosity cannot explain the full UV variability needed to match the z > 10 UV luminosity functions as estimated by Kravtsov & Belokurov (2024). However, the inclusion of the nebular emission can bring down the implied 'burstiness' of star formation.  \nOne important caveat is that our models assume an instantaneous burst of star formation rather than a realistic star formation history. Older stars can contribute to the UV luminosity while providing no additional ionizing photons that excite the nebular continuum. The general impact of including a realistic star formation history would be to decrease the fractional contribution of the nebular continuum in the UV. To better understand this effect, we perform a simple experiment and apply the C&B19 model with a maximum mass of 600 M ⊙ (our most optimistic model) to the realistic star formation histories of simulated galaxies from the SPHINX public data release (Katz et al. 2023). The stellar masses of these galaxies range between ≳ 10 7 M ⊙ to ≲ 10 10 M ⊙ , and all have 10 Myr-averaged SFRs > 0 . 3 M ⊙ yr -1 . We then account for the nebular contribution from the young stellar populations to the total SED following our fiducial photoionization model. Figure 4 shows a histogram of the change in UV magnitude from including the nebular continuum for > 1 , 400 star-forming galaxies from the SPHINX simulation between 4 . 6 ≤ z ≤ 10. While the typical magnitude increase is only ∼ 0 . 2 magnitudes, galaxies undergoing particularly strong bursts of star formation see much more substantial increases ( > 0 . 5 magnitudes) in UV luminosity due to the inclusion of the nebular continuum. Hence, during the burst, the UV magnitude is enhanced both due to the typical magnitude of the galaxy due to both the increase in SFR and the gaseous emission. We emphasize that the effect of  \nFig. 4.Histogram of the change in 1500 ˚ A UV magnitude of galaxies from the SPHINX simulation when accounting for the nebular contribution to the SED. In all cases, the nebular continuum increases the UV luminosity for galaxies with realistic star formation histories, in some extreme examples by more than 0.5 magnitudes.  \n<!-- image -->  \nthe nebular continuum is unavoidable unless the gas is at extreme densities or the escape fraction is very high.  \nModifying the upper IMF slope of the maximum stellar mass of the IMF can similarly increase the UV luminosity of a galaxy. We demonstrate this in Figure 5. In the top panel, we show the 1500 ˚ A luminosity of the Starburst99 models with different upper-mass IMF slopes as a function of stellar age, normalized by the model with an upper-mass slope of -2 . 3. A modest flattening of the slope from -2 . 3 to -1 . 7 increases the UV luminosity by a factor of five (1.75 magnitudes) at young stellar ages. Recall that under such circumstances, the nebular continuum contribution was not substantially increased and thus it is clearly the modification of the stellar component that leads to such a substantial increase in UV luminosity.  \nIn the bottom panel of Figure 5 we show the ratio of 1500 ˚ A luminosity of stars of different temperatures normalized to a 40,000 K star. Here we find that hotter stars tend to increase the UV luminosity; however, in contrast to the IMF, it is the nebular continuum that drives the substantial increase. At high gas densities, where twophoton emission is strongly suppressed, a 100,000 K star is not even twice as bright as a 40,000 K star. However, at low gas densities, when the nebular continuum is strong in the UV, the brightness can increase by a factor of five. Hence, we emphasize that there is a clear difference between simply assuming a more top-heavy IMF by modifying the upper-mass power-law slope and changing the maximum mass star in the IMF as the physical effects that lead to a luminosity increase are different. Nevertheless, the combination of both a top-heavy IMF and more massive stars can lead to substantial changes in the UV luminosity of early galaxies, especially when the nebular continuum is strong.  \nFig. 5.(Top) 1500 ˚ A UV luminosity as a function of stellar age for Starburst99 models using a density of 10 3 cm -3 with different upper-mass IMF slopes, normalized to the model with a slope of -2 . 3. (Bottom) 1500 ˚ A UV luminosity as a function of stellar temperature for zero-age main-sequence metal-poor massive star models, normalized to a star with a surface temperature of 40,000 K. All models include both the stellar and nebular continuum. In the case of IMF variations, the number of high mass stars drives the UV luminosity increase but in the case of stellar temperature, the nebular continuum is the primary factor that leads to an increase in UV luminosity.  \n<!-- image -->", '2.3.2. UV continuum slopes, β': 'The most common discussion on the implications of a strong nebular continuum is the impact on UV continuum slopes (e.g. Schaerer 2002; Bouwens et al. 2010; Dunlop et al. 2013; Cullen et al. 2024; Topping et al. 2024a). In particular, the nebular continuum acts to redden the observed spectrum. This is demonstrated in Figure 6 where we show the intrinsic (i.e. stellar only) and observed (stellar + nebular) UV slopes for all SPS models. Here we have measured UV slopes in the wavelength range 1400 ˚ A to 2600 ˚ A. For young stellar ages, the intrinsic UV slopes are all bluer than -3 . 0. After the nebular continuum is added to the spectrum, the observed slope reddens to ∼ -2 . 5. This is consistent with other works that have performed similar calculations (e.g. Schaerer 2002; Bouwens et al. 2010; Dunlop et al. 2013; Cullen et al. 2024; Topping et al. 2024a).  \nMaking the IMF more top-heavy only marginally changes the UV slope. The top panel of Figure 6 shows that as the IMF upper mass slope increases from -2 . 6  \n<!-- image -->  \nStellar Age [yr]  \nFig. 6.UV slope as a function of stellar age for various SPS models (top) or stellar temperature for the individual metal-poor star models (bottom). Solid lines represent the UV slope as measured from total (stellar+nebular) emission while dashed lines represent the UV slope of only the stellar component. The dotted red lines indicate a slope of -2 . 3 for reference in all panels. For the Starburst99 models, colors represent different upper-mass IMF slopes (see legend in top panel), while for the C&B19 models, colors indicate different upper-mass limits (see legend in the third panel). For the individual massive star models, colors indicate the gas density of the nebula.  \n<!-- image -->  \nto the extreme value of -0 . 5, β only increases by 0.15 at ages prior to the main sequence lifetime of massive stars. The reason the change is so mild is that the massive stars already dominate both the UV luminosity and the ionizing photon production. Adding more massive stars increases the contribution of these stars to the 1500 ˚ A luminosity, but this effect quickly saturates and the intrinsic β slope only becomes slightly bluer. At the same time, the nebular contribution is rising faster which reddens the slopes. Hence the net effect for a top-heavy IMF is that the UV continuum slope is redder than in the scenario with a more standard IMF (see also discussion in Cullen et al. 2024). We emphasize that these results are relatively independent of the chosen SPS model. Although the exact value of β can change slightly, the trends with IMF will be similar.  \nTo make β even redder than in models with a top-heavy IMF, hotter stars are needed. This is shown in the bottom panel of Figure 6 where UV slope is plotted as a function of stellar temperature. For typical ISM densities of 10 2 -10 3 cm -3 , β saturates at a slope of ∼ -2 . 0 consistent with earlier work that assumed Pop. III SEDs (e.g. Trussler et al. 2023). By increasing gas density, twophoton emission is suppressed and the free-bound emission dominates, which leads to extremely red UV slopes approaching -1. This is extremely important because we do not yet know the ISM conditions around Pop. III stars and selection criteria should not necessarily remove redder objects as the ISM could simply be at high density. Hence we encourage photometric selection criteria to allow for these exotic scenarios.  \nOne important caveat to these calculations is that we have computed the UV slopes without dust attenuation. Many high-redshift galaxies have Balmer decrements and UV slopes that are consistent with a dust-free scenario (e.g. Cullen et al. 2024; Sandles et al. 2023) and dust will only act to redden the SEDs. High quality spectra will be required to differentiate reddening due to dust versus the nebular continuum; however, it should be noted that β is not constant with wavelength when the nebular continuum is strong.  \nThe modification of β due to the nebular continuum has important implications for interpreting spectra. For example, Chisholm et al. (2022) recently proposed using the UV slope as an indicator for LyC escape. In their sample of galaxies, β is largely set by dust attenuation and the galaxies with less dust are more likely to be LyC leakers. Applying their model to our dust-free CLOUDY simulations where β can be as low at ∼ -2 . 55, one would predict an escape fraction of ∼ 17% at β = -2 . 55. Therefore, in the zero-dust limit, this relation should be revised or not applied (see also Choustikov et al. 2024) while at higher metallicities, when dust is present, the relation likely remains applicable.  \nAnother key result is that under normal nebular conditions and for typical SPS models, the UV slope should be bluer than approximately -2 . 3 for star-bursting galaxies 11 . This is indicated as the dotted red line on each  \n11 We emphasize that here we are discussing star-bursting galaxies. For older stellar populations, the UV slope will naturally redden. Based on the realistic star formation histories in the SPHINX simulation (Katz et al. 2023), even for galaxies with mass-weighted ages of 100 Myr, the UV slope measured from stellar+nebular emission remains ∼ -2 . 5.  \npanel of Figure 6 and represents the maximum slope that one obtains for standard SPS models. Redder slopes imply either dust or an increased nebular fraction (e.g. due to hotter/more massive stars), the former of which can be tested using the ratio of Balmer lines.', '2.3.3. Ionizing photon production efficiencies, ξ ion': 'The ionizing photon production efficiency ξ ion , is defined as  \nξ ion = Q/L ν, 1500 ˚ A erg -1 Hz , (2)  \nwhere Q is the number of ionizing photons produced per second and L ν, 1500 ˚ A is the 1500 ˚ A monochromatic luminosity. Note that throughout this work, we often refer to log 10 ( ξ ion ). As Q cannot be directly measured observationally, it is often inferred from emission lines (e.g. Balmer lines) using a conversion factor based on recombination theory (see Appendix B or e.g. Osterbrock & Ferland 2006).  \nBecause the nebular continuum can significantly impact the UV luminosity at 1500 ˚ A, measurements of ξ ion , which often normalize a Balmer line (e.g. H α ) to the 1500 ˚ A luminosity can also be impacted. More specifically, the ionizing photon production efficiency is intended to be a measurement that traces the intrinsic properties of the underlying stellar (or more generally photoionizing) population. However, if the nebular continuum represents a large fraction of the UV luminosity, these measurements will be systematically biased low.  \nThis is demonstrated in the left panel of Figure 7. For each SPS model, we adopt the fiducial CLOUDY model and measure the intrinsic ξ ion using the number of emitted ionizing photons ( Q ) by the intrinsic stellar continuum at 1500 ˚ A as well as the observed ξ ion where we add the nebular continuum to L UV . Independent of the chosen SPS model, the trend is nearly identical. For log 10 ( ξ ion / erg -1 Hz) ≲ 25 . 5 the intrinsic and observed ξ ion agree. This is a result of the fact that when ξ ion is low, there are not enough ionizing photons to create a significant nebular continuum. However, at log 10 ( ξ ion / erg -1 Hz) ≳ 25 . 5 the observed ξ ion deviates low compared to the intrinsic value. For the highest ξ ion values among all of our SPS models, this deviation can be nearly 1 dex. The primary conclusion here is that, while the different SPS models vary in terms of maximum intrinsic ξ ion , at high ionizing photon production efficiencies, the ξ ion value one would measure deviates from the true value in the exact same way in all cases.  \nThe reason why the underlying SPS models have very little impact on the relation between intrinsic and observed ξ ion is because the contribution of the nebular continuum to the 1500 ˚ A UV luminosity is primarily sensitive to the properties of the nebula. In the right panel of Figure 7 we adopt the hot star models of Larkin, Gerasimov & Burgasser (2023) 12 and vary the details of the CLOUDY models such that gas density is sampled in the range 10 2 -10 4 cm -3 , the log of the ionization parameter is varied from between -2 . 5 and -1 . 5, while  \nmetallicity ranges between 10 -3 -10 -1 Z ⊙ , all of which are conditions that are thought to be generally representative of the high-redshift Universe. The models show almost no sensitivity to ionization parameter (compare line styles), while there is a weak dependence on metallicity (compare line opacities). The primary impact of metallicity is that at lower Z , gas temperature is higher. At higher gas temperature, two-photon emission becomes stronger partially due to the departures from Case B (Raiter, Schaerer & Fosbury 2010). This effect is particularly noticeable when the metallicity of the nebula drops to 10 -3 Z ⊙ . Non-solar abundance patterns may modify where this effect matters, particularly if carbon or other strong coolants are removed from the gas.  \nThe largest impact on the ξ ion discrepancy arises with gas density (compare line colors). As already explained above, this is due to l -changing collisions being able to shift electrons between the 2 s and 2 p states which then preferentially decay as Ly α due to the large difference in Einstein coefficients. Nevertheless, in all cases, when log 10 ( ξ ion ) ≳ 25 . 5, the true intrinsic value must be higher that what is observed.  \nWe conclude from this experiment that under typical ISM conditions, if the gas in galaxies is being irradiated by a normal stellar population, log 10 ( ξ ion ) should not be observed to be ≳ 25 . 8. This is because the intrinsic log 10 ( ξ ion ) of even the optimistic SPS models peaks at ∼ 26 . 1 which corresponds to an observed log 10 ( ξ ion ) of 25 . 8. Measurements higher than this value would require any one of the following:  \n- · An incorrect dust correction. For dusty galaxies UV luminosity is suppressed more than the Balmer lines which would lead to an over-prediction of ξ ion in the case where dust acts as a screen.\n- · High gas densities such that the two-photon emission is suppressed and the true ξ ion can be measured.\n- · An exotic stellar population that intrinsically has a much higher ξ ion than that predicted by normal SPS models.\n- · Extreme nebular conditions (e.g. incredibly high temperatures) or geometries, or non-stellar ionizing sources.  \nWe discuss these options further below in the context of JWST galaxies.  \nFinally, we note that the bias in ξ ion measurements when the nebular continuum is strong is unlikely to impact inferences on reionization history. This is because in such calculations, ξ ion is multiplied by UV luminosity density (e.g. Robertson et al. 2013) and even though ξ ion is biased low, this is compensated by the fact that observed UV luminosity density is increased by the nebular continuum.', '2.3.4. UV turnovers masquerading as DLAs': 'A further implication of a strong nebular continuum is that as the ratio of the two-photon emission to the stellar continuum increases, a downturn in the UV emission just redward of Ly α begins to appear (see e.g. Cameron et al. 2023a). This downturn may be perceived to be a damped  \n<!-- image -->  \nFig. 7.(Left) Observed log 10 ( ξ ion ) as a function of intrinsic log 10 ( ξ ion ) for various SPS models for our fiducial photoionization model parameters. All models follow the same trend despite differences in the underlying spectra. (Right) Observed log 10 ( ξ ion ) as a function of intrinsic log 10 ( ξ ion ) for metal-poor, hot star models for various assumptions of gas density (colors), ionization parameters (line styles), and metallicities (line opacities). This demonstrates that the properties of the nebula dictate the observed ξ ion . Note that each intrinsic ξ ion value in the right panel corresponds to a single stellar temperature for the massive star models while in the left panel, the correspondence is with stellar age. In both panels, the diagonal dotted red line shows the one-to-one relation while the horizontal line represents log 10 ( ξ ion ) = 26 . 0.  \n<!-- image -->  \nLy α system (DLA). The deficit in UV luminosity near Ly α would raise its EW; however, at high-redshift, an increasingly neutral IGM may absorb emission close to rest-frame 1216 ˚ A. For this reason, high-redshift galaxies with significant two-photon emission may be perceived to have strong absorption from neutral hydrogen (i.e. a DLA).  \nTo demonstrate this, we fit the hot star SPS models in the wavelength interval (1250 ˚ A -2000 ˚ A) using a power law + DLA model. For this simple experiment, we assume that the galaxy is at z = 9, account for IGM absorption using the transmission curves from Garel et al. (2021), as well as the JWST/NIRSpec Prism/Clear line spread function 13 . In Figure 8 we show the inferred DLA column density as a function of stellar temperature. A clear trend emerges such that as stellar temperature increases, so does the inferred DLA column. At T ∼ 100 , 000 K, the inferred DLA column is ∼ 10 22 . 5 cm -2 , which is very high compared to known DLAs (e.g. Tanvir et al. 2019; Hu et al. 2023), and consistent with some of the extreme systems observed at high-redshift with JWST (Heintz et al. 2024c; Umeda et al. 2023).  \nWe emphasize two key points. First, the exact value of the inferred DLA column is highly subject to the assumed intrinsic SED. For our experiment, we have fit a power law to the continuum at wavelengths < 2000 ˚ A. However, the SED is not a perfect power law, noise can modify the inferred slope, and the wavelength baseline matters for the inferred UV slope. Flatter slopes lead to lower inferred column densities, while steeper slopes require more absorption. The exact shape of the underlying SED also contributes to the inferred DLA properties.  \nFig. 8.Strength of the UV downturn in terms of the inferred (masquerading) DLA column density as a function of stellar temperature. The UV downturn strength increases with stellar temperature as the two-photon emission becomes more visible.  \n<!-- image -->  \nHence Figure 8 should be taken more qualitatively than quantitatively (although in all cases the inferred DLA column will still be exceptionally high for high stellar temperatures). Second, the shape of a spectrum modulated by a DLA differs from the characteristic profile of two-photon emission. The DLA fits are thus not perfect representations of the true spectrum. This effect is likely diminished when noise is present and is subject to the redshift of the object as the spectral resolution of NIRSpec changes with wavelength. It should also be emphasized that strong nebular emission is not mutually exclusive with the presence of a DLA and both may occur simultaneously in the same galaxy. Nevertheless, a clear prediction of a nebula being irradiated by hot massive stars is a steep UV downturn that may masquerade as an extreme DLA.  \n<!-- image -->  \nFig. 9.(Left) 1500 ˚ A nebular fraction as a function of observed log 10 ( ξ ion ) for various SPS models for our fiducial photoionization model parameters. All models follow the same trend despite differences in the underlying spectra. (Right) Observed log 10 ( ξ ion ) as a function of intrinsic log 10 ( ξ ion ) for metal-poor, hot star models for various assumptions of gas density (colors), ionization parameters (line styles), and metallicities (line opacities).  \n<!-- image -->', '2.4. Inferring a high nebular contribution': 'As we have shown in Figure 2, for typical SPS models, the nebular continuum contribution does not exceed 40% unless more massive/hotter stars are included. Here we summarize the key characteristics of a galaxy that may have a high nebular contribution.  \n- 1. The presence of a Balmer jump. The only way to avoid a Balmer jump in systems with strong nebular continuum is to increase the temperature to extreme values, well beyond what are typically associated with H II regions. In this case, the strength of the jump becomes much smaller with respect to line emission and can be dominated by other continuum processes that may make detection/identification difficult.\n- 2. A UV slope redder than ∼ -2 . 7. The nebular continuum reddens the intrinsic stellar population and this reddening is dominated by the free-bound emission. In the extreme scenario where the continuum is dominated by two-photon emission, a lower limit of β ∼ -2 . 7 will be measured 14 . Assuming the galaxy is undergoing a starburst, UV slopes redder than ∼ -2 . 3 may require more extreme ionizing sources than those typically assumed in standard SPS models or dust attenuation.\n- 3. A high ξ ion as measured via a Balmer line. If ξ ion is measured to be high, ionizing photons are both being produced at a high rate and being efficiently absorbed by the gas (i.e. not leaking such that f esc ∼ 0). log 10 ( ξ ion ) values above ∼ 25 . 7 -25 . 8 may require more extreme ionizing sources than those typically assumed by standard SPS models.  \n14 Note that this value is slightly sensitive to how the wavelength range over which β is measured is set.  \n- 4. Strong UV downturns. The appearance of a UV downturn, which can easily be confused with a high column density DLA, scales very strongly with nebular contribution due to two-photon emission. When a UV downturn appears simultaneously with bright Ly α emission, this is an even cleaner signature of strong nebular emission because DLAs are optically thick to Ly α . The presence of Ly α is however not a necessary characteristic as both the IGM, dust, and geometry can absorb Ly α even if the nebular continuum is strong.\n- 5. High H α (or H β ) EWs. For instantaneous bursts of star formation, the H α EW can reach as high as ∼ 3 , 200 ˚ A across all of the different SPS models. Higher values require more extreme ionizing sources. We note that high H α EW is not required for extreme nebular emission as there are numerous scenarios where it can be lower than predicted for the instantaneous burst. For example, the presence of an older stellar population can impact the rest frame optical while having minimal contribution in the UV. Moreover, if the two-photon emission is particularly strong, the H α and H β EW would also be reduced compared to predictions from both nebular-only models that assume Case B recombination or SPS models where the temperature is set by the SPS model.  \nAny combination of the above spectral characteristics may indicate a galaxy is a candidate for having strong nebular continuum emission.  \nMore specifically, one can even quantify the exact nebular contribution if certain properties of the gas are known. Because nebular emission is only weakly sensitive to the shape of the input ionizing spectra and strongly sensitive to the properties of the nebular gas, if the underlying physical properties of the nebula are  \nknown, the nebular contribution at any wavelength 15 can be uniquely determined. Hence, all one needs to infer the nebular fraction is the observed ξ ion and a measure of gas density and temperature. We demonstrate this in the left panel of Figure 9 where we show the nebular contribution at 1500 ˚ A as a function of observed ξ ion (i.e. that which can be measured from the spectrum) for a metallicity of 1% Z ⊙ and a gas density of 100 cm -3 . The shape and normalization of the trend is relatively independent of the SPS model. The largest spread occurs between FSPS and Starburst99 where at an observed log 10 ( ξ ion ) of 25.71, the nebular contribution to the total spectrum differs by 6%.  \nIn the right panel of Figure 9 we show the 1500 ˚ A nebular fraction as a function of observed ξ ion for a selection of the metal-poor massive star models while varying gas density, ionization parameter, and metallicity. At typical ISM densities below 10 3 cm -3 the nebular fraction is only sensitive to metallicity (via gas temperature). Above this critical density, gas density becomes the dominant factor as discussed in the previous section. For all gas conditions, ionization parameter is a subdominant factor.  \nWe emphasize that there are important caveats to this analysis. First, the gas temperature in the photoionization models is determined via the gas metallicity and the underlying shape of the spectrum. Additional heating terms can modulate this relation such that lower observed ξ ion means a higher nebular fraction. A similar behavior can be achieved by assuming non-solar abundance patterns that decrease cooling (i.e. fixing the oxygen abundance but decreasing other metal species). This is because the two-photon emission is enhanced with respect to free-bound and Balmer emission in this scenario. Likewise, any other physics, e.g. collisional excitation, cooling flows, harder radiation, that increases two-photon emission without impacting the recombination lines will mean that f neb . is under-predicted for a given observed ξ ion .  \nThe photoionization models also assume that the nebula is ionization bounded. In the case of density bounded nebulae, the trends can differ. Finally, there are nuances to measuring ξ ion . A fixed conversion between H α (or another Balmer line) and ionizing luminosity is often assumed; however, this conversion depends on temperature. Measuring the Balmer line correctly is often difficult in JWST data where the spectrum can be very noisy in the continuum near H α . Finally, dust can significantly impact 1500 ˚ A luminosity and due to the large wavelength difference between the Balmer lines and 1500 ˚ A, an accurate dust correction is paramount.', '3. A SAMPLE OF HIGH-REDSHIFT JWST GALAXIES WITH STRONG NEBULAR CONTINUUM': "While the theory behind the nebular continuum is well understood, there are very few spectroscopically confirmed examples of high-redshift galaxies with a prominent nebular continuum contributions. The Balmer jump at rest-frame 3465 ˚ A is the most easily identifiable feature and this has only been reported for a select few high-redshift JWST galaxies (e.g. Cameron et al. 2023a; Welch et al. 2024), or in a stack of strong Ly α emit-  \nrs (Roberts-Borsani et al. 2024). The lack of reported Balmer jumps is perhaps surprising because SED fits to NIRCam photometry often predict strong nebular continuum emission (e.g. Endsley et al. 2023; Bradaˇc et al. 2024) for galaxies with extreme emission lines. For this reason, we have undertaken a survey of public JWST data to assemble a catalog of galaxies with strong nebular continuum emission.  \nMore specifically, we visually inspected all prism spectra in v2 of the DAWN JWST Archive (DJA, Heintz et al. 2024c) at z ≳ 2 . 5 for a possible spectral discontinuity at the location of the Balmer jump 16 . All spectra in this database were reduced using msaexp 17 and full details of the reduction can be found in Heintz et al. (2024c). No other sample selection criteria were applied except for the removal of galaxies with obvious artifacts in the 2D spectra that coincidentally overlapped with the location of the Balmer jump. This means no objects were removed for obvious broad features (that may indicate the presence of an AGN) or any other reason. Our broad selection criteria leads to an extreme diversity in the types of objects that are included. Figure 10 shows example spectra and RGB images of 20 galaxies in the data set. Our sample contains data from 14 different JWST programs as listed in Table 1, with the most coming from JADES GTO (IDs: 1180, 1181, 1210, 3215, PIs: Eisenstein, Leutzgendorf, Bunker et al. 2023a; D'Eugenio et al. 2024), RUBIES (ID: 4233, PI: de Graaff, de Graaff et al. 2024), and UNCOVER (ID: 2561, PI: Labb'e, Bezanson et al. 2022). We emphasize that because we adopt visual inspection, the sample is unlikely to be both complete and pure. The strength of spectral break that we are sensitive to is highly dependent on the signal-to-noise ratio of the continuum and we empirically find that most breaks in the sample have a > 20% reduction in flux between rest-frame 3500 ˚ A, and 4200 ˚ A.  \nThe final sample contains 58 objects in the redshift interval 2 . 5 ≲ z ≲ 9. Redshifts are calculated by fitting templates for emission lines to each prism spectrum or grating spectrum when available. The distribution of redshifts is shown in the top panel of Figure 11, with most of the objects falling in the redshift interval 4 ≲ z ≲ 6. The bias towards lower redshifts is likely a reflection of sample size and the ease of detecting the continuum at high enough signal-to-noise to identify the Balmer jump rather than strong nebular emission being more common at lower redshift. The data set spans more than four magnitudes in UV luminosity from M UV > -18 to M UV < -22 as shown in the bottom panel of Figure 11 18 .  \nFor each galaxy we measure various emission lines and continuum fluxes. We mask out strong emission lines and fit the remaining data with a seven parameter model using templates for the nebular continuum at various temperatures generated with PyNeb (Luridiana, Morisset & Shaw 2015) and for stellar spectra (Stanway & Eldridge 2018) including both young and old stellar populations. We use an MCMC sampler (Foreman-Mackey et al. 2013)  \n16 This possibly omits galaxies where the nebula is at such high temperatures that two-photon emission completely dominates over the Balmer jump. We leave such a search to future work.  \nFig. 10.A gallery of example Balmer jump galaxies. We show the observed spectrum of each galaxy as well as an RGB image combining JWST NIRCam filters F444W, F200W, and F150W in the red, green, and blue channels, respectively. Each image is 1' × 1'.  \n<!-- image -->  \n5000  \nFig. 11.Histograms of redshift (top), 1500 ˚ A UV magnitude (middle), and H α or [O III]+H β EW (bottom) for the 58 galaxies in our dataset. For comparison, we show the [O III]+H β EW distribution from the JWST primal database (Heintz et al. 2024b) of well-detected spectroscopically confirmed galaxies at z > 5 . 5. We have renormalized the JWST primal histogram by the ratio of Balmer jump galaxies to the Primal galaxies (i.e. 58/494).  \n<!-- image -->", 'TABLE 1': 'Table of program name, pi, program id, and the number of  \ngalaxies from each program that are part of our Balmer jump sample. Programs are sorted descending by the number of galaxies (and then alphabetically by program name).  \nto measure the posterior distribution on the contributions of each component the continuum as well as the gas temperature and ages of each stellar component. To measure emission line luminosities, we randomly sample the continuum posterior and subtract it from the observed spectrum. The residuals are then fit with Gaussians using specutils 19 , accounting for the error in the residuals. Line fluxes and their uncertainties are taken as the median and standard deviation of the Monte Carlo samples.  \nTypical rest-frame equivalent widths (EWs) for H α in our sample of galaxies is ≳ 1000 ˚ A. This is significantly higher than the 500 -800 ˚ A H α EWs found for the general population of high-redshift galaxies from Roberts-Borsani et al. (2024). A similar trend is seen for [O III]+H β EW where the typical galaxy in our sample has an EW a few times higher than that found in Roberts-Borsani et al. (2024) and similarly that estimated for photometric samples (e.g. Endsley et al. 2023). More quantitatively, comparing to the JWST Primal database (Heintz et al. 2024b), the median [O III]+H β EWis 2.3 × higher in our dataset compared to the typical high-redshift galaxy and the distributions are very different (see Figure 11). Some galaxies in our sample exhibit incredibly high EWs, approaching 5000 ˚ A in [O III]+H β .', '3.1. Star Formation Histories of Balmer Jump Galaxies': 'We have argued that the Balmer jump sample consists of galaxies with bursty star formation histories due to the visibility of the nebular continuum. To demonstrate this further, we have estimated the star formation history of each galaxy using SED fitting.  \nWe employ the Dense Basis code (Iyer & Gawiser 2017; Iyer et al. 2019), which reconstructs the star formation history of the galaxy using Gaussian Processes. The star  \nFig. 12.Normalized stacked star formation histories of the Balmer jump galaxies. The star formation history of each galaxy is normalized so that the SFR at the present time is 1. The black line and gray shaded region show the median star formation history and the 1 σ scatter.  \n<!-- image -->  \nformation history is split into five bins of equal stellar mass. We adopt a Jeffreys prior on the shape parameters. We first construct a sparse atlas using a wide range of final stellar mass, SFR, metallicity, and A V . Rather than fitting the spectrum directly (as we find that many SED fitting codes struggle to simultaneously fit the Balmer jump strength and the emission line properties), we measure photometry for all JWST NIRCam wide and medium bands from the spectrum and fit the photometry. This assumes that the spectra are properly flux calibrated as we do not rescale to NIRCam fluxes apart from what was applied at the time of reduction in the DJA, i.e. no additional slit loss corrections are applied. Uncertainties on the photometry are calculated by Monte Carlo resampling the observed spectrum from the error distribution on each pixel. Because the redshift is kept fixed in the fit to that measured from the spectrum, we consider only filters that are redward of (and do not include) Ly α . From the sparse atlas, we select the best fit model and create a refined atlas adopting priors that are 2 σ deviations around the best model. Our final best model comes from the refined atlas.  \nIn Figure 12, we show the normalized star formation history stack of all of the Balmer jump galaxies over the past 200 Myr. For each galaxy, we normalize the star formation history by the SFR at the current time. The histories are then interpolated onto the same time axis and we show the median and 1 σ deviations.  \nThe star formation histories appears to be nearly a delta function at the present time. The median history prefers almost no star formation beyond ∼ 35 Myr and the present day SFR of the median history is > 200 × that at 20 Myr. This exercise clearly demonstrates that the Balmer jump sample represents a unique class of galaxies undergoing extreme bursts of star formation, which potentially makes them analogues of the extreme highredshift galaxy population.', '3.2. The Characteristic Spectrum of Strong Nebular Emission': 'We begin by stacking the galaxies to increase the signal-to-noise. Each prism spectrum is de-redshifted to the rest-frame and we use the FluxConservingResampler from specutils to resample the spectra onto a common wavelength grid. Resampled spectra are then normalized to have a flux density of unity at 3200 ˚ A. Mean and median stacks are shown in Figure 13. Scatter in the median stack is calculated from 10,000 bootstrap resamples. Apart from exhibiting very strong emission lines, both in the rest-frame UV and optical, the most notable feature is the strong spectral discontinuity at 3645 ˚ A which demonstrates that these galaxies exhibit a very clear Balmer jump.  \nOnly minor differences are apparent between the mean and median stacks. Both exhibit remarkably little variation in UV flux as a function of wavelength between 1500 ˚ Aand 3000 ˚ A, indicating a spectral slope of β ∼ -2. The strength of the Balmer jump is also similar and both stacks exhibit a downturn in the UV, just redward of Ly α . Emission line fluxes tend to be slightly higher for the mean stack. This is particularly true for C IV λλ 1550 and Ly α , where we find that fewer than half of the Balmer jump galaxies are strong Ly α emitters.  \nThe stack exhibits strong rest-frame UV and optical emission lines, consistent with the individual EWs measured in our sample. There is also evidence for slightly rarer emission lines such as the combination of He II λ 4686 and the nearby [Ar IV ] λλ 4711,4740 doublet. The [O II ] quadruplet auroral line feature at rest-frame 7320-7330 ˚ A is clearly visible as was also recently seen in the lower-redshift stack of Strom et al. (2023). Numerous He I lines are present as well as weak [Ar III ] at 7136 ˚ A.', '3.3. The nebular contribution in Balmer jump galaxies': 'Now considering individual galaxy spectra, one of the key predictions of the SPS models is that the nebular contribution at 1500 ˚ A maximizes at ∼ 40% unless hotter/more massive stars are included. Having a sample of galaxies with bursty star formation histories and clear nebular continuum contributions in the optical, we can test whether this ∼ 40% upper limit holds for real galaxies.  \nInferring the nebular contribution is not straightforward unless the detailed properties of the nebula are known. Robust temperatures, densities, and ξ ion values are not available for the vast majority of galaxies in our sample and thus we resort to an alternative technique. In principle, the shape of the continuum encodes the contributions of stellar and nebular emission. However, the key confounding factor occurs in the UV where downturns can occur both due to neutral hydrogen absorption and two-photon emission. We thus fit the galaxies twice, once using our seven parameter MCMC fitting code and a second time allowing for a DLA in the fit. We then focus on sources where the nebular continuum is predicted to be strong independent of the model.', '3.3.1. Nebular dominated galaxies': "Based on this procedure, a few galaxies immediately appear to require an unusually strong nebular contribution. The first galaxy is 1210 13176 at z = 5 . 943 (also  \nFig. 14.Comparison of 1210 5217 at z = 4 . 890 with potential nebular-dominated galaxies 2189 7807 ( z = 5 . 387) and GS 9422 ( z = 5 . 946). All three show a clear Balmer jump and a very similar UV downturn. Spectral modeling predicts a very high nebular fraction in each of these objects.  \n<!-- image -->  \n-  \nFig. 13.Mean (cyan) or median (magenta) stack of all 58 galaxies in our sample. We highlight the spectral discontinuity at the location of the Balmer jump. The shaded region around the median stack shows the 5th-95th percentile distribution calculated from 10,000 bootstrap resamples. This demonstrates that no single galaxy has an extreme impact on the stack.  \n<!-- image -->  \n-  \nreferred to as NDG 9422 in Cameron et al. 2023a), the spectrum of which is shown as the blue line in Figure 14 and the spectral properties can be found in Table 2. Our template fitting indicates a 1 σ lower limit of 80% for the nebular contribution to the 1500 ˚ A luminosity. This galaxy is discussed at length in Cameron et al. (2023a) and the spectrum was first described in Saxena et al. (2024), so here we only highlight the key features.  \nWhat makes 1210 13176 interesting is that not only does it exhibit a very strong downturn in the UV, which if interpreted as a DLA would imply an extreme neutral column density of ∼ 10 23 cm -2 , it also has strong Ly α emission. If a such a high column density DLA is present, it remains unclear why such a large fraction ( > 25%, Saxena et al. 2024) of the Ly α escapes, especially considering that the galaxy is very compact with no visible companion that could otherwise account for the Ly α emission. Interestingly, Cameron et al. (2023a) demonstrated that the shape of the continuum is consistent with the expectations of the nebular continuum at a gas temperature similar to that measured from [O III ] λ 4363/[O III ] λ 5007. The log 10 ( ξ ion / erg -1 Hz) for this galaxy is 25.8, indicating that the nebular contribution in the UV could be very strong. While Cameron et al. (2023a) present numerous physical mechanisms to explain the shape of the spectrum, hot stars remain a natural explanation for the continuum shape of this galaxy.  \nMore recently, Schaerer et al. (2024) have argued, based on Case B recombination, that the H β equivalent width is too low for the spectrum to be fully nebular dominated. We argue that the discussion presented in Schaerer et al. (2024) does not necessarily apply to this object. It is clear from the Balmer decrement that the Case B assumption is not valid for this galaxy (e.g. Yanagisawa et al. 2024). Moreover, the temperature inferred from both the [O iii ] auroral line ratio and the Balmer jump is T e > 18 , 000 K, much higher than the typically assumed 10 4 K. At these temperatures, the expected Case B H β EW is 860 ˚ A (somewhat lower than the ∼ 1,000 ˚ A suggested in Schaerer et al. 2024, under the assumption of T e = 15 , 000 K).  \nMore importantly, as shown in Appendix A, the twophoton emission still contributes significantly in the restframe optical. Collisional excitation to the n = 2 state can significantly enhance the two-photon continuum, further lowering the observed H β EW. This has, in fact, been predicted for hot and metal-poor environments Raiter, Schaerer & Fosbury (2010); Mas-Ribas, Dijkstra &Forero-Romero (2016). The best-fit model of Cameron et al. (2023a) predicts such an enhancement which naturally explains both the observed H β and H α EW. Finally this system shows clear oxygen, carbon, and neon emission, which indicates that there would have been a previous generation of stars. This would further lower the Balmer EWs without impacting the fit in the UV. Thus, as presented in Cameron et al. (2023a), we stress that the spectrum of 1210 13176 is fully consistent with having a dominant contribution from the nebular continuum in the UV.  \nOther studies have suggested alternative explanations for the observed UV downturn and concomitant Lyα detection. Based on the fact that several other galaxies in  \nthe same field as 1210 13176 have DLAs, Heintz et al. (2024a); Terp et al. (2024) argue that a lower-redshift DLA at z ∼ 5 . 4 in a foreground galaxy cluster would offset the Ly α absorption trough, which could allow Ly α to leak and reconcile the differences in shape between the observed UV downturn and a DLA. However, this solution would imply an even larger gas column density. Likewise, Tacchella et al. (2024); Li et al. (2024) argue that an obscured black hole with an offset narrow line region and a large DLA column could provide a similar spectral shape. Because of the UV downturn and strong Ly α , any solution to this puzzle that invokes a DLA requires a 'fine-tuned' geometry. Therefore we argue that, for a single object these could be reasonable explanations; however, if more galaxies like 1210 13176 are identified, such models become less favorable.  \nInterestingly, a second galaxy in our sample, 2198 7807 at z = 5 . 387 from the Cycle 1 program: 'Quiescent or Dusty? Unveiling the Nature of Red Galaxies at z > 3' (Program ID: 2198; PIs: L. Barrufet & P. Oesch; Barrufet et al. 2024) shows a very similar UV downturn and strong Ly α emission as 1210 13176 (see the orange line in Figure 14. This galaxy is remarkably similar to 1210 13176 in UV magnitude, β slope, and ξ ion , and has similarly high emission line EWs (see Table 2). There is no indication of dust in the object and our 1 σ lower-limit on the 1500 ˚ A nebular fraction of 60%, much higher than that of typical SPS models (see Figure 2). As these two spectra are so similar, we refrain from a more detailed discussion since the same arguments presented for 1210 13176 in Cameron et al. (2023a) apply to 2198 7807. However, we briefly highlight a few key differences. The C IV λλ 1550 emission is much weaker in 2198 7807 and there is no evidence for He II λ 1640, which suggests the emission is even less likely to be powered by a black hole 20 . [O II ] λλ 3727 is also stronger in 2198 7807 compared to 1210 13176, suggesting it is less likely that the emitting H II regions are density bounded. Because 2198 7807 is at a redshift consistent with or lower than the presumed proto-cluster that was used in Heintz et al. (2024a); Terp et al. (2024) to explain both the UV downturn and the Ly α emission in 1210 13176, the same argument cannot apply to 2198 7807, unless there is another foreground cluster at even lower redshift. This would be a remarkable coincidence given that the shapes of the UV spectra are so similar. Finally, fine-tuned models that invoke strong DLAs (e.g. Tacchella et al. 2024; Li et al. 2024) are less appealing given that there is more than one spectra with an almost identical UV shape.  \nA third galaxy that is predicted to have a high 1500 ˚ A nebular fraction is 1210 5217 at z = 4 . 89 (green line in Figure 14). This galaxy was presented in Boyett et al. (2024) where they noted an extremely high log 10 ( ξ ion ) value of 25 . 9 ± 0 . 1 and no dust. Such a value is much higher than expected for standard SPS models when the nebular continuum is taken into account (see above). This value is also consistent with that measured in our work (26 . 04 +0 . 03 -0 . 03 ) to ∼ 1 σ . This UV luminosity of this galaxy is much fainter by two magnitudes compared to the previous two discussed in this section, but it has a  \n20 Although we emphasize that for 1210 13176, nebular diagnostics suggest that this object is not AGN dominated or at best a composite (Cameron et al. 2023a).  \nProperties of five galaxies predicted to have nebular contributions far exceeding that predicted from SPS models. We list the galaxy ID, redshift, 1500 ˚ A UV magnitude, UV slope, Balmer decrement, EWs of H α and [O III]+H β , H II region temperature predicted from the MCMC fitting of the spectra, ionizing photon production efficiency, 1 σ lower limit on the estimated nebular contribution at 1500 ˚ A, and finally the DLA gas column density or 1 σ upper limit. For each galaxy we list two values of ξ ion computed using either H α (top) or H β (bottom). The first three galaxies may host a population of massive stars while the bottom two spectra of different regions of a gravitationally lensed arc appear to be consistent with purely Case B recombination.  \nTABLE 2  \nvery similar UV slope of -1 . 99 +0 . 07 -0 . 09 . As this is much redder than that expected from SPS models for a single burst and given the lack of dust in this object, a high UV nebular fraction is a likely explanation. Moreover, this galaxy has extreme emission line equivalent widths. The equivalent widths of H α and [O III ] + H β are both ≫ 3 , 000 ˚ A, the highest H α EWamong the entire Balmer jump sample. No He II λ 1640 or λ 4686 is detected in the prism or gratings, despite a clear grating detection of O III λ 1666, and strong line ratios indicate that this source is not a candidate AGN (Scholtz et al. 2023; Curti et al. 2024).  \nWhile 1210 5217 is not a Ly α emitter, it also exhibits a strong UV downturn that has a remarkably similar shape to 1210 13176 and 2198 7807 (Figure 14). Even allowing for a DLA in our MCMC fit, the shape of the UV downturn prefers a nebular solution and hence we find a 1 σ lower limit of 63% for the nebular fraction at 1500 ˚ A. Like the other two galaxies, the gas temperature is predicted to be much higher than local H II regions, which could enhance the nebular contribution in the UV. We emphasize that a lack of Ly α emission does not preclude 1210 5217 from being nebular dominated. The IGM is not completely transparent at these redshifts (e.g. Inoue et al. 2014) and the actual emission profile of Ly α is highly viewing-angle dependent, even for galaxies without DLAs (e.g. Blaizot et al. 2023). The presence of massive stars is a natural explanation for the UV downturn and the relatively red β slope given the high ξ ion . If such stars are common at high redshift, we may expect more to look like 1210 5217 than either 1210 13176 or 2198 7807 due to the sight-line variability of Ly α and the fainter UV magnitude.", '3.3.2. No young stars in the shutter?': "One of the key features of the three galaxies in the previous section is that their Balmer decrements indicate that they are dust-free but their UV slopes are redder than that of typical SPS models. Redder slopes can be achieved if the intrinsic ionizing photon production efficiency is higher than the SPS models which enhances the nebular contribution at all wavelengths. Such is the case when hot/massive stars are present (see Section 2.2.3 and  \n-  \nalso Schaerer et al. 2024). However, even in the case of hot/massive stars, if the gas density is ≲ 10 4 cm -3 , the UV slope for a young stellar population cannot be redder than β ∼ -1 . 7. If a Balmer jump is present in the spectrum and the Balmer decrement indicates that there is no dust in the system, a β redder than ∼ -1 . 7 would either imply that the gas is either exceptionally dense, suppressing the two-photon continuum, or there is no underlying young stellar continuum and we are observing pure nebular emission in the UV. As we show in Appendix A, in the latter scenario, β from pure Case B recombination can approach values of ∼ -1 . 3.  \nThe latter scenario where there is little or no young stellar population can occur if the NIRSpec slit only covers the nebula. We find two example spectra where this may be occurring. In Figure 15 we show an image of a gravitationally lensed arc (The Cosmic Gummy Worm) being lensed by the Abell 2744 cluster. Two separate JWST programs (2561, PI: Labb'e and 2756, PI: Egami) placed shutters on different parts of the arc and the spectra are shown in the bottom two panels of Figure 15. This source has already been discussed extensively in the literature (Vanzella et al. 2022; Lin et al. 2023) and is thought to host proto-globular clusters. Most notably, Lin et al. (2023) provide a map of UV slope across the arc and find values ranging from -2 . 5 to -1 . 5 with the bluest values centered on the location of where there are thought to be young stars. Both previous analyses conclude that this system likely hosts a young, metal-poor stellar population with relatively little dust.  \nThe key features of these spectra (see Table 2) are that their Balmer decrements are consistent with Case B values indicating that little or no dust is present in the object and the log 10 ( ξ ion ) values are very high at 25 . 97 or 26 . 03. The UV slopes are very red with values of -1 . 58 and -1 . 33, which is much redder than can be generated by even the hot star models, unless extreme densities are assumed. Both spectra exhibit exceptionally high [O III ]+H β EWs > 3 , 000 ˚ A.  \nWhat is unique about our MCMC fits to these spectra is that not only do they require a very high nebular fraction at 1,500 ˚ A that is significantly higher than  \nwhat can be produced with typical SPS models, especially 2756 301 where the model prefers a fully nebular solution 21 , but the model also simultaneously infers a very strong DLA column between 10 22 . 50 -10 22 . 75 cm -2 . These two spectra are the only ones in our sample that exhibit such behavior. The continuum shapes of these spectra are consistent with being nebular continuum from Case B recombination with a DLA.  \nIn fact, this system is already known to host a DLA based on shorter wavelength data from MUSE that covers Ly β (Lin et al. 2023), with those authors inferring a DLA column of N H ∼ 10 21 . 8 cm -2 . The DLA column likely varies across the system which may explain why our measurement, which only covers part of the arc, finds a higher column density.  \nAs speculated by Lin et al. (2023), our results support a scenario in which the red portions of the arc are dominated by the nebular continuum. Spectra such as 2561 17467 and 2756 301 are extreme examples where gravitational lensing can help distinguish the gas from the underlying continuum source and hence the spectra may appear as being nebular dominated. More such examples exist in the DJA; however, none are as exemplary as those shown here. The high measured ξ ion values are thus not necessarily reflective of the underlying ionizing source population, but rather a manifestation of separating a continuum source from the nebula.  \nIt should be emphasized that this system is not a Ly α emitter. This is a clear example of a galaxy with a very high, spatially varying DLA column density, and yet the spectrum looks completely different from the objects above with strong UV downturns and Ly α emission. If the ISM of this object has any comparison to 1210 13176 and 2198 7807, it makes it even more peculiar that if a high column density DLA was the explanation for their spectrum, Ly α emission can leak significantly from those objects but not for 2561 17467 and 2756 301.", '4. IMPLICATIONS FOR REDSHIFT > 10 GALAXIES': "Both the theoretical and observational exercises presented here provide insight into the expected properties and physical nature of galaxies that are undergoing extreme bursts of star formation. Nearly all numerical simulations that attempt to resolve the ISM predict that star formation becomes increasingly bursty at high-redshift (e.g. Ma et al. 2018; Katz et al. 2023; Pallottini et al. 2022). Moreover, the leading explanation for the highredshift bright galaxy problem is that surveys are sampling the extreme 'tip-of-the-iceberg' of bursty star formation (e.g. Ren, Trenti & Mason 2019; Mason, Trenti & Treu 2023; Shen et al. 2023; Sun et al. 2023; Kravtsov & Belokurov 2024). For both of these reasons, the nebular continuum is thus expected to become increasingly important at z > 10.  \nIn this Section, we discuss some of the physical properties of known z > 10 galaxies in the context of the models we have presented for the nebular continuum which may help elucidate their physical origin. We emphasize that in order to explain the high-redshift bright galaxy problem, only a subset of the known objects would require  \n21 Note that the fit for 2756 301 requires a temperature extrapolation for the nebular continuum and there is a systematic uncertainty, particularly related to temperature associated with this fit.  \n<!-- image -->  \nFig. 15.(Top) RGB image of the Cosmic Gummy Worm combining JWST NIRCam filters F444W, F200W, and F150W in the red, green, and blue channels, respectively. Overlaid in magenta are the locations of the MSA shutters. (Bottom) Prism spectra of two regions of the gravitationally lensed arc at z = 3 . 99, both showing red continua and a Balmer jump.  \n<!-- image -->  \ndeviations from standard galaxy formation physics. The expectation is that many objects will be in agreement with standard SPS models of bursty galaxy formation, while others may require a different interpretation.", '4.1. Do the bright high-redshift galaxies have strong nebular continuum emission?': 'Except for extreme examples like GNz11 where the continuum is detected and a Balmer jump detection is not conclusive 22 , the signal-to-noise of the spectra of known z > 10 galaxies at rest-frame 3,645 ˚ A is too low to see a Balmer jump. Stacking all publicly available spectra in the DJA with z > 10 also leads to a too low signal to detect a Balmer jump. Thus we discuss these high-redshift galaxies using more indirect indicators of strong nebular emission.  \nUnfortunately, very few properties can be measured for most z > 10 galaxies apart from their UV magnitudes, continuum slopes, and an occasional emission line. Thus in Figure 16, we plot β as a function of M UV for all galaxies in the JWST Primal sample (Heintz et al. 2024b). The region shaded in cyan show the expected  \n22 In one pointing, a spectral discontinuity appears at the location of the Balmer jump. Since not all observations show this, the presence of a Balmer jump is inconclusive.  \nFig. 16.UV continuum slope, β , as a function of M UV for all z > 10 galaxies in the JWST Primal database (Heintz et al. 2024b) are shown as black points. The region shaded in cyan show the expected UV slopes from all of the different SPS models that we test. Slopes bluer than this (shaded in black) would require a non-zero escape fraction, while redder slopes (magenta), would indicate the presence of either hotter stars, the combination of hotter stars and high gas density, or dust. Red points represent the three galaxies in our Balmer jump sample that appear to have a very strong nebular contribution at 1500 ˚ A.  \n<!-- image -->  \nUV slopes from all of the different SPS models that we test. Slopes bluer than this (shaded in black) would require a non-zero escape fraction, while redder slopes (magenta), would indicate the presence of either hotter stars, the combination of hotter stars and high gas density, or dust.  \nMost of the z > 10 galaxies fall close to the cyan region indicating that vanilla nebular continuum presents a reasonable explanation. Recall that this does not rule out the presence of a flatter upper-mass slope, as long as the upper mass does not increase much above ∼ 100 -300 M ⊙ . Three galaxies have β < -2 . 55 which would imply a weaker nebular continuum contribution. This could be caused by a non-zero escape fraction, or if the massive stars output fewer ionizing photons than expected compared to the SPS models that we consider. Finally, 3/17 galaxies have β > -2 . 2, although within the uncertainties, two are consistent with the cyan region. The outlier is the z = 11 . 5 23 galaxy (CEERS2 588) that has β = -1 . 79 ± 0 . 13, more than 3 σ away from the vanilla prediction.  \nPerhaps the simplest explanation for this galaxy is that it is dusty. However, in the context of theoretical models, this explanation is unsatisfying because most pre-JWST predictions of the high-redshift bright galaxy number counts from full-box simulations catastrophically fail in reproducing observations (e.g. Finkelstein et al. 2023) except for those that adopt a top-heavy IMF (Cowley et al.  \n23 Note that this is the same galaxy as reported in Arrabal Haro et al. (2023) with MSA ID 10 (CEERS2 588). The Lyman break and possible detection of [O II] λλ 3727 provide different redshifts. Following the JWST primal database, we adopt the redshift from the break.  \n2018; Lu et al. 2024) or those with enhanced efficiency of star formation (e.g. McCaffrey et al. 2023; Katz et al. 2023). Since dust lowers the observed UV luminosity, if obscuration was the explanation for this galaxy, the intrinsic M UV would be much greater, exacerbating the tension between theoretical models and observations. In contrast, appealing to hotter stars with T ≳ 80 , 000 K and high gas densities of ∼ 10 4 cm -3 (or even cooler stars at higher gas density) would also naturally explain this object. Moreover, the inclusion of hotter stars would also help reconcile why the object is so bright.  \nHowever, as discussed previously, the hot/massive star solution implies other behavior. First, we might expect a small UV downturn. At the stellar temperatures considered here, the inferred DLA column density would be ≳ 10 21 . 8 cm -2 . Second, depending on the exact shape of the ionizing spectrum, one might expect high EW UV lines. The metal lines can easily be weakened by appealing to low metallicities; however, the key uncertainty is the strength of He II λ 1640. However, the regions of the spectrum at wavelengths below the He + ionizing threshold is arguably highly uncertain (e.g. Kewley, Nicholls & Sutherland 2019). For this reason, it is unclear how much flexibility there is in the He II λ 1640 predictions. Nevertheless, CEERS2 588 remains a very viable candidate for hosting hot, metal-poor stars at a high gas density. This same explanation can also be applied to the recently-discovered z > 14 candidate that has a UV slope of ∼ -2 . 2 ± 0 . 7 and very low EW UV lines (Carniani et al. 2024) and also GS-z12-0 from Curtis-Lake et al. (2023) at z = 12 . 6 with a measured β = -1 . 84 ± 0 . 19.  \nIn summary, all but two z > 10 galaxies are consistent with having strong nebular continuum emission. Most galaxies do not appear to require hot massive stars in very high density gas, but with limited information, we cannot rule out massive stars with normal ISM densities, a top-heavy IMF, or contributions from an extended star formation history.', '4.2. UV downturns and β slopes at z > 10 ?': "Akey signature of the presence of a strong nebular continuum is the appearance of downturns in the spectrum in the UV, just redward of Ly α , or more generally deviations from a power-law slope produced by the stellar continuum. Recently, Heintz et al. (2024b) defined the 'Ly α damping parameter' that encapsulates UV downturns when the value is strongly positive. While Ly α emission can impact this metric, at z ≳ 10, the IGM is likely to absorb the vast majority of Ly α emission and thus, the difference between the measured damping that is expected from the IGM provides a clear probe of the spectral deviation from a pure power-law.  \nHeintz et al. (2024b) show that the fraction of galaxies with perceived absorption that would imply a neutral column density in excess of 10 21 cm -2 mildly increases from ∼ 60% at z ∼ 6 to 65% -90% at z > 8. Moreover, Umeda et al. (2023) show that all of their z > 10 galaxies have UV downturns corresponding to absorption by column densities > 10 22 cm -2 . Hence, there is tentative evidence that UV downturns, in particular strong ones, are more common at high-redshift.  \nIn Figure 17, we plot the Ly α damping parameter ( D Ly α ) as a function of UV slope for the various SPS models (left) and the individual metal-poor star models  \n<!-- image -->  \nFig. 17.Ly α damping parameter, D Ly α , as a function of UV slope compared to various SPS models (left) and for hot metal-poor star models at various gas densities (right). Black data points with error bars represent individual z > 10 from the JWST Primal database (Heintz et al. 2024b). Red points represent the three galaxies in our Balmer jump sample that appear to have a very strong nebular contribution at 1500 ˚ A.  \n<!-- image -->  \nat various ISM densities in the range 10 2 -10 5 cm -3 (right). Note that for this calculation, we have recalculated β to match the definition used in Heintz et al. (2024b), although this results in only a marginal change. We also apply IGM absorption assuming the highest redshift curve from Garel et al. (2021).  \nWhile many of the objects are consistent with the SPS models, two z > 10 galaxies seem to be outliers for having high D Ly α ∼ 90 ˚ A and a relatively flat UV slope, slightly bluer than -2 . 0, consistent with the galaxies in our Balmer jump sample that seem to exhibit large nebular contribution at 1500 ˚ A 24 . The standard SPS models, even with IGM absorption, cannot reproduce these galaxies and hence very dense DLAs would be required as was the interpretation in Heintz et al. (2024b); Umeda et al. (2023). However, the right panel of Figure 17 shows that these galaxies are very consistent with the hot star models and a low-density ISM. Since our CLOUDY calculations were stopped when the electron fraction reached 1%, no additional absorption is required as the perceived damping is entirely due to H I two-photon emission. These galaxies are a z = 12 . 5 galaxy from JADES (Curtis-Lake et al. 2023) and a z = 11 . 4 galaxy from CEERS (Arrabal Haro et al. 2023).  \nAnother cluster of high-redshift galaxies appears to have much lower D Ly α and redder β than predicted by the SPS models. As we show in the right panel of Figure 17, these can be reconciled with increasing gas density up to 10 5 cm -3 (although not necessarily hot massive stars). In this case, two-photon emission is suppressed which reduces the amount of perceived absorption. This can also be reconciled with Ly α emission, though this is rare at these redshifts (although c.f. Bunker et al. 2023b).  \nFinally, there is one object (1210 14220) that has  \nD Ly α ∼ 80 ˚ A and β ∼ -3 which cannot be explained by any of our models. This would imply a dust-free DLA as to not redden the β , but the DLA must be far enough from the continuum source to keep the nebular emission low enough as to not redden the UV slope. The physical nature of this particular galaxy thus remains uncertain.", '4.3. Considerations for Finding Population III Stars at High-Redshift': 'The stellar spectra of hot massive stars are extremely blue (e.g. Schaerer 2003) and without a nebular contribution, we expect Pop. III galaxies to have β < -3 . 0. However, the nebular emission is clearly non-negligible and the search for Pop. III galaxies has primarily consisted of looking for blue ( β ≲ -2) galaxies that show strong He II λ 1640 emission (e.g. Trussler et al. 2023). However, as discussed in Cameron et al. (2023a) and shown again here (see also Trussler et al. 2023; Zackrisson, Inoue & Jensen 2013; Inoue 2011), the two-photon emission is necessarily dominant if the stellar population consists of stars with T ≳ 100 , 000 K. Hence if strong He II emission is found but no UV downturn is present (e.g. Saxena et al. 2020), either the ISM of the galaxy is at very high density (which suppresses the two-photon emission), or the galaxy is not hosting hot massive stars. The latter does not rule out the galaxy from being Pop. III; however, this would imply that not all Pop. III stars are extremely hot and might populate a lower range of masses, in agreement with some numerical simulations (e.g. Stacy, Bromm & Lee 2016). In the case of highdensities, the strength of the UV downturn is correlated with the UV slope such that weaker downturns imply redder spectra. We re-emphasize here that the ISM density surrounding Pop. III stars is unknown. ISM densities surrounding massive stars are typically regulated by a combination of photoheating, radiation pressure, and stellar winds. However in the Pop. III scenario, winds are much weaker due to the lower opacities of the stellar at-  \ne (e.g. Vink, de Koter & Lamers 2001), radiation pressure is caused by UV absorption and Ly α scattering rather than IR photons scattering on dust (e.g. Kimm et al. 2018), and finally, the ISM is at much higher pressure due to inefficient cooling (e.g. Omukai et al. 2005), which lowers the pressure gradient between the H II region and the PDR. For these reasons, one may reasonably expect that the gas densities around Pop. III stars may be considerably higher than in the local ISM. We therefore stress that selection criteria for high-redshift Pop. III galaxies should not rule out systems with red β slopes approaching values of -1.', '5. CONCLUSIONS': "Understanding the physical nature of galaxies at cosmic dawn remains a primary goal of modern cosmology and galaxy formation. Arguably one of the biggest surprises of early high-redshift JWST observations is the number density of both photometrically selected and spectroscopically confirmed objects at z > 10 as their number counts are significantly higher that expected from many pre-JWST models (e.g. Finkelstein et al. 2023; Harikane et al. 2024; Leung et al. 2023; Chemerynska et al. 2023). One of the leading explanations for this anomaly is UV variability driven by extremely bursty star formation histories (e.g. Ren, Trenti & Mason 2019; Mason, Trenti & Treu 2023; Shen et al. 2023; Sun et al. 2023; Kravtsov & Belokurov 2024). While most numerical simulations with a resolved ISM predict that star formation histories at early epochs are much burstier than in the low-redshift Universe (e.g. Ma et al. 2018; Katz et al. 2023), there remains debate on whether this is enough to reconcile theory with observations (Pallottini & Ferrara 2023; Gelli, Mason & Hayward 2024). Hence a detailed study of galaxies with extremely bursty star formation histories is warranted.  \nOne method for probing galaxies with extremely bursty star formation histories is strong nebular emission, and specifically, the visibility of the nebular continuum. More generally, the nebular continuum can significantly modify the observed spectrum of a galaxy, by reddening its spectral slope (e.g. Bouwens et al. 2010; Dunlop et al. 2013; Cullen et al. 2024; Topping et al. 2024a), contributing to the 1500 ˚ A luminosity of the galaxy, and causing a downturn in the UV spectrum, just redward of Ly α (e.g. Fosbury et al. 2003; Raiter, Schaerer & Fosbury 2010; Cameron et al. 2023a). For this reason, we began with a theoretical exercise, adopting five different stellar population synthesis models to quantify how changes in the underlying stellar population synthesis model, IMF slope, stellar temperature, and nebular conditions impact the observed spectra of galaxies under the assumption of an extreme burst of star formation. The conclusions from our theoretical exercise can be summarized as follows:  \n- 1. 1500 ˚ A UV luminosities can increase by up to ∼ 0 . 7 magnitudes in the presence of extreme bursts of star formation for published stellar population synthesis models at young stellar ages ( ≲ 5 Myr). This result holds even when considering realistic star formation histories from state-of-the-art numerical simulations. If common at high-redshift, the nebular continuum contribution at 1500 ˚ A has the potential to significantly reduce the amount of\n- 'burstiness' required to reproduce the measured abundance of bright high-redshift galaxies.\n- 2. Assuming a top-heavy IMF only marginally increases the nebular contribution to the 1500 ˚ A UV luminosity and only slightly reddens the SED as long as the upper-mass limit of the IMF remains fixed. In contrast, assuming hotter (more massive) stars can redden UV slopes to values of β > -2 . 0. When the gas density is allowed to increase above the critical density of H I two-photon emission, β slopes around hot, massive stars can approach ∼ -1. Hence, searches for Population III galaxies should not rule out such red objects.\n- 3. When the nebular continuum is strong, twophoton emission can masquerade as DLA absorption with apparent column densities reaching up to ∼ 10 23 cm -2 . When UV downturns are observed at high-redshift, the DLA solution should be compared with predictions from strong nebular continuum emission.\n- 4. Because ξ ion measurements adopt the observed 1500 ˚ A UV luminosity, observationally inferred ξ ion values deviate from intrinsic ξ ion at values > 10 25 . 5 erg -1 Hz. Assuming normal nebular conditions, published SPS models predict that measured values of ξ ion should never be ≳ 10 25 . 81 erg -1 Hz. When the effects of dust can be rules out, values above this limit would imply either 1) non-standard ISM conditions, 2) an intrinsic ξ ion ≳ 10 26 . 2 erg -1 Hz which would point to objects like hot, massive stars or an AGN, or 3) that the spectrum does not fully cover both the continuum source and the nebular emission.  \nWe then proceeded to study bursty star formation observationally by compiling a sample of 58 galaxies from the Dawn JWST Archive (DJA, Heintz et al. 2024c) in the redshift interval 2 . 5 < z < 9 that have clear Balmer jumps in their R ∼ 100 JWST NIRSpec prism spectra. This sample increases the number of reported highredshift JWST Balmer jumps by a factor of ∼ 30 × , confirming some expectations from SED fitting (e.g. Endsley et al. 2023). The conclusions from our observational exercise can be summarized as follows:  \n- 1. SED fitting using non-parametric star formation histories indicates that nearly all sources in the sample are well described by an extreme recent burst of star formation, consistent with expectations for when the Balmer jump becomes visible.\n- 2. Five spectra in the sample are consistent with having a nebular dominated spectrum - significantly more nebular emission is preferred in the continuum fit to the observed spectrum compared to expectations from SPS models. Two of these spectra come from the same strongly gravitationally lensed galaxy at z ∼ 4 and have red UV slopes, but have Balmer decrements consistent with Case B recombination expectations. We argue that in these spectra, the MSA slit preferentially covers the nebular gas and not the underlying ionizing source due to  \nthe high effective spatial resolution enabled by lensing. In contrast, three other spectra have β ∼ -2, two of which have Ly α emission in addition to a Balmer jump, and the third has an H α EW of ∼ 3 , 300 ˚ A. All three have a remarkably similar shape for their UV downturn that is consistent with being two-photon emission from H I . The leading explanation is that these galaxies host hot, massive stars as the similarities across multiple galaxies seem to disfavor fine-tuned geometric scenarios that require a DLA.  \nFinally, we discussed some of the recently detected z > 10 galaxies in the context of our nebular continuum models. While many galaxies are consistent with expectations for ordinary bursty star formation combined with nebular emission, some seem to deviate by either having excessive downturns in the UV and redder UV slopes. Some of these galaxies are consistent with hot massive star models, which may help explain some of the highredshift bright galaxy problem.  \nBecause the nebular continuum is sensitive to both ISM density and temperature, a sample such as this provides a unique opportunity to constrain the properties of the high-redshift ISM using emission that has different sensitivities to temperature and density compared to existing strong-line methods. Ideally one could constrain the abundance of Balmer jump galaxies at high-redshift;  \nhowever, the visual selection and ill-constrained selection functions of the various observational programs used in this work currently inhibits accurate estimates for this quantity. Topics such as this will be explored in future work.", 'ACKNOWLEDGMENTS': "HK thanks Andrey Kravtsov for insightful comments and thoughtful discussions. We sincerely thank the PIs and Co-Is of the JWST programs where spectral data was made publicly available on the DJA. We refer interested readers to the following papers for survey descriptions regarding the spectral data: Bunker et al. (2023a); D'Eugenio et al. (2024); Bezanson et al. (2022); Barrufet et al. (2024); de Graaff et al. (2024); Finkelstein et al. (2024); Glazebrook et al. (2024); Pierel et al. (2024); Siebert et al. (2024); Maseda et al. (2024). 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. These observations are associated with programs listed in Table 1. AJC and AS acknowledge funding from the 'FirstGalaxies' Advanced Grant from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 789056).", 'REFERENCES': "- Abel T., Bryan G. L., Norman M. L., 2002, Science, 295, 93\n- Arrabal Haro P. et al., 2023, Nature, 622, 707 Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481 Barrufet L. et al., 2024, arXiv e-prints, arXiv:2404.08052 Bate M. R., 2023, MNRAS, 519, 688 Bezanson R. et al., 2022, arXiv e-prints, arXiv:2212.04026 Blaizot J. et al., 2023, MNRAS, 523, 3749 Bouwens R. J. et al., 2010, ApJ, 708, L69 Bouwens R. J., Smit R., Labb'e I., Franx M., Caruana J., Oesch P., Stefanon M., Rasappu N., 2016, ApJ, 831, 176 Boyett K. et al., 2024, arXiv e-prints, arXiv:2401.16934 Bradaˇc M. et al., 2024, ApJ, 961, L21 Bromm V., Coppi P. S., Larson R. B., 2002, ApJ, 564, 23 Bunker A. J. et al., 2023a, arXiv e-prints, arXiv:2306.02467 Bunker A. J. et al., 2023b, A&A, 677, A88 Cameron A. J., Katz H., Rey M. P., 2023, MNRAS, 522, L89 Cameron A. J., Katz H., Witten C., Saxena A., Laporte N., Bunker A. J., 2023a, arXiv e-prints, arXiv:2311.02051 Cameron A. J. et al., 2023b, A&A, 677, A115 Carniani S. et al., 2024, arXiv e-prints, arXiv:2405.18485 Castellano M. et al., 2024, arXiv e-prints, arXiv:2403.10238 Chatzikos M. et al., 2023, Rev. Mexicana Astron. Astrofis., 59, 327 Chemerynska I. et al., 2023, arXiv e-prints, arXiv:2312.05030 Chisholm J. et al., 2022, MNRAS, 517, 5104 Chon S., Ono H., Omukai K., Schneider R., 2022, MNRAS, 514, 4639 Choustikov N. et al., 2024, MNRAS, 529, 3751 Conroy C., Gunn J. E., 2010, ApJ, 712, 833 Conroy C., Gunn J. E., White M., 2009, ApJ, 699, 486 Cowley W. I., Baugh C. M., Cole S., Frenk C. S., Lacey C. G., 2018, MNRAS, 474, 2352 Cueto E. R., Hutter A., Dayal P., Gottlober S., Heintz K. E., Mason C., Trebitsch M., Yepes G., 2024, A&A, 686, A138 Cullen F. et al., 2024, MNRAS, 531, 997 Cullen F. et al., 2021, MNRAS, 505, 903 Curti M. et al., 2024, A&A, 684, A75\n- Curtis-Lake E. et al., 2023, Nature Astronomy, 7, 622\n- Dayal P., Ferrara A., Saro A., 2010, MNRAS, 402, 1449 de Graaff A. et al., 2024, arXiv e-prints, arXiv:2404.05683 Dekel A., Sarkar K. C., Birnboim Y., Mandelker N., Li Z., 2023, MNRAS, 523, 3201 D'Eugenio F. et al., 2024, arXiv e-prints, arXiv:2404.06531 Dijkstra M., 2009, ApJ, 690, 82 Dunlop J. S. et al., 2013, MNRAS, 432, 3520 Endsley R. et al., 2023, arXiv e-prints, arXiv:2306.05295 Fardal M. A., Katz N., Gardner J. P., Hernquist L., Weinberg D. H., Dav'e R., 2001, ApJ, 562, 605 Faucher-Gigu'ere C.-A., Kereˇs D., Dijkstra M., Hernquist L., Zaldarriaga M., 2010, ApJ, 725, 633 Ferrara A., Pallottini A., Dayal P., 2023, MNRAS, 522, 3986 Finkelstein S. L. et al., 2023, arXiv e-prints, arXiv:2311.04279 Finkelstein S. L. et al., 2024, ApJ, 969, L2 Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306 Fosbury R. A. E. et al., 2003, ApJ, 596, 797 Gardner J. P. et al., 2006, Space Sci. Rev., 123, 485 Garel T., Blaizot J., Rosdahl J., Michel-Dansac L., Haehnelt M. G., Katz H., Kimm T., Verhamme A., 2021, MNRAS, 504, 1902 Gelli V., Mason C., Hayward C. C., 2024, arXiv e-prints, arXiv:2405.13108 Glazebrook K. et al., 2024, Nature, 628, 277 Goerdt T., Dekel A., Sternberg A., Ceverino D., Teyssier R., Primack J. R., 2010, MNRAS, 407, 613 Guseva N. G., Izotov Y. I., Papaderos P., Fricke K. J., 2007, A&A, 464, 885 Guseva N. G., Izotov Y. I., Thuan T. X., 2006, ApJ, 644, 890 Harikane Y., Nakajima K., Ouchi M., Umeda H., Isobe Y., Ono Y., Xu Y., Zhang Y., 2024, ApJ, 960, 56 Heintz K. E. et al., 2024a, arXiv e-prints, arXiv:2407.06287 Heintz K. E. et al., 2024b, arXiv e-prints, arXiv:2404.02211 Heintz K. E. et al., 2024c, Science, 384, 890 Hirano S., Hosokawa T., Yoshida N., Umeda H., Omukai K., Chiaki G., Yorke H. W., 2014, ApJ, 781, 60 Hu W. et al., 2023, ApJ, 956, 39  \nInoue A. K., 2011, MNRAS, 415, 2920  \n- Inoue A. K., Shimizu I., Iwata I., Tanaka M., 2014, MNRAS, 442, 1805\n- Isobe Y., Ouchi M., Nakajima K., Harikane Y., Ono Y., Xu Y., Zhang Y., Umeda H., 2023, ApJ, 956, 139\n- Iyer K., Gawiser E., 2017, ApJ, 838, 127\n- Iyer K. G., Gawiser E., Faber S. M., Ferguson H. C., Kartaltepe J., Koekemoer A. M., Pacifici C., Somerville R. S., 2019, ApJ, 879, 116\n- Izotov Y. I., Guseva N. G., Thuan T. X., 2011, ApJ, 728, 161 Jin Y., Kewley L. J., Sutherland R. S., 2022, ApJ, 934, L8 Katz H. et al., 2023, The Open Journal of Astrophysics, 6, 44 Kereˇs D., Katz N., Weinberg D. H., Dav'e R., 2005, MNRAS, 363, 2\n- Kewley L. J., Nicholls D. C., Sutherland R. S., 2019, ARA&A, 57, 511\n- Kimm T., Haehnelt M., Blaizot J., Katz H., Michel-Dansac L., Garel T., Rosdahl J., Teyssier R., 2018, MNRAS, 475, 4617 Kobayashi C., Karakas A. I., Lugaro M., 2020, ApJ, 900, 179 Kravtsov A., Belokurov V., 2024, arXiv e-prints, arXiv:2405.04578 Larkin M. M., Gerasimov R., Burgasser A. J., 2023, AJ, 165, 2 Laseter I. H. et al., 2024, A&A, 681, A70  \nLeitherer C. et al., 1999, ApJS, 123, 3  \n- Leung G. C. K. et al., 2023, ApJ, 954, L46\n- Li Y., Leja J., Johnson B. D., Tacchella S., Naidu R. P., 2024, ApJ, 969, L5  \nLin X. et al., 2023, ApJ, 944, L59  \n- Lu S., Frenk C. S., Bose S., Lacey C. G., Cole S., Baugh C. M., Helly J. C., 2024, arXiv e-prints, arXiv:2406.02672\n- Luridiana V., Morisset C., Shaw R. A., 2015, A&A, 573, A42 Ma X. et al., 2018, MNRAS, 478, 1694\n- Mas-Ribas L., Dijkstra M., Forero-Romero J. E., 2016, ApJ, 833, 65\n- Maseda M. V. et al., 2020, MNRAS, 493, 5120\n- Maseda M. V. et al., 2024, arXiv e-prints, arXiv:2403.05506  \nMason C. A., Trenti M., Treu T., 2023, MNRAS, 521, 497  \nMatthee J., Sobral D., Best P., Khostovan A. A., Oteo I.,  \nBouwens R., Rottgering H., 2017, MNRAS, 465, 3637  \nMcCaffrey J., Hardin S., Wise J. H., Regan J. A., 2023, The  \nOpen Journal of Astrophysics, 6, 47  \n- Morishita T. et al., 2024, arXiv e-prints, arXiv:2402.14084 Omukai K., Tsuribe T., Schneider R., Ferrara A., 2005, ApJ, 626, 627\n- Osterbrock D. E., Ferland G. J., 2006, Astrophysics of gaseous nebulae and active galactic nuclei\n- Pallottini A., Ferrara A., 2023, A&A, 677, L4\n- Pallottini A. et al., 2022, MNRAS, 513, 5621\n- Panagia N., 2002, arXiv e-prints, astro  \nPecaut M. J., Mamajek E. E., 2013, ApJS, 208, 9  \n- Pequignot D., Petitjean P., Boisson C., 1991, A&A, 251, 680\n- Pierel J. D. R. et al., 2024, arXiv e-prints, arXiv:2406.05089\n- Raiter A., Schaerer D., Fosbury R. A. E., 2010, A&A, 523, A64\n- Ren K., Trenti M., Mason C. A., 2019, ApJ, 878, 114", 'SPECTRAL PROPERTIES OF THE NEBULAR CONTINUUM': 'Our primary analysis has focused on spectral properties of the nebular continuum combined with the transmitted stellar continuum. However, as we have shown in Section 3, a subset of galaxies appear to have spectra with potentially higher nebular fractions than predicted by the various SPS models. Here we discuss the spectral properties of purely nebular spectra. For all calculations we use PyNeb (Luridiana, Morisset & Shaw 2015) assuming Case B recombination, accounting for hydrogen and singly ionized helium with He I /H I = 0 . 08. A gas density of 10 2 cm -3 is used for all models.', 'Nebular-only UV slopes': "In Figure 18 we show the UV slope, β , as a function of temperature for the nebular continuum. As above, β is measured in the rest-frame wavelength interval 1400 ˚ A ≤ λ ≤ 2600 ˚ A. For typical H II region temperatures of T ≳ 10 4 K, β is redder than -1 . 6, reaching a peak value of ∼ -1 . 25. For gas temperatures below ∼ 18 , 000 K, β is a direct measure of gas temperature. However, above this temperature threshold, the temperature corresponding to UV slope becomes double valued. Other indicators would need to be used to elucidate this ambiguity.  \nAt higher gas densities, where two-photon emission is suppressed, β can further redden to strongly positive values but we do not explore this further as this is unlikely for most H II regions.  \nRoberts-Borsani G. et al., 2024, arXiv e-prints, arXiv:2403.07103 Robertson B. E. et al., 2013, ApJ, 768, 71 Rosdahl J., Blaizot J., 2012, MNRAS, 423, 344 Sanders R. L., Shapley A. E., Topping M. W., Reddy N. A., Brammer G. B., 2023, ApJ, 955, 54 Sandles L. et al., 2023, arXiv e-prints, arXiv:2306.03931 Saxena A. et al., 2024, A&A, 684, A84 Saxena A. et al., 2020, A&A, 636, A47 Schaerer D., 2002, A&A, 382, 28 Schaerer D., 2003, A&A, 397, 527 Schaerer D., Guibert J., Marques-Chaves R., Martins F., 2024, arXiv e-prints, arXiv:2407.12122 Scholtz J. et al., 2023, arXiv e-prints, arXiv:2311.18731 Senchyna P., Plat A., Stark D. P., Rudie G. C., Berg D., Charlot S., James B. L., Mingozzi M., 2024, ApJ, 966, 92 Shen X., Vogelsberger M., Boylan-Kolchin M., Tacchella S., Kannan R., 2023, MNRAS, 525, 3254 Siebert M. R. et al., 2024, arXiv e-prints, arXiv:2406.05076 Smith A., Ma X., Bromm V., Finkelstein S. L., Hopkins P. F., Faucher-Gigu'ere C.-A., Kereˇs D., 2019, MNRAS, 484, 39 Stacy A., Bromm V., Lee A. T., 2016, MNRAS, 462, 1307 Stanway E. R., Eldridge J. J., 2018, MNRAS, 479, 75 Strom A. L. et al., 2023, ApJ, 958, L11 Sun G., Faucher-Gigu'ere C.-A., Hayward C. C., Shen X., Wetzel A., Cochrane R. K., 2023, ApJ, 955, L35 Tacchella S. et al., 2024, arXiv e-prints, arXiv:2404.02194 Tacchella S. et al., 2022, MNRAS, 513, 2904 Tanvir N. R. et al., 2019, MNRAS, 483, 5380 Terp C., Heintz K. E., Watson D., Brammer G., Carnall A., Witstok J., Smit R., Vejlgaard S., 2024, arXiv e-prints, arXiv:2404.06543 Topping M. W. et al., 2024a, MNRAS, 529, 4087 Topping M. W. et al., 2024b, arXiv e-prints, arXiv:2407.19009 Topping M. W. et al., 2024c, MNRAS, 529, 3301 Trussler J. A. A. et al., 2023, MNRAS, 525, 5328 Umeda H., Ouchi M., Nakajima K., Harikane Y., Ono Y., Xu Y., Isobe Y., Zhang Y., 2023, arXiv e-prints, arXiv:2306.00487 Vanzella E. et al., 2022, ApJ, 940, L53 Vink J. S., de Koter A., Lamers H. J. G. L. M., 2001, A&A, 369, 574 Wang B. et al., 2023, ApJ, 957, L34 Welch B. et al., 2024, arXiv e-prints, arXiv:2405.06631 Willott C. J. et al., 2023, arXiv e-prints, arXiv:2311.12234 Yajima H., Li Y., Zhu Q., Abel T., 2012, MNRAS, 424, 884 Yanagisawa H. et al., 2024, arXiv e-prints, arXiv:2403.20118 Yung L. Y. A., Somerville R. S., Finkelstein S. L., Wilkins S. M., Gardner J. P., 2024, MNRAS, 527, 5929 Zackrisson E., Inoue A. K., Jensen H., 2013, ApJ, 777, 39  \n- Zavala J. A. et al., 2024, arXiv e-prints, arXiv:2403.10491  \n<!-- image -->  \n<!-- image -->  \nFig. 18.(Left) UV slope, β as a function of gas temperature for the nebular continuum assuming Case B recombination. (Center) Equivalent widths of H α and H β as a function of temperature. Solid and dashed lines show the results with and without two-photon emission, respectively. (Right) log 10 ( ξ ion ) as a function of gas temperature for the nebular continuum.  \n<!-- image -->", 'Balmer Equivalent Widths': 'In the middle panel of Figure 18 we show the equivalent widths (EWs) of H α and H β as a function of temperature. These quantities are strong functions of temperature such that at higher temperatures, the predicted EWs drop significantly. For example, at a temperature of 10 , 000 K, the H β EW is nearly 1400 ˚ A, but at 30 , 000 K, it drops to 658 ˚ A. One of the key features of the Balmer line EWs is that they are also sensitive to the presence of two-photon emission. This is perhaps counter-intuitive as the two-photon emission falls off steeply at wavelengths longer than ∼ 1450 ˚ A (in f λ ). The dashed lines in the middle panel of Figure 18 shows the Balmer line EWs when the two-photon emission has been removed and they are clearly significantly higher in the absence of two-photon emission.  \nIt is important to consider that these calculations all assume Case B recombination. However, as discussed in Raiter, Schaerer & Fosbury (2010); Mas-Ribas, Dijkstra & Forero-Romero (2016), in low-metallicity environments, Case B departures are expected because the gas is hotter and collisional excitation can preferentially enhance both the Ly α emission and the two-photon emission. The Balmer lines are not enhanced nearly as much due to their higher energies and thus lower collisional excitation coefficients. The effect of these Case B departures is that the Balmer EWs would then further decrease in low-metallicity environments.', 'Ionizing photon production efficiencies': 'In the right panel of Figure 18 we show the observed ξ ion that would be inferred from H α and L UV , 1500 ˚ A as a function of gas temperature for the nebular continuum. This quantity again varies as a function of temperature such that the measured ξ ion decreases with increasing gas temperature. This reflects the decrease in line emissivity with respect to the continuum at high-temperatures. Note that the inferred ξ ion values from a nebular-only spectrum are not too dissimilar from the intrinsic values of the SPS models. We again emphasize that any Case B departures or cooling radiation would systematically decrease the observed ξ ion from the Case B values presented here.', 'FITTING FORMULAE FOR ξ ION MEASUREMENTS': 'Measuring the ionizing photon production efficiency, ξ ion requires converting between an observed emission line (e.g. a Balmer line such as H α or H β ) and the total number of ionizing photons. Typically a single value is used for this conversion that assumes Case B recombination and a fixed temperature of, for example, 10 4 K (e.g. Bouwens et al. 2016; Matthee et al. 2017; Maseda et al. 2020). Under these assumptions the value adopted is ∼ 7 . 3 × 10 11 erg -1 (Osterbrock & Ferland 2006). However, as we measure the full posterior distribution on electron temperature from our galaxy fits, we can adopt a temperature-dependent conversion rate. For convenience, we provide fitting formulae to convert H α and H β luminosity into ionizing photon luminosity such that:  \nQ [s -1 ] = 10 -0 . 007358 log 10 ( T ) 2 +0 . 1729 log 10 ( T )+11 . 289 [erg -1 ] × L H α [erg s -1 ] (B1)  \nand  \nQ [s -1 ] = 10 0 . 0309 log 10 ( T ) 2 -0 . 212 log 10 ( T )+12 . 679 [erg -1 ] × L H β [erg s -1 ] . (B2)  \nThese equations assume that the escape fraction is zero and are based on atomic data from Pequignot, Petitjean & Boisson (1991). Note that the ratio of the two conversion factors provides the H α /H β ratio as a function of temperature. Because H II regions in the early Universe tend to be hotter than their low-redshift counterparts, likely due to lower metallicity (e.g. Laseter et al. 2024; Morishita et al. 2024), the effect of using a temperature-dependent conversion factor is that the measured ξ ion values can be elevated by ≳ 10%.  \nOne final effect to consider is that these conversions assume that all Lyman continuum photons emitted by a source population interact with hydrogen. In practice other species such as helium (and possibly dust) will be present which  \nreduces the ionizing photon budget for hydrogen by ∼ 10% (e.g. Tacchella et al. 2022). In this situation, the Q measured from a Balmer line will be an underestimate of the true value.  \nThis paper was built using the Open Journal of Astrophysics L A T E X template. The OJA is a journal which provides fast and easy peer review for new papers in the astro-ph section of the arXiv, making the reviewing process simpler for authors and referees alike. Learn more at http://astro.theoj.org .'}