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

ABSTRACT

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 }~\sim 27~{{\ \rm per\ cent}}$|⁠) is also consistent with strong nebular two-photon continuum emission. Modelling the UV turnover with a DLA requires extreme column densities (⁠|$N_{\rm HI}\,\,\gt\,\, 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 ionizing source has ‘turned off’ more recently than the recombination time-scale (⁠|$\sim$|1000 yr), hence we may be catching the object at a very specific time. Alternatively, hot stars with |$T_{\rm eff}\gtrsim 10^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 |$\gtrsim 50~{\rm M}_\odot$| stars relative to |$\sim 5\!-\!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.

galaxies: ISM, galaxies: starburst, galaxies: star formation

1 INTRODUCTION

The unprecedented sensitivity of JWST in the near-infrared has revolutionized 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 |$z\gtrsim 5$|⁠, providing new insights in to the conditions of the interstellar media (ISM) of these galaxies (Mascia et al. 2023; Matthee et al. 2023; Sanders et al. 2023; Cameron et al. 2023b; 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. 2024; Umeda et al. 2023), assembly of massive galaxies in the early Universe (Carnall et al. 2023; Glazebrook et al. 2024), and have also provided the many of the highest spectroscopic redshift confirmations to date (Arrabal Haro et al. 2023; Curtis-Lake et al. 2023). Spectral 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. 2024). 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 |$\lambda _{\rm 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 (⁠|$\xi _{\text{i}on}$|⁠). They have been detected in some highly star-forming galaxies at low-redshift (Peimbert & Costero 1969; Guseva, Izotov & Thuan 2006; Guseva et al. 2007). Many numerical simulations predict Balmer jumps to be common at high redshift (e.g. Wilkins et al. 2024; Katz et al. 2023a), 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 |$\xi _{\rm 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 |$\sim 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 |$\lambda _{\rm rest}\approx 1430$| Å (Schaerer 2002; Raiter, Schaerer & Fosbury 2010; Zackrisson et al. 2011; Trussler et al. 2023). Observations of the Lynx arc, a strongly lensed extreme emission system at |$z=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 |$T_{\rm eff}\sim 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). None the less, the identification of these objects, should they exist, offers powerful insight into the properties of systems with extreme ionizing spectra.

In this paper, we identify one object, JADES-GS+53.12175-27.79763 (GS-NDG-9422 hereafter), at |$z=5.943$|⁠, that exhibits high equivalent-width emission lines as well as a strong spectral discontinuity near |$\lambda _{\rm 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 |$\alpha$| absorption (DLA) system (e.g. Heintz et al. 2024), 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-line-based 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 ionizing source. Section 6 outlines the broader implications of the scenarios presented. We then summarize our findings in Section 7. Throughout this paper, we adopt the cosmology of Planck Collaboration XIII (2016) with |$H_0=67.31$| km s|$^{-1}$| Mpc|$^{-1}$| and |$\Omega _{\rm m}=0.315$|⁠, which gives a luminosity distance to GS-NDG-9422 of |$D_l=58.5$| Gpc.

2 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 h integration in Prism/CLEAR and 7 h 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 (Bunker et al. 2023; Eisenstein 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 Fig. 1, we show photometry from the publicly released photometric catalogues (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 arcsec 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 |$\alpha$| emission, |$F444W$| and |$F460M$|⁠. In each case, the flux from imaging is 14 and 17 per cent higher, respectively, suggesting the flux of H |$\alpha$| may be marginally underestimated in the Prism/CLEAR flux calibration. We note that the |$F460M$| flux changes by |$\lt 2$| % when measured within apertures with radii of 0.10 or 0.25 arcsec, suggesting that this offset is not driven by spatial variations in H |$\alpha$|⁠.

Figure 1.

Top: 2D Prism/CLEAR spectrum of GS-NDG-9422. Middle: 1D Prism/CLEAR spectrum of GS-NDG-9422 shown in |$f_\nu$| (upper middle) and |$f_\lambda$| (lower middle). Open blue squares show 0.15 arcsec radius aperture photometry from JWST/NIRCam medium- and broad-band imaging, while filled orange squares show the predicted photometry obtained by convolving the Prism spectrum with the NIRCam filter transmission profiles. Bottom left: Three-colour image (⁠|$F090W$|⁠, |$F200W$|⁠, |$F444W$|⁠) 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] |$\lambda$|5007 in the G395M spectrum. Bottom right: Zoom-in of the region surrounding H |$\alpha$| in the G395H spectrum.

2.1 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.

Table 1.

Emission line flux measurements for GS-NDG-9422 across each observed spectrum. Fluxes are reported in units of |$\times 10^{-19}$| erg s|$^{-1}$| cm|$^{-2}$|⁠.

	Prism/CLEAR	G140M	G235M	G395M	G395H
Ly |$\alpha$|	|$110.1\pm 2.8$|	|$101.4\pm 4.3$|	...	...	...
N v|$\lambda$|1239		|$\lt 5.4$|	...	...	...
N v|$\lambda$|1243		|$\lt 4.9$|	...	...	...
N iv] |$\lambda$|1483		|$\lt 2.4$|	...	...	...
N iv] |$\lambda$|1486		|$\lt 3.1$|	...	...	...
C iv|$\lambda$|1548		|$21.5\pm 1.1^\dagger$|	...	...	...
C iv|$\lambda$|1550		|$17.6\pm 1.1^\dagger$|	...	...	...
C iv|$\lambda \lambda$|1549	|$37.8\pm 0.9$|
He ii|$\lambda$|1640		|$5.1\pm 0.9$|	...	...	...
O iii] |$\lambda$|1660		|$\lt 2.5^\dagger$|	...	...	...
O iii] |$\lambda$|1666		|$4.9\pm 1.0^\dagger$|	...	...	...
He ii + O iii]	|$15.1\pm 1.1$|
N iii] |$\lambda \lambda$|1750		|$\lt 1.9$|	...	...	...
C iii] |$\lambda \lambda$|1909	|$13.6\pm 0.9$|	|$8.0\pm 1.0$|	...	...	...
[O ii] |$\lambda \lambda$|3727	|$2.2\pm 0.2$|	...	|$3.0\pm 0.6$|	...	...
H11	...	...	|$\lt 1.0$|	...	...
H10	...	...	|$\lt 0.8$|	...	...
H9	...	...	|$1.2\pm 0.3$|	...	...
[Ne iii] |$\lambda$|3869		...	|$6.2\pm 0.5$|	...	...
He i|$\lambda$|3889		...	|$1.6\pm 0.4$|	...	...
[Ne iii] + He i	|$7.6\pm 0.4$|
H |$\delta$|	|$4.3\pm 0.5$|	...	...	...	...
H |$\gamma$|	|$8.3\pm 0.4^\dagger$|	...	|$7.5\pm 0.5$|	|$8.1\pm 0.5$|	|$5.9\pm 0.6$|
[O iii] |$\lambda$|4363	|$2.8\pm 0.4^\dagger$|	...	|$2.6\pm 0.7$|	|$3.4\pm 0.5$|	|$2.8\pm 0.5$|
He i|$\lambda$|4471	|$0.8\pm 0.1$|	...	|$\lt 3.1$|	|$\lt 14.6$|	|$2.1\pm 0.5$|
He ii|$\lambda$|4686	|$1.1\pm 0.2^\dagger$|	...	...	|$\lt 1.5$|	|$\lt 2.4$|
[Ar iv] |$\lambda 4711^\ddagger$|	|$0.8\pm 0.2^\dagger$|	...	...	|$1.5\pm 0.4$|	|$\lt 1.3$|
[Ar iv] |$\lambda$|4740	|$0.4\pm 0.2^\dagger$|	...	...	|$\lt 0.8$|	|$\lt 1.2$|
H |$\beta$|	|$17.2\pm 0.3$|	...	...	|$19.4\pm 0.5$|	|$19.6\pm 0.5$|
[O iii] |$\lambda$|4959	|$32.3\pm 0.4$|	...	...	|$36.3\pm 0.7$|	|$36.4\pm 0.6$|
[O iii] |$\lambda$|5007	|$97.1\pm 0.8$|	...	...	|$102.2\pm 1.1$|	|$98.7\pm 1.5$|
He i|$\lambda$|5875	|$2.2\pm 0.1$|	...	...	|$1.3\pm 0.4$|	...
[O i] |$\lambda$|6300	|$\lt 0.5$|	...	...	|$\lt 1.7$|	|$\lt 1.3$|
H |$\alpha$|	|$45.5\pm 0.6$|	...	...	...	|$52.5\pm 0.9$|
[N ii] |$\lambda$|6583	...	...	...	...	|$\lt 1.4$|
[S ii] |$\lambda$|6716	|$\lt 0.5$|	...	...	...	|$\lt 1.6$|
[S ii] |$\lambda$|6731	|$\lt 0.5$|	...	...	...	|$\lt 1.7$|
He i|$\lambda$|7065	|$1.2\pm 0.1$|	...	...	|$\lt 2.9$|	|$\lt 2.1$|

Notes.|$^\dagger$| Fit simultaneously with neighbouring line and fluxes reported separately.

|$^\ddagger$| Blended with He i|$\lambda$|4713.

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 |$\lt 200$| km s|$^{-1}$| (Table 2).

Table 2.

Rest-frame equivalent widths of prominent lines measured from the Prism/CLEAR spectrum, and line widths measured from |$R\sim 2700$| G395H/F290LP grating observations. We also reported the equivalent width of H |$\alpha$| measured from |$F460M$| and |$F430M$| photometry.

Line	EW|$_0$| (Å)	EW|$_0$| (Å)	FWHM (km s|$^{-1}$|⁠)
	Prism/CLEAR	Imaging	G395H/F290LP
Ly-|$\alpha$|	|$328\pm 25$|
H|$\beta$|	|$383\pm 29$|		|$184 \pm 5$|
[O iii] |$\lambda$|5007	|$2275\pm 159$|		|$177 \pm 3$|
H|$\alpha$|	|$1784\pm 340$|	|$2195\pm 400$|	|$161 \pm 3$|

Fluxes derived from different observations are mildly discrepant in some cases. Notably, H |$\beta$|⁠, [O iii] |$\lambda$|4959, [O iii] |$\lambda$|5007 and H |$\alpha$| 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 Fig. 1 between the Prism spectrum and NIRCam photometry, with the possible exception of H |$\alpha$|⁠. We note that in low-resolution data, the continuum level for some emission lines can be difficult to determine, especially for Ly |$\alpha$|⁠, He ii  + O iii] and [O ii] |$\lambda \lambda$|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 |$\alpha$| flux between the imaging and spectroscopy, we also derive EW|$_0$|(H |$\alpha$|⁠) directly from the imaging. Comparing the measured flux in |$F460M$| (⁠|$196\pm 5$| nJy), shown in Fig. 1 to be clearly elevated due to contamination from H |$\alpha$|⁠, with |$F430M$| (⁠|$25\pm 4$| nJy), which captures clean continuum between He i|$\lambda$|5875 and H |$\alpha$|⁠, we find the imaging implies EW|$_0$|(H |$\alpha) = 2195\pm 400$| Å. Within the uncertainty, the measured difference between EW|$_0$|(H |$\alpha)_{\rm Prism}$| and EW|$_0$|(H |$\alpha)_{\rm Imaging}$|⁠, Table 2) is consistent with the discrepancy noted above between the observed |$F460M$| flux, and the predicted |$F460M$| 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.

3 DERIVATION OF PHYSICAL CONDITIONS

We now explore the physical conditions of the gas in GS-NDG-9422 and the basic properties of the ionizing source, based on the measured emission line fluxes. Throughout this section, we make use of pyneb (Luridiana, Morisset & Shaw 2015) using atomic data from chianti (version 10.0.2; Dere et al. 1997; Del Zanna et al. 2021).

3.1 Diagnostic diagrams

To explore the properties of the ionizing source, in Fig. 2, we look at nebular line ratio diagnostic diagrams, comparing measured line ratios from GS-NDG-9422 with photoionization model predictions for star formation and active galactic nuclei. For star forming models, we adopt the predictions of Gutkin, Charlot & Bruzual (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 |$z\sim 6$| stellar populations, while nebular abundance patterns allow for some variation in C/O and N/O ratios. Star forming models with |$Z=0.001$| (the grid value closest to the gas-phase oxygen abundance derived for GS-NDG-9422; see Section 3.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, Charlot & Gutkin (2016) are shown in red (⁠|$Z=0.001$|⁠) and pink (all other Z).

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 |$Z=0.001$| (⁠|$Z\approx 0.07 Z_\odot$|⁠; |$12+\log ({\rm O}/{\rm H})\approx 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] |$\lambda \lambda$|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|$\lambda$|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 |$z\sim 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.

The non-detection of [S ii] |$\lambda$||$\lambda$|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 left panel of Fig. 2) and in a region that is difficult to reconcile with emission powered by an AGN. Non-detections of [N ii] |$\lambda$|6583 and [O i] |$\lambda$|6300 paint a similar picture in the [N ii]-BPT diagram (Baldwin, Phillips & Terlevich 1981) and the [O i]-VO87 diagram, although these are not shown.

The top right panel of Fig. 2 shows that the He ii|$\lambda$|4686/H |$\beta$| ratio of GS-NDG-9422 exceeds that predicted by any of the star-forming models in Gutkin et al. (2016), especially those at |$Z=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, Fragos & Izotov 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|$\lambda$|4686/H |$\beta$| 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 stars, and/or a top-heavy stellar initial mass function (Kehrig et al. 2015, 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|$\lambda$|4686/H |$\beta$| ratio. However, we consider the measured value of this line ratio to be inconclusive.

In the bottom left panel of Fig. 2, we see that the O iii] |$\lambda \lambda$|1660, 1666/He ii|$\lambda$|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 |$z\sim 2-4$| He ii-selected 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 photoionization models struggle to reproduce a number of the observed emission line ratios, especially those involving low-ionization 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 ionizing source in GS-NDG-9422 is most likely of stellar origin.

3.2 Balmer decrements and dust extinction

Balmer decrements H |$\delta$|/H |$\beta$| and H |$\gamma$|/H |$\beta$| from the Prism and H9/H |$\beta$| from the grating are consistent with Case B values, indicating that there is no significant dust reddening in GS-NDG-9422 (Fig. 3). Measured H |$\alpha$|/H |$\beta$| ratios are lower than theoretically predicted for Case B at |$T=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 |$T_e\gtrsim 2\times 10^4$| K, possibly indicative of a very hot nebula. As noted in Section 2.1, H |$\alpha$| may be underestimated in the Prism which could lead to the slightly lower observed ratio. H |$\gamma$|/H |$\beta$| 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 |$\gamma$|/H |$\beta =0.3\pm 0.03$|⁠. However, we note that the continuum is undetected in the high-resolution grating, which contributes uncertainty to the measured ratio.

Figure 3.

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 |$T_e=10^4$| K and |$n_e=100$| cm|$^{-3}$|⁠. Coloured points show theoretical ratios for a range of temperatures (⁠|$10^4\le T_e \le 2.5 \times 10^4$| K; colour) and densities (⁠|$n_e=10$|⁠, 100, 1000 cm|$^{-3}$|⁠; marker size). The measured H |$\alpha$|/H |$\beta$| is somewhat lower than the theoretical value, but we note that H |$\alpha$| may be underestimated in the prism spectrum (Table 2).

Adopting a luminosity distance of |$D_l=58.5$| Gpc, we derive |$L_{\rm H\alpha }=1.86\times 10^{42}$| erg s|$^{-1}$|⁠. Assuming no dust, a metallicity of |$Z=0.1 Z_\odot$|⁠, and a typical IMF, this would correspond to a star formation rate of |$\sim 3.2$| M|$_\odot$| yr|$^{-1}$| (Eldridge et al. 2017). Under the assumption of no dust, Case B recombination, |$n_e=100$| cm|$^{-3}$| and |$T_e=1.8\times 10^4$| K, we derive a Ly |$\alpha$| escape fraction of |$f_{\rm esc,Ly\alpha }=0.29\pm 0.01$| from the Ly |$\alpha$|/H |$\alpha$| ratio, or |$f_{\rm esc,Ly\alpha }=0.27\pm 0.01$|⁠, measured from the Ly |$\alpha$|/H |$\beta$| ratio, the latter of which may be more reliable owing to the H |$\alpha$| flux uncertainty discussed above.

3.3 Electron temperature

The temperature-sensitive [O iii] |$\lambda$|4363/|$\lambda$|5007 ratio can be measured from each of the Prism, G395M and G395H observations, yielding three consistent, independent temperature measurements (⁠|$T_e = 1.83\pm 0.15$|⁠, |$1.99\pm 0.18$|⁠, and |$1.81\pm 0.18\times 10^4$| K, respectively). The temperature derived from the medium-resolution O iii] |$\lambda$|1666/[O iii] |$\lambda$|5007 ratio is somewhat lower (⁠|$T_e=1.70^{+0.05}_{-0.06}\times 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 (⁠|$T_e=1.83\times 10^4$| K) implies |$\lambda \lambda$|1660,1666/|$\lambda 5007=0.08$|⁠, which gives |$f_{\lambda \lambda 1660,1666}=8.2\pm 0.1\times 10^{-19}$| erg s|$^{-1}$| cm|$^{-2}$| based on the [O iii] |$\lambda$|5007 Prism flux, suggesting that O iii] contributes |$\sim$|55 per cent 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|$\lambda$|1640 flux (⁠|$6.9\pm 1.1\times 10^{-19}$| erg s|$^{-1}$| cm|$^{-2}$|⁠) gives a He ii|$\lambda$|1640/|$\lambda$|4686 ratio of |$6.3\pm 1.5$| which is consistent with the theoretical value of 7.19 assuming Case B at |$T_e = 1.8\times 10^4$| K and |$n_e=100$| cm|$^{-3}$|⁠. Note that, this further supports the conclusion of a lack of dust in this system since the He ii|$\lambda$|1640 would be subject to extremely high extinction, relative to the |$\lambda$|4686 line. In systems with a strong nebular continuum component, |$T_e$| 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 |$T_e$|(H|$^+$|⁠) implied by the Balmer jump is consistent with |$T_e$|(O|$^{++}$|⁠) measured from the [O iii] auroral line ratio.

3.4 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] |$\lambda \lambda$|4711, 4740 doublet in the Prism spectrum. These lines are partially blended with each other and with He ii|$\lambda$|4686, while [Ar iv] |$\lambda$|4711 is also completely blended with He i|$\lambda$|4713. Our three-component fit to this complex yields [Ar iv] |$\lambda$|4711/|$\lambda 4740=1.6\pm 0.8$| after subtracting the predicted He i|$\lambda$|4713 contribution (assuming |$\lambda$|4713/|$\lambda$|4471 = 0.15), consistent with the low-density limit (⁠|$n_e\lesssim 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 l-changing collisions at higher densities. The presence of the two-photon continuum will be discussed in Section 4.

3.5 Chemical abundances

We derive chemical abundances for GS-NDG-9422 adopting |$T_e=1.83\times 10^4$| K and a density of |$n_e=100$| cm|$^{-3}$| following the procedure in Cameron et al. (2023a). Given the apparent agreement between |$T_e({\rm H}^+)$| and |$T_e({\rm O}^{++})$|⁠, we assume a constant temperature for all ionization zones. We derive |$12+\log ({\rm O}/ {\rm H})=7.59\pm 0.01$| from the [O ii], [O iii], and H |$\beta$| lines, obtaining a consistent result with both Prism and grating line fluxes. We derive the carbon abundance from the C iii] |$\lambda \lambda$|1907, 1909 / [O iii] |$\lambda$|5007 ratio measured from the Prism, assuming no dust, finding |$\log ({\rm C}/{\rm O})=-0.73\pm 0.03$| after applying the ionization correction factor (ICF) presented in Amayo, Delgado-Inglada & Stasińska (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] |$\lambda \lambda$|1750/[O iii] |$\lambda$|5007 limit measured from the Prism, and the N iii] |$\lambda \lambda$|1750/O iii] |$\lambda$|1666 limit measured from the grating. These yield 3|$\sigma$| upper limits of |$\log ({\rm N}^{++}/{\rm O}^{++})\lt -0.85$| and |$\lt -1.01$| respectively. We adopt the former as our preferred limit.

We detect several helium recombination lines from both the singly and doubly ionized states. Prism measurements of He i|$\lambda$|4471/H |$\beta$| and He i|$\lambda$|5875/H |$\beta$| yield consistent measurements of |${\rm He}^{+}/{\rm H}^{+}=0.10\pm 0.005$| and |$0.10\pm 0.01$|⁠, respectively. Deriving a |${\rm He}^{++}/{\rm H}^{+}$| abundance using the He ii|$\lambda$|4686 line results in only a |$6\pm 1$| % contribution to the total helium abundance, giving |${\rm He}^{++}/{\rm H}^{+}=0.005\pm 0.001$|⁠. Together, this implies a total helium abundance of |${\rm He}/{\rm H}=0.11\pm 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, Torres-Peimbert & Ruiz 1992).

The values obtained from these measurements are summarized in Table 3.

Table 3.

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

Property	Derived value
|$T_e({\rm O}^{++})$|	|$18,300\pm 1,500$| K
|$n_e$|	|$\lesssim 10^3$| cm|$^{-3}$|
|$12+\log ({\rm O}^{++}/{\rm H}^+)$|	|$7.58\pm 0.01$|
|$12+\log ({\rm O}^{+}/{\rm H}^+)$|	|$5.79^{+0.04}_{-0.03}$|
|$12+\log ({\rm O}/{\rm H})$|	|$7.59\pm 0.01$|
|$\log ({\rm N}^{++}/{\rm O}^{++})$|	|$\lt -0.85$|
|$\log ({\rm C}^{++}/{\rm O}^{++})$|	|$-0.76\pm 0.03$|
|$\log ({\rm C}/{\rm O})$|	|$-0.73\pm 0.03$|
|$\log ({\rm Ne}/{\rm O})$|	|$-0.85^{+0.03}_{-0.04}$|
|${\rm He}^{+}/{\rm H}^+$|	|$0.10\pm 0.01$|
|${\rm He}^{++}/{\rm H}^+$|	|$0.005\pm 0.001$|
|${\rm He}/{\rm H}$|	|$0.11\pm 0.01$|
|$\log ({\rm O}^{++}/{\rm O}^+)$|	|$1.79^{+0.03}_{-0.04}$|
|${\rm He}^{++}/{\rm He}^+$|	|$0.06\pm 0.01$|

3.6 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|$_{\odot }$| and a high-mass slope of |$-2.35$|⁠. We consider a population with an age of 3 Myr and a metallicity of 0.1 |$Z_{\odot }$|⁠, 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 per cent or when the neutral column density reaches |$10^{18.7}\ {\rm cm^{-2}}$|⁠, consistent with the minimal observed Ly |$\alpha$| emission offset (⁠|$\lesssim 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] |$\lambda$|4363/|$\lambda$|5007, until we reproduce the line strengths of [O iii] |$\lambda$|5007, [O ii] |$\lambda \lambda$|3727, C iii] |$\lambda \lambda$|1909, H |$\alpha$|⁠, and H |$\gamma$| with respect to H |$\beta$|⁠. 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 |$\alpha$|⁠. Note that, by focusing on line ratios, the resulting continuum is a prediction of the model. The intrinsic spectrum of our best fit model (⁠|$n=10^{3}\ {\rm cm^{-3}}$|⁠, |$\log _{10}(U)=1.2$|⁠, |$\log _{10}(Z_{\rm O}/Z_{\odot })=-1.1$|⁠, and |$\log _{10}(Z_{\rm C}/Z_{\odot })=-1.4$|⁠) is shown as the magenta line in Fig. 4.

Figure 4.

Prism spectrum of GS-NDG-9422 (black) compared with the best-fitting models using a standard SSP (magenta) accounting for both the |$z=6$| IGM opacity (dashed yellow) and a DLA with column density of |$10^{23.1}\ {\rm cm^{-2}}$| (dotted cyan).

As can be seen in Fig. 4, the strengths of the emission lines are well reproduced by our best-fitting 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¹, 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, 2018; Schaerer et al. 2019; Saxena et al. 2020; Senchyna et al. 2020).

4 CONTINUUM SHAPE

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

4.1 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 Fig. 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] |$\lambda \lambda$|3726, 3729 and [Ne iii] |$\lambda$|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 |$f_\nu$|⁠) over the range |$\lambda _{\rm rest}\gt 1930$| Å with a two-part linear function, broken at the Balmer limit (⁠|$\lambda _{\rm rest}=3645$| Å). This results in a clear detection of a spectral jump with |$15.0\pm 0.9$| nJy in the observed frame (second panel Fig. 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. 2024) and in simulations (Wilkins et al. 2024; Katz et al. 2023a).

As shown in Fig. 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.

4.2 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 |$\lambda _{\rm obs}\approx 1$|μm (⁠|$\lambda _{\rm rest}\approx 1430$| Å) (Fig. 1). While our best-fitting 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 |$\lambda _{\rm rest}\lesssim 1430$| Å. Here, we explore the origin of this discrepancy.

4.2.1 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 |$\alpha$| absorption (DLA) systems (Heintz et al. 2024).

Beginning with the IGM, we apply the |$z=6$| IGM transmission curves from Garel et al. (2021), to our best-fitting photoionization model and we find that IGM damping is unable to produce the observed UV turnover from this BPASS model (yellow dashed line in Figs 4 and 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. (2024) for targets at much higher redshift. This is not particularly surprising because if the IGM were responsible, many more galaxies at |$z\gtrsim 6$| would show very strong UV turnovers.

Figure 5.

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

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-fitting BPASS model through slabs of neutral hydrogen with increasing column densities, finding that column densities of |$N_{\rm H}\sim 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 Fig. 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 |$\alpha$| emission. The middle panel of Fig. 5 shows that this can, in principle, be reconciled by invoking a DLA with 30 per cent 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 |$\alpha$| emission (Heintz et al. 2024; Umeda et al. 2023). In principle, one could shift the DLA to a lower redshift which would allow for some Ly |$\alpha$| emission to escape as the emitted Ly |$\alpha$| 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 per cent 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.

Other geometries may exist (apart from the picket fence or lower-redshift 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 |$\alpha$| escape fraction of |$\sim 27~{{\ \rm per\ cent}}$| which again seems unlikely. For this reason, we consider other alternatives for the origin of the UV turnover.

4.2.2 Two-photon continuum emission

The nebular continuum consists of three components: free–bound, which gives rise to the spectral discontinuity at |$\lambda _{\rm rest}\sim 3645$| Å, free–free, which is typically subdominant compared to free–bound and only impacts |$\lambda _{\rm rest}\gtrsim 3000$| Å, and two-photon emission (Fig. 6). Two-photon emission arises from transitions from |$2s\rightarrow 1s$| in neutral hydrogen, resulting in the emission of two photons whose energies sum to that of Ly |$\alpha$|⁠. The distribution of photons emitted via this process is symmetric around 2431 Å (⁠|$\nu = \frac{1}{2}\nu _{\rm Ly\alpha }$|⁠) when expressed in terms of number of photons per second per frequency interval. However, when expressed in units |$f_\lambda$|⁠, it takes the form of a broad asymmetric peak which turns over at |$\lambda _{\rm rest}\approx 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 ionizing 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.

Figure 6.

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 (⁠|$\lambda _{\rm rest}\approx 3645$| Å). Two-photon continuum has a fixed shape with a turnover at |$\lambda _{\rm rest}\approx 1430$| Å, but is usually subdominant compared to continuum emission of the ionizing source. Free–free emission is usually sub-dominant at these wavelengths. The relative contribution of two-photon continuum is highly dependent on nebular conditions.

To test whether the continuum of GS-NDG-9422 is consistent with being primarily nebular, we consider a model where the spectrum consists of:

• Two-photon emission

• Free–bound & Free–free emission

• Young stars with ages |$\lt 10^{7.5}$| yr

• Old stars with ages |$\gt 10^{7.5}$| yr

The shape of the spontaneous two-photon continuum is fixed², so we only vary its overall normalization. The combination of free–bound 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 Fig. 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 |$\lambda _{\rm rest}\lt 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 |$\chi ^2$| of |$\chi ^2_{\nu }=0.80$| indicating that this 50th percentile model is a very good fit to the continuum of GS-NDG-9422.

Figure 7.

Top: Best-fitting 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-fitting model. The blue line shows the overall best fit. Bottom: Posterior distribution on gas temperature predicted by the continuum fitting. The median and 1|$\sigma$| values (red lines) are in good agreement with the nebular temperature measured from the [O iii] |$\lambda$|4363/|$\lambda$|5007 ratio (Table 3; blue).

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 Fig. 7, we show the posterior distribution on gas temperature predicted by the continuum fitting (black histogram) compared to that measured using the [O iii] |$\lambda$|4363 auroral line (cyan), finding that the two are consistent within |$1\sigma$| 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. 2006; Guseva et al. 2007).

An even stronger test is to compare the observed H |$\beta$| equivalent width versus that predicted by the continuum-fitting model. The H |$\beta$| emission should primarily arise from recombination, the same as the free-bound continuum. Thus, in a nebular-dominated scenario, EW(H |$\beta$|⁠) should only depend on gas temperature and the relative contribution of free–bound and free–free to two-photon emission. The top panel of Fig. 8 shows that the |$1\sigma$| uncertainty on the observed EW(H |$\beta$|⁠) significantly overlaps the |$1\sigma$| spread in the posterior distribution of predicted EW(H |$\beta$|⁠). This again demonstrates that the information in the continuum is sufficient to explain many of the properties of the emission lines.

Figure 8.

Top: Posterior distribution on H|$\beta$| equivalent width predicted by the continuum fitting. The median and 1|$\sigma$| values (red lines) are consistent with that measured from the Prism spectrum (blue). Bottom: Same as top panel, but for H |$\alpha$|⁠. Here, the median and 1|$\sigma$| 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.

We can apply the same test for H |$\alpha$| emission. Since the continuum was only fit up to rest-frame 5800 Å, predicting the H |$\alpha$| equivalent width requires extrapolating the model. In Fig. 8, we compare the posterior distribution on predicted EW(H |$\alpha$|⁠) to that measured from the prism as well as that inferred from imaging. As discussed above, the EW|$_0$|(H |$\beta$|⁠) from the imaging data is somewhat higher than that measured from the prism (although the |$1\sigma$| contours overlap). We find that our predicted values fall high compared to the prism measurement, but the value from the imaging falls on top of our |$1\sigma$| confidence interval. Hence, our model is formally very consistent with the imaging data and consistent within |$2\sigma$| of the prism. Because we have not fit the continuum near H |$\alpha$| 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~Z_{\odot }$| is much higher than that predicted for the IGM at |$z\sim 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 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 |$\alpha$|⁠.

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?

5 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.

5.1 Is the ionizing 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. 2006; Guseva 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 ionizing 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|$\lambda$|4686/H |$\beta$| is |$\gt 6\times$| higher in the presumed QSO light echo compared to the galaxy studied here (Fig. 2). He ii|$\lambda$|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 (Fig. 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 ionizing source would be detected and we would be observing only the remnant nebular emission. However, the H|$^+$| recombination time-scale, assuming Case B recombination at a temperature of 18 300 K, is |$\sim 5\,000$| yr for a density of |$10^2\ {\rm cm^{-3}}$| or |$\sim 500$| yr for |$10^3\ {\rm 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 |$\lt 10^4$| yr (Eilers, Hennawi & Davies 2018; Eilers et al. 2021). While such QSOs are a rare subset of the general population, these lifetimes are broadly consistent with the recombination time-scale. 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 time-scale of |$\sim 300$| yr for a density of |$10^2\ {\rm cm^{-3}}$| or |$\sim 30$| yr for |$10^3\ {\rm 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.

5.2 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 ionizing photon production efficiency (⁠|$\xi _{\rm ion}$|⁠) than standard SSP models, necessitating hotter blackbody temperatures. We explore this by running cloudy simulations with input blackbody SEDs, with a set-up 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] |$\lambda$|4363/|$\lambda$|5007) and systematically vary the density and blackbody temperature (Fig. 9). A weak UV turnover begins to appear at blackbody temperatures |$T_{\rm BB}\gtrsim 65\,000$| K, which is hotter than a typical O star. A strong UV turnover requires at least |$T_{\rm BB}\gtrsim 90\,000$| K, much hotter than massive O stars (up to |$\sim$|50 000 K; Walborn et al. 2004; Evans et al. 2011; Bressan et al. 2012).

Figure 9.

Top: Normalized spectrum of hydrogen-only gas with |$n=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 |$T=10^5$| K at different gas densities (as given in the colour bar). The magenta line represents a density of |$10^3\ {\rm cm^{-3}}$|⁠.

The two-photon continuum is also sensitive to gas density because, at high densities, l-changing collisions will suppress the two-photon emission relative to Ly |$\alpha$|⁠. This is seen in the bottom panel of Fig. 9 where we vary gas density for a fixed blackbody temperature of 100 000 K. The UV turnover is strongly suppressed at |$n\gtrsim 10^3\ {\rm 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 |$\alpha$|⁠). Furthermore, there exist significant differences in atomic data predictions for the strength of l-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] |$\lambda$|5007/[O ii] |$\lambda \lambda$|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. Optimization of this model (top panel of Fig. 10) demonstrates that a blackbody model with |$T=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.

Figure 10.

Top: Spectrum of GS-NDG-9422 (black) compared with the best-fitting blackbody model (magenta) with |$T_{\rm BB}=10^{5.05}$| K and a He ii leakage fraction of 0.22 (Table 4; Section 5.2). Middle: Spectrum of GS-NDG-9422 (black) compared with photoionziation 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-fitting blackbody SED. All models can also reproduce the observed spectrum well (dotted lines). See Section 5.2 for details.

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 (Fig. 11). The weak He ii|$\lambda$|1640 and He ii|$\lambda$|4686 in GS-NDG-9422 places tight constraints on the ionizing spectrum at |$E\gt 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|$^+$|-ionizing photons (⁠|$\lambda \lt 228$| Å; e.g. Smith, Norris & Crowther 2002). To explore this, we run models varying the leakage fraction of photons with |$E\gt 4$| Rydberg from 0  per cent (i.e. no He|$^+$|-ionizing photons) to 100 per cent (unattenuated blackbody). We conclude that |$70~{{\ \rm per\ cent}}\!-\!75~{{\ \rm per\ cent}}$| 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-fitting model. The parameters for our optimized model are reported in Table 4.

Figure 11.

Zoom in on the region around He ii|$\lambda$|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-fitting model from Fig. 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|$\lambda$|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|$\lambda$|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  per cent He ii leakage model shown in green.

Table 4.

Parameters for the optimized blackbody model.

Parameter	Value
n	|$10^3\ {\rm cm^{-3}}$|
|$T_{\rm BB}$|	|$10^{5.05}$| K
log U	0.0
|$T_{\rm gas}$|	|$10^{4.3}$| K
|$Z_{\rm O}/Z_{\odot }$|	|$-1.1$|
|$Z_{\rm C}/Z_{\odot }$|	|$-1.6$|
O32	37.1
He ii leakage	0.25

We note that other sources, including active galactic nuclei (AGNs) or high-mass X-ray binaries (HMXBs), can have significant high-energy photon outputs, and have been invoked to explain peculiar emission line ratios seen at high-redshift (Maiolino et al. 2024; 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 overprediction 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 (⁠|$T\gt 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 |$\xi _{\rm 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-fitting blackbody model with |$T_{\rm 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)³ 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 |$Z=0.07~Z_{\odot }$|⁠, similar to that measured for GS-NDG-9422. We identify model 13–10 as most similar to our optimized blackbody, for which |$T=100\,000$| K and the luminosity is fixed at |$L=10^{5.3}L_{\odot }$| (Fig. 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 |$z\sim 6$|⁠. Abundances of carbon, nitrogen, and oxygen are assumed to have undergone significant CNO burning (see Todt et al. 2015 for details).

Figure 12.

Best-fitting blackbody SED model with He-rich atmosphere (black solid) compared to hot star model 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 |$Z=0.07~Z_\odot$| 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 ionizing 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-ionizing flux. The bottom right panel shows individual massive Pop. III stars from Larkin, Gerasimov & Burgasser (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 (Fig. 10), 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.

Adopting parameters from the optimized 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 optimization is performed and the resulting spectrum is shown as the dashed orange line in Fig. 10, nearly identical to the optimized 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|$\lambda$|1640 and He ii|$\lambda$|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, 2023). Compared to our best-fitting 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 (Fig. 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 (Fig. 12 bottom left) indicates these produce too many He|$^{+}$| ionizing photons, while predictions for even less top-heavy IMFs begin to fall short of the required |$\xi _{\rm ion}$|⁠. In contrast, some individual Pop. III star models from Larkin et al. (2023) with effective temperatures of |$\sim 97\,000$| K and masses of 85 to 108 M|$_{\odot }$| reasonably reproduce the required ionizing SED (Fig. 12 bottom right). We strongly emphasize that the measured metallicity of |$12+\log _{10}({\rm O/H})=7.59\pm 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 and 99 000 K also results in a very good fit to the spectrum of GS-NDG-9422 (green dashed lines in Fig. 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 (Fig. 10; see also schematic in Fig. 13). Wolf–Rayet stars with |$Z=0.07~Z_\odot$|⁠, 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.

Figure 13.

Schematic of how the hot star scenario gives rise to a nebular-dominated spectrum in a simple spherical nebula. The incident stellar spectrum (⁠|$T_{\rm eff}=10^{5.05}$| K) is shown in blue. The portion of this transmitted through the nebula is shown in the darker shade, while the lighter shade shows the portion that is absorbed to photoionization 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 |$T_{\rm eff}\gtrsim 10^{5}$| K emit significantly more ionizing 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 |$\alpha$| (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 primarily contribute flux at longer wavelengths.

6 DISCUSSION

6.1 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, DeGioia-Eastwood & Waterhouse 2001). Masses of |$\ge 37\ {\rm {\rm M}_{\odot }}$| have been estimated for Wolf–Rayet stars in the Small Magellanic Cloud (SMC) with |$T_{\rm eff}\gtrsim 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 |$\sim 100~{\rm {\rm M}_{\odot }}$|⁠. While our search may not be exhaustive, we conclude from Fig. 12 that the nebular-dominated spectrum of GS-NDG-9422 is consistent with ionization powered by low-metallicity massive stars (⁠|$\gtrsim 50~{\rm {\rm M}_{\odot }}$|⁠), perhaps in the Wolf–Rayet phase.

If we adopt a progenitor mass of |$50\!-\!100~{\rm {\rm M}_{\odot }}$|⁠, then, assuming a typical IMF, we expect to form one such star per every |$\sim 1300\!-\!3800~{\rm {\rm M}_{\odot }}$| of stellar mass. To explore whether our model is consistent with this, we repeat our photoionization modelling with our best-fitting 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 low-metallicity 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|$_{\odot }$|⁠) 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

$$\begin{eqnarray} M_{\rm SSP}=\frac{U_{\rm SSP}}{U_{\rm star}}\frac{Q_{\rm star}}{Q_{\rm SSP}}, \end{eqnarray}$$

(1)

where |$Q_{\rm star}=10^{49.36}\ \gamma \ {\rm s^{-1}}$| for the Wolf–Rayet star model or |$10^{49.98}\ \gamma \ {\rm s^{-1}}$| for the low-metallicity massive star model with |$T_{\rm eff}=97\,352$| K. |$Q_{\rm SSP}=10^{46.73}$| or |$10^{46.03}\ \gamma \ {\rm s^{-1}{\rm M}_{\odot }^{-1}}$| at 3 and 5 Myr, respectively. In both cases, we find the maximal allowable contribution to be |$\frac{U_{\rm SSP}}{U_{\rm star}}\le 6.3~{{\ \rm per\ cent}}$|⁠. Therefore, we can place upper limits of the BPASS SSP mass of |$M_{\rm SSP}\le 112$| or 137 M|$_{\odot }$| per one low-metallicity massive star or Wolf–Rayet star, respectively. This is clearly discrepeant from the values of |$3800$| and |$1300~{\rm {\rm M}_{\odot }}$| implied from a typical IMF. In other words, there is a |$35\times$| excess in the number of massive stars in the case of our low-metallicity massive star model, or a |$9.5\times$| 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 ionizing photon production rate, dominated in our model by stars more massive than |$\gtrsim ~50$| M|$_{\odot }$|⁠, and the continuum flux at |$\sim 1500$| Å, which in standard young stellar population models is dominated by stars with |${\rm M}\gtrsim 5~{\rm M}_{\odot }$| (e.g. Byler et al. 2017). It is therefore only sensitive to the high-mass end (⁠|${\rm M}\gtrsim 5~{\rm M}_{\odot }$|⁠) of the IMF, and is completely insensitive to the impact of stars with |${\rm M}\lesssim 5~{\rm M}_{\odot }$|⁠.

Taken at face value, this calculation implies that the ratio of |$\gtrsim 50$| M|$_{\odot }$| to |$\sim 5\!-\!50$| M|$_{\odot }$| stars in GS-NDG-9422 is of order |$\sim 10\times$| 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 |$\xi _{\rm ion}$|⁠, the excess drops as |$Q_{\rm 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 characterization 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 characterization 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, Omukai & Schneider 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, Covey & Meyer 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 top-heavy 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, van Dokkum & Villaume 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. 2024), both of which are sensitive to the mass-to-light ratios. We note that these measurements are sensitive to the ratios of stars with |$M\lesssim 0.5 {\rm M}_{\odot }$| and |$M\sim 0.5\!-\!1\,\, {\rm M}_{\odot }$|⁠, 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 preceded 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 |$z\sim 2-3$| dusty starburst galaxies have been interpreted as arising from a top-heavy 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).

6.2 Total cluster mass and number of hot stars

The hydrogen ionizing photon luminosity, Q, can be approximated from the H |$\beta$| flux as

$$\begin{eqnarray} Q \approx \frac{4\pi D_l^2 I_{\rm H\beta }}{(1-f_{\rm esc})h\nu _{\rm H\beta }}\frac{\alpha _{B}}{\alpha _{\rm B}^{\rm eff}}, \end{eqnarray}$$

(2)

where |$D_l$| is the luminosity distance, |$I_{\rm H\beta }$| is the measured H |$\beta$| flux, h is Planck’s constant, |$\nu _{\rm H\beta }$| is the frequency of the H |$\beta$| transition, |$f_{\rm esc}$| is the escape fraction of ionizing photons, |$\alpha _{B}$| is the total case B recombination rate, and |$\alpha _{\rm B}^{\rm eff}$| is the effective H |$\beta$| recombination rate. To calculate a luminosity distance, we adopt a cosmology with |$H_0=67.31$| km s|$^{-1}$| Mpc|$^{-1}$| and |$\Omega _{\rm m}=0.315$| (Planck Collaboration XIII 2016). This gives |$D_l=58.5$| Gpc, implying an H |$\beta$| luminosity of |$L_{\rm H\beta }=7.0\times 10^{41}$| erg s|$^{-1}$|⁠. Using the escape fraction derived from our best fitting blackbody photoionization model of 7.3 per cent and recombination rates evaluated at 20 000 K (Osterbrock & Ferland 2006), we estimate a hydrogen ionising photon luminosity of |$Q=1.65\times 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 |$10^{5.3} L_\odot$|⁠, which results in a hydrogen ionizing photon luminosity of |$Q_{\rm WR}=10^{49.36}$| s|$^{-1}$|⁠. However, known Wolf–Rayet stars in the most comparable local environments typically have significantly higher luminosities of |$\sim 10^{6.2} L_\odot$| (Shenar et al. 2016), implying that |$Q_{\rm WR}\approx 10^{50.3}$| s|$^{-1}$| might be a more realistic value. Based on the H |$\beta$| luminosity, this would imply |$\sim 10\,000$| of the |$Z=0.07~Z_\odot$| Wolf–Rayet stars would be needed to power the spectrum of GS-NDG-9422. In the alternative model, |$\sim 17\,000$| very metal-poor stars of |$\sim 100$| M|$_{\rm \odot }$| with |$T_{\rm 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}\ {\rm {\rm M}_{\odot }}$| for the Wolf–Rayet star model, or |$10^{6.55}\ {\rm {\rm M}_{\odot }}$| for the low-metallicity 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.

6.3 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 Fig. 14, we highlight two objects that both show Ly |$\alpha$| 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 |$z=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 |$\alpha$| 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 (⁠|$T_{\rm eff}\gtrsim 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.

Figure 14.

Zoom in on the rest-frame UV region for GS-NDG-9422 (black) compared with the spectra of the |$z=3.35$| Lynx arc (top) and A2744-NDG-ZD4 at |$z=7.88$| (bottom). We have normalized the spectra and shifted them all to the rest frame.

We identify a second galaxy, A2744-NDG-ZD4 at |$z=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 |$\alpha$| escape and a UV turnover (Fig. 14 bottom panel). The spectrum of this object, retrieved from Dawn JWST Archive (DJA)⁴ which had been reduced using the custom-made pipeline msaexp v.0.6.12⁵ (see Heintz et al. 2024 for details), received only 2.3 h 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.

7 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 |$z=5.943$|⁠. The identification of a spectral discontinuity of |$15\pm 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] |$\lambda$|5007, H |$\beta$|⁠, and H |$\alpha$|⁠. We measure high nebular temperatures (⁠|$(1.83\pm 0.15)\times 10^4$| K), low densities (⁠|$n_e\lesssim 10^3$| cm|$^{-3}$|⁠), and low metallicity [|$12+\log ({\rm O}/{\rm H})=7.59\pm 0.01$|⁠; |$\sim$|8 per cent (O/H)|$_\odot$|] 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 |$\lambda _{\rm rest}\approx 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 Fig. 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 |$\alpha$| escape fraction, and the UV turnover redward of Ly |$\alpha$| .

The primary challenge faced by models explaining this feature with a DLA is that the magnitude of the turnover implies unprecedented column densities of |$N_{\rm HI}\gtrsim 10^{23}$| cm|$^{-2}$|⁠, higher than any known DLA. Furthermore, the observation of narrow, high EW Ly |$\alpha$| emission with a high-escape fraction (⁠|$\sim 27~{{\ \rm per\ cent}}$|⁠) requires some mechanism by which this emission can reach the observer. This can be reconciled if Ly |$\alpha$| 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 |$\alpha$| might be allowed to reach the observer if the DLA has only partial coverage (Scenario #2) or 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.

We demonstrated in Fig. 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 ionizing source?

One possibility is that no flux from the ionizing source is observed, either due to a spatial offset whereby the ionizing source falls outside the slit (Scenario #4) or because the ionizing 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 (Fig. 1). The predicted recombination time-scales in Scenario #5 (⁠|$\sim$|500–5 000 yr for H|$^+$| and |$\sim$|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 time-scale of |$\lesssim 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 ionizing source is indeed present in the slit, but its ionizing 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 Fig. 13). This is a generic prediction of all photoionization models that include a significant population of hot stars with |$T_{\rm eff}\gtrsim 10^5$| K (e.g. Schaerer 2002; Raiter et al. 2010; Inoue 2011; Zackrisson et al. 2011; Trussler et al. 2023). We identify several model stellar SEDs that can attain these temperatures, showing that models with ionizing spectra dominated by |$\sim$|7  per cent |$Z_\odot$| Wolf–Rayet stars and/or very low-metallicity massive stars are able to reproduce the peculiar shape of the observed continuum. The critical aspect to this scenario is that it requires a top-heavy IMF, with the implied ratio of |$\gtrsim$|50 M|$_\odot$| to |$\sim$|5–50 M|$_\odot$| 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. None the less, 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 ionizing source on time-scales 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 |$z=3.35$| and A2744-NDG-ZD4 at |$z=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.

ACKNOWLEDGEMENTS

We thank Bob Fosbury for useful discussions and providing us with the LRIS spectrum of the Lynx arc. We thank Paul Crowther for providing insightful and constructive comments. We thank Thibault Garel for providing us with the simulated IGM attenuated curves from the SPHINX simulation. We thank Anna Feltre and Stéphane Charlot for sharing their latest AGN photoionization model grids. AJC, AS, and AJB, 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). CW thanks the Science and Technology Facilities Council (STFC) for a PhD studentship, funded by UKRI grant 2602262. NL acknowledges support from the Kavli foundation.

DATA AVAILABILITY

Reduced spectra, imaging, and measured photometry of GS-NDG-9422 were made publicly available as part of JADES Data Release 1 (Bunker et al. 2023; Rieke et al. 2023) and can be found at https://archive.stsci.edu/hlsp/jades. Photoionization models presented in this work will be shared on reasonable request.

Footnotes

1

For example, the He ii|$\lambda$|1640/H |$\beta$| ratio is measured to be |$0.26\pm 0.05$| whereas the photoionization models predict a value of 0.03.

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 |$f_\lambda$|⁠. However, we found this effect to be completely negligible (⁠|$\lt 1$| Å) in any of the modelling considered in this work.

3

https://www.astro.physik.uni-potsdam.de/~wrh/PoWR/powrgrid1.php

4

https://dawn-cph.github.io/dja/index.html

5

https://zenodo.org/record/8319596

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 |$\alpha$|⁠. 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 two-photon 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 Fig. 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}\ {\rm cm^{-2}}$|⁠, which is slightly lower than reported in Section 3, but very consistent within |$1\sigma$|⁠. In either case, the DLA column required is exceptionally high compared to known DLAs.

Figure A1.

Top: Best-fitting 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.

Interestingly, despite the fact that Ly |$\alpha$| was not included in the fitting, the MCMC tightly constrains the amount of leakage to be |$\sim 35~{{\ \rm per\ cent}}$|⁠. 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] |$\lambda$|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, 2007). Moreover, in this model, the predicted HĪ|$\beta$| EW is 767  Å|$\pm 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\,\, {\rm {\rm M}_{\odot }}$| and disc 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 Fig. B1, we show our XRB model for 25 M|$_{\odot }$| black holes optimised to reproduce the continuum shape of GS-NDG-9422. This model significantly over-predicts the strength of the He ii|$\lambda$|1640 and |$\lambda$|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.

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\ {\rm {\rm M}_{\odot }}$| (magenta). The solid and dashed magenta lines indicate models with different accretion disc radii.