Nebular dominated galaxies: insights into the stellar initial mass function at high redshift

We identify a low-metallicity ($12+\log({\rm O}/{\rm H})=7.59$) Ly$\alpha$-emitting galaxy at $z=5.943$ with evidence of a strong Balmer jump, arising from nebular continuum. While Balmer jumps are sometimes observed in low-redshift star-forming galaxies, this galaxy also exhibits a steep turnover in the UV continuum. Such turnovers are typically attributed to absorption by a damped Ly$\alpha$ system (DLA); however, the shape of the turnover and the high observed Ly$\alpha$ escape fraction ($f_{\rm esc,Ly\alpha}~\sim27\%$) is also consistent with strong nebular two-photon continuum emission. Modelling the UV turnover with a DLA requires extreme column densities ($N_{\rm HI}>10^{23}$ cm$^{-2}$), and simultaneously explaining the high $f_{\rm esc,Ly\alpha}$ requires a fine-tuned geometry. In contrast, modelling the spectrum as primarily nebular provides a good fit to both the continuum and emission lines, motivating scenarios in which (a) we are observing only nebular emission or (b) the ionizing source is powering extreme nebular emission that outshines the stellar emission. The nebular-only scenario could arise if the ionising source has `turned off' more recently than the recombination timescale ($\sim$1,000 yr), hence we may be catching the object at a very specific time. Alternatively, hot stars with $T_{\rm eff}\gtrsim10^5$ K (e.g. Wolf-Rayet or low-metallicity massive stars) produce enough ionizing photons such that the two-photon emission becomes visible. While several stellar SEDs from the literature fit the observed spectrum well, the hot-star scenario requires that the number of $\gtrsim50~{\rm M}_\odot$ stars relative to $\sim5-50~{\rm M}_\odot$ stars is significantly higher than predicted by typical stellar initial mass functions (IMFs). The identification of more galaxies with similar spectra may provide evidence for a top-heavy IMF at high redshift.


INTRODUCTION
The unprecedented sensitivity of JWST in the near-infrared has revolutionised our ability to study the rest-frame ultraviolet to optical spectra of high-redshift galaxies.JWST spectroscopy has already unveiled large samples of emission line galaxies at  ≳ 5, providing new insights in to the conditions of the interstellar media (ISM) of these galaxies (Matthee et al. 2023;Cameron et al. 2023b;Mascia et al. 2023;Sanders et al. 2023;Boyett et al. 2024;Roberts-Borsani et al. 2024).Furthermore, the sensitivity of JWST/NIRSpec has even enabled high fidelity spectroscopic continuum measurements, which have provided insights into high-redshift neutral gas (Heintz et al. 2023;Umeda et al. 2023), assembly of massive galaxies in the early Universe (Carnall et al. 2023;Glazebrook et al. 2023), and have also provided the many of the highest spectroscopic redshift confirmations to date (Curtis-Lake et al. 2023;Arrabal Haro et al. 2023).Spectral ★ E-mail: alex.cameron@physics.ox.ac.uk energy distribution (SED) fitting of many photometric data sets has indicated a need for strong nebular emission, including predictions for strong Balmer jumps (Endsley et al. 2023;Topping et al. 2023).The strong nebular contribution implied by these fits motivates more detailed studies examining the contribution of the nebular continuum to the integrated spectra of high-redshift galaxies.
The nebular continuum is comprised of three main components: (1) free-bound emission, (2) free-free emission, and (3) two-photon emission.Free-bound emission is the component that is most commonly observed to make a significant contribution to the integrated optical spectra of galaxies, presenting in the form of the 'Balmer jump', a discontinuity in the spectrum at  rest = 3645 Å, corresponding to the Balmer limit.This feature arises due to electron recombinations with ionized hydrogen to the first excited state.Balmer jumps only appear in the spectra of galaxies that contain young stellar populations with high ionizing photon production efficiencies ( i ).They have been detected in some highly star-forming galaxies at lowredshift (Peimbert & Costero 1969;Guseva et al. 2006Guseva et al. , 2007)).Many numerical simulations predict Balmer jumps to be common at high redshift (e.g.Katz et al. 2023a;Wilkins et al. 2023), in line with the early SED fitting results outlined above.In contrast, predictions based on typical stellar population models indicate that the other nebular continuum components, free-free and two-photon emission, are generally expected to be subdominant compared to the stellar and free-bound contributions in the integrated spectra of galaxies.
The predicted subdominance of two-photon emission arises because the  ion values of typical low-metallicity stellar populations are not high enough to drive nebular continuum that overcomes the steep UV stellar continuum slopes (Leitherer et al. 1999).In contrast, very hot stars, with blackbody temperatures of ∼ 100, 000 K, produce enough ionizing photons to power nebular continuum emission that outshines their stellar UV emission, and are predicted to exhibit a strong two-photon continuum bump, peaking at  rest ≈ 1430 Å (Schaerer 2002;Raiter et al. 2010;Zackrisson et al. 2011;Trussler et al. 2023).Observations of the Lynx arc, a strongly-lensed extreme emission system at  = 3.357, have suggested that the UV continuum of this system has a strong contribution from two-photon continuum, leading to the suggestion of the presence of hot stars with  eff ∼ 80, 000 K (Fosbury et al. 2003).But this has so far remained an isolated candidate for this type of system, with other examples of two-photon continuum being the domain of unusual nebular-only systems (e.g.Lintott et al. 2009).Nonetheless, the identification of these objects, should they exist, offers powerful insight into the properties of systems with extreme ionising spectra.
In this paper, we identify one object, JADES-GS+53.12175-27.79763(GS-NDG-9422 hereafter), at  = 5.943, that exhibits high equivalent-width emission lines as well as a strong spectral discontinuity near  rest = 3645 Å that we interpret as a Balmer jump.The UV continuum exhibits a strong turnover which could be indicative of either a strong damped Lyman- absorption (DLA) system (e.g.Heintz et al. 2023), or the two-photon continuum introduced above.We present a detailed analysis of this system, outlining a number of scenarios that can potentially explain the physical origin of this intriguing spectrum.
The paper is structured as follows: Section 2 outlines the details of data used in this work.In Section 3 we present emission-linebased analysis of the nebular conditions.We perform fitting to the continuum shape in Section 4, comparing the results of DLA and nebular models.Section 5 then considers the nebular case in more detail, exploring the necessary properties of the ionising source.Section 6 outlines the broader implications of the scenarios presented.We then summarise our findings in Section 7. Throughout this paper we adopt the cosmology of Planck Collaboration et al. (2016) with  0 = 67.31km s −1 Mpc −1 and Ω m = 0.315, which gives a luminosity distance to GS-NDG-9422 of   = 58.5 Gpc.

DATA
JWST/NIRSpec spectroscopy of GS-NDG-9422 was taken as part of the JADES survey (PID: 1210, PI: Luetzendorf) in five spectral modes, receiving 28 hours integration in Prism/CLEAR and 7 hours integration in each of G140M/F070LP, G235M/F170LP, G395M/F290LP and G395H/F290LP.We use the reduced spectra released as part of the JADES Public Data Release (Eisenstein et al. 2023;Bunker et al. 2023).GS-NDG-9422 falls within the JWST/NIRCam footprint of JADES (PID: 1880, PI: Eisenstein) as well as the JWST Extragalactic Medium-band Survey (JEMS; PID: 1963; PI: Williams; Williams et al. 2023), combining to provide imaging in 14 wide-and medium-band filters.In Figure 1 we show photometry from the publicly released photometric catalogs (Rieke et al. 2023) compared to the observed Prism/CLEAR spectrum.Blue squares in the second-top panel show aperture photometry measured within a 0.15" radius.Orange squares show the predicted values by convolving the observed NIRSpec spectrum with the NIRCam filter transmission profiles.We find there is good agreement across the full spectral range, suggesting the flux calibration of the Prism/CLEAR spectrum is robust.The exception to this is the two filters covering H emission, 444 and 460.In each case, the flux from imaging is 14 % and 17 % higher, respectively, suggesting the flux of H may be marginally underestimated in the Prism/CLEAR flux calibration.We note that the 460 flux changes by < 2 % when measured within apertures with radii of 0.10" or 0.25", suggesting that this offset is not driven by spatial variations in H.

Emission line fitting
Where possible, emission lines were fit with a single component Gaussian profile with the continuum modelled as a 1D spline.In cases where lines are sufficiently blended, we fit the whole complex with one component.In some cases, partially blended lines are fit simultaneously with neighbouring lines and fluxes reported separately.These are marked in Table 1.We fit all identifiable lines and report upper limits for notable undetected lines.
Line fluxes from higher resolution grating spectra of GS-NDG-9422 were measured independently.Emission line widths are only spectrally resolved in the high-resolution G395H grating.These show no evidence of a broad component and are well fit with a single component with velocity dispersions < 200 km s −1 (Table 2).
Fluxes derived from different observations are mildly discrepant in some cases.Notably, H, [Oiii] 4959, [Oiii] 5007 and H lines exhibit higher fluxes in the grating modes.This behaviour is reported in Bunker et al. (2023) who suggest that the Prism flux calibration is more reliable.This conclusion is supported by the good agreement observed in Figure 1 between the Prism spectrum and NIRCam photometry, with the possible exception of H.We note that in low-resolution data, the continuum level for some emission lines can be difficult to determine, especially for Ly, He ii + O iii] and [O ii] 3727, which introduces uncertainty into the emission line flux.The clear detection of the continuum across almost the entire Prism coverage, allows us to derive equivalent widths directly from this spectrum (Table 2).Given the noted discrepancy on the H flux between the imaging and spectroscopy, we also derive EW 0 (H) directly from the imaging.Comparing the measured flux in 460 (196 ± 5 nJy), shown in Figure 1 to be clearly elevated due to contamination from H, with 430 (25 ± 4 nJy), which captures clean continuum between He i 5875 and H, we find the imaging implies EW 0 (H) = 2195 ± 400 Å.Within the uncertainty, the measured difference between EW 0 (H) Prism and EW 0 (H) Imaging (Table 2) is consistent with the discrepancy noted above between the observed 460 flux, and the predicted 460 flux obtained by convolving the Prism spectrum with the NIRCam filter profile.Throughout our analysis, we adopt the Prism fluxes where possible.However, we ensure that conclusions presented in this work are also consistent with measured grating flux ratios.

DERIVATION OF PHYSICAL CONDITIONS
We now explore the physical conditions of the gas in GS-NDG-9422 and the basic properties of the ionising source, based on the measured emission line fluxes.Throughout this section, we make

Diagnostic diagrams
To explore the properties of the ionising source, in Figure 2, we look at nebular line ratio diagnostic diagrams, comparing measured line ratios from GS-NDG-9422 with photoionisation model predictions for star formation and active galactic nuclei.For star forming models, we adopt the predictions of Gutkin et al. (2016).These use input stellar SEDs based on plane-parallel non-local thermal equilibrium models calculated with Tlusty (e.g.Lanz & Hubeny 2007).Stellar abundance patterns are assumed to follow scaled-solar, which may not be representative of  ∼ 6 stellar populations, while nebular abundance patterns allow for some variation in C/O and N/O ratios.Star forming models with  = 0.001 (the grid value closest to the gasphase oxygen abundance derived for GS-NDG-9422; see Section 3.   2016) across a large range of metallicities.Navy (star-forming) and red (AGN) points show the subsets of these model grids that have  = 0.001 ( ≈ 0.07 ⊙ ; 12 + log(O/H) ≈ 7.54) which is the closest grid value to that measured for GS-NDG-9422 (see Section 3.5).Top left: The strong upper limit on [S ii] 6716, 6731 positions GS-NDG-9422 well below the theoretical maximum starburst limit from Kewley et al. (2001) (black line).Top right: The weak detection of He ii 4686 is also below the maximum starburst limit from Shirazi & Brinchmann (2012) (black line).Strong He ii-emitting star-forming galaxies from Shirazi & Brinchmann (2012) and Senchyna et al. (2017) are also shown as orange and green diamonds respectively.We also show measurements from Hanny's Voorwerp, suggested to be a quasar light echo (Lintott et al. 2009; purple line-connected points), which we discuss in Section 5.1.Points of increasing size indicate larger offset from the quasar ranging from 13 to 31 kpc.Bottom left: GS-NDG-9422 lies beyond the AGN region defined by Mingozzi et al. (2024) (black lines), and is more coincident with the star-forming models than AGN models.He ii-emitting star-forming galaxies at  ∼ 2.5 − 4 from Saxena et al. (2020) are shown for comparison as orange squares.Bottom right: GS-NDG-9422 lies at the tip of the parameter space covered by AGN models, while it is well within the bounds of that covered by star-forming models.left panel of Figure 2) and in a region that is difficult to reconcile with emission powered by an AGN.Non-detections of [N ii] 6583 and [O i] 6300 paint a similar picture in the [N ii]-BPT diagram (Baldwin et al. 1981) and the [O i]-VO87 diagram, although these are not shown.
The top right panel of Figure 2 shows that the He ii 4686 / H ratio of GS-NDG-9422 exceeds that predicted by any of the star-forming models in Gutkin et al. (2016), especially those at  = 0.001.However, it is well known that, at low metallicity, star-forming galaxies are often observed with He ii emission exceeding that which can be powered by standard stellar population models (Shirazi & Brinchmann 2012;Kehrig et al. 2015;Senchyna et al. 2017;Schaerer et al. 2019;Saxena et al. 2020).He ii-selected star-forming galaxies from Shirazi & Brinchmann (2012) and Senchyna et al. (2017) are shown as orange and green diamonds respectively.We find that while GS-NDG-9422 is at the upper end of the He ii 4686 / H distribution from these works, it does not exceed the very highest ratios, and also falls below the He ii maximum-starburst demarcation presented in Shirazi & Brinchmann (2012).The origin of strong He ii emission at low metallicity remains an unsolved problem with possible solutions including revised stellar wind properties, X-ray binaries, very hot binary stellar evolution products, rotating massive low-metallicity Green squares give the measured ratios from the Prism, while diamonds give ratios measured from gratings, where available.Solid black lines give the theoretical ratio for case B recombination at   = 10 4 K and   = 100 cm −3 .Coloured points show theoretical ratios for a range of temperatures (10 4 ≤   ≤ 2.5 × 10 4 K; colour) and densities (  = 10, 100, 1000 cm −3 ; marker size).The measured H/H is somewhat lower than the theoretical value, but we note that H may be underestimated in the prism spectrum (Table 2).stars, and/or a top-heavy stellar initial mass function (Kehrig et al. 2015(Kehrig et al. , 2018;;Schaerer et al. 2019;Senchyna et al. 2020;Olivier et al. 2022).We return to this question in Section 5. We note that GS-NDG-9422 was included in a selection of narrow-line AGN in Scholtz et al. (2023) on the basis of the He ii 4686 / H ratio.However, we consider the measured value of this line ratio to be inconclusive.
In the bottom left panel of Figure 2, we see that the Oiii] 1660, 1666 / He ii 1640 ratio measured for GS-NDG-9422 exceeds the AGN locus defined by Mingozzi et al. (2024) and has a value which cannot be reproduced by the Feltre et al. (2016) AGN models.The measured ratio is more in line with that measured in the  ∼ 2−4 He iiselected star-forming galaxies from Saxena et al. (2020).Meanwhile the (C iv + C iii]) / He ii ratio in the bottom right panel is only reproduced by the very tip of the AGN model parameter space, being more readily reproduced by the star-forming models.
In summary, AGN photoionisation models struggle to reproduce a number of the observed emission line ratios, especially those involving low-ionisation emission lines.In contrast, GS-NDG-9422 resides within regions of line-ratio space consistent with emission powered by stars.We therefore conclude that the ionising source in GS-NDG-9422 is most likely of stellar origin.

Balmer decrements and dust extinction
Balmer decrements H/H and H/H from the Prism and H9/H from the grating are consistent with Case B values, indicating that there is no significant dust reddening in GS-NDG-9422 (Figure 3).Measured H/H ratios are lower than theoretically predicted for Case B at  = 10 4 K. Note, this is not suggestive of dust reddening, which would act in the opposite direction.At higher temperatures, the theoretical ratio decreases, and our G395H/F290LP measurement is consistent with theoretical ratios with   ≳ 2 × 10 4 K, possibly indicative of a very hot nebula.As noted in Section 2.1, H may be underestimated in the Prism which could lead to the slightly lower observed ratio.H/H measured from the medium-resolution grating is marginally below the theoretical value.In isolation, this could suggest non-zero dust reddening; however, this evidence is outweighed by the other measured ratios.The high-resolution grating returns a much lower H/H = 0.3 ± 0.03.However, we note that the continuum is undetected in the high-resolution grating, which contributes uncertainty to the measured ratio.
Adopting a luminosity distance of   = 58.5 Gpc, we derive  H = 1.86 × 10 42 erg s −1 .Assuming no dust, a metallicity of  = 0.1 ⊙ , and a typical IMF, this would correspond to a star formation rate of ∼ 3.2 M ⊙ yr −1 (Eldridge et al. 2017).Under the assumption of no dust, Case B recombination,   = 100 cm −3 and   = 1.8×10 4 K, we derive a Ly escape fraction of  esc,Ly = 0.29±0.01from the Ly/H ratio, or  esc,Ly = 0.27 ± 0.01, measured from the Ly/H ratio, the latter of which may be more reliable owing to the H flux uncertainty discussed above.

Electron temperature
The temperature-sensitive [O iii] 4363/5007 ratio can be measured from each of the Prism, G395M and G395H observations, yielding three consistent, independent temperature measurements (  = 1.83±0.15,1.99±0.18,and 1.81±0.18×10 4K, respectively).The temperature derived from the medium-resolution O iii] 1666 / [O iii] 5007 ratio is somewhat lower (  = 1.70 +0.05 −0.06 × 10 4 K).However, the He ii+O iii] flux measured from the medium-resolution G140M grating is significantly lower than that of the Prism.Instead, the measured temperature from above (  = 1.83 × 10 4 K) implies 1660,1666/5007 = 0.08, which gives  1660,1666 = 8.2 ± 0.1 × 10 −19 erg s −1 cm −2 based on the [O iii] 5007 Prism flux, suggesting that O iii] contributes ∼55 % of the measured Prism He ii+O iii] blend.This is consistent with the O iii] / He ii ratio measured in G140M.The implied He ii 1640 flux (6.9 ± 1.1 × 10 −19 erg s −1 cm −2 ) gives a He ii 1640/4686 ratio of 6.3 ± 1.5 which is consistent with the theoretical value of 7.19 assuming Case B at   = 1.8 × 10 4 K and   = 100 cm −3 .Note that this further supports the conclusion of a lack of dust in this system since the He ii 1640 would be subject to extremely high extinction, relative to the 4686 line.In systems with a strong nebular continuum component,   can also be constrained from the magnitude of the Balmer jump (e.g.Guseva et al. 2006;Pérez-Montero 2017).We return to this in Section 4 where we show that   (H + ) implied by the Balmer jump is consistent with   (O ++ ) measured from the [O iii] auroral line ratio.

Electron density
The density-sensitive doublets of C iii] and [O ii] are not resolved in our observations, while N iv] and [S ii] are not detected.We report a marginal detection of the [Ar iv] 4711, 4740 doublet in the Prism spectrum.These lines are partially blended with each other and with He ii 4686, while [Ar iv] 4711 is also completely blended with He i 4713.Our three-component fit to this complex yields [Ar iv] 4711/4740 = 1.6 ± 0.8 after subtracting the predicted He i 4713 contribution (assuming 4713/4471 = 0.15), consistent with the low-density limit (  ≲ 10 3 cm −3 ; Kewley et al. 2019).A consistent density constraint arises if the UV continuum turnover in GS-NDG-9422 is driven by two-photon nebular continuum, since the feature strongly suppressed by -changing collisions at higher densities.The presence of the two-photon continuum will be discussed in Section 4.He ++ /He + 0.06 ± 0.01

Chemical abundances
We derive chemical abundances for GS-NDG-9422 adopting   = 1.83 × 10 4 K and a density of   = 100 cm −3 following the procedure in Cameron et al. (2023a).Given the apparent agreement between   (H + ) and   (O ++ ), we assume a constant temperature for all ionisation zones.We derive 12 and H lines, obtaining a consistent result with both Prism and grating line fluxes.We derive the carbon abundance from the C iii] 1907, 1909 / [O iii] 5007 ratio measured from the Prism, assuming no dust, finding log(C/O) = −0.73 ± 0.03 after applying the ionisation correction factor (ICF) presented in Amayo et al. (2021).However, we note the significant emission from C iv in the spectrum.The ICF may not be representative of the extreme conditions in GS-NDG-9422 (e.g.Berg et al. 2019), and the quoted C/O may represent a lower limit.Since no nitrogen lines are detected in GS-NDG-9422, we can only place upper limits on the nitrogen abundances.The low O + abundance implies that the N + abundance is likely negligible.We instead consider the N ++ abundance using N iii] 1750 / [O iii] 5007 limit measured from the Prism, and the N iii] 1750 / O iii] 1666 limit measured from the grating.These yield 3- upper limits of log(N ++ /O ++ ) < −0.85 and < −1.01 respectively.We adopt the former as our preferred limit.
We detect several helium recombination lines from both the singly and doubly ionised states.Prism measurements of He i 4471 / H and He i 5875 / H yield consistent measurements of He + /H + = 0.10 ± 0.005 and 0.10 ± 0.01 respectively.Deriving a He ++ /H + abundance using the He ii 4686 line results in only a 6 ± 1 % contribution to the total helium abundance, giving He ++ /H + = 0.005±0.001.Together, this implies a total helium abundance of He/H = 0.11 ± 0.01, higher than typical values observed in low-metallicity systems (Matsumoto et al. 2022), but consistent with some massive globular clusters (Piotto et al. 2007).However, we caution that our derived helium abundance does not account for collisional emission or self-absorption, which could be significant (e.g.Peimbert et al. 1992).
The values obtained from these measurements are summarised in Table 3.

Photoionization Modelling of the Emission Lines
The diagnostic diagrams, electron temperature, and line widths all point to a scenario where GS-NDG-9422 is powered by emission from young stars.We explore the feasibility of reproducing the emission lines with young stellar populations with photoionization models using CLOUDY v23 (Ferland et al. 2017).We assume that the intrinsic spectrum is powered by a standard young SSP model, adopting BPASS v2.2.1 SSP models including binary stellar evolution (Eldridge et al. 2017) with an IMF maximum mass of 300 M ⊙ and a high-mass slope of −2.35.We consider a population with an age of 3 Myr and a metallicity of 0.1 ⊙ , consistent with that measured for the gas from the spectrum.We note that the abundance patterns assumed in these models are scaled solar, which may not be representative of the stellar populations forming at this redshift.We adopt a spherical geometry with an inner radius of 0.1 pc.The calculation is stopped at an electron fraction of 1% or when the neutral column density reaches 10 18.7 cm −2 , consistent with the minimal observed Ly emission offset (≲ 100 km s −1 ) (Verhamme et al. 2015).We then vary gas density, ionization parameter, metallicity (assuming solar abundance patterns from Grevesse et al. 2010), and carbon abundance, fixing the gas temperature to that measured from [O iii] 4363/5007, until we reproduce the line strengths of , H, and H with respect to H.Fe is assumed to be heavily depleted or to have not yet been produced.We also assume resonant lines from low-ionization species have the same escape fraction as that measured for Ly.Note that by focusing on line ratios, the resulting continuum is a prediction of the model.The intrinsic spectrum of our best fit model ( = 10 3 cm −3 , log 10 () = 1.2, log 10 ( O / ⊙ ) = −1.1, and log 10 ( C / ⊙ ) = −1.4) is shown as the magenta line in Figure 4.
As can be seen in Figure 4, the strengths of the emission lines are well reproduced by our best fit CLOUDY model, highlighting the fact that relatively metal-poor, young stellar populations are able to explain the emission lines seen in this galaxy.There remains a discrepancy between the observed He ii emission and that predicted by the photoionization models 1 , but this is consistent with numerous low-redshift metal-poor galaxies that show anomalous He ii emission with no signatures of a black hole or an AGN (Kehrig et al. 2015(Kehrig et al. , 2018;;Schaerer et al. 2019;Senchyna et al. 2020;Saxena et al. 2020).

CONTINUUM SHAPE
We have shown in the previous section that the emission lines are consistent with ionisation by young metal-poor stellar populations.In this section, we explore the possible origin of the continuum emission in GS-NDG-9422.

Balmer jump
Low-metallicity, young stellar populations can have ionizing photon production efficiencies that are high enough that the free-bound emission from recombining hydrogen can outshine the stellar continuum at optical wavelengths, leading to the observation of a Balmer jump (e.g.Byler et al. 2017).As seen in Figure 1, already from the broad-and medium-band photometry alone, there is evidence of a spectral discontinuity at the location of the Balmer limit.The Balmer jump is not typically observed as a sharp spectral feature since, at realistic spectral resolutions, blending of the tail of the Balmer series and strong emission lines like [O ii] 3726, 3729 and [Ne iii] 3869 can occur (e.g.Schirmer 2016).We perform an empirical fit for this discontinuity by masking out regions contaminated by strong emission lines and fitting the spectrum (in units of   ) over the range  rest > 1930 Å with a two-part linear function, broken at the Balmer limit ( rest = 3645 Å).This results in a clear detection of a spectral jump with 15.0 ± 0.9 nJy in the observed frame (second panel Figure 1; dotted line), demonstrating the strong contribution of nebular continuum to the spectrum of GS-NDG-9422.Such features have been observed in the spectra of metal-poor star-forming galaxies at low-redshift (Guseva et al. 2006(Guseva et al. , 2007) ) and at high-redshift (Roberts-Borsani et al. 2024), while they have also been widely predicted at high-redshift in SED modelling of photometric data (Endsley et al. 2023;Topping et al. 2023) and in simulations (Katz et al. 2023a;Wilkins et al. 2023).
As shown in Figure 4 the spectral discontinuity at the location of the Balmer jump is well reproduced by our best fit photoionization model.It is important to emphasize here that the presence of spectral discontinuity in the continuum is a prediction of the photoionization model which was purely designed to reproduce the emission lines rather than the shape of the continuum.

UV continuum turnover
Even more striking than the spectral discontinuity at the location of the Balmer jump is the presence of a steep turnover in the UV continuum of GS-NDG-9422 at  obs ≈ 1 m ( rest ≈ 1430 Å) (Figure 1).While our best-fit photoionization model is successful in reproducing the line emission and the Balmer jump, where the model fails is that it significantly over-predicts the emission at  rest ≲ 1430 Å.Here we explore the origin of this discrepancy.

Absorption from neutral hydrogen?
Similar UV turnovers are often observed as a result of absorption from foreground neutral hydrogen -either the neutral intergalactic medium (IGM) (Miralda-Escudé 1998) or Damped Lyman- absorption (DLA) systems (Heintz et al. 2023).
Beginning with the IGM, we apply the  = 6 IGM transmission curves from Garel et al. (2021), to our best-fit photoionization model and we find that IGM damping is unable to produce the observed UV turnover from this BPASS model (yellow dashed line in Figures 4 & 5).Indeed, the observed turnover in GS-NDG-9422 requires absorption far exceeding the maximal IGM damping wing, calculated following the formalism presented in Miralda-Escudé (1998) and Barkana & Loeb (2001).Similar results were found in Heintz et al. (2023) for targets at much higher redshift.This is not particularly surprising because if the IGM were responsible, many more galaxies at  ≳ 6 would show very strong UV turnovers.
We next consider the presence of a DLA.We model damping due to the presence of a DLA with CLOUDY by calculating the transmission of the best-fit BPASS model through slabs of neutral hydrogen with increasing column densities, finding that column densities of  H ∼ 10 23 cm −2 are needed to reproduce the magnitude of the turnover.Such a value is higher than any previously reported galaxy-scale DLA system (e.g.Tanvir et al. 2019;Umeda et al. 2023; see Figure 5).
The primary issue with naively assuming a DLA is that such high column densities imply zero transmission at 1216 Å, conflicting with the strong observed Ly emission.The middle panel of Figure 5 shows that this can, in principle, be reconciled by invoking a DLA with 30 % leakage.However, the plausibility of such an extreme column density with a low covering fraction is unclear, given the fact that other known DLA systems with very high gas columns do not show Ly emission (Umeda et al. 2023;Heintz et al. 2023).In principle, one could shift the DLA to a lower redshift which would allow for some Ly emission to escape as the emitted Ly will be redshifted out of resonance in the reference frame of the DLA; however, this would require much higher column densities than 10 23 cm −2 in order to reproduce the UV turnover.Furthermore, at such high column densities, the gas can self-shield from the local-radiation field, and the core would be expected to be fully molecular.It is not clear whether such high gas columns can be maintained without the gas collapsing and forming stars (Schaye 2001).Moreover, GS-NDG-9422 shows no signatures of dust extinction based on the Balmer decrement, the He II ratio, and the fact that the photoionization model can easily reproduce the observed UV slope without assuming dust (Section 3).If the DLA was present within GS-NDG-9422 at a metallicity of nearly 10% solar, one would need to explain why the dust has not formed or could not survive in a thick neutral cloud.Given the fine-tuned requirements, we consider this picket-fence scenario of a thick DLA with optically-thin channels highly unlikely.Free-bound emission can give rise to a spectral discontinuity at the Balmer limit ( rest ≈ 3645 Å).Two-photon continuum has a fixed shape with a turnover at  rest ≈ 1430 Å, but is usually subdominant compared to continuum emission of the ionising source.Free-free emission is usually sub-dominant at these wavelengths.The relative contribution of two-photon continuum is highly dependent on nebular conditions.
Other geometries may exist (apart from the picket fence or lowerredshift DLA) that could possibly explain GS-NDG-9422.For example, one could consider a foreground DLA and background or spatially offset clouds where emission could be reflected.These scenarios all must be reconciled with the very high Ly escape fraction of ∼ 27% which again seems unlikely.For this reason, we consider other alternatives for the origin of the UV turnover.

Two-photon continuum emission
The nebular continuum consists of three components: free-bound, which gives rise to the spectral discontinuity at  rest ∼ 3645 Å, freefree, which is typically subdominant compared to free-bound and only impacts  rest ≳ 3000 Å, and two-photon emission (Figure 6).Two-photon emission arises from transitions from 2 → 1 in neutral hydrogen, resulting in the emission of two photons whose energies sum to that of Ly.The distribution of photons emitted via this process is symmetric around 2431 Å ( = 1 2  Ly ) when expressed in terms of number of photons per second per frequency interval.However, when expressed in units   , it takes the form of a broad asymmetric peak which turns over at  rest ≈ 1430 Å, remarkably close to the wavelength of the observed UV turnover in GS-NDG-9422.The two-photon continuum is typically subdominant compared to the stellar continuum, but has been predicted to be observable in systems with extremely high ionising photon production efficiency (Fosbury et al. 2003;Raiter et al. 2010).Here we consider the possibility that the observed UV turnover is the two-photon continuum.
To test whether the continuum of GS-NDG-9422 is consistent with being primarily nebular, we consider a model where the spectrum consists of: (i) Two-photon emission (ii) Free-bound & Free-free emission (iii) Young stars with ages < 10 7.5 yrs (iv) Old stars with ages > 10 7.5 yrs The shape of the spontaneous two-photon continuum is fixed 2 , so we only vary its overall normalization.The combination of freebound and free-free emission has a shape that is sensitive to gas temperature and thus we vary both gas temperature and the normalization of this component.Finally the shapes of the stellar spectra are sensitive to age, so we consider the normalizations and ages of each stellar population.In total, the model has seven free parameters.As above, we assume that the stellar component follows BPASS v2.2.1, while the nebular continuum components are computed with PYNEB (Luridiana et al. 2015).Posterior distributions for each parameter are computed using an MCMC (Foreman-Mackey et al. 2013).
In the top panel of Figure 7 we show the model corresponding to the 50th percentile distribution of each of the seven parameters (blue line) as well as the points used for the fit (red) and the observed spectrum (green).The MCMC prefers a model where the continuum at  rest < 5800 Å is dominated by the nebular continuum.Indeed the two-photon continuum is able to reproduce of the shape of the UV downturn while free-bound emission peaks high enough to create the observed spectral discontinuity near the Balmer jump.Applying frequentest statistics, we find a reduced  2 of  2  = 0.80 indicating that this 50th percentile model is a very good fit to the continuum of GS-NDG-9422.
Although this nebular-dominated scenario provides a good fit to the continuum, it remains an open question of whether the model is consistent with the observed emission lines.This can be verified in two ways.Because the free-free and free-bound emission are sensitive to gas temperature, we can check whether the temperature predicted from the continuum is consistent with that measured from the oxygen emission lines.In the bottom panel of Figure 7, we show the posterior distribution on gas temperature predicted by the continuum fitting (black histogram) compared to that measured using the [O III] 4363 auroral line (cyan), finding that the two are consistent within 1 uncertainty.While formally the O iii temperature does not have to be the same as that of the H II gas, empirical measurements from low-redshift galaxies suggest that the two temperatures are often very similar (e.g.Guseva et al. 2006Guseva et al. , 2007)).
An even stronger test is to compare the observed H equivalent width versus that predicted by the continuum-fitting model.The H emission should primarily arise from recombination, the same as the free-bound continuum.Thus, in a nebular-dominated scenario, EW(H) should only depend on gas temperature and the relative contribution of free-bound and free-free to two-photon emission.The top panel of Figure 8 shows that the 1 uncertainty on the observed EW(H) significantly overlaps the 1 spread in the posterior distribution of predicted EW(H).This again demonstrates that the information in the continuum is sufficient to explain many of the properties of the emission lines.
We can apply the same test for H emission.Since the continuum was only fit up to rest-frame 5800 Å, predicting the H equivalent width requires extrapolating the model.In Figure 8 we compare the posterior distribution on predicted EW(H) to that measured from the prism as well as that inferred from imaging.As discussed above, the EW 0 (H) from the imaging data is somewhat higher than that measured from the prism (although the 1 contours overlap).We 2 A significant contribution of induced two-photon emission can alter the width of the probability distribution of emitted photons (e.g.Chluba & Sunyaev 2006), which can in turn shift the wavelength at which the two-photon continuum turns over when expressed in   .However, we found this effect to be completely negligible (< 1 Å) in any of the modelling considered in this work.3; blue).find that our predicted values fall high compared to the prism measurement, but the value from the imaging falls on top of our 1 confidence interval.Hence our model is formally very consistent with the imaging data and consistent within 2 of the prism.Because we have not fit the continuum near H the model could be missing a contribution from older stars, which can increase at these wavelengths without impacting our current fit.Since GS-NDG-9422 has metals, these must have originated somewhere.A metallicity of 0.1  ⊙ is much higher than that predicted for the IGM at  ∼ 6 (e.g.Madau & Dickinson 2014) and thus it is unlikely the system was enriched externally.While it is possible that we are witnessing the illumination of the immediate enrichment from the current population of massive stars, the abundance patterns are such that it is The median and 1 values (red lines) are consistent with that measured from the Prism spectrum (blue).Bottom: Same as top panel, but for H.Here, the median and 1 values (red lines) are somewhat higher than that measured from the Prism spectrum (blue), but in good agreement with the equivalent width implied by the medium-band imaging.
likely that there might be an underlying population of stars that is no longer UV-bright.Therefore our current model remains flexible enough to accommodate the scenario of an older population of stars that contributes to the continuum at wavelengths near H.
Nevertheless, given that our model prefers that the UV and optical part of the spectrum arises primarily from nebular emission, it poses the question: what drives this behaviour?

WHAT COULD BE DRIVING TWO-PHOTON EMISSION IN GS-NDG-9422?
Under the assumption that the continuum of GS-NDG-9422 is nebular-dominated, we consider the numerous scenarios that may give rise to the observed spectrum.The magenta line represents a blackbody temperature of 100,000 K. Bottom: Normalized spectrum of hydrogen-only gas irradiated by a blackbody with  = 10 5 K at different gas densities (as given in the colour bar).The magenta line represents a density of 10 3 cm −3 .

Is the ionising source present in the spectrum?
Although Balmer jumps have commonly been observed in the spectra of young, low-metallicity star-forming galaxies (e.g.Guseva et al. 2006Guseva et al. , 2007)), the steep UV slopes associated with the spectra of these stellar populations means that the nebular contribution is completely sub-dominant in the FUV where the two-photon continuum peaks.However, one can envision a scenario where the nebular emission is offset from the location of ionising source.If the slit were to contain only nebular gas, the two-photon emission could dominate the observed spectrum.
One example of this is Hanny's Voorwerp, which has been suggested to be quasar light echo (Lintott et al. 2009).In this case, both the free-bound and two-photon emission are expected to be strong.However, there are three key differences between Hanny's Voorwerp and GS-NDG-9422: (1) He ii 4686/H is > 6× higher in the presumed QSO light echo compared to the galaxy studied here (Figure 2).He ii 1640 is also much stronger in Hanny's Voorwerp (Keel et al. 2012), (2) while the two-photon continuum is likely important in Hanny's Voorwerp, another unidentified source has been postulated to explain the excess UV emission.In other words, the spectrum is not fully nebular dominated in the UV and may consist of additional scattered light, and (3) the spatial extent of Hanny's Voorwerp is tens of kpc.Conversely, GS-NDG-9422 is very compact, and it is well-centred within the NIRSpec micro-shutter (Figure 1).While other QSO light echoes have been detected (e.g.Schirmer et al. 2016), the spectra and morphologies do not seem to be consistent with GS-NDG-9422, leading us to disfavour this scenario.
An alternative explanation is that the ionization source has flickered on and off and we are capturing the system just after it has shut off.In this case, no continuum from the ionising source would be  4; Section 5.2).Middle: Spectrum of GS-NDG-9422 (black) compared with photoionsiation models powered by various individual massive Pop.III stars with different effective temperatures.Bottom Spectrum of GS-NDG-9422 (black) compared with the photoionisation model powered by a model Wolf-Rayet SED with the most similar spectrum to the best-fit blackbody SED.All models can also reproduce the observed spectrum well (dotted lines).See Section 5.2 for details.  10. Green to blue lines show decrease He ii leakage fraction applied to the blackbody, modelling the effect of a He-rich atmosphere.Models with a high leakage fraction clearly overpredict the He ii 4686 line, although we note this feature is blended with nearby [Ar iv] lines.We note a slight flux excess blue-ward of He ii 4686 which, from inspection of the 2D spectrum, appears to be driven by a noisy pixel.This excess flux is still below the level of the excess emission implied by the 100 % He ii leakage model shown in green.
detected and we would be observing only the remnant nebular emission.However, the H + recombination timescale, assuming Case B recombination at a temperature of 18,300 K, is ∼ 5, 000 yr for a density of 10 2 cm −3 or ∼ 500 yr for 10 3 cm −3 .This is much shorter than the main-sequence lifetimes of massive stars.Thus, in the case of flickering, it is highly unlikely that stars are powering the emission.Some QSO proximity zones show evidence for QSO lifetimes of < 10 4 yr (Eilers et al. 2018(Eilers et al. , 2021)).While such QSOs are a rare subset of the general population, these lifetimes are broadly consistent with the recombination timescale.Because these QSOs are detected via their near zone, the QSOs are currently "on" and thus the UV continuum is still dominated by the AGN and not the nebular emission.Evidence for AGN fading on such short time-scales has also been observed locally (e.g.French et al. 2023), but the spectra of these objects, in particular the low ionization state lines is inconsistent with GS-NDG-9422.Moreover, we have only considered the hydrogen recombination time scale.If we consider He ++ under the same assumptions, we arrive at a recombination timescale of ∼ 300 yr for a density of 10 2 cm −3 or ∼ 30 yr for 10 3 cm −3 .The fact that we observe He ii emission is thus difficult to reconcile with this 'flickering' scenario, and imply that we would be catching the source at a very specific moment time.Given the very special timing required, the identification of other objects with spectra like GS-NDG-9422 would seemingly disfavour this scenario.We discuss the identification of similar candidates in Section 6.3.

Is GS-NDG-9422 powered by hot stars?
Assuming that the ionizing source remains luminous and is present within the slit, in order for both the two-photon and free-bound continua to dominate over the stellar spectrum in the rest-frame UV and optical, the source population must have a much larger ion-ising photon production efficiency ( ion ) than standard SSP models, necessitating hotter blackbody temperatures.We explore this by running CLOUDY simulations with input blackbody SEDs, with a setup closely based on that described in Section 4.2.1.To gain insight into the requirements for the two-photon continuum to dominate over the ionizing SED, we initialize hydrogen-only gas at constant gas temperature (measured from [O iii] 4363/5007) and systematically vary the density and blackbody temperature (Figure 9).A weak UV turnover begins to appear at blackbody temperatures  BB ≳ 65, 000 K, which is hotter than a typical O star.A strong UV turnover requires at least  BB ≳ 90, 000 K, much hotter than massive O stars (up to ∼50,000 K; Walborn et al. 2004;Evans et al. 2011;Bressan et al. 2012).
The two-photon continuum is also sensitive to gas density because, at high densities, -changing collisions will suppress the two-photon emission relative to Ly.This is seen in the bottom panel of Figure 9 where we vary gas density for a fixed blackbody temperature of 100, 000 K. The UV turnover is strongly suppressed at  ≳ 10 3 cm −3 .We emphasize that the details of this calculation are sensitive to the chosen column density at which the model is truncated (which in our case is set by the velocity offset of Ly).Furthermore, there exist significant differences in atomic data predictions for the strength of -changing collisions as a function of temperature (Guzmán et al. 2017).
Under the constraints determined above, we adopt an empirical approach to determine the possible underlying spectrum of GS-NDG-9422.We assume that the ionizing spectrum, to first order, can be modelled as a blackbody.To optimize the blackbody model fit to GS-NDG-9422, we begin with the parameters of the best fit BPASS model (Section 4.2.1) and update ionization parameter, gas density, and blackbody temperature to simultaneously reproduce the emission line ratios and continuum shape.We allow for both density and ionization bounded nebulae by adding an additional stopping criterion to reproduce the measured [O iii] 5007/[O ii] 3727 (O32) ratio.This stopping criterion often supersedes the neutral gas column stopping criteria used above.The O32 ratio and the gas temperature are allowed to vary within their observational uncertainties.Optimisation of this model (top panel of Figure 10) demonstrates that a blackbody model with  = 10 5.05 K can provide a good fit to both the shape of the continuum and the majority of lower-ionization state emission line ratios in GS-NDG-9422.
However, where our simple blackbody model fails is that it strongly overpredicts the flux of He ii lines compared to those observed in the spectrum (Figure 11).The weak He ii 1640 and He ii 4686 in GS-NDG-9422 places tight constraints on the ionizing spectrum at  > 4 Rydberg.The opacity of helium in the stellar atmosphere is well known to play an important role in mediating the flux of stellar populations at these energies, suppressing the flux of He +ionising photons ( < 228 Å; e.g.Smith et al. 2002).To explore this, we run models varying the leakage fraction of photons with  > 4 Rydberg from 0% (i.e.no He + -ionizing photons) to 100% (unattenuated blackbody).We conclude that 70% − 75% of the He +ionizing photons emitted by the blackbody must be extinguished to match the observed He ii emission, with a leakage fraction of 0.25 in our best-fit model.The parameters for our optimized model are reported in Table 4.
We note that other sources, including active galactic nuclei (AGN) or high-mass X-ray binaries (HMXBs), can have significant highenergy photon outputs, and have been invoked to explain peculiar emission line ratios seen at high-redshift (Maiolino et al. 2023;Katz et al. 2023b).However, the weak He ii emission disfavours power-law SEDs that extend past the He + -ionizing edge.We exclude HMXBs due to the strong over-prediction of He ii emission in these models (see Appendix B), while the presence of an AGN is discussed above.Finally, we note that spectral fitting of some low-redshift 'high-redshift analogues' has suggested a similar need for significant ionizing contribution from a hot ( > 80, 000 K) blackbody (Olivier et al. 2022), however the key difference in GS-NDG-9422 presented here is the dominance of the two-photon continuum.This necessitates an extremely high  ion that cannot be produced if a substantial fraction of the ionizing photons are emitted by typical OB stars.
We now turn to the question of what sort of stars have sufficiently high surface temperatures to reproduce our best-fit blackbody model with  BB = 10 5.05 K.
Wolf-Rayet stars not only exhibit extremely high surface temperatures, but can also have helium atmospheres that provide the necessary opacity to reduce their He + ionizing output (Crowther 2007).We explore models from grid of the PoWR Wolf-Rayet models (Todt et al. 2015) 3 .Given the gas-phase metallicity measured for GS-NDG-9422, we consider the WNL-H40 grid of the PoWR Wolf-Rayet models (Todt et al. 2015) with  = 0.07  ⊙ , similar to that measured for GS-NDG-9422.We identify model 13-10 as most similar to our optimised blackbody, for which  = 100, 000 K and the luminosity is fixed at  = 10 5.3  ⊙ (Figure 12 top left).We note that this model assumes iron group abundances to be scaled solar, which may not be representative of stars forming at  ∼ 6. Abundances of carbon, nitrogen, and oxygen are assumed to have undergone significant CNO burning (see Todt et al. 2015 for details).
Adopting parameters from the optimised blackbody model but removing the He ii leakage parameter, we replaced the blackbody SED with this theoretical low-metallicity Wolf-Rayet star spectrum.No further optimisation is performed and the resulting spectrum is shown as the dashed orange line in Figure 10, nearly identical to the optimised blackbody model.
Although Wolf-Rayet stars are typically associated with broad emission lines due to their high-velocity stellar winds, the absence of these features in GS-NDG-9422 could simply be the result of the extremely bright nebular component outshining these wind features, as predicted by our photoionization models.Furthermore, in the absence of iron, which is the dominant source of stellar atmospheric opacity, wind speeds can drop below 500 km s −1 , even at solar oxygen abundance (Gräfener & Hamann 2008).Hence, high-redshift Wolf-Rayet dominated galaxies may not show broad He ii 1640 and He ii 4686 lines (Gräfener & Vink 2015).
Stars stripped in binaries can also exhibit the required surface temperatures above 10 5 K (Götberg et al. 2018(Götberg et al. , 2023)).Compared to our best-fit extinguished blackbody, we find that stripped star SEDs (Götberg et al. 2018) produce too few He + ionizing photons, and the He ii emission is underpredicted by these stars (Figure 12 top right).
Massive metal-free Population III stars are predicted to have sufficiently high temperatures to power a two-photon dominated spectrum (Schaerer 2002;Trussler et al. 2023).While the IMF of Population III stars is highly uncertain (Klessen & Glover 2023), comparing extremely top-heavy (Pop.III.1) and moderately top-heavy (Pop.III.2) Yggdrasil models (Zackrisson et al. 2011) with our extinguished blackbody (Figure 12 bottom left) indicates these produce too many He + ionizing photons, while predictions for even less topheavy IMFs begin to fall short of the required  ion .In contrast, some individual Pop.III star models from Larkin et al. (2023) with effective temperatures of ∼ 97, 000 K and masses of 85 M ⊙ to 108 M ⊙ reasonably reproduce the required ionizing SED (Figure 12 bottom right).We strongly emphasize that the measured metallicity of 12 + log 10 (O/H) = 7.59 ± 0.01 suggests the presence of Pop.III stars is highly unlikely, unless we are witnessing the immediate enrichment and illumination of metals produced by primordial stars.Hence, we are not advocating that GS-NDG-9422 hosts a population of primordial stars.Nevertheless, stellar atmosphere models have large uncertainties at low metallicity, and few such models exist in the literature, so we consider these Pop.III star models in our analysis under the assumption that the atmospheres may be representative of hot massive stars at very low metallicity.Repeating the exercise described for our Wolf-Rayet models, this time replacing the blackbody SED with low-metallicity massive star models from Larkin et al. (2023) with effective temperatures between 95,000 K and 99,000 K also results in a very good fit to the spectrum of GS-NDG-9422 (green dashed lines in Figure 10).In summary, we identify three classes of model stellar SED that have the surface temperatures required to reproduce the observed two-photon continuum turnover in GS-NDG-9422 (Figure 10; see also schematic in Figure 13).Wolf-Rayet stars with  = 0.07  ⊙ , equal to the measured nebular metallicity, and hot massive stars with low-metallicity atmospheres provide a remarkably good fit, while existing model SEDs from stripped stars only fall short in having too strong suppression of He + -ionizing flux.While a perfect match to the observed spectrum of GS-NDG-9422 will require fine-tuning of the photoionization model, the fact that the continuum shape and the vast majority of the emission lines can be reproduced under very simple assumptions is encouraging that stars such as those considered in this section may be present in GS-NDG-9422.

Implications for the stellar initial mass function if GS-NDG-9422 hosts a population of hot stars
We now consider what our modelling implies for the mass distribution of the stellar population in GS-NDG-9422.The progenitor masses of Wolf-Rayet stars at the metallicity of GS-NDG-9422 are not well constrained (Massey et al. 2001).Masses of ≥ 37 M ⊙ have been estimated for Wolf-Rayet stars in the Small Magellanic Cloud (SMC) with  eff ≳ 100, 000 K (Hainich et al. 2015), implying even higher progenitor masses.Meanwhile, the low-metallicity massive star models have somewhat higher progenitor masses of ∼ 100 M ⊙ .While our search may not be exhaustive, we conclude from Figure 12 that the nebular-dominated spectrum of GS-NDG-9422 is consistent with ionization powered by low-metallicity massive stars (≳ 50 M ⊙ ), perhaps in the Wolf-Rayet phase.
If we adopt a progenitor mass of 50 − 100 M ⊙ , then, assuming a typical IMF, we expect to form one such star per every .Schematic of how the hot star scenario gives rise to a nebular-dominated spectrum in a simple spherical nebula.The incident stellar spectrum ( eff = 10 5.05 K) is shown in blue.The portion of this transmitted through the nebula is shown in the darker shade, the lighter shade shows the portion that is absorbed to photoionisation of the gas.The red shows the resulting nebular spectrum, with the black showing the sum of nebular and transmitted stellar, which is what an observer will see.For a fixed 1500 Å flux density, hot stars with  eff ≳ 10 5 K emit significantly more ionising photons (light blue) than a typical stellar population, which in turn allows them to power much stronger nebular emission (red), which ultimately outshines their own spectrum at wavelengths longer than Ly (1216 Å).Note that a significant population of old stars can be accommodated within this scenario without affecting the nebular dominance at UV wavelengths, since these old stars primarly contribute flux at longer wavelengths.
∼ 1, 300 − 3, 800 M ⊙ of stellar mass.To explore whether our model is consistent with this, we repeat our photoionization modelling with our best-fit massive star models while simultaneously adding a second component from a BPASS model.We assume the BPASS stellar population has the same metallicity as the gas.In the case of lowmetallicity massive stars, we assume the BPASS model has an age of 3 Myr (the approximate lifetime of massive stars), while for the Wolf-Rayet models, we assume the progenitor stars are lower mass (50 M ⊙ ) and live for slightly longer (5 Myr).We then progressively increase the ionization parameter of the second population until the UV turnover from the two-photon continuum begins to weaken.We calculate the stellar mass of a BPASS SSP that can be present per hot star as where  star = 10 49.36  s −1 for the Wolf-Rayet star model or 10 49.98  s −1 for the low-metallicity massive star model with  eff = 97, 352 K.  SSP = 10 46.73 or 10 46.03  s −1 M −1 ⊙ at 3 Myr and 5 Myr, respectively.In both cases, we find the maximal allowable contribution to be  SSP  star ≤ 6.3%.Therefore, we can place upper limits of the BPASS SSP mass of  SSP ≤ 112 or 137 M ⊙ per one low-metallicity massive star or Wolf-Rayet star, respectively.This is clearly discrepeant from the values of 3, 800 M ⊙ and 1, 300 M ⊙ implied from a typical IMF.In other words, there is a 35× excess in the number of massive stars in the case of our low-metallicity massive star model, or a 9.5× excess for our Wolf-Rayet model.Either scenario implies that the IMF in GS-NDG-9422 is very top-heavy.
We note that the measurement presented here is sensitive to the IMF by way of constraining the ratio between the ionising photon production rate, dominated in our model by stars more massive than ≳ 50 M ⊙ , and the continuum flux at ∼ 1500 Å, which in standard young stellar population models is dominated by stars with M ≳ 5 M ⊙ (e.g.Byler et 2017).It is, therefore, only sensitive to the high-mass end (M ≳ 5 M ⊙ ) of the IMF, and is completely insensitive to the impact of stars with M ≲ 5 M ⊙ .
Taken at face value, this calculation implies that the ratio of ≳ 50 M ⊙ to ∼ 5 − 50 M ⊙ stars in GS-NDG-9422 is of order ∼10× higher than that predicted by a typical IMF.However, the results of this calculation are very sensitive to the chosen underlying SSP, the assumed age, the progenitor masses of Wolf-Rayet stars, and the mass of the low-metallicity massive stars.For SSPs that have intrinsically lower  ion , the excess drops as  SSP appears in the denominator of Equation 1.If the ages at which the Wolf-Rayet stars evolve off the main sequence is longer than what we have assumed here, the excess will also decrease.More generally, we highlight that the origin of such hot stars in high-redshift environments and the modelling of their atmospheres is highly uncertain and motivates detailed characterisation of the effects of binary stripping and stellar wind mass loss, particularly at low metallicity (Götberg et al. 2019;Vink 2022).Thus, while we have demonstrated qualitatively that explaining the observed spectrum without variations to IMF is extremely difficult, a detailed characterisation of the high-mass end of the IMF in GS-NDG-9422 will require more advanced understanding of stellar evolutionary processes in these environments.
Understanding the shape, upper-mass, and lower-mass cutoff of the IMF, and whether these evolve with initial conditions, is critical to the interpretation of nearly all extragalactic observables (Hopkins 2018).Theoretical works have predicted the IMF to get progressively more top-heavy for low-metallicity gas at high pressure (Chon et al. 2021(Chon et al. , 2022;;Sneppen et al. 2022), while others indicate that increased CMB temperature can cause the IMF to become more bottom-light at high-metallicity (Bate 2023).
While many studies of local field stars and young clusters have found no strong evidence for variation in the IMF of these systems (e.g.Bastian et al. 2010 and references therein), a number of lines of evidence have emerged suggesting the IMF may vary in some environments.Top-heavy IMFs have been derived in some local massive young star clusters (Kalari et al. 2018;Schneider et al. 2018), while some models of globular clusters have invoked a topheavy IMF at early times to explain how gas expulsion proceeded from the young cluster (Marks et al. 2012).
Spectral line indices in the spectra of early-type galaxies have been widely demonstrated to show that the IMF is bottom-heavy in these systems (van Dokkum & Conroy 2010;Spiniello et al. 2014;Martín-Navarro et al. 2015;Conroy et al. 2017;Maksymowicz-Maciata et al. 2024).Similar conclusions have been drawn from gravitational lensing studies (Treu et al. 2010;Smith 2020) and dynamical modelling (Cappellari et al. 2012;Poci et al. 2022;Lu et al. 2023), both of which are senstive to the mass-to-light ratios.We note that these measurements are sensitive to the ratios of stars with  ≲ 0.5 ⊙ and  ∼ 0.5 − 1 ⊙ , which our constraint is insensitive to.Furthermore, we note that chemical evolution modelling of massive ellipticals has shown that a time-invariant bottom-heavy IMF cannot explain the observed [Mg/Fe] and [Fe/H] abundance ratios in these systems with these studies instead suggesting that this bottom-heavy phase may have been preceeded by a short-lived top-heavy phase (Vazdekis et al. 1997;Weidner et al. 2013;Jeřábková et al. 2018;De Masi et al. 2019).
Abundances of CNO elements and their isotopes are expected to be highly sensitive to variations in the high-mass end of the IMF (Romano 2022).Variations in the 13 C/ 18 O ratio observed in  ∼ 2−3 dusty starburst galaxies have been interpreted as arising from a topheavy IMF (Zhang et al. 2018), while a sub-solar C/O abundance ratios arising from a top-heavy IMF has been suggested to explain anomalous IR emission line ratios at high-redshift (Katz et al. 2022).Furthermore, a rapid starburst with a top-heavy IMF has been invoked as a possible explanation the enhanced N/O abundance ratio observed in GN-z11 (Bekki & Tsujimoto 2023).
Top-heavy IMFs with low mass-to-light ratios have also been widely invoked as a possible cause of the surprising abundance of UV-bright galaxies observed at high redshift (e.g.Finkelstein et al. 2023;Harikane et al. 2024;Yung et al. 2024).

Total cluster mass and number of hot stars
The hydrogen ionising photon luminosity, , can be approximated from the H flux as: where   is the luminosity distance,  H is the measured H flux, ℎ is Planck's constant,  H is the frequency of the H transition,  esc is the escape fraction of ionising photons,   is the total case B recombination rate, and  eff B is the effective H recombination rate.To calculate a luminosity distance, we adopt a cosmology with  0 = 67.31km s −1 Mpc −1 and Ω m = 0.315 (Planck Collaboration et al. 2016).This gives   = 58.5 Gpc, implying an H luminosity of  H = 7.0 × 10 41 erg s −1 .Using the escape fraction derived from our best fitting blackbody photoionization model of 7.3% and recombination rates evaluated at 20,000 K (Osterbrock & Ferland Nebular dominated galaxies 17 2006), we estimate a hydrogen ionising photon luminosity of  = 1.65 × 10 54 s −1 .This value allows us to estimate the mass of the star clusters that have formed in GS-NDG-9422.The Wolf-Rayet star models presented in Section 6.1 have a fixed luminosity of 105.3  ⊙ , which results in a hydrogen ionising photon luminosity of  WR = 10 49.36 s −1 .However, known Wolf-Rayet stars in the most comparable local environments typically have significantly higher luminosities of ∼ 10 6.2  ⊙ (Shenar et al. 2016), implying that  WR ≈ 10 50.3 s −1 might be a more realistic value.Based on the H luminosity, this would imply ∼ 10, 000 of the  = 0.07  ⊙ Wolf-Rayet stars would be needed to power the spectrum of GS-NDG-9422.In the alternative model, ∼ 17, 000 very metal-poor stars of ∼ 100 M ⊙ with  eff = 97, 000 K would be required.Based on the calculations presented in Section 6.1, this implies a maximum star cluster mass of 10 6.22 M ⊙ for the Wolf-Rayet star model, or 10 6.55 M ⊙ for the lowmetallicity massive star model.We caution that these mass estimates should be treated only as a guide since they are very sensitive to the adopted hot star SED as well as the same assumptions outlined in Section 6.1.We also note that the quoted value includes only the mass of recently formed population contributing significant UV luminosity.A considerable mass of older stellar populations could also be accommodated within these models.

Are there other galaxies like GS-NDG-9422?
While the physics driving the abnormal continuum shape of GS-NDG-9422 remains uncertain, it is important to understand whether GS-NDG-9422 is unique, or whether other objects show similar spectral features.Identifying larger samples of objects that are similar to GS-NDG-9422 will help rule out scenarios that require fine-tuning.
While we have not performed an exhaustive or systematic search, in Figure 14 we highlight two objects that both show Ly and a UV turnover, similar to GS-NDG-9422.The first is the Lynx arc (Fosbury et al. 2003), a gravitationally-lensed, extreme emission system identified at  = 3.35.The emission lines of the Lynx arc differ considerably from GS-NDG-9422 in that He ii is weaker and there is strong N iv] emission which is not observed in GS-NDG-9422.Nevertheless the UV downturn and Ly are remarkably similar.Indeed, analysis in Fosbury et al. (2003) came to similar conclusions as presented here for GS-NDG-9422 -namely, that hot stars ( eff ≳ 80, 000 K) are powering the emission in this galaxy which leads to the visibility of the two-photon continuum.Hence, it is clear that the shape of the continuum in GS-NDG-9422 is not unique among the galaxy population.
We identify a second galaxy, A2744-NDG-ZD4 at  = 7.88, observed as part of the Ultra-deep NIRCam and NIRSpec Observations Before the Epoch of Reionization (UNCOVER) program (ID: 2561, PI: Labbe; Bezanson et al. 2022) which also appears to show Ly escape and a UV turnover (Figure 14 bottom panel).The spectrum of this object, retrieved from Dawn JWST Archive (DJA)4 which had been reduced using the custom-made pipeline MsaExp v.0.6.12 5 (see Heintz et al. 2023 for details), received only 2.3 hours integration.Thus, the continuum is not well detected in the rest-frame optical and we can not determine whether it shows a Balmer jump.Deeper data will be required to conclusively determine whether A2744-NDG-ZD4 is truly of the same nature as the Lynx arc and GS-NDG-9422.

CONCLUSIONS
In this paper, we have presented a discussion of the possible physical origin of the observed continuum and line emission in GS-NDG-9422, a galaxy at  = 5.943.The identification of a spectral discontinuity of 15 ± 0.9 nJy at the location of the Balmer limit is clearly indicative of very strong nebular emission, as are the high equivalent widths observed for [O iii] 5007, H, and H.We measure high nebular temperatures ((1.83 ± 0.15) × 10 4 K), low densities (  ≲ 10 3 cm −3 ), and low metallicity (12 + log(O/H) = 7.59 ± 0.01; ∼8 % (O/H) ⊙ ) in this system.We show that, considering only the emission line ratios, this system is highly consistent with emission powered by a typical population of young, metal-poor stars, with no strong evidence for the presence of AGN activity.
What is more difficult to explain is the observed shape of the continuum.In particular, the strong UV turnover observed at  rest ≈ 1430 Å.The shape of this feature is a remarkably consistent with the hydrogen two-photon continuum, but could also be explained by a very high column density of neutral hydrogen.We consider several scenarios for how each of these effects could self-consistently give rise to the observed spectrum of GS-NDG-9422, all of which are summarized in Figure 15.
The primary challenge faced by models explaining this feature with a DLA is that the magnitude of the turnover implies unprecedented column densities of  HI ≳ 10 23 cm −2 , higher than any known DLA.Furthermore, the observation of narrow, high EW Ly emission with a high escape fraction (∼ 27%) requires some mechanism by which this emission can reach the observer.This can be reconciled if Ly arises from a reflection spectrum (Scenario #1), although known examples of this are orders of magnitude larger in physical scale than GS-NDG-9422.Alternatively, the Ly might be allowed to reach the observer if the DLA has only partial coverage (Scenario #2) or

5.
6. is at lower redshift (Scenario #3); however, each of these scenarios imply even higher neutral gas column densities which means their plausibility is questionable.Why such high gas columns would not be fully molecular and why at this moderate metallicity there are no signatures of dust attenuation with so much neutral gas remain open questions in the DLA scenario.For these reasons, we do not favour the DLA solution.

3.
We demonstrated in Figure 7 that modelling the continuum of GS-NDG-9422 with strong two-photon, free-bound and free-free continuum plus a subdominant young stellar and old stellar component provides a very good fit to the spectrum, with the predicted nebular temperature and equivalent widths of Balmer emission highly consistent with what is measured from the emission lines.The main challenge for explaining the UV turnover with two-photon continuum emission is therefore simply: how can the system power such strong nebular emission without the rest-UV continuum becoming dominated by the spectrum of the ionising source?
One possibility is that no flux from the ionising source is observed, either due to a spatial offset whereby the ionising source falls outside the slit (Scenario #4) or because the ionising source has recently turned off and we are observing emission from its relic (Scenario #5).The former is disfavoured by the compact morphology of GS-NDG-9422, which is well-centred in the slit (Figure 1).The predicted recombination timescales in Scenario #5 (∼500 -5,000 yr for H + and ∼30-300 yr for He ++ ) are incredibly short meaning we would need to have caught this object at a very specific moment in time, although this could perhaps be explained with an AGN fading on a timescale of ≲ 10 4 yr (e.g.French et al. 2023).We consider this model to be possible, but the identification of larger samples of similar may rule this scenario out.
The final scenario we consider is that the ionising source is indeed present in the slit, but its ionising photon production efficiency is so large that it powers sufficiently strong nebular emission that the two-photon continuum provides the dominant contribution to the rest-UV continuum (Scenario #6 and Figure 13).This is a generic prediction of all photoionization models that include a significant population of hot stars with  eff ≳ 10 5 K (e.g.Schaerer 2002;Raiter et al. 2010;Zackrisson et al. 2011;Inoue 2011;Trussler et al. 2023).We identify several model stellar SEDs that can attain these temperatures, showing that models with ionising spectra dominated by ∼7 %  ⊙ Wolf-Rayet stars and/or very low-metallicity massive stars are able to reproduce the puculiar shape of the observed continuum.The critical aspect to this scenario is that it requires a top-heavy IMF, with the implied ratio of ≳50 M ⊙ to ∼5-50 M ⊙ stars significantly enhanced relative to typical stellar populations.We note that the precise magnitude of this excess is highly dependent on the assumptions around the evolution and atmospheric properties of these hot, massive, low-metallicity stars, and this work motivates a pursuing improved models of these stars.Nonetheless, we note that although such an IMF is unprecedented in nearby star clusters, it is broadly consistent with other suggestions of enhanced massive star formation in extreme systems at early times (see discussion in Section 6.1).
All three classes of solutions we have proposed to explain the spectrum of GS-NDG-9422 pose interesting physical questions about high-redshift galaxy formation.Either we have discovered the highest DLA column density for a star-forming galaxy, we have witnessed the flickering of a powerful ionising source on timescales much shorter than typical QSO lifetimes, or the ionizing source represents an exotic population of stars with high surface temperature.The similarities with the Lynx Arc at  = 3.35 and A2744-NDG-ZD4 at  = 7.88, indicate that GS-NDG-9422 may not be unique.This motivates a more systematic search of such objects to unravel the origin of their unique spectra.

Nebular dominated galaxies 21 APPENDIX A: EMPIRICAL DLA COLUMN DENSITY MEASUREMENTS
In Section 3, we optimized a CLOUDY photoionization model in order to reproduce the emission lines of GS-NDG-9422, which then provided an estimate of the continuum.We then applied various DLA columns with leakage to this model in order to reproduce the UV turnover and the escape of Ly.However, the exact column density needed is degenerate with the underlying shape of the continuum.In this section, we repeat the experiment trying only to match the shape of the continuum.
Following the approach in Section 4, we run an MCMC to fit the continuum of GS-NDG-9422.In this case, we do not allow the twophoton continuum to have a free normalization (and hence we remove this parameter), and we replace it with a variable DLA column density and leakage.The model thus has eight free parameters (compared to the seven used in Section 4).
In Figure A1 we show the fit corresponding to the 50th percentile model in the posterior distributions of each parameter as well as the marginalized posterior distributions on H I column density and leakage.Indeed we find values close to 10 23.3 cm −2 , which is slightly lower than reported in Section 3, but very consistent within 1.In either case, the DLA column required is exceptionally high compared to known DLAs.
Interestingly, despite the fact that Ly was not included in the fitting, the MCMC tightly constrains the amount of leakage to be ∼ 35%.This is because the shape of the transmission curve for the DLA is inconsistent with the shape of the UV downturn in GS-NDG-9422, while in contrast, it is much more consistent with that of the two-photon continuum.Nevertheless, it is encouraging to see that the MCMC requires the same leakage as discussed in Section 3.
While this DLA model is able to reproduce the shape of the continuum, we find two inconsistencies.First, the 50th percentile temperature of 22,369 K is much higher than that measured by the [O iii] 4363 auroral line.As stated above, these two temperatures do not necessarily have to agree, but are typically observed to be close (Guseva et al. 2006(Guseva et al. , 2007)).Moreover, in this model, the predicted H EW is 767 Å±40 Å which is very inconsistent with what is observed.This is in stark contrast to the model where the UV is dominated by two-photon emission.Hence the results in this section confirm and strengthen our conclusions that the spectrum of GS-NDG-9422 is more consistent with being primarily nebular in origin.

APPENDIX B: MODELS WITH XRBS
During the course of the modelling presented in Section 5, we also considered the possibility that the emission is powered by X-ray binaries.Following Senchyna et al. (2020) and Katz et al. (2023b), we model the X-ray sources using a modified blackbody spectrum (Mitsuda et al. 1984).We assume black hole masses in the range 6 − 25 M ⊙ and disk radii between 10 3 − 10 4 gravitational radii.In order for a UV turnover to appear, the ionizing photon output from the XRBs must dominate over the stellar population so that the nebular continuum can outshine the stars, implying a very high ratio of X-ray luminosity to star formation rate.In Figure B1, we show our XRB model for 25 M ⊙ black holes optimised to reproduce the continuum shape of GS-NDG-9422.This model significantly over-predicts the strength of the He ii 1640 and 4686 lines.Varying the black hole masses does not resolve this issue.We therefore conclude that XRBs are not the dominant ionizing source in GS-NDG-9422.This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
Figure 1.Top: 2D Prism/CLEAR spectrum of GS-NDG-9422.Middle: 1D Prism/CLEAR spectrum of GS-NDG-9422 shown in   (upper middle) and   (lower middle).Blue squares show 0.15" radius aperture photometry from JWST/NIRCam medium-and broad-band imaging, while orange squares show the predicted photometry obtained by convolving the Prism spectrum with the NIRCam filter transmission profiles.Bottom left: Three-colour image (090, 200, 444) of GS-NDG-9422 showing the positioning of the three NIRSpec micro-shutters across the three nod positions.Bottom centre: Zoom-in of the region surrounding [O iii] 5007 in the G395M spectrum.Bottom right: Zoom-in of the region surrounding H in the G395H spectrum.
5) are shown in dark blue, while all other star-forming models are shown in light blue.Model predictions for active galactic nuclei from Feltre et al. (2016) are shown in red ( = 0.001) and pink (all other ).The non-detection of [S ii] 6716, 6731 places GS-NDG-9422 firmly below the Kewley et al. (2001) maximum-starburst limit of the classical [S ii]-VO87 diagram (Veilleux & Osterbrock 1987; top

Figure 2 .
Figure 2. GS-NDG-9422 plotted onto several line ratio diagnostic diagrams.Light blue points show model predictions for star-forming regions from Gutkin et al. (2016), while pink points show AGN model predictions from Feltre et al. (2016) across a large range of metallicities.Navy (star-forming) and red (AGN) points show the subsets of these model grids that have  = 0.001 ( ≈ 0.07 ⊙ ; 12 + log(O/H) ≈ 7.54) which is the closest grid value to that measured for GS-NDG-9422 (see Section 3.5).Top left: The strong upper limit on [S ii] 6716, 6731 positions GS-NDG-9422 well below the theoretical maximum starburst limit from Kewley et al. (2001) (black line).Top right: The weak detection of He ii 4686 is also below the maximum starburst limit from Shirazi & Brinchmann (2012) (black line).Strong He ii-emitting star-forming galaxies from Shirazi & Brinchmann (2012) and Senchyna et al. (2017) are also shown as orange and green diamonds respectively.We also show measurements from Hanny's Voorwerp, suggested to be a quasar light echo(Lintott et al. 2009; purple line-connected points), which we discuss in Section 5.1.Points of increasing size indicate larger offset from the quasar ranging from 13 to 31 kpc.Bottom left: GS-NDG-9422 lies beyond the AGN region defined byMingozzi et al. (2024) (black lines), and is more coincident with the star-forming models than AGN models.He ii-emitting star-forming galaxies at  ∼ 2.5 − 4 fromSaxena et al. (2020) are shown for comparison as orange squares.Bottom right: GS-NDG-9422 lies at the tip of the parameter space covered by AGN models, while it is well within the bounds of that covered by star-forming models.

Figure 3 .
Figure3.Observed Balmer decrements compared to theoretical values at different temperatures and densities.Green squares give the measured ratios from the Prism, while diamonds give ratios measured from gratings, where available.Solid black lines give the theoretical ratio for case B recombination at   = 10 4 K and   = 100 cm −3 .Coloured points show theoretical ratios for a range of temperatures (10 4 ≤   ≤ 2.5 × 10 4 K; colour) and densities (  = 10, 100, 1000 cm −3 ; marker size).The measured H/H is somewhat lower than the theoretical value, but we note that H may be underestimated in the prism spectrum (Table2).

Figure 4 .
Figure 4. Prism spectrum of GS-NDG-9422 (black) compared with the best fit models using a standard SSP (magenta) accounting for both the  = 6 IGM opacity (dashed yellow) and a DLA with column density of 10 23.1 cm −2 (dotted cyan).

Figure 5 .
Figure 5. Top: Zoom-in on UV turnover for models plotted in Figure 4. Middle: Fitting to GS-NDG-9422 using a DLA model with a 70 % covering fraction allows Ly escape, but implies an extremely high column density.Bottom: Comparison of implied column density with known DLAs.

Figure 6 .
Figure6.Schematic of the main nebular continuum components across the rest UV to optical range.Free-bound emission can give rise to a spectral discontinuity at the Balmer limit ( rest ≈ 3645 Å).Two-photon continuum has a fixed shape with a turnover at  rest ≈ 1430 Å, but is usually subdominant compared to continuum emission of the ionising source.Free-free emission is usually sub-dominant at these wavelengths.The relative contribution of two-photon continuum is highly dependent on nebular conditions.

Figure 7 .
Figure 7. Top: Best fit continuum model from fitting described in Section 4.2.2.Continuum points using in the fitting are shown in red, while the full spectrum is shown in green.Coloured lines show the continuum components arising from young stars (pink), old stars (light blue), free-bound + free-free nebular (brown), and two-photon (yellow) in the best fit model.The blue line shows the overall best fit.Bottom: Posterior distribution on gas temperature predicted by the continuum fitting.The median and 1 values (red lines) are in good agreement with the nebular temperature measured from the [O iii] 4363/5007 ratio (Table3; blue).

Figure 8 .
Figure8.Top: Posterior distribution on H equivalent width predicted by the continuum fitting.The median and 1 values (red lines) are consistent with that measured from the Prism spectrum (blue).Bottom: Same as top panel, but for H.Here, the median and 1 values (red lines) are somewhat higher than that measured from the Prism spectrum (blue), but in good agreement with the equivalent width implied by the medium-band imaging.

Figure 9 .
Figure 9. Top: Normalized spectrum of hydrogen-only gas with  = 10 3 cm −3 irradiated by blackbodies of different temperatures (as given in the colour bar).The magenta line represents a blackbody temperature of 100,000 K. Bottom: Normalized spectrum of hydrogen-only gas irradiated by a blackbody with  = 10 5 K at different gas densities (as given in the colour bar).The magenta line represents a density of 10 3 cm −3 .

Figure 10 .
Figure10.Top: Spectrum of GS-NDG-9422 (black) compared with the best fit blackbody model (magenta) with  BB = 10 5.05 K and a He ii leakage fraction of 0.22 (Table4; Section 5.2).Middle: Spectrum of GS-NDG-9422 (black) compared with photoionsiation models powered by various individual massive Pop.III stars with different effective temperatures.Bottom Spectrum of GS-NDG-9422 (black) compared with the photoionisation model powered by a model Wolf-Rayet SED with the most similar spectrum to the best-fit blackbody SED.All models can also reproduce the observed spectrum well (dotted lines).See Section 5.2 for details.

Figure 11 .
Figure11.Zoom in on the region around He ii 4686 demonstrating the need to reduce the He ii ionizing flux in the blackbody models.The spectrum of GS-NDG-9422 is shown in black.The magenta line shows the best-fit model from Figure10.Green to blue lines show decrease He ii leakage fraction applied to the blackbody, modelling the effect of a He-rich atmosphere.Models with a high leakage fraction clearly overpredict the He ii 4686 line, although we note this feature is blended with nearby [Ar iv] lines.We note a slight flux excess blue-ward of He ii 4686 which, from inspection of the 2D spectrum, appears to be driven by a noisy pixel.This excess flux is still below the level of the excess emission implied by the 100 % He ii leakage model shown in green.

Figure 12 .
Figure 12.Best fit blackbody SED model with He-rich atmosphere (black solid) compared to hot star SEDs from the literature.The dotted black line shows the blackbody spectrum prior to extinguishing the He ii ionizing radiation.The top left panel shows select  = 0.07  ⊙ Wolf-Rayet stars from Todt et al. (2015), selected to match the measured nebular metallicity.These show remarkable resemblance to our model blackbody.The top right panel shows stripped star models from Götberg et al. (2018).These attain sufficiently high surface temperatures, but the atmospheres in existing models exhibit too much suppression of the He ii ionising continuum.The bottom left panel shows Yggdrasil Pop.III models (Zackrisson et al. 2011) with two different assumed Pop.III IMFs (very top heavy, Pop.III.1; blue) and (moderately top heavy, Pop.III.2; magenta).These models strongly overpredict the He ii-ionising flux.The bottom right panel shows individual massive Pop.III stars from Larkin et al. (2023) with different effective temperatures as shown in the colour bar.Some of these models exhibit the required SED properties to reproduce the observed spectrum of GS-NDG-9422 (Figure10), although we note the discrepancy between these zero-metallicity models and the measured nebular metallicity.Hence, we are not claiming the presence of Pop.III stars in GS-NDG-9422.

Figure 14 .
Figure14.Zoom in on the rest-frame UV region for GS-NDG-9422 (black) compared with the spectra of the  = 3.35 Lynx arc (top) and A2744-NDG-ZD4 at  = 7.88 (bottom).We have normalised the spectra and shifted them all to the rest frame.
N HI ≳10 23 cm -2 blocks the continuum source leading to the UV downturn, but Lya is reflected from offset clouds A DLA with N HI ≫10 23 cm -2 partially blocks the continuum source leading to both the UV downturn and Lya escape A DLA with N HI ≫10 23 cm -2 is present at lower redshift and fully blocks the source but some Lya passes through as it is shifted out of resonance The NIRSpec slit only covers the nebula and not the stellar continuum.The UV downturn is caused by twophoton emission An AGN flickers off in ≲ 1,000 years so that no continuum remains.The UV downturn is caused by twophoton emission Sources with T eff ≳ 10 5 K illuminate a relatively low column density cloud.Lya escapes with minimal scattering and the rest-frame UV and optical are dominated by nebular continuum as the stellar emission is shifted to shorter wavelengths 1.

Figure 15 .
Figure 15.Summary of the three classes of scenarios that can explain the spectrum on GS-NDG-9422.For each scenario, we show various cartoons for how the physics may manifest and give rise to a steep UV slope, a high Ly escape fraction, and the UV turnover redward of Ly.

Figure A1 .
Figure A1.Top: Best-fit continuum model from fitting described in Appendix A. Young stellar component + DLA is shown in pink, nebular component with a fixed low two-photon continuum contribution is shown in brown, while the old stellar component is shown in light blue.Middle: Resulting posterior distribution of neutral hydrogen column density.Bottom: Resulting posterior distribution of leakage fraction.

MFigure B1 .
Figure B1.Spectrum of GS-NDG-9422 (black) compared with photoionization models that include high-mass X-ray binaries with a black hole mass of 25 M ⊙ (magenta).The solid and dashed magenta lines indicate models with different accretion disk radii.

Table 2 .
Rest-frame equivalent widths of prominent lines measured from the Prism/CLEAR spectrum, and line widths measured from  ∼ 2700 G395H/F290LP grating observations.We also reported the equivalent width of H measured from 460 and 430 photometry.

Table 3 .
Physical properties and chemical abundances derived for GS-NDG-9422.

Table 4 .
Parameters for the optimized blackbody model.