First insights into the ISM at z > 8 with JWST: possible physical implications of a high [O iii] λ4363/[O iii] λ5007

ABSTRACT

We present a detailed analysis of the rest-frame optical emission line ratios for three spectroscopically confirmed galaxies at z > 7.5. The galaxies were identified in the James Webb Space Telescope (JWST) Early Release Observations field SMACS J0723.3 − 7327. By quantitatively comparing Balmer and oxygen line ratios of these galaxies with various low-redshift ‘analogue’ populations (e.g. Green Peas, Blueberries, etc.), we show that no single analogue population captures the diversity of line ratios of all three galaxies observed at z > 7.5. We find that S06355 at z = 7.67 and S10612 at z = 7.66 are similar to local Green Peas and Blueberries. In contrast, S04590 at z = 8.50 appears to be significantly different from the other two galaxies, most resembling extremely low-metallicity systems in the local Universe. Perhaps the most striking spectral feature in S04590 is the curiously high [O iii] λ4363/[O iii] λ5007 ratio (RO3) of 0.048 (or 0.055 when dust-corrected), implying either extremely high electron temperatures, >3 × 10⁴ K, or gas densities >10⁴ cm⁻³. Observed line ratios indicate that this galaxy is unlikely to host an AGN. Using photoionization modelling, we show that the inclusion of high-mass X-ray binaries or a high cosmic ray background in addition to a young, low-metallicity stellar population can provide the additional heating necessary to explain the observed high RO3 while remaining consistent with other observed line ratios. Our models represent a first step at accurately characterizing the dominant sources of photoionization and heating at very high redshifts, demonstrating that non-thermal processes may become important as we probe deeper into the Epoch of Reionization.

1 INTRODUCTION

The first stars, galaxies, and black holes to form at Cosmic Dawn directly impact all subsequent galaxy formation in the Universe. The stars in the first galaxies synthesize and distribute the first significant amounts of heavy elements, which modifies how gas can cool and form stars (e.g. Omukai et al. 2005), while their ionizing photons heat and pressure smooth the intergalactic medium, impacting how galaxies can accrete gas from cosmic filaments (e.g. Okamoto, Gao & Theuns 2008; Katz et al. 2020). Understanding the interstellar medium (ISM) within the first few generations of galaxies, and both the efficiency and properties of early star formation, is crucial for constraining the physics of the early Universe.

Rest-frame optical and UV emission line ratios provide robust diagnostics of the ISM conditions of galaxies, and the underlying photoionization or other energetic mechanisms that shape galaxy spectra. Emission line ratios have been widely used at lower redshifts (z ≲ 3) to distinguish photoionization due to star-formation from active galactic nuclei (AGN) activity and other ionization sources (e.g. BPT diagnostics; Baldwin, Phillips & Terlevich 1981; Kewley et al. 2001) as well as to constrain the underlying physical conditions (e.g. density, metallicity, ionization parameter) within the ISM of galaxies (e.g. Kewley, Nicholls & Sutherland 2019). However, to date, the lack of sensitive space-based near-infrared spectroscopy has meant that detailed emission line studies of the ISM within the Epoch of Reionization (z ≳ 6) have only been possible in the rest-frame far-infrared (most notably with ALMA), where few systems have been observed (e.g. Carniani et al. 2020). As a result, many studies have instead placed constraints on the ISM conditions of galaxies within the Epoch of Reionization using various samples of local galaxies deemed to be analogues of high-redshift galaxies (e.g. ‘Green Peas’; Cardamone et al. 2009; Yang et al. 2017a, ‘Blueberries’; Yang et al. 2017b, low-metallicity galaxies; Izotov et al. 2019, etc.).

The recent successful launch and commissioning of the James Webb Space Telescope (JWST; Rigby et al. 2022) is set to change this landscape. In particular, the Near-Infrared Spectrograph (NIRSpec; Ferruit et al. 2022; Jakobsen et al. 2022) enables deep spectroscopy to be obtained across the full ∼1–5 |$\mu$|m wavelength range, allowing us to study rest-optical emission line ratios of z ≳ 6 galaxies with unprecedented sensitivity. Already, NIRSpec spectroscopic follow up of more than 40 new high-redshift galaxy candidates identified from NIRCam Early Release Observations (ERO; Program ID 2736, PI: K. Pontoppidan) of the SMACS J0723.3-7327 lensing cluster has revealed prominent detections of numerous rest-frame optical emission lines (Pontoppidan et al. 2022), enabling for the first time: (1) comparisons between the ISM conditions of ‘high-redshift analogues’ and true z ≳ 6 galaxies; and (2) direct constraints on the ISM conditions of z ≳ 6 galaxies from rest-optical emission lines.

To constrain the formation and evolution of some of the earliest structures in the Universe, it is vital to compare the observed spectroscopic properties with predictions from models. This includes both photoionization models of the ISM (e.g. Gutkin, Charlot & Bruzual 2016; Nakajima & Maiolino 2022), or numerical simulations that attempt to model the relevant physics from first principles (e.g. Katz et al. 2019). Such comparisons will allow us to accurately infer the dominant sources of ionizing photons and the important cooling and heating processes in high-redshift galaxies. To this end, in this work, we investigate the rest-frame optical emission lines in three galaxies at z > 7.5 in the SMACS0723 field and compare their spectroscopic properties with various low-redshift galaxy populations and infer the physics that is governing the ISM in the early Universe.

The layout of this paper is as follows: we discuss the JWST NIRSpec ERO data used in this study, along with reduction steps, properties of the three z > 7.5 galaxies and emission line fitting methodology in Section 2. We discuss insights obtained from emission line ratios and compare the observations of galaxies at z > 7.5 with those in the local Universe in Section 3. In Section 4, we explore the physics that can generate the curiously high [O iii] λ4363/[O iii] λ5007 ratio observed at z = 8.5 using a suite of photoionization models. Our discussion and conclusions are presented in Section 5.

2 DATA AND EMISSION LINE FITTING

We use data from JWST NIRSpec Micro-Shutter Assembly (MSA) observations of the SMACSJ0723.3-7327 galaxy cluster taken on 2022 Jun 30, as part of the public ERO data release (Program ID 2736, PI: K. Pontoppidan, Pontoppidan et al. 2022). The 0.2 arcsec × 0.46 arcsec MSA shutters are configurable and can be used to acquire spectra of hundreds of sources in a single field-of-view (Ferruit et al. 2022). The spectra were taken in two filter/disperser combinations: F170LP/G235M and F290LP/G395M which combine to span from 1.66–5.10 μm with a resolving power of R ∼ 1000. The MSA configuration involved placing slitlets comprised of a column of at least three consecutive open shutters in a single column on each target. Two 8753-s observations were obtained for this field (Obs 007 and Obs 008 hereafter), each of which involved a three-point nodding pattern to enable local background subtraction, for a total of ∼17 507 s.

Our analysis in this paper is performed on the publicly available¹ extracted 1D spectra presented in Curti et al. (2022). These resulted from processing Stage 1 outputs from the JWST data calibration pipeline with the GTO pipeline (NIRSpec GTO collaboration, in preparation). These data were reduced assuming a point-like source profile, although from NIRCam imaging we see that the sources are marginally resolved. While this assumption may introduce some uncertainty on the absolute flux calibration, the analysis presented in this paper is based only on emission line flux ratios. We estimate that the assumed source profile affects these flux ratios at the |$\lt 10{{\ \rm per\ cent}}$| level. Background subtraction was performed locally by subtracting the background flux as measured by neighbouring ‘empty’ shutters. The reader is directed to Curti et al. (2022) for specific details about the calibration of these data. Finally, we note that, although we used the same 1D extracted spectra, the emission line fitting presented here was done independently of that in Curti et al. (2022).

This study focuses on three objects confirmed to lie at z > 7.5: S04590 at z = 8.50 (RA: 110.89533, Dec: −73.44916), S06355 at z = 7.67 (RA: 110.84452, Dec: −73.43508), and S10612 at z = 7.66 (RA: 110.89533, Dec: −73.43454). These objects are all gravitationally lensed by the foreground cluster; however, crucially for this study, the achromatic nature of gravitational lensing means that the emission line ratios are magnification-independent, and therefore, the results presented here are not affected by magnification at all. Extracted 1D spectra of all three objects along with the identified emission lines are shown in Fig. 1.

We run emission line fitting for the following rest-optical emission lines, identified as being detected in the spectra: [O ii] λλ3726, 3729, [Ne iii] λ3869, H δ, H γ, [O iii] λ4363, H β, [O iii] λ4959, and [O iii] λ5007. We note that all of the emission lines used in this study are measured in the F290LP/G395M filter/disperser setting, so there is no uncertainty introduced by comparing the relative calibration of these two settings. We find the emission line profiles can be well fit by a single Gaussian function and use the best fitting parameters to infer line fluxes and ratios.

To obtain the best-fitting Gaussian, we model the continuum (where detected) or noise level around each emission line as a linear fit and measure its median and standard deviation. Best-fitting Gaussians (red) for the main emission lines investigated in this study are shown in Fig. 1. Although we do not force the width of the Gaussians to be consistent across all emission lines fitted, we find good agreement in the measured full-width at half-maxima of the best-fitting Gaussians (≈200–250 km s⁻¹, consistent with being unresolved at R ∼ 1000). We find that our ‘local’ estimates of the continuum/noise that we use to fit Gaussians (blue dashed line) are broadly consistent with the noise estimates in the spectra returned by the data reduction pipeline (grey solid line) in line with expectations from weak to no continuum detection. A notable exception is the [O iii] λ4363 in S06355 (middle panel in Fig. 1), where two noise peaks appear to straddle the [O iii] λ4363 peak.

We find secure detections of H β, [O iii] λ4959, and [O iii] λ5007 across all three galaxies. The [O ii] |$\lambda \lambda \, 3726,3729$| doublet is unresolved, and is weak with respect to [O iii] λ5007, as evidenced by the relatively high [O iii]/[O ii] ratios (see Table 1) for all three galaxies. Interestingly, the [O iii] λ4363 auroral line is securely detected in S04590 and S10612 with S/N ≈ 5.3 and ≈6.3, respectively, with a lower S/N ≈ 2.8 detection in S06355. We also securely detect H δ and H γ in the galaxies, in addition to the [Ne iii] λ3869 and λ3967 lines, although the [Ne iii] λ3967 line is blended with H ϵ.

We then correct the emission line fluxes for the effects of dust extinction based on the observed Balmer decrement H γ/H β, applying the appropriate dust correction following the Small Magellanic Cloud extinction curve from Gordon et al. (2003), assuming an intrinsic value of 0.468 corresponding to n_e = 100 cm⁻³ and T = 10⁴ K (e.g. Osterbrock & Ferland 2006). For S04590 we find that dust correction of with A_V = 0.64 ± 0.20 is required to achieve a Balmer decrement of 0.468,² in line (albeit slightly lower) with what was reported by Curti et al. (2022). This dust extinction value is consistent with the low-metallicity, low-stellar mass galaxies at high-redshift in the SPHINX²⁰ simulation (Katz et al. 2022a). We note that even without a dust correction, the RO3 value for this particular galaxy is much higher than is typically observed. We find that the other two galaxies require no dust correction, i.e. A_V < 0.1.

Despite employing slightly different emission line measurement techniques, the emission line ratios we measure are in agreement with those reported by Curti et al. (2022). The emission line ratios and their uncertainties employed in this study are given in Table 1.

3 INSIGHTS FROM REST-FRAME OPTICAL EMISSION LINES

We start our analysis of the observed emission line ratios of the three z > 7.5 galaxies by comparing them with various populations of lower redshift galaxies. Our goal is to address the question: how similar are reionization-epoch galaxies to low-redshift galaxies and to local analogues of high-redshift galaxies?

In Fig. 2 we show various strong-line diagnostics: [O iii] λ4363/[O iii] λ5007 (RO3) versus [O iii] λ5007/[O ii] λλ3726, 3729 (O32; top left); RO3 versus [Ne iii] λ3869/[O ii] λλ3726, 3729 (bottom left); RO3 versus [O iii] λ5007/H β (top right); and [O iii] λλ3726, 3729/H δ versus [Ne iii] λ3869/H δ (bottom right). In all panels, we show emission line ratio measurements from SDSS³ for galaxies and quasars, split based on the Kauffmann et al. (2003) BPT criteria, selecting only those sources that have ≥8σ detection of [O iii] λ4363. For additional comparison, where data is available, we show the locations of low-redshift Seyfert 2 nuclei (Armah et al. 2021), low redshift analogues in the form of Green Peas (Yang et al. 2017a), Blueberries (Yang et al. 2017b), and extremely low-metallicity galaxies (12 + log (O/H) ≲ 7.5; Izotov et al. 2019), as well as extreme emission-line galaxies (Amorín et al. 2015).

Interestingly, we find that there is no single low-redshift reference sample that can perfectly explain the properties of all three z > 7.5 galaxies (when selecting only those that exhibit [O iii] λ4363). Within error bars, S06355 at z = 7.67 (lime green star) is qualitatively consistent with Green Peas, Blueberries, and the edge of the distribution of SDSS galaxies.

Based on oxygen emission lines and H β, S10612 at z = 7.66 (cyan star) is consistent with Blueberries, but seems to have higher O32 and RO3 compared to the Green Peas, and significantly higher [O iii] λ5007/H β compared to the extremely low-metallicity galaxies from Izotov et al. (2019). This galaxy also shows exceptionally strong [Ne iii] λ3869, falling above the AGN demarcation line from Pérez-Montero et al. (2007) when comparing [Ne iii] λ3869/H δ versus [O ii] λλ3726, 3729/H δ. We show below that such emission line ratios can also be driven by strong shocks, so this galaxy is not conclusively an AGN, nor has a signature of higher excitation energy emission lines such as He ii λ4686.

Similar to S06355, S04590 at z = 8.5 (magenta star) shows very high O32, RO3, and [Ne iii] λ3869/[O iii] λλ3726, 3729. This galaxy resides at the very edges of the distributions of Blueberries and extremely low-metallicity galaxies. The dust-corrected RO3 value of S04590, at 0.059, is among the highest currently known in the local Universe. Unlike S06355, S04590 is characterized by very high Balmer emission compared to the emission from metal lines, with [O iii] λ5007/H β lying below the typical values in SDSS. Furthermore, both [Ne iii] λ3869/H δ and [O ii] λλ3726, 3729/H δ are low compared to the typical SDSS galaxy, as can be seen in the bottom right-hand panel of Fig. 2.

To quantify the similarity between these three high-redshift galaxies and each low-redshift reference population, we train an Isolation Forest (Liu, Ting & Zhou 2008) on the three emission line ratios of O32, RO3, and [O iii] λ5007/H β. We compute the probability that each z > 7.5 galaxy is an ‘inlier’ for each reference galaxy population by randomly sampling the three flux ratios 10 000 times across the error distribution of the reference populations and using the Isolation Forest to classify⁴ the galaxy.

In Table 2 we report the fraction of these 10 000 samples that are classified as inliers for each reference population. Consistent with our qualitative analysis above, we see that S06355 at z = 7.67 is marginally consistent with SDSS galaxies, Green Peas, and consistent with the spread in properties exhibited by Blueberries, while S10612 at z = 7.66 is only consistent with Green Peas and Blueberries. In contrast, S04590 at z = 8.50 is most similar to the extremely low metallicity galaxies of Izotov et al. (2019) and sits at the very edge of the distribution of Blueberries. This quantitatively demonstrates that these three high-redshift galaxies, despite lying at similar redshifts, are diverse in nature, and multiple different low-redshift analogue samples are needed to realistically encapsulate their properties (cf. Schaerer et al. 2022; Rhoads et al. 2022).

Since all three galaxies have high O32 and RO3 compared to SDSS galaxies, the ISM in these high-redshift systems likely have high ionization parameters, lower metallicities, and/or higher H ii region electron temperatures and densities compared to that in the local Universe (see also Trump et al. 2022). However, with this initial analysis, using the dust-corrected values, two key features stand out: the very high RO3 in S04590 and the very high [Ne iii] λ3869/H δ and O32 in S10612. The latter can likely be explained by a high ionization parameter and low metallicity; however, the former requires more investigation, which we discuss further below.

4 CURIOUSLY HIGH [O iii] |$\mathbf {\lambda 4363}$|/[O iii] |$\mathbf {\lambda 5007}$| RATIO AT REDSHIFT 8.5

Qualitatively, the spectral feature that is immediately interesting is the strength of the [O iii] λ4363 line with respect to [O iii] λ5007 in S04590 at z = 8.5. For galaxies in SDSS with measured [O iii] λ4363, we find log₁₀(RO3) = −1.91 ± 0.3 (i.e. 0.012). The value is larger for low-metallicity galaxies (−1.41 or 0.039; Izotov et al. 2019), however, S04590 has a log₁₀(RO3) of −1.32 (−1.26 when corrected for dust), which is a 2σ outlier compared to the SDSS population with detected [O iii] λ4363.

Fig. 3 shows contours of constant RO3 as a function of electron temperature and density, demonstrating the physical conditions required to give rise to this high RO3 value. High RO3 values can either be driven by high electron temperatures (T_e) or high electron densities (n_e). The change in level populations with temperature drives the behaviour of the ratio with temperature, while the dependence on density is caused by the fact that the critical density of the λ5007 line is lower than that of the λ4363 line (e.g. Osterbrock & Ferland 2006).

In order to obtain the dust-corrected value of RO3 (∼0.055) as is seen in S04590, one either requires temperatures ∼3 × 10⁴ K, or densities >10⁴ cm⁻³. Such a temperature is higher than what is typically seen in galactic H ii regions (e.g. Deharveng et al. 2000) where T_e ∼ 7000–14000 K and n_e ∼ 100 cm⁻³ (Osterbrock & Ferland 2006), and those computed from photoionization models for typical stellar spectral energy distribution (SED) choices (e.g. Xiao, Stanway & Eldridge 2018), while such electron densities are an order of magnitude greater than those seen in star-forming galaxies at z ∼ 2, where H ii region densities can be directly inferred (e.g. Sanders et al. 2016). To demonstrate this further, in the top row of Fig. 4, we show the results of a set of cloudy models (Ferland et al. 2017) where we have varied the stellar age, ionization parameter, and metallicity of the gas and stellar population that is illuminating the cloud (see Appendix  A for details on the models). Only the models with the lowest metallicities (10⁻²Z_⊙) and the highest ionization parameters (10^−1.5) begin to approach the values of RO3 and O32 seen in S04590; however, these models cannot simultaneously reproduce the strength of [O iii] λ5007/H β due to the low metallicity. Hence, under relatively standard assumptions (see Appendix  B for an extended set of models) on H ii region properties, it is difficult to reconcile what is driving the emission line ratios observed in S04590.

Before investigating the physical conditions required to drive such a high RO3, we reiterate the fact that JWST/NIRSpec is a new instrument and the observations underpinning this work are among the first taken at such high redshifts. We note that in an independent analysis of these spectra, Curti et al. (2022) reported a slightly lower value of RO3 = 0.045 ± 0.007 for S04590 (before dust correction). While this value is lower than what we find here (0.048), we emphasize that the models described above are still unable to simultaneously reproduce those RO3 and O32 ratios reported in Curti et al. (2022). If further observations continue to find such high RO3 ratios in high-redshift galaxies, there must be some additional physics beyond star formation that is responsible. We consider four possible solutions below:

• The presence of an AGN

• Strong shocks

• X-ray heating

• Cosmic ray heating

4.1 The presence of an AGN?

The observation of high RO3 is reminiscent of the temperature problem in low-redshift Seyfert galaxies (e.g. Koski & Osterbrock 1976; Osterbrock 1978; Ferland & Netzer 1983; Dopita & Sutherland 1995; Nagao, Murayama & Taniguchi 2001; Baskin & Laor 2005; Binette et al. 2022). The narrow line region (NLR) around an AGN can reach the densities and temperatures that are perhaps required to explain the high RO3 in S04590 (e.g. as was shown in the photoionization models of Baskin & Laor 2005). The line widths of S04590 in the spectrum remain spectrally unresolved at R ∼ 1000, resulting in an upper limit of <250 km s⁻¹, placing them on the low-end of what one might expect for an NLR. Nevertheless, we first consider an AGN as a potential solution for the high RO3 ratio.

Not only is S04590 inconsistent with SDSS QSOs and low-redshift Seyfert 2 nuclei on all diagnostics considered, it falls well below the AGN-galaxy demarcation proposed by (Pérez-Montero et al. 2007) on the [O ii] λλ3726, 3729/H δ versus [Ne iii] λ3869/H δ diagnostic. Without a measurement of [N ii] λ6583, Hα, we cannot distinguish between star formation and AGN activity using the BPT diagram. Furthermore high excitation lines such as He ii lines,⁵ or any of the higher ionization state Fe lines (these are likely too weak to be detected) are absent from the spectrum, making it hard to conclusively determine whether this galaxy hosts an AGN. However, none of the available evidence suggests that it does host an AGN. For these reasons, we consider alternative solutions to explain the spectrum of this object.

4.2 Strong shocks

Fast, radiative shocks can also increase the temperature of the gas well above 10⁴ K (Allen et al. 2008; Dopita & Sutherland 2017), potentially in the regime needed to explain the high RO3 in S04590. To determine whether such shocks can plausibly reproduce the emission line ratios we see at z > 7.5, we compare line ratios from the MAPPINGS models of Allen et al. (2008) to the three high-redshift galaxies in Fig. C1 across the four diagnostics presented above.

As before, we find that S04590 remains a clear outlier even when compared with these shock models. Therefore, we conclude that shocks alone are unlikely to explain the emission line ratios seen in this particular galaxy. We highlight that shocks can explain the [Ne iii] λ3869/H δ and [O ii] λλ3726, 3729/H δ ratios seen in S10612, suggesting that, although this galaxy lies above the Pérez-Montero et al. (2007) AGN demarcation line in this diagnostic plot (lower right-hand panel of Fig. 2), the line ratios can be explained without invoking the presence of an AGN.

4.3 X-ray binaries

High-energy photons from X-ray sources (including high-mass X-ray binaries or HMXBs) can potentially provide additional heating in the ISM (e.g. Lebouteiller et al. 2017; Lehmer et al. 2022), above what one would expect from photoionization due to stars alone. Observational data has shown that the X-ray luminosity per unit star formation rate increases significantly at low metallicities (e.g. Mapelli et al. 2010; Basu-Zych et al. 2013; Saxena et al. 2021). Furthermore, HMXBs have been postulated as an important ingredient to explain He ii emission in low-metallicity star-forming galaxies (e.g. Schaerer, Fragos & Izotov 2019; Senchyna et al. 2020; Saxena et al. 2020; Umeda et al. 2022). Since faint galaxies at z > 6 are predicted to have low metallicities, HMXBs could be an important source of heating that is often missing from typical H ii region photoionization or shock models. S04590 is predicted to have a much lower metallicity (12 + log(O/H) = 6.99 ± 0.11) than S06355 and S10612 (12 + log(O/H) = 8.24 ± 0.07, and 7.73 ± 0.12, respectively; Curti et al. 2022), which means that HMXBs are likely to have a larger impact on S04590 compared to the other two galaxies.

To assess the impact of HMXBs, we have expanded our cloudy model grid, varying the ratio of X-ray luminosity⁶ to star formation rate (L_X/SFR) in the range 10³⁸–10⁴² erg s⁻¹/M_⊙ yr⁻¹, where observational constraints out to z ∼ 5 currently exist (Saxena et al. 2021). Details of these models are presented in Appendix  A. In the middle row of Fig. 4 we show three strong-line diagnostics along with a subset⁷ of the results from our cloudy models.

The models with weak X-ray luminosities are nearly identical (on the diagnostics that we consider) to the runs with only stars. Once the X-ray luminosity is increased, the diagnostics fan away from the photoionization models that only include stars to high RO3 due to the additional heating from X-rays. As before, the models with lower metallicity tend to result in higher RO3 values due to the fact that cooling becomes less efficient with decreasing metallicity. Across all three diagnostics, the cloudy models with HMXBs can provide a reasonable explanation for the observed emission line ratios if high L_X/SFR > 10⁴¹ erg s⁻¹/M_⊙ yr⁻¹ is assumed. There remains some discrepancy in [Ne iii] λ3869/[O ii] λλ3726, 3729, but we highlight how all three diagnostics and in particular, [O iii] λ5007/H β, are significantly better at reproducing the observations compared to models with stars only. In summary, the inclusion of HMXBs with high L_X/SFR as predicted by extrapolating observed trends found z ∼ 3–5 (Saxena et al. 2021) combined with low metallicity can plausibly reproduce the emission line ratios seen in S04590. Further tuning of the photoionization models (e.g. higher ionization parameters, changing the mass of the HMXB, non-solar abundance ratios) may improve the agreement. However, although physically motivated, we emphasize that this solution may not be unique because there are other heating mechanisms that can increase ISM temperature.

4.4 High cosmic ray densities

Ferland & Mushotzky (1984) considered cosmic rays as an additional heating mechanism to explain the observed RO3 ratios in low-redshift AGN. Assuming a simple model where galaxy volume scales with the size of the Universe (as (1 + z)⁻³),⁸ it is plausible that the cosmic ray densities for high-redshift galaxies with similar star formation rates as the Milky Way may be ∼1000 times higher than what is inferred from the galactic neighbourhood. The cosmic ray density should further scale with star formation rate if the cosmic rays are sourced by supernova explosions. Papadopoulos et al. (2011) showed that in cosmic ray-dominated regions where the cosmic ray density is 10³–10⁴ times higher than the local neighbourhood, temperatures can easily rise to 5–10 × 10⁴ K, which could help explain the high temperature needed to reproduce the observed RO3 in S04590. Such scenarios are predicted to operate in local ultraluminous infrared galaxies or high-redshift submillimeter galaxies. With the additional pressure and temperature, the stellar initial mass function is likely boosted in cosmic ray-dominated regions due to the increased Jeans mass, which aligns with certain predictions of a more top-heavy initial mass function at high-redshift (e.g. Katz et al. 2022b).

To assess the impact of cosmic rays, we have further extended our cloudy grid by systematically increasing the strength of the cosmic ray background. The bottom row of Fig. 4 shows the results of these cloudy models compared to the three z > 7.5 JWST galaxies. Here we can see fundamentally different behaviour compared to the other two grids. As the cosmic ray background is increased, RO3 increases due to the additional heating from the cosmic rays. However, initially O32, [Ne iii] λ3869/[O ii] λλ3726, 3729, and [O iii]/H β drop as well due to the increased temperature in the outer regions of the gas cloud allowing for more [O ii] λλ3726, 3729 and H β emission. The models with 10⁻²Z_⊙, U = 10^−1.5, and a cosmic ray background intensity that is 100 × greater than that of the Milky Way can plausibly explain both the O32, RO3, and [Ne iii] λ3869/[O ii] λλ3726, 3729, but in this model, such low metallicity combined with the additional H β emission means that [O iii]/H β is systematically underpredicted.

Interestingly, we can reproduce the behaviour of S04590 using cosmic rays in two ways. For each choice of ionization parameter and metallicity, there is a particular cosmic ray background intensity where the cloud switches from being ionization bounded to being density bounded. In this regime, O32 increases (see e.g. Nakajima & Ouchi 2014) as the O ii region in the nebula is destroyed. For the same reason, [Ne iii] λ3869/[O ii] λλ3726, 3729 also increases. Finally, at sufficiently hot temperatures, [O iii]/H β also increases because the recombination emissivity for H β drops with temperature. If the cosmic ray background is very fine-tuned, the models that are right at the edge of the density and ionization bounded regimes can conspire to reproduce all of the flux ratios shown in Fig. 4. Due to the fine tuning required, this is not our preferred solution, although we highlight that in this unique regime, cosmic rays will be the dominant driver of Lyman continuum leakage. To our knowledge, a cosmic ray-dominated ISM has not been considered as a mechanism for causing ionizing photons to leak from galaxies.

More generally, if we increase the ionization parameter to U = 10^−0.5, we can fit all the emission line ratios of S04590 shown in Fig. 4. This solution still requires a cosmic ray background intensity that is more than 100 × that of the Milky Way; however, such a strong background is not unrealistic due to the small sizes and bursty nature of high-redshift galaxies as predicted by numerical simulations (e.g. Rosdahl et al. 2018).

5 CONCLUSIONS

We have presented an analysis on the spectra of the first three z > 7.5 galaxies observed with the JWST/NIRSpec Micro-Shutter Assembly in the SMACS0723 field, realased as part of the Early Release Observations public data release.

We have compared emission line ratios from these high-redshift galaxies to various low-redshift galaxy (SDSS), QSO (SDSS), and potential high-redshift analogue (Amorín et al. 2015; Yang et al. 2017a, b; Izotov et al. 2019) populations. We find that the properties of the high-redshift galaxies are diverse and not a single analogue population well represents the flux ratios of all three systems across various diagnostics. Nevertheless, studying only oxygen and Balmer lines, the two galaxies (S06355 and S10612) at z ∼ 7.6 are consistent with the low-redshift analogues, while S04590 at z = 8.5 seems to reside at the very edge of the low-redshift analogue distribution.

S04590 is characterized by an extremely high [O iii] λ4363/[O iii] λ5007 ratio, which is at the upper limit of what is seen in very low metallicity galaxies in the local Universe (Izotov et al. 2019). We consider the possibility that additional physics beyond low-metallicity star formation is required to explain the nature of this source.

We show that neither AGN nor shock models can simultaneously reproduce the flux ratios observed in S04590. However, if we assume that high-mass X-ray binaries or cosmic rays contribute significantly to the heating of the ISM, we can reasonably reproduce many of the properties of S04590. In such a scenario, either the L_X/SFR is very high and consistent with extrapolated predictions from z = 3–5 (e.g. Saxena et al. 2021), or the cosmic ray background intensity is orders of magnitude greater than what is observed in the Milky Way, in line with expectations of decreasing galaxy size at high-redshift or generally more concentrated star formation.

Our explanations of the observed flux ratios using cosmic rays or HMXBs are unlikely to be unique. In principle, any mechanism that can enhance the ionizing photon count and/or drive additional heating could potentially enhance RO3.

Due to the early nature of these observations, we remain cautious about over-interpreting the spectra, particularly the high RO3. Nevertheless, based on the available evidence, this work favours additional heating in the ISM, beyond photoheating from stars and photoelectric heating from dust grains to explain the high RO3. Cosmic rays or HMXBs could provide viable explanations for the high RO3 measured in S04590, indicating that we may need to revise our models of the ISM in the early Universe and how photons leak galaxies and contribute to reionization. JWST/NIRSpec observations of larger samples of z ≳ 8 galaxies will quantify whether or not these additional heating mechanisms are important in the high-redshift ISM.

ACKNOWLEDGEMENTS

We thank the referee for their detailed comments that improved the manuscript. AS acknowledges financial support from European Research Council (ERC) Advanced Grant FP7/669253. AJC 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). SA acknowledges funding from grant PID2021-127718NB-I00 by the Spanish Ministry of Science and Innovation/State Agency of Research (MCIN/AEI). RB acknowledges support from an STFC Ernest Rutherford Fellowship [grant number ST/T003596/1]. ECL acknowledges support of an STFC Webb Fellowship (ST/W001438/1). JW gratefully acknowledges support from the ERC Advanced Grant 695671, ‘QUENCH’, and the Fondation MERAC.

DATA AVAILABILITY

The 1D spectra used in this work are publicly available at https://doi.org/10.5281/zenodo.6940561. Any other data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

REFERENCES

APPENDIX A: DESCRIPTION OF cloudy MODELS

We have employed cloudy v17 (Ferland et al. 2017) to calculate emission line fluxes for simplified photoionization models accounting for stellar radiation, X-ray binaries, and a cosmic ray background. For all models, we assume a 1 kpc thick shell with a constant density of 100 H cm⁻³. In our fiducial set of models the inner radius of the slab is fixed to 10³⁰ cm for a plane-parallel model and the ionization parameter is varied between 10^−3.5 and 10^−1.5. The shape of the stellar SED is adopted from BPASSv2.0 (Stanway, Eldridge & Becker 2016; Eldridge & Stanway 2016). We assume an instantaneous burst and vary the age of the stellar population between 10⁶ and 10⁸ yr. We vary the metallicity of the gas and stellar population between 10⁻²Z_⊙ and 2Z_⊙. Metal abundances ratios are assumed to follow Asplund et al. (2009). Dust grains are included using the GRAINS ISM command and the dust-to-gas-mass ratio is assumed to scale with metallicity. For all models, we account for metal depletion on to dust grains. Furthermore, we have turned off molecules.

To account for X-ray binaries, in addition to the stellar SED and CMB, we have added a modified blackbody spectrum (Mitsuda et al. 1984) that we assume captures the SED of XRBs. We fix the black hole mass to 25 M_⊙. Following Senchyna et al. (2020), we assume the disk has a maximum radius of 10⁴ cm and compute the spectra using ARES, which is then combined with the stellar SED (Mirocha et al. 2012). Following observational constraints at high redshift (Saxena et al. 2021), we vary L_X/SFR in the range 10³⁸ − 10⁴² erg s⁻¹/M_⊙ yr⁻¹, where L_X is the X-ray luminosity in the range 2–10 keV.

To assess the impact of cosmic rays, we have expanded our cloudy grid by removing the HMXB SED and replacing it with a cosmic ray background (using the COSMIC RAYS BACKGROUND command). We vary the cosmic ray background intensities from 1 to 10^4.8 times the mean background strength of the Milky Way.⁹ For the cosmic ray runs, we also explore ionization parameters up to 10^−0.5 at a metallicity of 10⁻¹Z_⊙.

For all models, the calculation is stopped when the temperature reaches 4 000 K or when the outer radius of 1 kpc is reached. Note that these models are not designed to reproduce any particular galaxy population. Our goal is to qualitatively assess how HMXBs or cosmic rays may impact the relevant emission line ratios.

APPENDIX B: BPASS MODELS WITH HIGHER DENSITY AND IONIZATION PARAMETER

In Fig. B1, we show cloudy model grids assuming a BPASSv2.0 SED a single burst with an age of 5 Myr extending up to densities of 10⁵ cm⁻³ and ionization parameters up to 10^−0.5. We show that for any combination of these parameters, the emission line properties of S04590 are not reproduced.

APPENDIX C: COMPARISON WITH SHOCK MODELS

In Fig. C1, we compare the three z > 7.5 galaxies with the shock models from Allen et al. (2008) across the four primary strong line diagnostics analyzed in this work.