Nebular Line Emission During the Epoch of Reionization

TL;DR

This paper predicts nebular emission line properties of early galaxies during the epoch of reionization using simulations and photoionisation models, aiding future observations with the Webb Telescope.

Contribution

It combines the Bluetides simulation with photoionisation modeling to forecast nebular emission lines at high redshifts, filling a gap before Webb Telescope data becomes available.

Findings

Good agreement with existing observational constraints on equivalent widths.

Poorer agreement with luminosity function constraints.

Predictions will assist in interpreting Webb Telescope observations.

Abstract

Nebular emission lines associated with galactic HII regions carry information about both physical properties of the ionised gas and the source of ionising photons as well as providing the opportunity of measuring accurate redshifts and thus distances once a cosmological model is assumed. While nebular line emission has been extensively studied at lower redshift there are currently only few constraints within the epoch of reionisation (EoR, $z>6$), chiefly due to the lack of sensitive near-IR spectrographs. However, this will soon change with the arrival of the Webb Telescope providing sensitive near-IR spectroscopy covering the rest-frame UV and optical emission of galaxies in the EoR. In anticipation of Webb we combine the large cosmological hydrodynamical simulation Bluetides with photoionisation modelling to predict the nebular emission line properties of galaxies at $z=8\to 13$. We find good agreement with the, albeit limited, existing direct and indirect observational constraints on equivalent widths though poorer agreement with luminosity function constraints.

Equations

$$U_S\propto {\left(Q{ϵ}^2{n}_H\right)}^{1/3}$$

$$U_S=U_{S,\mathrm{ref}}{\left(\frac{Q}{Q_{\mathrm{ref}}}\right)}^{1/3}.$$

$${\tau }_{\lambda }={\tau }_V\times (\lambda /550\mathrm{nm}{)}^{-1}.$$

Full text

††pubyear: 2019††pagerange: Nebular Line Emission During the Epoch of Reionization–LABEL:lastpage

1 Introduction

Massive stars and active galactic nuclei (AGN) are often intense sources of Lyman-continuum (LyC, or hydrogen ionising) photons resulting in the formation of regions of ionised gas in their surroundings (e.g. Hii regions). The emission from these regions carries information about the physical conditions in the interstellar medium (ISM) as well as the source of the ionising photons themselves. Key properties that can be constrained include the star formation rate (e.g. Kennicutt & Evans, 2012; Wilkins et al., 2019), gas metallicity (e.g. Tremonti et al., 2004), temperature, density, dust content (e.g. Reddy et al., 2015), and the presence of an AGN (e.g. Baldwin et al., 1981). Nebular line emission also enables the accurate measurement of redshifts, and thus distances once a cosmological model is assumed.

While there has been extensive progress in observing line emission at low (e.g. Brinchmann et al., 2004) and intermediate (e.g. Steidel et al., 1996; Shapley et al., 2003) redshifts there are few direct constraints at high-redshift. This is predominantly due to lack of sensitive near-IR spectrographs, particularly at $>2\mu$m where the rest-frame optical lines lie at $z>4$, and the comparative lack of strong lines, other than Lyman-$\alpha$, in the rest-frame UV. The small number of detections at high-redshift come overwhelmingly from Lyman-$\alpha$ (e.g. Stark et al., 2010; Pentericci et al., 2011; Stark et al., 2011; Caruana et al., 2012; Stark et al., 2013; Finkelstein et al., 2013; Caruana et al., 2014; Stark et al., 2017) though there has now been a handful of detections of the [Civ]$\lambda 1548$ and [Ciii],Ciii]$\lambda 1907,1090$ lines (Stark et al., 2015a, b, 2017). The presence of extremely strong optical lines can also be indirectly inferred from their impact on broadband photometry (e.g. Schaerer & de Barros, 2010; Stark et al., 2013; Wilkins et al., 2013; Smit et al., 2014; Wilkins et al., 2016b; De Barros et al., 2019) yielding constraints now available up to $z\approx 8$ (De Barros et al., 2019).

While existing observational constraints in the EoR are limited this will soon change with the arrival of the Webb Telescope. Webb’s NIRSpec instrument will provide deep near-IR single slit, multi-object, and integral field spectroscopy from $\sim 0.7-5\mu$m, while the NIRISS and NIRCam instruments will, together, provide wide field slitless spectroscopy over a similar range. This is sufficient to encompass all the strong optical lines to $z\sim 7$ with [Oii] potentially accessible to $z\sim 12$. Webb’s mid-infrared instrument (MIRI) will provide mid-IR single slit, and slitless spectroscopy at $\lambda>5\mu$m, albeit at much lower sensitivity and thus will likely only detect line emission for the brightest sources in the EoR.

The existing direct and indirect constraints and the nearing prospect of Webb motivates us to produce predictions for the nebular emission line properties of galaxies in the EoR. In this paper, we combine the large BlueTides hydrodynamical simulation with photoionisation modelling to predict the nebular line properties of galaxies in the EoR, specifically ($z=8\to 13$). As part of this paper we also explore some of the photon-ionisation modelling assumptions including the choice of stellar population synthesis (SPS) model and initial mass function (IMF). These predictions can be used to optimise the design of Webb surveys prior to launch, targeting emission lines in the EoR. The observation of these lines will also provide a powerful constraint on the physics incorporated into galaxy formation models.

This work builds upon other recent efforts to model nebular lines using both simple analytical models (e.g. Charlot & Longhetti, 2001; Schaerer & de Barros, 2009; Inoue et al., 2014; Gutkin et al., 2016; Feltre et al., 2016; Nakajima et al., 2018) and the modelling in semi-analytical galaxy formation models (e.g. Orsi et al., 2014) and hydrodynamical simulations (e.g. Shimizu et al., 2016; Hirschmann et al., 2017; Hirschmann et al., 2019).

This article is structured as follows: in Section 2 we describe the BlueTides simulation and our methodology for calculating spectral energy distribution including nebular emission (§2.2). In Section 3 we present our predictions. In this section we also explore the impact of some modelling assumptions (§3.2), including the choice of stellar population synthesis (SPS) model (§3.2.1), initial mass function (IMF, §3.2.2), photoionisation modelling parameters, including the impact of dust (§3.2.3). In Section 3.4 we compare our predictions to existing observational constraints. In Section 4 we present our conclusions.

2 Modelling Nebular Emission in BLUETIDES

2.1 The BlueTides Simulation

The BlueTides simulation111http://bluetides-project.org/ (see Feng et al., 2015, 2016, for description of the simulation physics) is an extremely large cosmological hydrodynamical simulation designed to study the early phase of galaxy formation and evolution with a specific focus on the formation of the massive galaxies. BlueTides phase 1 evolved a $(400/h\approx 577)^{3}\,{\rm cMpc^{3}}$ cube to $z=8$ using $2\,\times\,7040^{3}$ particles assuming the cosmological parameters from the Wilkinson Microwave Anisotropy Probe ninth year data release (Hinshaw et al., 2013). The dark matter particle and gas particle initial masses are $1.2\times 10^{7}h^{-1}\ {\rm M_{\odot}}$ and $2.36\times 10^{6}h^{-1}\ {\rm M_{\odot}}$ respectively. The gravitational softening length is $\epsilon_{\rm grav}=1.5h^{-1}\ {\rm kpc}$. This is sufficient to allow us to easily resolve galaxies to $M_{\rm halo}\approx 10^{9}\,{\rm M_{\odot}}$ though in this work we adopt a more conservative approach only focussing on galaxies with stellar masses $>10^{8}\,{\rm M_{\odot}}$ which contain at least approximately 100 star particles in order to star formation history samplings effects (see Appendix B for exploration of sampling effects). The properties of galaxies in the simulation are extensively described in Feng et al. (2015, 2016); Wilkins et al. (2016a); Wilkins et al. (2016b); Waters et al. (2016a, b); Di Matteo et al. (2017); Wilkins et al. (2017, 2018). While BlueTides contains super-massive black holes, and even a small number of AGN dominated sources at $z=8$, in this work we focus on the line emission arising solely from gas excited by stellar sources.

2.1.1 Ages and Metallicities of Galaxies in BLUETIDES

As emission line luminosities and equivalent widths are predominantly driven by galaxy star formation and metal enrichment histories it is useful to explore the average ages and metallicities predicted by BlueTides. The mean stellar age and mean metallicity of young ($<10\,{\rm Myr}$) stars are shown as a function of stellar mass for a range of redshifts in Figure 1. These correlations were previously discussed in more detail in Wilkins et al. (2017) while a more detailed analysis of the joint star formation and metal enrichment history is presented in Fairhurst et al. in-prep. In short, the mean stellar age appears to show little dependence on mass but evolves strongly with redshift while the mean metallicity of young stars shows a power-law dependence ($Z\propto M^{0.4}$) on stellar mass but little evolution with redshift.

2.2 Spectral Energy Distribution Modelling

We model the spectral energy distributions (SEDs) of galaxies in BlueTides as the sum of the SEDs of each star particle identified as belonging to each galaxy.

We begin by associating each star particle with a pure stellar SED according to its age and metallicity. To obtain this SED we interpolate publicly available grids produced by stellar population synthesis (SPS) models. By default we make the following modelling choices: we assume the bpass v2.2.1 SPS model (Stanway & Eldridge, 2018; Eldridge et al., 2017)222https://bpass.auckland.ac.nz and a modified version of the Salpeter IMF containing a flattened ($\alpha=-1.3$) power-law at low-masses ($m<0.5\,{\rm M_{\odot}}$). This IMF produces very similar ($<0.05\,{\rm dex}$) results to the assumption of a Chabrier (2003) IMF but permits a fairer comparison with the alternative IMFs available for bpass. In §3.2 we explore the impact of these, and other assumptions.

2.2.1 Nebular emission

Using the intrinsic stellar SED, and assuming no escape of LyC photons, we use the cloudy photoionisation code (Ferland et al., 2017)333https://www.nublado.org to associate each young ($t<10\,{\rm Myr}$) star particle with an individual Hii region or birth cloud. The metallicity of this region is assumed to be identical to the star particle itself and we adopt the same interstellar abundances and dust depletion factors as described in Gutkin et al. (2016).

To model the nebular emission associated with a stellar population we adopt a similar approach to Charlot & Longhetti (2001) (see also Gutkin et al., 2016; Feltre et al., 2016). Like these works we choose characterise our photoionisation modelling using the density of hydrogen ($n_{H}$) and ionisation parameter at the Stromgren radius $U_{S}$. This is defined as,

[TABLE]

where $\epsilon$ is the effective gas filling factor.

We differ from previous approaches by parameterising models for the ionising spectrum in terms of an ionisation parameter defined at a reference age ($t=1\,{\rm Myr}$) and metallicity ($Z=0.02$) - $U_{S,{\rm ref}}$. Because of this the actual ionisation parameter passed to cloudy depends on the ionising photon production rate relative to the reference value, i.e.

[TABLE]

This ensures that the assumed geometry of the Hii region, encoded in the $\epsilon^{2}n_{H}$ term, is fixed for different metallicities/ages. By default we assume $\log_{10}(U_{S,{\rm ref}})=-2$ and $\log_{10}(n_{H}/{\rm cm^{-3}})=2.5$.

In our modelling we include the effect of dust grains which can not only boost certain lines van Hoof et al. (see 2004); Nakajima et al. (see 2018) but also provide an additional source of attenuation. Specifically we include cloudy’s default implementation of Orion-type graphite and silicate grains but scale the abundances to match those assumed for carbon and silicon in the Hii region. The impact of the inclusion of grains is discussed in more detail in §3.2.3.

In our calculations we assume the default cloudy stopping temperature (4000K) which is suitable for UV/optical recombination lines.

2.2.2 Modelling attenuation by dust in the ISM

BlueTides includes a simple model to account for dust in the wider intervening ISM. This ISM dust component is modelled using a simple scheme which links the smoothed metal density integrated along a single consistent line-of-sight to each star particle within each galaxy to the dust optical depth in the $V$-band ($550\,{\rm nm}$). Attenuation at other wavelengths is determined using a simple attenuation curve of the of the form,

[TABLE]

This model has a single free parameter which effectively links the surface density of metals to the optical depth. This parameter is tuned to recover the shape of the of the observed $z=8$ far-UV luminosity function. For a full description see Wilkins et al. (2017).

3 Predictions for BLUETIDES

Using the methodology outlined above we calculate the luminosities and equivalent widths of twelve prominent rest-frame UV and optical single lines or doublets (see Fig. 2 for a list of lines and their accessibility at high-redshift to Webb) for all galaxies in BlueTides at $z=8-13$ with $M_{\star}>10^{8}\,{\rm M_{\odot}}$. Equivalent widths (EWs) are calculated simply by dividing the line luminosities by the underlying continuum emission (which includes contributions from both transmitted starlight and nebular continuum emission).

3.1 Line Luminosities and Equivalent Width Distributions

Detailed diagnostic plots for each of the 12 calculated single lines or doublets are presented in Appendix C. Tabulated results for all 12 lines are also all available in electronic form at https://github.com/stephenmwilkins/BluetidesEmissionLines_Public. Fig. 3 provides a summary showing the median equivalent widths of all 12 lines at $z=8$ in bins of stellar mass.

A more detailed summary, concentrating on only 6 lines or doublets444Here we have combined both the [Neiii] and [Oiii] doublets into a single quantity., is presented in Fig. 4. Here we show predictions for both the luminosity function, and EW, $L/M_{\star}$, and $L/L_{\rm FUV}$ distributions as a function of redshift, stellar mass, and FUV luminosity.

The luminosity function (Fig. 4, row 1) of each line broadly follows a similar trend to the UV luminosity function: intrinsically the LF is approximated by a single power-law; the inclusion of dust however causes a strong break at high luminosities. Like the UV LF the line luminosity function evolves strongly with redshift, increasing by a factor $\approx 1000$ from $z=13\to 8$. At fixed stellar mass equivalent widths (Fig. 4, row 2) mostly increase to higher-redshift. This predominantly reflects that higher-redshift galaxies are generally younger (see §2.1.1) and thus their SED has a larger contribution from the most massive (LyC producing) stars. This also results in a wavelength dependence with the EWs of bluer lines evolving less strongly with redshift. The trend of line EW with stellar mass is more complex due to the correlation of stellar mass with metallicity (§2.1.1). For example, for this reason the EW of the hydrogen recombination lines drops at higher stellar mass while that of the [Oii]$3726,3729$ line peaks at $M^{\star}\sim 10^{9.5}\,{\rm M_{\odot}}$ (at $z=8$). The trend of EW with UV luminosity (Fig. 4, row 3) shows a similar trend with redshift. However, the trend with UV luminosity is less pronounced compared to with stellar mass due to the weaker correlation between observed UV luminosity and metallicity. The specific line luminosity ($L/M_{\star}$) (Fig. 4, row 4) shows a clear increase to higher-redshift, again this is driven by the fact that at higher-redshift the average age of the stellar populations are typically younger and thus produce more ionising photons per unit stellar mass. There is also a strong trend with stellar mass. While some of this is affected by metallicity it is predominantly dominated by the effect of dust. In contrast to the other quantities the ratio of line luminosity to FUV luminosity (Fig. 4, row 5) shows little evolution with redshift. This is because both ionising and FUV photons are dominated by the most massive stars. There is however a strong FUV luminosity dependence due to increased effect of dust.

3.2 Impact of Modelling Assumptions

These predictions depend not only simulation physics but the additional modelling assumptions made in §2.2. In this section we explore the impact of some of these assumptions.

3.2.1 Stellar Population Synthesis Model

The production rate of LyC photons and the shape of the ionising spectrum (or hardness) predicted for a given simple stellar population is sensitive to the a range of stellar evolution and atmosphere modelling assumptions and thus choice of stellar population synthesis (SPS) model (see §A.1.2 for more details). Fig. 5 shows the impact of changing the assumed SPS model on the line luminosities and equivalent widths predicted by BlueTides. Adopting the previous release (2.2) of bpass produces only relatively small changes ($<0.1\,{\rm dex}$) to the predicted line luminosities and equivalent widths. On the other hand adopting the Pegase.2 SPS model (and $m_{\rm up}=100\,{\rm M_{\odot}}$) produces a significant decrease in the luminosities and equivalent widths relative to our default model. For most lines luminosities drop by $\sim 0.5$ dex while equivalent widths drop by $\sim 0.3$. While some of this decrease can be attributed to the small upper-mass cutoff of the IMF most of the effect is attributed to wider modelling differences between Pegase.2 and bpass, in particular the inclusion of binary stars in the latter.

3.2.2 Initial Mass Function

The production rate and hardness of LyC photons are also affected by the choice of initial mass function (see §A.1.3 for more details). Fig. 5 shows the impact of assuming an alternative IMF on the average predicted line luminosities and equivalent widths as a function of redshift. Unsurprisingly, increasing the fraction of high-mass stars, either by extending the high-mass cutoff or flattening the high-mass slope results in increased line luminosities and equivalent widths. As continuum luminosities are also increased the impact on equivalent widths is smaller than on line luminosities themselves. It is also important to note that flattening the IMF in this way would also increase the far-UV continuum luminosities of galaxies breaking the otherwise good agreement with observations (see Wilkins et al., 2017). This could however be ameliorated by having more aggressive dust attenuation. The effect of steepening the slope produces the opposite effect. Steepening the IMF to this extent will also significantly decrease the far-UV continuum luminosities again breaking the good agreement with observational constraints. In this case the good agreement can not be recovered without more drastic changes to the simulation physics.

3.2.3 Photoionisation Modelling

The luminosity of each line is also sensitive the parameters encapsulating the geometry, density, excitation, and dust content/composition of the Hii region. The impact of the ionisation parameter, which effectively encodes the geometry of the region, and the hydrogen density are discussed in §A.2.1. As demonstrated in Fig. 17 the choice of these parameters can have a significant ($>\pm 0.5\ {\rm dex}$) impact on the luminosities and equivalent widths of metal lines.

As noted previously by default we include dust-grain physics. Within our model framework Fig. 18 shows the impact of grains on the emergent line luminosities as a function of metallicity for a constant burst of star formation while in Fig. 5 we show the effect of turning off grain physics on our overall results. Removing grains results in a boost of $\sim 0.05$ dex to both luminosities and equivalent widths.

3.3 Comparison to other models

Like this work, Shimizu et al. (2016) (S16) model the UV/optical line emission of galaxies in the EoR by combining photonionsation modelling with a cosmological hydrodynamical simulation albeit with several key differences, including the base simulation physics, choice of SPS model and IMF, photoionisation model, and wider dust model. Overall we find good agreement between our predictions and S16.

3.4 Comparison with existing observational constraints

As noted in the introduction there remain relativel few constraints (direct or otherwise) on optical/UV line emission at very high-redshift.

The majority of direct constraints come from observations of Lyman-$\alpha$ (e.g. Stark et al., 2010; Pentericci et al., 2011; Stark et al., 2011; Caruana et al., 2012; Stark et al., 2013; Finkelstein et al., 2013; Caruana et al., 2014; Stark et al., 2017). However, the Lyman-$\alpha$ line is resonantly scattered by the ISM/IGM significantly complicating its modelling (see Smith et al., 2019). For this reason we have omitted a detailed comparison with Lyman-$\alpha$ observations. We do however nevertheless make predictions for the Lyman-$\alpha$ properties including dust attenuation but not scattering by the ISM/IGM. These predictions are presented in Appendix C.

Recently Stark et al. (2015a) and Stark et al. (2017) have obtained constraints on the [Ciii],Ciii] doublet at $z=6-8$. These constraints are shown in Fig. 6 alongside predictions from BlueTides. The two faintest objects (A383-5.2, GN-108036) have EWs statistically consistent with the BlueTides predictions assuming our default modelling choices. However, were we to alternatively assume the Pegase.2 SPS model (see Fig. 5) our predictions would lie below both these observations, albeit without strong statistical significance given the small number of objects and large measurement uncertainties. The brightest object (EGS-zs8-1 Stark et al., 2017) lies $\sim 0.35$ dex above the median prediction for the same FUV luminosity. However, as this object is very bright this raises the possibility of a contribution from an AGN, which would raise the EW (see e.g. Nakajima et al., 2018). While BlueTides includes AGN their contribution is not included in this work but instead deferred to a future study.

Indirect constraints on the strength of the strongest optical lines are possible through the effect of these lines on broad-band photometry (e.g. Schaerer & de Barros, 2010; Stark et al., 2013; Wilkins et al., 2013; Smit et al., 2014; Wilkins et al., 2016b; De Barros et al., 2019). De Barros et al. (2019) recently combined Hubble and Spitzer observations probing the rest-frame UV and optical to constrain the prominent H$\beta$ and [Oiii]$\lambda$4959,5007 lines at $z\approx 8$. As shown in Fig. 7 the H$\beta$ $+$ [Oiii] EW distribution measured by De Barros et al. (2019) has an almost identical median to that predicted by BlueTides for our default assumptions. However we do fail to explain the tail of very-high ($>2000{\rm\AA}$) and low EW sources. As noted in A.2.1 the [Oiii]$\lambda$4959,5007 lines are sensitive to the choice of ionisation parameter. If we instead of a single reference ionisation parameter we chose a distribution this would naturally produce extended tails. It is also worth noting that if instead we adopted the Pegase.2 SPS model our predictions would fall below these observational constraints (though this could be ameliorated by assuming a more high-mass biased IMF). Similarly, adopting a model in which line emission is more strongly attenuated by dust would only be consistent by also changing the IMF.

Unfortunately this good agreement is not seen in the luminosity function, as shown in Fig. 8. The observed luminosity function is systematically offset to higher luminosities ($\sim 0.4\,{\rm dex}$) or higher space densities ($\sim 0.7\,{\rm dex}$) than that predicted by BlueTides. The cause of this discrepancy appears to lie in the relation between the combined line luminosity and the far-UV luminosity, which is used by De Barros et al. (2019) to convert the observed FUV luminosity function to a line luminosity function. The individual values measured by De Barros et al. (2019) for this are shown in Fig. 9 and are compared to the values predicted by BlueTides. The measured values are on average $\sim 0.4\,{\rm dex}$ higher than predicted by BlueTides. As many of the measured values are above the intrinsic expectation (see Fig. 10) one interpretation of this discrepancy is that De Barros et al. (2019) measure higher dust attenuations than predicted by BlueTides. It is also possible that differences between the measured and predicted values of the metallicity, age, ionisation parameter, and hydrogen density can have an effect. Given the limited observational constraints such differences may not be surprising considering the range of degeneracies present. As direct emission line measurements become available from the Webb Telescope and other upcoming facilities the cause of this discrepancy should become clearer.

4 Conclusions

Using the large cosmological hydrodynamical simulation BlueTides combined with photoionisation modelling we have made detailed predictions for the luminosities (including luminosity function) and equivalent widths of twelve prominent rest-frame UV and optical emission lines for galaxies with $M^{\star}>10^{8}\,{\rm M_{\odot}}$ across the EoR ($8<z<13$). As part of this analysis we also explored the impact of various modelling assumptions including the choice of stellar population synthesis model, initial mass function, and photoionisation modelling, finding that these can have a significant impact.

At present there are few observational constraints on line emission available in the EoR with only a handful of direct constraints on the [Ciii],Ciii] doublet along with indirect constraints on H$\beta$ and [Oiii]$\lambda$4959,5007 based on Spitzer photometry (De Barros et al., 2019). Overall the agreement with these observations is mixed with good agreement with the H$\beta$ + [Oiii]$\lambda$4959,5007 equivalent width distribution but with the observationally inferred line luminosity function offset to higher luminosities or space densities. One possible explanation to this discrepancy is that De Barros et al. (2019) measured higher FUV dust attenuation than predicted by BlueTides.

With the arrival of the Webb Telescope it will be possible to obtain direct measurements of individual line luminosities and equivalent widths for a large range of galaxies at $z\sim 8$ and beyond. Combined with other observational constraints this will allow us to test not only the base simulation but also the assumptions involved in modelling the nebular emission.

Acknowledgements

We acknowledge support from NSF ACI-1036211, NSF AST-1009781, and NASA ATP grants NNX17AK56G and 80NSSC18K1015. The BlueTides simulation was run on facilities at the National Center for Supercomputing Applications. TDM acknowledges support from Shimmins and Lyle Fellowships at the University of Melbourne.

Appendix A The Impact of Photoionisation Modelling Assumptions

In this section we describe, within the context of our nebular emission model, the impact of various assumptions on the hydrogen ionising (Lyman-continuum, LyC) photon production and subsequent line emission.

A.1 The Production of Ionising Photons

Young, massive, hot stars produce LyC photons. The photons can be reprocessed by surrounding gas (and dust) into nebular continuum and line emission. Consequently the production rate of these photons by a stellar population is sensitive to the joint distribution of stellar ages and metallicities. This is demonstrated in Fig. 11, where we the LyC production rate $\dot{n}_{\rm LyC}$ as a function of age for a range of different metallicities assuming our default choices of SPS model and IMF is shown. The production rate drops rapidly after the first few million years at higher ages and metallicities, declining by $\sim 10-100$ as the population ages from $1\,\to 10\,{\rm Myr}$ and then again by a factor of $\sim 100$ from $t=10\,\to 100\,{\rm Myr}$. At young ages ($<20$ Myr) the lowest metallicity populations can produce up-to $10$ times as many LyC photons, though at later times this trend reverses. The overall difference in the production of LyC photons as a function of metallicity is summarised in Fig. 12 where we show the total number of LyC photons produced by an SSP from $t=0\to 10^{7}\,{\rm yr}$ and $t=0\to 10^{8}\,{\rm yr}$. The lowest metallicity modelled ($Z=10^{-5}$) SSP considered produces approximately double the number LyC photons over its lifetime compared to $Z=0.01$.

A.1.1 The Ionising Photon Hardness

More complex atoms have a range of potential ionisation states each excited by photons of different energies. For example, helium can be singly ionised by photons with $E_{\gamma}>24.6\,{\rm eV}$ and doubly ionised by those with $E_{\gamma}>54.4\,{\rm eV}$. For this reason it is useful to also consider the ionising photon hardness, essentially a ratio of the number of more energetic photons to $\dot{n}_{\rm LyC}$. The left panel of Fig. 13 shows the hardness the LyC by comparing the number of LyC and Oii ionising ($>35.1\,{\rm eV}$) photons as a function of age for two metallicities $Z\in\{0.02,0.004\}$. At the youngest ages the higher metallicity stellar population produces significantly fewer ($\sim 1$ dex) harder photons. For older ($>10\,{\rm Myr}$) populations the hardness is similar. The impact of this will be line ratios that vary as a function of the age of the ionising stellar population.

A.1.2 The Effect of SPS Model Choice

The number of LyC photons and the shape of the ionising spectrum predicted for a given stellar population is also sensitive to the a range of stellar evolution and atmosphere modelling assumptions and thus choice of stellar population synthesis (SPS) model (see also Wilkins et al., 2019). The middle panel of Fig. 14 shows a comparison between $\dot{n}_{\rm LyC}$ for different SPS models/versions; these include the three most recent versions of bpass (v2.2.1, v2.2, v2.1) and the Pegase.2 model (Fioc & Rocca-Volmerange, 1997). This analysis reveals relatively small differences between the different bpass versions but larger differences between bpass and Pegase.2. This difference is particularly acute at ages $>5$ Myr where the LyC production rate predicted by Pegase.2 drops of much more rapidly than in bpass. The left panel of Fig. 14 shows the difference in the hardness between the default model and bpass; again, the most notable feature is the difference at $>5$ Myr.

A.1.3 The Effect of the Choice of IMF

Both the production rate and hardness are also affected by the choice of IMF. The right hand panel of Fig. 14 shows the production rate relative to our default model for several different high-mass slopes $\alpha\in\{2.0,2.35,2.7\}$ and for a lower ($100\,{\rm M_{\odot}}$) high-mass cut-off. Assuming a shallower slope ($\alpha=2.0$) yields more around double the number of LyC photons overall with the enhancement decreasing with age. Assuming a steeper slope has the opposite effect albeit with a slightly larger magnitude. Adopting a lower high-mass cut-off reduces the number of LyC photons produced at the youngest ages, overall leading to around $\sim 30-50\%$ less LyC photons produced by the SSP over its lifetime, depending on the metallicity.

A.2 Photoionisation Modelling

Using the modelling procedure described above we now make predictions for line luminosity and equivalent widths. We concentrate here on 12 prominent UV and optical lines. In making these predictions we assume a constant star formation history with fixed metallicity.

Fig. 15 shows the predicted line luminosities (per unit stellar mass) and equivalent widths (EWs) for a range of prominent rest-frame UV and optical emission lines as a function of metallicity. In both cases we assume continuous star formation for 10 Myr. The luminosity of the hydrogen lines largely track the change in the LyC production rate with metallicity with the luminosity dropping by $\sim 0.5$ dex over the metallicity range considered. The non-hydrogen lines exhibit more complicated behaviour with an increase to $Z\sim 10^{-2.5}$ before declining to higher metallicities. The rapid increase broadly reflects the increasing abundance of each element in the ISM while the drop at high metallicities reflects the decline in the number of suitably energetic photons. The metallicity dependence of the EW of each line exhibits a similar behaviour, albeit often with reduced magnitude.

The equivalent width of any line is also sensitive the star formation history of the stellar population. Fig. 16 shows how the equivalent width of H$\alpha$ varies with the duration of continuous star formation for a range of metallicities. The equivalent width declines from $\sim 1000-2000{\rm\AA}$ after to 10 Myr to $\sim 150-300{\rm\AA}$ after 1 Gyr.

A.2.1 The Effect of Photoionisation Modelling Assumptions

In addition to the choice of SPS model and IMF the strength of lines are also senstive to the parameters encapsulating the geometry, density, and excitation of the H ii region in addition to the presence of dust grains. Fig. 17 demonstrates the effect on changing both the reference ionisation parameter $U_{S,{\rm ref}}$ and hydrogen density $n_{H}$ on the strengths of the same 12 prominent lines considered previously. While the hydrogen lines are largely insensitive to these choices many of the other lines, and in particular line ratios, are strongly sensitive with the effect dependent on the metallicity. The incusion of dust-grains in the H ii region not only provide an additional source of LyC photon attenuation but also play a role in photoelectric heating of the gas. Fig. 18 shows the predicted line luminosities when grains are omitted from the model. Omitting grains generally results in stronger lines with the effect being strongly metallicity dependent.

Appendix B Star Formation History Sampling Effects

The Lyman continuum photon production rate is a strong function of the age, and to a lesser extent, metallicity of the stellar population. As the star formation history of each galaxy is sampled, at the lowest masses considered, by a small ($\sim 100$) number of individual star particles, this raises the possibility that the predicted line properties differs from the truth because of SFH sampling affects. To gauge the impact of this effect we re-sample the average star formation history of galaxies at $z=8$ using different numbers of fixed mass star particles ($n=10^{2}-10^{4}$); corresponding roughly to stellar masses of $10^{8}-10^{10}\,{\rm M_{\odot}}$. The result of this analysis is shown in Fig. 19. This test reveals that there is no significant bias in the average (median or mean) of the predicted line luminosity (in this case H$\alpha$), even at the lowest masses considered in this study. However, at low-masses there is some scatter ($\approx 0.1$ dex for $n=300$ particles / $M_{\star}\approx 2.5\times 10^{8}\,{\rm M_{\odot}}$). Because of the steepness of the luminosity function this will have the effect of flattening the LF at the lowest-luminosities.

Appendix C Detailed Predictions for Individual Lines

In Figs. 20-23 we show the EW and luminosity as a function of stellar mass and observed UV luminosity for each line in addition to the luminosity function. Tabulated values for each line and redshift are available at https://github.com/stephenmwilkins/BluetidesEmissionLines_Public.