© ESO 2019

1. Introduction

Ionizing photons, those emitted at wavelengths shortwards of the hydrogen Lyman Limit at 912 Å, have played a key role in shaping the evolution of the intergalactic medium and the galaxies which form and interact within it. In the local Universe, the emission of such photons originates from a mixture of sources. Amongst these, two are most prominent: sources powered by accretion onto compact objects and supermassive black holes, and the emission from the photospheres of hot, massive stars. In the distant Universe, massive active galactic nuclei (AGN) are rarer. Their number density drops far more rapidly than that of star forming galaxies and hence star formation dominates the ionizing photon background (Robertson et al. 2015). Lyman continuum (hydrogen-ionizing) photons escaping from star forming galaxies were responsible for the process of ionizing the intergalactic medium at redshifts of z ≳ 6, and also have an important impact on the thermal conditions within, and emission spectra observed from, the galaxies themselves. As a result there has been strong interest in recent years both in the ways in which ionizing photons leak from their origin galaxies, and in how these photons are produced in the first place (e.g. Rosdahl et al. 2018; Stanway et al. 2016; Wilkins et al. 2016; Ma et al. 2015, 2016; Nestor et al. 2013).

An important development in this field has been the development of population synthesis models which account for the effects of binary interactions on stellar evolution pathways. When these are combined with stellar atmosphere models describing the spectral energy distribution of individual stars, the resulting spectral synthesis yields predictions that can be compared to the observed spectral properties of stellar populations (Eldridge & Stanway 2009, 2012; Stanway et al. 2016). This has led to an increasing recognition of the importance of massive stars in young stellar populations for reionization, twinned with an appreciation of the uncertainties in their evolution. In particular, both rotation (Topping & Shull 2015) and stellar multiplicity (Sana et al. 2012; Eldridge et al. 2017) have a strong effect on massive star evolution, and the impact of stellar binaries or multiples may be still more important at low metallicities (Moe et al. 2018). Binary star evolution can both boost the ionizing photon flux from a stellar population (Stanway et al. 2016) and extend its lifetime (Stanway et al. 2014; Xiao et al. 2018). Crucially this latter property may allow for increased ionizing photon escape fractions, as supernovae from the most massive stars blow channels through the circumstellar and circumgalactic medium, allowing ionizing photon leakage (Rosdahl et al. 2018).

Improved spectroscopy has made analysis of the ionizing spectra of distant galaxies possible in recent years, primarily through study of those ionizing photons reprocessed by nebular gas in the interstellar environment. With the detection of multiple, high signal-to-noise absorption and emission features in galaxies at z ∼ 2 − 4 (e.g. Du et al. 2018; Steidel et al. 2016; Maseda et al. 2017), and even higher redshifts in rare cases (e.g. Stark et al. 2017), it is now possible to construct photoionization models which yield insight into the physical conditions in these star-forming systems. While there is always some degeneracy between assumed gas conditions and the spectrum of the ionizing potential, a number of distant systems show evidence for the presence of a very hard ionizing radiation field. This is observed both in (rest-frame) optical line ratios such as [O III λ5007 Å]/[O II λ3727 Å] or [O III λ5007 Å]/H β λ4681 Å, and directly in the strength of nebular emission lines including C III] λ 1909 Å (which requires a 48 eV ionization potential) and He II λ 1640 Å (at 54.4 eV) in the ultraviolet (Du et al. 2018; Rigby et al. 2015; Erb et al. 2018, 2010; Berg et al. 2018) Interestingly, the class of low redshift extreme emission line galaxies and Lyman-continuum leakers identified as analogues to galaxies in the distant Universe also show hard ionizing spectra (e.g. Izotov et al. 2018).

Given the association of these hard spectra with intensely star forming galaxies and the absence of clear indications for AGN, a natural question to ask is whether the photospheric emission from stars can provide sufficient ionizing photons, with appropriate energies. The early indications have indicated mixed answers to this question. While stellar populations including interacting binaries have performed well in reproducing the mean properties of galaxies in stacks of high redshift spectra (Steidel et al. 2016), and certainly produce a harder radiation field than a single star population (Eldridge & Stanway 2012; Eldridge et al. 2017; Xiao et al. 2018) and sustain this for significantly longer timescales (Stanway et al. 2014, 2016), reproducing the highest ionization potential lines remains a challenge. In particular, He II emission (in either the ultraviolet λ1640 Å or optical λ4686 Å lines) appears to be consistently underestimated given default stellar population synthesis model assumptions (e.g. Erb et al. 2018; Berg et al. 2018).

Evolutionary population and spectral synthesis - the complex process of creating a model spectrum for a stellar system of known age, mass and metallicity based on a combination of individual stellar models - involves a number of approximations and assumptions. At its simplest, it requires a grid of stellar evolution models, a grid of stellar atmosphere models to predict observable properties of a star at given temperature and surface gravity, and a function which provides the weighting with which to combine models of different initial mass. More complex cases may include evolutionary tracks which include the effects of stellar rotation, Roche-lobe overflow, common envelope and other binary interactions, and expanded grids of atmosphere models which account for surface composition and stellar winds. For the short-lived massive stars that dominate the ionizing spectrum, the binary fraction is close to unity and the effects of binary interactions can be dramatic (see e.g. De Donder & Vanbeveren 2004; Sana et al. 2012; De Marco & Izzard 2017; Eldridge et al. 2017). Including these pathways, given an assumed set of mass-dependent binary parameter (period, mass ratio) distributions, is thus essential to accurately predict the ionizing spectrum of high redshift galaxies.

However a more fundamental assumption than the binary parameter assumptions is the underlying initial mass function (IMF). This describes the probability distribution of stellar masses (in the binary case, of primary star masses) at zero age in a starburst, before stellar evolution has any significant effect. For many years, the assumption of a “universal” mass function has been used to simplify modelling of galaxies. This is usually described as a power law function in the form dN/dM ∝ M^α where N is the number of stars in a population at a given mass M and the function extends between 1 and 100 M_⊙. The most commonly used prescription, based on that of Salpeter (1955), uses α = −2.35. However, it has become increasingly clear in recent years that this oversimplification is untenable. As the extensive reviews in Hopkins (2018) and Kroupa & Jerabkova (2019) discuss in detail, not only is it necessary to break the power law to avoid overpredicting sub-Solar mass stars, but the power law slope at high masses is also highly uncertain. There is evidence for local variations in the IMF and potentially also some metallicity dependence. Studies of the massive (M > 8 M_⊙) and very massive (M > 100 M_⊙) star population in the 30 Doradus region have suggested that the IMF in the M > 30 M_⊙ regime may be as shallow as $\alpha =-{1.9}_{-0.26}^{+0.37}$ (Schneider et al. 2018a,b). Indirect measurements of the IMF in the extragalactic regime also suggest a shallow IMF may be appropriate (e.g. Wilkins et al. 2008), and similar shallow functions have been seen in individual local massive galaxies, while other studies have suggested steeper, “bottom heavy” IMF in elliptical galaxies (albeit in the low stellar mass regime, Parikh et al. 2018; Vaughan et al. 2018)

Determining the effects of the plausible range of IMF variation based on local observations is crucial to deciding whether an additional source of hard ionizing photons is necessary to explain the observed results from intense starbursts. A shallow IMF slope will lead to a population that is “top-heavy”, i.e. has more massive stars per unit total starburst mass than expected from a “Salpeter”-like slope. A similar, although subtly different, uncertainty arises from the upper mass cut-off in a stellar population model. While most IMFs do not consider stars with M > 100 M_⊙, identification of the very massive star population in the Large Magellanic Cloud (Crowther et al. 2016) has led to suggestions that the IMF may extend as far as 300 M_⊙ in intense starburst regions. Given the dependence of the hard ultraviolet spectrum of a population on its hottest stars, the inclusion or neglect of even a small population of very massive stars can have a dramatic effect on the photons capable of ionizing not only hydrogen but also helium. It is thus both timely and necessary to explore the photon flux generated above both 13.6 eV and 54.4 eV by young stellar populations, as a function not only of metallicity but also of the assumed IMF in a binary stellar population synthesis model.

Is it also instructive to consider the ionizing spectrum at still higher energies. The far-ultraviolet O VI λ1032,1038 Å doublet (with an ionization potential of 138 eV) is frequently seen in absorption in the halos around low redshift star forming galaxies, and has thus been identified both as sensitive to specific star formation rate in galaxies and as a potential tracer of reservoirs of warm ionized gas in the circumgalactic environment (Tumlinson et al. 2011). In this case the primary ionization source is the extragalactic ionizing background, which is largely powered in turn by emission from AGN, while the association with star formation is provided by the necessity for large scale starburst-driven outflows to account for the density, velocity profile and composition of the absorbing gas (e.g. Muzahid et al. 2015). Indeed ionization powered by accretion onto AGN, whether indirectly through the extragalactic background, or directly through short-duty-cycle bursts of activity within the galaxies in question, can likely account for many low redshift O VI absorption features (e.g. Tripp et al. 2008; Oppenheimer et al. 2018; Nelson et al. 2018).

However there are signs that this interpretation may not apply universally. Analysis of O VI absorbers in the haloes of L* galaxies suggests that ionization by the extragalactic background is strongly disfavoured, and either additional local sources of photoionization or mechanisms such as collisional or shock excitation are required (Werk et al. 2016; Turner et al. 2015). Observations at high redshift, where the extragalactic background is softer due to evolution in the AGN luminosity function, again imply additional sources of hard photons or collisional mechanisms (Chisholm et al. 2018), while there are also now a handful of examples of intense starbursts which show the ultraviolet O VI lines in emission (e.g. Hayes et al. 2016).

There is general consensus that photoionization models driven by emission from hot stars are unable to produce sufficient photons at these energies to explain observations. Nonetheless, the hints of hard ionizing spectra seen in distant objects suggest that we should not rule out this possibility without exploring the limits of the high mass stellar population given plausible IMF models. While observations of the O VI doublet remain challenging at both high and low redshift, it nonetheless provides a useful metric for the hardness of the ionizing spectrum and so will also be considered here.

In this paper, we use a custom set of binary population and spectral synthesis (BPASS) models to explore the effects of IMF assumptions on the hard ionizing spectra produced by stellar population synthesis as a function of metallicity. In Sect. 2 we introduce our methodology. In Sect. 3 we explore the effects of assumed IMF slope and the upper mass cut-off in the massive star population. In Sect. 4 we introduce a second break in the IMF to shallow slopes at very high mass. In Sect. 5 we attempt to quantify the effects of total starburst mass on ionizing photon flux. Finally in Sect. 6 we discuss the implications of this study and summarise our conclusions in Sect. 7.

2. Method

2.1. Population synthesis

The BPASS v2.2 data release (Stanway & Eldridge 2018) includes calculation of results for nine IMF prescriptions. Of these, one is an unbroken power law and two follow the low mass cut-off prescription of Chabrier (2003) with different upper mass cut-offs. The remainder comprise three pairs of broken power law models with an IMF slope of −1.3 between 0.1 and 0.5 M_⊙ and −2.0, −2.35 or −2.7 between 0.5 M_⊙ and an upper mass cut-off of either 100 M_⊙ or 300 M_⊙. Calculating such a range of models already represents a substantial investment of resources, but nonetheless this sparse model grid is inadequate to fully explore the effects of IMF on the ionizing spectrum.

In order to make such an exploration tractable we modify the BPASS population synthesis to truncate populations at an age of 1 Gyr. This removes the utility of the models as predictors of old stellar populations and astrophysical transient rates, but does substantially speed up the calculation of models with different input IMF parameters. We also limit our analysis to a subset of five metallicities, namely Z = 10⁻⁵, 10⁻⁴, 0.001, 0.010 and 0.020, where Z_⊙ = 0.020 is the metal mass fraction of our Sun. Apart from this, we make no changes to the underlying BPASS v2.2 assumptions (e.g. stellar models, atmosphere models, mass-dependent binary parameters, mass loss rates, mass transfer prescription or synthesis procedure).

For the purposes of this analysis we adopt a twice-broken power law for the IMF, of the form:

$$\frac{\mathrm{d}N(M)}{\mathrm{d}M}=A_{\mathrm{l}}{\displaystyle {\int }_{M_0}^{M_{\mathrm{l}}}{M}^{{\alpha }_{\mathrm{l}}}}\mathrm{d}M+A_{\mathrm{m}}{\displaystyle {\int }_{M_{\mathrm{l}}}^{M_{\mathrm{m}}}{M}^{{\alpha }_{\mathrm{m}}}}\mathrm{d}M+A_{\mathrm{u}}{\displaystyle {\int }_{M_{\mathrm{m}}}^{M_{\mathrm{u}}}{M}^{{\alpha }_{\mathrm{u}}}}\mathrm{d}M,$$

where α_l is the low mass slope between M₀ and lower break mass M_l, α_m is the intermediate mass slope valid between M_l and intermediate break mass M_m and α_u is the upper mass slope which is effective between M_m the upper cut-off mass limit M_u. The constants A_l, A_m and A_u ensure the function is continuous and normalised to yield a total population of 10⁶ M_⊙. In the case where α_m = α_u the break mass M_m becomes arbitrary and the form collapses to a simple broken power law.

In Table 1 we summarise both the standard BPASS IMFs and the parameter ranges explored in this work, according to this formulation. In all cases we fix the lower mass limit and power law slope. The stellar binary fraction implemented in BPASS v2.2 is mass-dependent (based on that of Moe & Di Stefano 2017). Hence a consequence of varying the IMF is that the mass fraction of single stars versus binaries (primary and secondary models) is also IMF dependent. In Fig. 1 we show how the single star fraction varies as a function of both IMF slope and upper mass cut-off. Note that this is independent of metallicity in the current implementation.

Fig. 1. Effect of IMF slope on binary star fraction in a population of total mass 106 M⊙, given the stellar mass-dependent binary fractions implemented in BPASS v2.2. The lower region between the x-axis and the first pair of lines indicates the fraction of the total stellar mass for a given IMF slope αu which is assigned to single star models, the central region the mass fraction in primary star models and the upper region the mass fraction assigned to their less massive secondary companions. Regions between each pair of lines (with hatching in two styles) switch between stellar model types depending on the IMF. Models assume a single power law break (αm = αl), The upper cut-off mass, Mu, is shown by line style (100 M⊙ - solid; 300 - dot-dash). The curves for Mu =150 & 200 M⊙ lie between these curves.

2.2. Ionizing spectrum parameterisation

We define the photon flux generated by each model as:

$$Q(X)={\displaystyle {\int }_{1\mathrm{\AA }}^{{\lambda }_c}\frac{\lambda L_{\lambda }}{hc}}\mathrm{d}\lambda ,$$

where the critical wavelength λ_c = 912 Å for Q(H I), 227.9 Å for Q(He II) and 89.79 Å for Q(O VI), equivalent to the relevant ionization potentials of 13.6, 54.42 and 138.12 eV respectively. The resultant fluxes are quoted in units of emitted photons per second.

These are calculated at each time step, for each output spectrum generated by our stellar population synthesis models. We consider photon fluxes rather than directly studying emission or absorption line profiles in the spectra for two reasons:

Firstly, the flux in spectral features associated with these transitions (e.g. Lyα λ 1216 Å, Hα λ 6563 Å, He II λ 1640 Å, He II λ 4686 Å, O VI λ 1035 AA) is likely to be modified by absorption and re-emission in the nebular gas surrounding young star forming regions. Interpretation of line strengths thus requires an extensive set of assumptions regarding the conditions, covering fraction and composition of the nebular gas screen.

Secondly, we wish to avoid possible uncertainties introduced into stellar absorption lines by the modelling of stellar atmospheres for rare subgroups within the stellar population. BPASS implements a custom grid of OB star atmospheres that are required by binary interaction products (Eldridge et al. 2017) and also the Potsdam Wolf-Rayet (PoWR¹, Sander et al. 2015) atmospheres for stripped helium stars. The PoWR models are extrapolated to lower metallicities than they were created at, which may lead to overprediction in certain metal lines, and modification of the detailed strengths of lines in hydrogen and helium. Due to technical constraints, the current version of BPASS also does not make use of the helium dwarf atmospheres generated by Götberg et al. (2017, 2018) and these may have a substantial impact on the strength of some emission lines.

As a result, we opt to report directly the ionizing flux generated by a starburst rather than studying the detailed outputs that result from our spectral synthesis.

3. Effects of IMF slope and upper mass limit

As shown in Table 1, we expand the standard BPASS IMF grid with power law models with single breaks at 0.5 M_⊙ and α_m = α_u in the range −1.8 to −2.7 at intervals of 0.1. We calculate each of these with four possible upper mass cut-offs: 100, 150, 200 and 300 M_⊙, and for each of the output spectra calculate the time evolution of the ionizing photon flux from an instantaneous burst. As Fig. 2 demonstrates, Q(H I) is largely independent of assumed IMF upper mass limit, and mildly sensitive to the upper mass slope, which changes the relative fraction of low mass stars with respect to their massive siblings.

Fig. 2. Time evolution of photon production rates in H I (black), He II (blue) and O VI (yellow-green) at five metallicities, from an instantaneous burst of total mass 106 M⊙ occurring at an age of zero. At each metallicity three single-break power law (αu = αm = −1.8, −2.35, −2.70) IMF models are shown at two upper mass limits (Mu = 100, 300 M⊙), with line style indicating different slopes and line colour different upper mass cut-offs as indicated in the key. The rest of the models lie within the bounds of those shown at fixed age.

By contrast, Q(He II) shows sensitivity to M_u, particularly at low metallcity and log(age/years) = 6.3−6.5. At these ages very massive stars are approaching the end of their life, swelling to become very luminous giants and (when stripped by winds or binary interactions) Wolf-Rayet stars. The presence or absence of these stars is strongly dependent on the upper mass cut-off.

The photon flux shortwards of 89.8 Å (138 eV) capable of powering O VI emission and absorption, Q(O VI), is also sensitive to the presence or absence of these massive stars, both in terms of the peak flux attained by the population and its timing,as well as to metallicity. Stellar populations including very massive stars produce higher Q(O VI) and at substantially later ages than those which lack very massive stars. At high metallicities the strong stellar winds exhibited by very massive stars prevent them from ever becoming sufficiently hot or luminous to produce a significant Q(O VI) flux.

In Fig. 3 we consider in more detail the effect of varying metallicity, IMF slope and upper mass limit in each photon flux. Here we show the peak flux obtained by the population in a given line, which is usually seen in a single time bin at log(age/years) = 6.0−6.5, although we caution that the precise peak time may vary model by model. The strongest influence on hydrogen ionizing flux, Q(H I) appears to be the IMF slope, α_u. Over the reasonable range of slopes explored here, peak Q(H I) varies smoothly by ∼ 1.2−1.3 dex at each metallicity considered. By contrast, peak Q(H I) only varies by ∼ 0.3−0.4 dex when the slope is fixed and the upper mass limit is altered. Overall, the range of peak fluxes is confined to a range of log(Q(H I)) = 51.9–53.7 for our 10⁶ M_⊙ populations, narrower than that seen in Q(He II) and Q(O VI).

Fig. 3. Peak ionizing photon fluxes reached by an ageing instantaneous starburst of total mass 106 M⊙ as a function of IMF slope (assuming αm = αu). In each panel, metallicity is indicated by colour while upper mass cut-off Mu is indicated by line style (100 M⊙ = solid, 300 M⊙ = dot-dash).

While the metallicity dependence of peak values of Q(He II) is perhaps the most visible feature of Fig. 3 (middle panel), this flux also shows more dependence on IMF slope (∼ 1.2–1.6 dex at fixed M_u) than upper mass cut-off (∼ 0.6–1.2 dex at fixed α_u). The photon flux from the highest metallicity models (Z = 0.020) is both very low and relatively stable to IMF variations.

The IMF dependence of peak Q(O VI) flux is also governed by metallicity. At the lowest metallicities considered here, IMF slope has a stronger effect than cut-off mass (1.6 dex compared to 0.8 dex). However at intermediate metallicities, around Z = 0.001, the two effects are comparable (∼ 1.6 dex each). At the higher metallicities, Z > 0.010, the ionizing flux is both very low and not dependent on IMF. In these cases, as Fig. 2 shows, Q(O VI) is rising at late times and is generated by relatively low mass helium dwarfs, formed through binary interactions and largely independent of the very massive star population (see e.g. Götberg et al. 2018).

While models with upper mass limits M_u = 100 and 300 M_⊙ provide the bounding cases at fixed power law slope and metallicity, we also consider two intermediate upper mass cut-offs at M_u = 150 and 200 M_⊙. We show the variation in peak fluxes as a function of upper mass limit and metallicity for a fixed IMF slope (α_m = α_u = −2.35) in Fig. 4. The peak H I ionizing flux increases smoothly with upper mass cut-off at all metallicities considered here. At most metallicities, the same is true of peak He II flux and O VI flux but anomalous behaviour emerges in our models with Z ∼ 0.001. At this metallicity (5% of Solar) there is a very strong dependence of hard ionizing flux on upper mass limit which results from a combination of physical effects. These models have sufficient metal line opacities in their stellar atmospheres to drive mass loss through stellar winds, but at rates sufficiently low to allow angular momentum to be retained and rotational mixing to strongly influence the evolution of stars. In BPASS models, this chemically homogenous evolution is deemed to occur if significant mass is accreted by stars with Z < 0.004 and M > 20 M_⊙. The result is that where the upper mass cut-off is as low as 100 M_⊙ the majority of massive stars quickly drop below the BPASS mass limit for onset of chemically homogenous evolution and very few rotationally mixed stars are formed in our models. By contrast, when M_u = 300 M_⊙ a reasonable amount of mass loss can occur through winds, while still leaving stars with sufficient mass for rotational mixing to influence their evolution. Since chemically homogenous stars tend to burn hotter and for longer than a non-rotating star of comparable mass, these stars dominate emission at wavelengths capable of ionizing He II. The result is a strong dependence on M_u in the hard ionizing flux at Z = 0.001 relative to those at higher metallicities (where stellar winds always dominate and fewer stars are stripped) and those at the lowest metallicities. At Z ≤ 10⁻⁴, stellar winds are very weak due to the low metal ion opacities in the stellar atmospheres, and very little of the initial stellar mass is lost to winds. These stars are more likely than those at higher metallicities to be rotationally mixed, but those at the lowest metallicity also tend to become inflated to larger stellar radii which increases the probability that their atmospheres may become stripped by a binary companion (and so again create a hot, blue spectrum). The detailed balance of these different physical mechanisms for producing hot stars leads to a somewhat counterintuitive inversion in the metallicity trends of Q(He II) and Q(O VI) with the Z = 0.001 models producing more hard ionizing photons than those at Z = 10⁻⁴.

Fig. 4. Peak ionizing fluxes, colour coded by metallicity, as in Fig. 3, but now for single-break power law initial mass functions with a varying upper mass limit, Mu. We fix αm = αl = −2.35. The high metallicity fluxes in the He II and O VI ionizing photons are too low to show here, but follow a similar trend.

4. Introducing a second break

Studies of high mass star forming regions in the local Universe, notably the 30 Doradus complex in the Large Magellanic Cloud, have suggested that the IMF for massive stars may be substantially shallower than the Salpeter-like IMF that appears common for near-Solar mass stars (Schneider et al. 2018a,b). A possible interpretation of these findings is that the power-law IMF may show a second break in the high- to very-high-mass regime.

To explore the effects of such an IMF, we calculate a limited set of models with two power law breaks. We fix M₀ = 0.5 M_⊙, M_u =300 M_⊙, α_l = −1.3, α_m = −2.35 and α_u = −2.0 for consistency with constraints from local studies, while allowing the intermediate break mass, M_m, to vary. Again, we look at the peak photon fluxes attained by the population as a function of metallicity, and break mass, as shown in Fig. 5.

Unsurprisingly, the effects of introducing a second break on peak photon flux are relatively minor. The model fluxes lie between those of single-break power law IMF with slopes of α_m = α_u = −2.0 and −2.35, tending towards one or the other depending on the intermediate break mass, M_m selected. There is no clear sign of distinctive behaviour of double-break models at the metallicities considered, and in the absence of other probes these would likely be interpreted as single break models with a slope between the two limits if the appropriate photon flux ratios were observed.

Fig. 5. Peak ionizing fluxes, as in Fig. 3, but now for power law IMF with two breaks. We fix M0 = 0.5 M⊙, Mu=300 M⊙, αl = −1.3, αm = −2.35 and αu = −2.0, while allowing the intermediate break mass, Mm, to vary. The high metallicity fluxes in the He II and O VI ionizing photons are too low to show here, but follow the same trend.

5. Effects of total stellar mass

BPASS, in common with most other stellar population synthesis models, assumes that the star forming population is sufficiently massive that the IMF is fully sampled. An obvious problem with this assumption arises in the case of small starbursts; if the total mass of a newborn stellar population is only 1000 M_⊙, the probability of forming a single 300 M_⊙ star is very low. The fully sampled statistical IMF approach is equivalent to assigning fractions of the flux of such very massive stars to the composite spectrum, rather than considering a binary probability (1 or 0) that the star contributes. Since even a single individual massive star can dominate the hard ionizing radiation from a population, and the spectrum of that star is strongly mass dependent, this presents a challenge to modelling.

An alternate approach is stochastic sampling - drawing individual stars from an underlying probability distribution until the desired total mass is reached (Cerviño 2013). Because this is fundamentally probabilistic, many iterations must be performed for each IMF and each total mass under consideration, in order to form a complete picture. While models which do this exist and have been used to explore the properties of the ultraviolet continuum (e.g. Cerviño et al. 2003; Pflamm-Altenburg et al. 2009; Villaverde et al. 2010; Fumagalli et al. 2011; Fouesneau et al. 2012; Vidal-García et al. 2017; Paalvast & Brinchmann 2017; Ashworth et al. 2018), these have been limited to treatment of single star populations. A detailed stellar evolution-based code such as BPASS is ill-suited for such a procedure. Reading each stellar model in turn, matching it with an atmosphere model at each timestep and calculating its contribution to the composite spectrum in a given time bin is a slow process which does not lend itself to many iterations or large grids of parameters. Experimentation with an earlier BPASS version suggested that the effects of binary interactions broaden the range of evolutionary pathways at each initial mass and allow for stars to become more massive via mass transfer and mergers that mitigate some of the stochastic sampling effects (Eldridge 2012).

Here we consider an intermediate approach between simple IMF weighting and full stochastic sampling. For a selected metallicity (Z = 0.001), upper IMF slope (α_m = α_u = −2.35) and nominal mass cut-off (M_u = 300 M_⊙), we consider hybrid IMFs in which the weighting distribution for individual models is generated assuming full sampling but then modified. Given a range of assumed total stellar masses ranging from 100 to 10⁸ M_⊙, we remove individual models from the spectral synthesis if their contribution to a stellar population of that initial mass is less than 0.5 stars. We then readjust the overall scaling to account for the diminished mass of stars contributing, giving the final results as the peak photon flux per stellar mass formed in the initial starburst.

Fig. 6. Peak ionizing photon fluxes per Solar mass of stars, reached by an ageing instantaneous starburst at Z = 0.001 as a function of initial total mass assuming a total-mass-dependent cut-off in the high mass stellar population as described in Sect. 5.

In the case of a single star population, in which the probability of a model occurring only depends on the stellar initial mass, this approach is equivalent to an upper mass cut-off, M_u, that scales with the total starburst size. As a binary population synthesis code, BPASS includes a large number of models at each initial mass, with weightings which also scale with the (primary mass-dependent) period and mass ratio distributions. As a result, while the highest mass stars will always be excluded from a small population, some intermediate mass primaries will be represented by a subset of period and mass ratio models, while less frequently occurring models will be omitted.

We note that this analysis does not account for the full stochastic scatter in the population. While the probability of each high mass stellar model occurring in a small population may be low, there are likely to be a few individual stars in the high mass range which this prescription will not capture. In rare cases, a low mass population will host a large number of massive stars; in the majority of cases it will not. Similarly it is possible for a high mass starburst to lack any high mass stars when stochastic sampling is taken into acount. Thus the models calculated here are indicative of the scale of the uncertainties rather than their spread.

The results of this analysis are shown in Fig. 6. As the figure indicates, the ionizing flux from a population does not become a simple product of total stellar mass until M_total ≥ 10⁶ M_⊙, at which point the IMF is fully populated. The Q(He II) and Q(O VI) flux from BPASS models show plateaus at intermediate masses. This is likely a modelling artifact resulting from the relatively coarse grid of BPASS stellar models at high masses, and emphasizing the importance of individual very massive stars in producing these lines. The transition between the ionizing photon flux from a population hosting no very massive stars to one hosting even a single 300 M_⊙ star is abrupt and dramatic in this implementation. As a result, it seems likely that stochastic sampling of the IMF may be important for interpretation of high ionization potential spectral features, even in the galaxy-wide starbursts seen in the distant Universe. This has been seen before in work exploring spectral synthesis with single star evolution models (e.g. Vidal-García et al. 2017; Paalvast & Brinchmann 2017) but the analysis here suggests that the result holds true for binary models, despite the increased range of evolutionary pathways: the ionizing flux of the most massive star is simply too dominant in contributing the hard ionizing photons under consideration.

To test the limits of stochastic sampling, we generate three additional population synthesis models. In the first, only 300 M_⊙ primary stars and their binary companions (given by the standard BPASS v2.2 binary parameter distribution) are included in a population of total mass 10⁶ M_⊙, in the second, only 30 M_⊙ stars (and their companions) and in the third, only 3 M_⊙ stars (single and primary, with their binary companions). The ionizing fluxes from these three models as a function of metallicity are shown in Fig. 7. The short lifetimes of the 300 M_⊙ stars and their rapid evolution leads to extreme variation on their hard ionizing flux in short timescales, even while the Q(H I) hydrogen ionizing flux declines slowly and smoothly. Stars at 30 M_⊙ have hydrogen ionizing fluxes almost an order of magnitude lower than the very massive stars, but this is sustained over a much longer lifetime, as is Q(He II). The only photons capable of ionizing O VI from this population arise from the stripped binary products present at relatively late times. Finally, the 3 M_⊙ population produces very little flux shortwards of the Lyman limit, and negligible emission in Q(O VI). All stellar populations show much lower ionizing fluxes at Solar metallicity (Z = 0.020) than at lower metallicities. The lowest ratio in Q(H I)/Q(He II) attained by the 300 M_⊙ population at Z = 10⁻⁵ is 1.3 dex.

6. Discussion

The range of variation in peak fluxes, given our grid of “reasonable” IMF parameters is shown in Fig. 8 at each metallicity, where the values are shown relative to the BPASS default IMF model, which is a power law with a single break at M_l = 0.5 M_⊙, α_m = α_l = −2.35 and M_u = 300 M_⊙. As the figure makes clear, the hydrogen ionizing flux Q(H I) is relatively predictable in its response to IMF variations, as the effect of the IMF is almost independent of metallicity - in other words, if a value is calculated for one IMF and metallicity, the photon flux given different assumptions can be approximated by application of a simple scaling factor. By contrast, the hard ionizing flux is far more sensitive to variations in IMF and metallicity. In all cases, the default, Salpeter-like upper mass slope, model lies in the centre of the range of α_l variations, but our default cut-off M_u = 300 M_⊙ at the top end of the upper mass limit variations. At moderate metallicities (Z = 10⁻⁴ − 0.010), it is the upper mass cut-off that has the strongest effect on hard ionizing flux; very massive stars with M > 100 M_⊙ are critical for the production of photons sufficiently strong to power He II or O VI spectral features. The introduction of a second break to a shallower slope at high masses can also boost hard photon production by the same mechanism: increasing the number of very massive stars in the population.

Fig. 7. Ionizing photon flux time evolution for stellar populations comprised entirely of single stars or primaries at a given mass (3, 30 or 300 M⊙) and appropriate binary companions. Three metallicities are shown for each mass, and each energy threshold. For Q(O VI) the emission from low mass and high metallicity stars is negligible.

While the peak fluxes attained by a population may be critical for interpreting some of the hardest observed spectra in both the distant and local Universe, for the majority of star forming regions, this is attained only for a very brief interval in the starbursts' evolution. A more generally useful quantity may be the time-averaged photon flux over the first 20 Myr of stellar evolution (i.e. at log(age/years) < 7.3). As Fig. 2 demonstrated, this interval is sufficient to capture the lifetimes of the most massive stars, and hence the interval in which hard fluxes might be expected. In Fig. 9 we show the equivalent variation in this quantity. As might be expected, it shows very similar trends to those seen in the peak fluxes, but with the variation mitigated by the time average, yielding somewhat narrower ranges.

The origin of the short-lived peaks in hard photon production is worth consideration. These are traditionally associated with the evolution of Wolf-Rayet stars. The stripped envelopes and evolved cores of these stars allow them to expand and become luminous while remaining hot and blue. As Fig. 10 (left panel) demonstrates, the timescales for formation and evolution of these stars are in good agreement with the Q(O VI) flux in our models at low metallicity, suggesting that these might be a dominant source of hard ultraviolet emission. It is notable that the dominant contribution at early times is, in fact, from the WNH population of partially-stripped stars which still show measurable surface hydrogen fractions, rather than the most stripped (WN or WC) Wolf-Rayet stars.

Fig. 8. Range of variation in peak ionizing photon fluxes across the grid of models considered here. Range is shown separately for each metallicity, first for single-break IMFs with Mu = 300 M⊙ and varying αm = αl, then for models with fixed αm = αl = −2.35 and varying Mu, and finally for double break models with αm = −2.35, αl = −2.00 and Mu = 300 M⊙ but varying Mm. All ranges are shown relative to the BPASS default IMF model, which is a single break power law with αm = αl = −2.35 and Mu = 300 M⊙.

Fig. 9. As in Fig. 8, but now considering the mean photon flux averaged over the first 20 Myr of the starbursts’ evolution.

However, an examination of the right-hand panel of Fig. 10 suggests that care should be taken in such interpretation. A high metallicity starburst will generate as many, or even more, Wolf-Rayet stars as one that is near-primordial, while showing very little hard ionizing flux. In fact, at Z = 0.020 (Solar metallicity), the entirety of the high energy flux arises from the more numerous but far less luminous helium dwarf population generated by binary interactions of intermediate mass stars.

Fig. 10. Population of stars with a surface hydrogen mass fraction of Xsurf < 0.4, i.e. helium stars, as a function of time assuming an initial population of 106 M⊙, compared to the photon flux Q(O VI). Those with very low hydrogen (Xsurf < 0.001, classic Wolf-Rayet stars) are subdivided by the presence (WC) or absence (WN) of carbon. The population is further subdivided into high luminosity (solid lines) and low luminosity (dot-dash) populations. We show results at the highest and lowest metallicities considered here.

Fig. 11. Herzsprung–Russell diagram for stellar populations at a log(age/years) = 6.4, shown at high and low metallicity. Red regions indicate stars with a surface hydrogen mass fraction of Xsurf > 0.4, blue those with Xsurf < 0.001, and green helium stars with an intermediate surface hydrogen fraction. The same density scaling has been applied to both metallicities. Although the number of hot stars is comparable at Z = 0.020 to Z = 10−5, each star is substantially less luminous than its low metallicity counterpart.

However, while Solar-metallicity Wolf-Rayet stars are more numerous, they are also substantially less luminous than their counterparts at low metallicities. The weak stellar winds associated with low photospheric metal fractions lead to low metallicity stars retaining much of their mass through their evolution which leads to higher core temperatures to achieve hydrostatic equilibrium and thus higher stellar luminosities. For stars that accrete material through binary mass transfer, the weaker wind mass loss also allows these stars to retain more angular momentum and so they experience strong rotational mixing that causes chemically homogenous evolution to occur. These stars have high surface temperatures for their entire hydrogen-burning lifetimes and contribute significantly to the ionizing flux. The factors all result is much higher luminosities.

Figure 11 illustrates this with a traditional Hertzsprung–Russell diagram, indicating the luminosity and temperature distribution of stars at the hard emission peak age of log(age/years) = 6.4, as a function of their surface hydrogen mass fraction. The fraction of partially-stripped stars is higher, and the individual stars typically ∼ 0.5 dex more luminous, at Z = 10⁻⁵ than at Z = 0.020. At the same time, the metal abundance in the photospheres is lower. As a result, the spectrum of such starburst populations may be expected to produce substantially weaker traditional “Wolf-Rayet” spectral features (dominated by wind-broadened helium and carbon lines), while still being significant sources of hard ionizing photons. It is unclear whether the Wolf-Rayet atmosphere models currently in use for spectral synthesis fully capture this behaviour at low metallicities.

As we have demonstrated, IMF variations can substantially affect the ionizing flux. It is interesting to note, however, that the variations in the ionizing spectrum are somewhat mitigated by the fact that the entire spectrum is responding in a similar way to IMF. As Fig. 12 illustrates, even the most extreme ionizing populations in the grid we have calculated reach peak Q(He II)/Q(H I) ratios (i.e. spectral hardness measures) which only vary by a few tenths of a dex from those of the default BPASS IMF at the same metallicity.

Fig. 12. Q(H I)/Q(He II) ratios attained by our most extreme models (typically those with αm = αu = −1.8 and Mu = 300 M⊙) compared to those of the BPASS default model ( αm = αu = −2.35 and Mu = 300 M⊙) and observational constraints from the literature. Left-hand panel: time evolution of the photon flux ratio at two metallicities, illustrating the brief lifetime of the hardest ratios. Right hand panel: both the hardest ratios attained at a given metallicity and the mean ratio over the first 20 Myr of stellar evolution are shown. In both panels, we compare to observational data as indicated in the labels and in Sect. 6. The filled region in the left hand panel (and error bars on the relevant point in the right hand panel) indicates the range of values measured for different regions in the resolved galaxy SBS 0335-052E (D = 54 Mpc, Kehrig et al. 2018).

This presents a challenge. The increasing sample of ionizing galaxies in the literature showing evidence for a hard ionizing photon spectrum motivated this study and constraints on the spectral properties of these sources are also shown in Fig. 12. We have restricted our analysis to a representatitive sample of extreme sources with measurements of both the H_β and He II 1640 Å or 4686 Å recombination lines. While other lines, including [C IV] 1540, C III] 1909 and O III] 1665 are also indicative of hard ionizing fields, they are also highly sensitive to nebular gas conditions complicating interpretation of line strength. By contrast, the recombination lines of hydrogen and helium are only weakly sensitive to electron density and temperature, such that,

$$\frac{Q(\mathrm{HI})}{Q(\mathrm{HeII})}=\frac{F_{\mathrm{H}}{\lambda }_{\mathrm{H}}}{F_{\mathrm{He}}{\lambda }_{\mathrm{He}}}\frac{{\alpha }_{\mathrm{B}}^{\mathrm{H}}}{{\alpha }_{\mathrm{eff}}^{\mathrm{H}}}\frac{{\alpha }_{\mathrm{eff}}^{\mathrm{He}}}{{\alpha }_{\mathrm{B}}^{\mathrm{He}}},$$

where F_X and λ_X are the (extinction and lensing corrected) flux and wavelength of the lines used, ${\alpha }_{\mathrm{B}}^{\mathrm{X}}$ are the case B recombination coefficients and ${\alpha }_{\mathrm{eff}}^{\mathrm{X}}$ are the effective recombination coefficients relevant to each specific line. For the observational data in Fig. 12 we calculate ${\alpha }_{\mathrm{B}}^{\mathrm{X}}$ and ${\alpha }_{\mathrm{eff}}^{\mathrm{X}}$ for N_e =10³ cm⁻³ and T_e = 15,000 K, using a simple linear interpolation between the values given by Osterbrock & Ferland (2006) for each line. The electron density and temperature of extreme ionizing sources appear to vary from source to source, from region to region and even between different tracers from species to species within a source. The electron density and temperature assumed here is motivated by observations of SL2S J021737-051329 at z = 1.8 (Berg et al. 2018) and by the mean properties of stacked z ∼2.3 Lyman break galaxies (Steidel et al. 2016). They may be slightly higher than are appropriate for Q2343-BX418 (z = 2.3, Erb et al. 2010) and SGAS J105039.6+001730 (z = 3.6, Bayliss et al. 2014), but the inferred values for Q(H I)/Q(He II) are only weakly dependent on this. For SBS 0335-052E (D = 54 Mpc, Kehrig et al. 2018) we use the values provided by the authors in their table 1 for both the range of ratios in spatially resolved regions (shown as a shaded region) and the integrated value for the whole galaxy (shown by a line). Where the observational uncertainties on line flux are small, we apply a minimum uncertainty in the photon flux ratio of 0.2 dex, accounting for the possible effects of extinction, lensing corrections and uncertainty on the recombination coefficients. None of these sources show any significant evidence for AGN activity.

As the figure makes clear, even our most extreme physically-motivated IMF stellar population synthesis models are unable to reproduce the high line ratios (and hence inferred photon flux ratios) observed in some of the most extreme known systems at either high or low redshift. Although the highest flux ratios reached by our models may be consistent with the observational data for SBS 0335-052E and the z ∼ 2.3 Lyman break galaxy stack, these will only be attained for brief periods, or at the very lowest metallicities we consider. If the mean ratio observed over a 20 Myr period is considered instead, the requirement for exceptionally low metallicities (inconsistent with the observed nebular gas) becomes still stronger. For other objects at z ∼ 2 – 4, BPASS stellar spectrum models are simply unable to reproduce the line ratios, regardless of assumed IMF. Indeed, even a population comprised entirely of 300 M_⊙ stars and their binary companions never drops below Q(H I)/Q(He II) = 1.3, given our stellar evolution and atmosphere prescriptions. Hence a stochastic IMF sampling heavy in very massive stars cannot reproduce the observational data. This strongly implies that additional sources of hard ionizing photons need to be considered when analysing highly ionized galaxies, either in the local Universe or at high redshift. It is likely that the source of these hard ionizing photons arises from accretion around compact objects. The open question is whether it arises from acretion onto neutron stars and stellar mass black holes in X-ray binaries (Mirabel et al. 2011; Justham & Schawinski 2012) or supermassive black holes at the centre of galaxies acting as AGN Netzer (2015), neither of which are currently included within the BPASS models.

Fig. 13. Uncertainty in the extreme ultraviolet spectrum of a 1 Myr starburst as predicted by a range of stellar population synthesis codes, given the same two metallicities and (αu = αm = −2.35) IMF, but different stellar evolution and atmosphere input models. High resolution models have been slightly smoothed for display purposes. Model references: BPASS (Eldridge et al. 2017), BC03 (Bruzual & Charlot 2003), Starburst99 (Leitherer et al. 1999, 2011), Yggdrasil (Zackrisson et al. 2011). The 2016 version of the BC03 models has not been fully described in the refereed literature but some details are to be found in Chevallard & Charlot (2016) and Wofford et al. (2016).

This analysis has necessarily focussed on the output of just one stellar evolution and population synthesis code, the assumptions and methodology of which have been clearly described in the literature. However given the range of such models in the literature, it is appropriate to consider whether an alternate choice of stellar population synthesis model might yield a different conclusion. All model sets face two key questions: which stellar evolution tracks to employ, and what stellar atmosphere models to match them to. The former question includes uncertainties such as the impact of binary evolution or rotation, while the latter incorporates our uncertainty regarding radiative transfer through stellar envelopes as well as the difficulty of testing or calibrating predictions for the extreme ultraviolet (which is effectively unobservable due to line-of-sight absorption by interstellar hydrogen). The degree of uncertainty in the extreme ultraviolet emission of a population even at a main sequence age of 1 Myr is demonstrated in Fig. 13. As the figure demonstrates BPASS models produce substantially more ionizing flux than the older Starburst99 (Geneva, non-rotating, Leitherer et al. 1999, 2011) and BC03 (galaxev, Bruzual & Charlot 2003) population synthesis model at Solar metallicity. The 2016 release of the BC03 models presents a mixed picture. Both when using a BaSeL stellar atmosphere grid, and when using the semi-empirical MILES stellar spectra grid, the BC03-2016 model spectra extend to shorter wavelengths than BPASS at Solar metallicity, but the flux normalisation in the ionizing regime appears to be strongly dependent on the stellar atmosphere library in use, with the MILES library producing a calibration similar to BPASS just below the Lyman limit but failing to match the flux predicted by both BPASS and Starburst99 between 912 and 1216 Å, and exceeding both shortwards of 500 Å. The origin of the very blue spectrum at Solar metallicity in the BC03-2016 models is unclear (the details of implementation in this code have not appeared in the literature) but is likely associated with a boost in the number or luminosity of Wolf-Rayet spectra incorporated by the population synthesis. The same trends are seen at lower metallicities. At Z = 10⁻⁴ the extreme ultraviolet spectrum is unconstrained by observations, and the implementation favoured by BC03-2016 is significantly bluer than that of BPASS for the same input IMF. The only other SPS model which is available at this very low metallicity is the Yggdrasil model of (Zackrisson et al. 2011) which has Z = 4 × 10⁻⁴. This is very similar to the Z = 10⁻⁴ model of BC03-2016 longwards of the Lyman limit, but more similar to BPASS below it and deficient in photons with λ < 500 Å.

Fig. 14. Effect of binary stellar pathways on the evolution of the Q(H I)/Q(He II) ratios attained by stellar population synthesis models at low metallicities. The tracks show the ionizing photon flux ratio with time for the BPASS and BC03-2016 models at Z = 1 × 10−4 and the Yggdrasil model at Z = 4 × 10−4. The region spanned by the observational constraints in Fig. 13 is shaded. While the single-star evolution models begin life bluer at these metallicities, binary pathways generate hard ionizing flux at later ages.

However, as Fig. 14 demonstrates, the behaviour of the SPS models at effectively zero age may not be the key diagnostic in interpreting the spectra of high redshift stellar populations. Both the low metallicity BC03-2016 models and the Yggdrasil models evolve rapidly, with strong Q(He II) emission fading away within 2–3 Myr after the stars are formed. By contrast, binary evolution pathways generate more hot stars at later ages, and as a result the BPASS models maintain strong Q(He II) emission relative to Q(H I) over timescales of 30 Myr or longer. While all the models fail to match the values inferred from the observational constraints, timescale arguments suggest that we are unlikely to be catching so many systems at ages <3 Myr, somewhat favouring the binary evolution models.

7. Conclusions

Our main conclusions can be summarised as follows:

1. Young starbursts at low metallicity can generate the high energy photons (shortwards of 228 Å and 89.8 Å) required to power emission or absorption features in the He II and O VI lines.

2. The photon flux at these energies is 2.5 orders of magnitude lower than that shortwards of 912 Å for Q(He II), and five orders of magnitude lower for Q(O VI), even when at its peak and at very low metallicities.

3. Hard ionizing radiation, particularly Q(O VI), is associated with the most massive stars in the population (M > 100 M_⊙), and so is short-lived and boosted by shallower power law IMF or higher mass IMF cut-offs.

4. The probability of these very massive stars occurring in a starburst with a low total mass is ≪ 1, meaning the fluxes in Q(He II) and Q(O VI) are highly susceptible to stochastic variation, even in galaxy-wide starbursts.

5. Even the most extreme IMF in a grid motivated by observations of stellar masses in the local Universe are unable to reproduce the low Q(H I)/Q(He II) ratios observed in extreme ionizing sources in either the local or distant Universe.

The analysis presented here demonstrates the range of hard photon emission rates that result from plausible assumptions regarding the stellar IMF. We caution that detailed interpretation of features in individual galaxy spectra is inevitably going to be subject to uncertainties in the IMF of its contributing starbursts. We remind the community that the IMF is fundamentally a statistical construct, and that stellar population synthesis models are most effective when considering entire galaxy populations rather than individual objects. Nonetheless, it seems likely that for some of the hardest ionization potentials inferred for galaxies in the distant Universe, IMF variations are inadequate to explain the observed emission, and additional sources of short-wavelength photons beyond those from stellar photospheres are required.

Acknowledgments

BPASS would not be possible without the computational resources of the University of Auckland’s NeSI Pan Cluster and the University of Warwick’s Scientific Computing Research Technology Platform (SCRTP). New Zealand’s national facilities are provided by the NZ eScience Infrastructure and funded jointly by NeSI’s collaborator institutions and through the Ministry of Business, Innovation & Employment’s Research Infrastructure programme. We thank both the University of Warwick and the University of Auckland for travel support. ERS thanks John Chisholm and Steve Wilkins for useful discussions. We also thank Polychronis Papaderos and the other organisers of the “Escape of Lyman radiation from Galactic Labyrinths” conference (OAC, Crete, 2018) for the stimulus that motivated this work.

References

All Tables

All Figures

Fig. 1. Effect of IMF slope on binary star fraction in a population of total mass 106 M⊙, given the stellar mass-dependent binary fractions implemented in BPASS v2.2. The lower region between the x-axis and the first pair of lines indicates the fraction of the total stellar mass for a given IMF slope αu which is assigned to single star models, the central region the mass fraction in primary star models and the upper region the mass fraction assigned to their less massive secondary companions. Regions between each pair of lines (with hatching in two styles) switch between stellar model types depending on the IMF. Models assume a single power law break (αm = αl), The upper cut-off mass, Mu, is shown by line style (100 M⊙ - solid; 300 - dot-dash). The curves for Mu =150 & 200 M⊙ lie between these curves.

In the text

Fig. 2. Time evolution of photon production rates in H I (black), He II (blue) and O VI (yellow-green) at five metallicities, from an instantaneous burst of total mass 106 M⊙ occurring at an age of zero. At each metallicity three single-break power law (αu = αm = −1.8, −2.35, −2.70) IMF models are shown at two upper mass limits (Mu = 100, 300 M⊙), with line style indicating different slopes and line colour different upper mass cut-offs as indicated in the key. The rest of the models lie within the bounds of those shown at fixed age.

In the text

	Fig. 3. Peak ionizing photon fluxes reached by an ageing instantaneous starburst of total mass 106 M⊙ as a function of IMF slope (assuming αm = αu). In each panel, metallicity is indicated by colour while upper mass cut-off Mu is indicated by line style (100 M⊙ = solid, 300 M⊙ = dot-dash).
In the text

	Fig. 4. Peak ionizing fluxes, colour coded by metallicity, as in Fig. 3, but now for single-break power law initial mass functions with a varying upper mass limit, Mu. We fix αm = αl = −2.35. The high metallicity fluxes in the He II and O VI ionizing photons are too low to show here, but follow a similar trend.
In the text

	Fig. 5. Peak ionizing fluxes, as in Fig. 3, but now for power law IMF with two breaks. We fix M0 = 0.5 M⊙, Mu=300 M⊙, αl = −1.3, αm = −2.35 and αu = −2.0, while allowing the intermediate break mass, Mm, to vary. The high metallicity fluxes in the He II and O VI ionizing photons are too low to show here, but follow the same trend.
In the text

	Fig. 6. Peak ionizing photon fluxes per Solar mass of stars, reached by an ageing instantaneous starburst at Z = 0.001 as a function of initial total mass assuming a total-mass-dependent cut-off in the high mass stellar population as described in Sect. 5.
In the text

	Fig. 7. Ionizing photon flux time evolution for stellar populations comprised entirely of single stars or primaries at a given mass (3, 30 or 300 M⊙) and appropriate binary companions. Three metallicities are shown for each mass, and each energy threshold. For Q(O VI) the emission from low mass and high metallicity stars is negligible.
In the text

Fig. 8. Range of variation in peak ionizing photon fluxes across the grid of models considered here. Range is shown separately for each metallicity, first for single-break IMFs with Mu = 300 M⊙ and varying αm = αl, then for models with fixed αm = αl = −2.35 and varying Mu, and finally for double break models with αm = −2.35, αl = −2.00 and Mu = 300 M⊙ but varying Mm. All ranges are shown relative to the BPASS default IMF model, which is a single break power law with αm = αl = −2.35 and Mu = 300 M⊙.

In the text

	Fig. 9. As in Fig. 8, but now considering the mean photon flux averaged over the first 20 Myr of the starbursts’ evolution.
In the text

Fig. 10. Population of stars with a surface hydrogen mass fraction of Xsurf < 0.4, i.e. helium stars, as a function of time assuming an initial population of 106 M⊙, compared to the photon flux Q(O VI). Those with very low hydrogen (Xsurf < 0.001, classic Wolf-Rayet stars) are subdivided by the presence (WC) or absence (WN) of carbon. The population is further subdivided into high luminosity (solid lines) and low luminosity (dot-dash) populations. We show results at the highest and lowest metallicities considered here.

In the text

Fig. 11. Herzsprung–Russell diagram for stellar populations at a log(age/years) = 6.4, shown at high and low metallicity. Red regions indicate stars with a surface hydrogen mass fraction of Xsurf > 0.4, blue those with Xsurf < 0.001, and green helium stars with an intermediate surface hydrogen fraction. The same density scaling has been applied to both metallicities. Although the number of hot stars is comparable at Z = 0.020 to Z = 10−5, each star is substantially less luminous than its low metallicity counterpart.

In the text

Fig. 12. Q(H I)/Q(He II) ratios attained by our most extreme models (typically those with αm = αu = −1.8 and Mu = 300 M⊙) compared to those of the BPASS default model ( αm = αu = −2.35 and Mu = 300 M⊙) and observational constraints from the literature. Left-hand panel: time evolution of the photon flux ratio at two metallicities, illustrating the brief lifetime of the hardest ratios. Right hand panel: both the hardest ratios attained at a given metallicity and the mean ratio over the first 20 Myr of stellar evolution are shown. In both panels, we compare to observational data as indicated in the labels and in Sect. 6. The filled region in the left hand panel (and error bars on the relevant point in the right hand panel) indicates the range of values measured for different regions in the resolved galaxy SBS 0335-052E (D = 54 Mpc, Kehrig et al. 2018).

In the text

Fig. 13. Uncertainty in the extreme ultraviolet spectrum of a 1 Myr starburst as predicted by a range of stellar population synthesis codes, given the same two metallicities and (αu = αm = −2.35) IMF, but different stellar evolution and atmosphere input models. High resolution models have been slightly smoothed for display purposes. Model references: BPASS (Eldridge et al. 2017), BC03 (Bruzual & Charlot 2003), Starburst99 (Leitherer et al. 1999, 2011), Yggdrasil (Zackrisson et al. 2011). The 2016 version of the BC03 models has not been fully described in the refereed literature but some details are to be found in Chevallard & Charlot (2016) and Wofford et al. (2016).

In the text

Fig. 14. Effect of binary stellar pathways on the evolution of the Q(H I)/Q(He II) ratios attained by stellar population synthesis models at low metallicities. The tracks show the ionizing photon flux ratio with time for the BPASS and BC03-2016 models at Z = 1 × 10−4 and the Yggdrasil model at Z = 4 × 10−4. The region spanned by the observational constraints in Fig. 13 is shaded. While the single-star evolution models begin life bluer at these metallicities, binary pathways generate hard ionizing flux at later ages.

In the text