© Y. Sibony et al. 2022

1. Introduction

The stars we know from our local Universe formed from giant molecular clouds, which have experienced chemical enrichment (or pollution, depending on one’s point of view) from previous stellar generations. As a result all the stars we have ever observed contain at least traces of metals. Since all the known mechanisms that produce elements heavier than beryllium are stellar (or of stellar origin, e.g. supernovae or neutron star mergers), there must have been a series of generations of stars, each one forming from the ejecta of the previous one, being more and more metal-rich as time went on. Logically then, stars without any metals – the so-called Population III (Pop III) stars – while not having been observed yet, must have existed as the first generation of stars, forming from primordial Big Bang nucleosynthesis material around a redshift of z ∼ 30 − 10 (e.g. Barkana & Loeb 2001).

Because the Universe was very different then than it is now (chemically, but also in terms of average density and temperature, as well as the intensity of magnetic fields), the processes of Pop III star formation (e.g. Abel et al. 2000; Bromm 2013) and stellar evolution (see e.g. Heger & Woosley 2010; Yoon et al. 2012; Murphy et al. 2021) were also very different from what we can observe in the local Universe (for a review article concerning Pop III stars, see Bromm 2017). A number of studies on the theoretical properties of these stars have been conducted in recent decades, concluding–among many other things–that Pop III stars should be more massive and their initial mass function (IMF) more top-heavy than current stars (Greif et al. 2011; Stacy & Bromm 2013; Susa et al. 2014; Wollenberg et al. 2020). Furthermore, for the same mass, a Pop III star should exhibit a larger effective temperature and luminosity than its metal-enriched counterpart (Ekström et al. 2008; Yoon et al. 2012; Murphy et al. 2021). Also, Pop III stars may be fast rotators because they rapidly accrete gas with high angular momentum, which is typical in hot, thick Pop III star-forming disks with inefficient cooling (Stacy & Bromm 2013). Rotational mixing renders stars more luminous and hotter during their main sequence (MS) phase, thus enhancing their capacity to produce ionising photons. This implies that their production rate of ionising photons could be very large; indeed Pop III stars are expected to be one of the causes (along with active galactic nuclei and Population II stars around z ∼ 20 − 15, although the latter is expected to exhibit much softer spectra than their metal-free counterparts, see Schaerer 2003) of the early reionisation period (z ∼ 30 − 15) which followed the dark ages of the Universe.

Very fast rotation may induce such strong mixing inside stars so that they follow a homogeneous chemical evolution, representing an extreme limit of what single star evolution can achieve in terms of ionising power. In order to place upper limits on the ionising flux emitted by the first stars, one can consider the case where all Pop III stars would follow a chemically homogeneous evolution (CHE) path. Indeed, stars evolving under this constraint of being perfectly mixed at all times would be even more compact and would thus have an even hotter surface. Furthermore, their bluewards (instead of the classical redwards) evolution during the MS, coupled with their larger fuel supply compared to standard models, would make chemically homogeneous stars the largest emitters of integrated ionising radiation over the population’s lifetime. We propose an analytical model including chemically homogeneous stellar evolution and ionising radiation in order to quantify the aforementioned upper limit. We then derive upper limits on the reionisation of the intergalactic medium (IGM) by Pop III stars. We compare our results to those obtained with classical non-rotating Pop III models computed with the Geneva stellar evolution code (GENEC).

This article is organised as follows. In Sect. 2 we detail the analytical model for homogeneous evolution. In Sect. 3 we present the ionisation model, based on the self-similar champagne flow solution (Shu et al. 2002). Section 4 presents our results and compares them to those obtained with the standard Pop III models. Section 5 discusses some limitations of the present approach. Finally, we summarise our main findings in Sect. 6.

2. Chemically homogeneous stellar models

When the matter processed by nucleosynthesis in the stellar core is mixed with the whole stellar mass on timescales much shorter than the nuclear timescale, stars have a chemically homogeneous evolution. This type of evolution has been advocated for explaining young blue stragglers in young open star clusters (Maeder 1987), the evolution of long soft Gamma Ray Bursts (Yoon & Langer 2005), and the evolution leading to the formation of close binary black holes (Mandel & de Mink 2016; de Mink & Mandel 2016; Marchant et al. 2017; du Buisson et al. 2020). At the moment, the observational evidence for stars following such an evolution remains scarce. However, as discussed in Martins et al. (2009), Marchant et al. (2017), some early-type massive stars appear to be very good candidates. The physical mechanism responsible for the strong mixing is probably linked to the hydro- or magnetohydrodynamical instabilities triggered by a fast rotation. For now however, the physics of rotation still remain a debated question and different hypotheses lead to different results (see e.g. the discussion in Song et al. 2016). In the present work, we do not focus on the physical cause of homogeneous evolution but rather on the consequences such an evolution would have on cosmic reionisation. This is why we use here a very simple approach with an analytical model for such stars.

The quantities we aim to derive from our analytical model are the star’s luminosity and effective temperature. Such models have been made in the literature with polytropes (Beech & Mitalas 1989). While this could also be a physically relevant approach for advanced phases of evolution, we focus here on the MS phase. We do not consider mass loss as the stars do not reach the Eddington luminosity, and at primordial chemical composition radiatively line-driven winds should be negligible.

A polytrope model is characterised by a global relation between pressure and density:

$$\begin{array}{c}\hfill P(r)=K\rho {(r)}^{\frac{n+1}{n}}.\end{array}$$(1)

Here, we consider the special case n = 3 (the so-called Eddington Standard model). In this case, the Eddington quartic relation provides a relation between the mass m_⋆ of the star, its mean molecular weight μ (constant throughout the star by virtue of chemical homogeneity), and its ratio of gas to total pressure, β:

$$\begin{array}{c}\hfill \frac{1-\beta }{{\beta }^4}=A{\mathrm{\mu }}^4{\left(\frac{m_{\star }}{M_{\odot }}\right)}^2,\end{array}$$(2)

where A = 0.002986 is a constant. Equation (2) has four solutions for β, two of which are complex, one is a negative real number, and one is a positive real number. Hence for a given stellar mass (we do not consider the case of mass loss), at any stage of its MS phase (since μ is directly linked to the hydrogen mass fraction X by μ = (0.75 + 1.25X)⁻¹ for a Pop III star), we can compute the value of β.

Given the three quantities m_⋆, X, and β, we use Eq. (5.47) of Prialnik (2009) to compute the luminosity:

$$\begin{array}{c}\hfill \frac{L}{L_{\odot }}=0.003\frac{4\pi cG{M}_{\odot }}{{\kappa }_s{L}_{\odot }}{\mathrm{\mu }}^4{\beta }^4{\left(\frac{m_{\star }}{M_{\odot }}\right)}^3,\end{array}$$(3)

where κ_s is the surface electron scattering opacity assuming that the medium is completely ionised:

$$\begin{array}{c}\hfill {\kappa }_s=0.2(1+X){\mathrm{cm}}^2/\mathrm{g}.\end{array}$$(4)

We note that the mass-luminosity relation simplifies to L ∝ m_* for very massive stars (m_* ≳ 100 M_⊙; Bromm et al. 2001). Our next quantity of interest is the star’s effective temperature. Here we use the following relation between the central and effective temperature of a homogeneous n = 3 polytrope from Hoyle & Fowler (1963):

$$\begin{array}{c}\hfill \left(\frac{T_{\mathrm{eff}}}{{10}^4\mathrm{K}}\right)=1.010{(1-\beta )}^{1/8}{(1+X)}^{-1/4}\left(\frac{T_c}{{10}^6\mathrm{K}}\right),\end{array}$$(5)

where T_c is the star’s central temperature.

The system of Eqs. (3) and (5) is not fully constrained by the sole choice of a stellar mass, which determines only the luminosity on the zero-age main sequence (ZAMS). We also need to fix an initial value for the radius, or equivalently, effective temperature, central temperature, central density, or mean density. To that end, we extract the central temperature on the ZAMS from the Pop III models of Murphy et al. (2021). This initial condition constrains the initial effective temperature and thus the star’s initial radius. We then evolve the central temperature using the following relation (which is just a consequence of the perfect gas equation of state and of the hydrostatic equilibrium equation):

$$\begin{array}{c}\hfill T_c\propto \frac{\beta \mathrm{\mu }{m}_{\star }}{R_{\star }},\end{array}$$(6)

though it turns out that the change in T_c over the MS for the CHE models ends up being negligible.

The method goes as follows: for the masses at which Murphy et al. (2021) computed Pop III models (complemented by models at 180, 250, and 300 M_⊙ computed with the same method), we vary X from its initial primordial value of X = 0.75 down to X = 0 (as hydrogen is burned by the star), and compute the evolution of β, μ, and κ_s. This allows us to compute the evolution of the luminosity with Eq. (3) and of the central temperature with Eq. (6). Then, using Eq. (5), we can get to the effective temperature. As a result, the stars’ evolution in the HR diagram can be plotted, and their radius computed from L and T_eff.

In our model, the mass fraction of hydrogen X serves as a proxy for time. However, in the halo ionisation model that we shall introduce further on, we need the actual time evolution of the ionising radiation. In order to go from X to t, we use the equivalence of mass and energy and the fact that hydrogen burning generates 0.007 grams of energy for each gram of hydrogen consumed. Thus in the time dt, the mass m_⋆ dX of hydrogen is consumed, generating a luminosity

$$\begin{array}{c}\hfill L(t)=0.007\frac{\mathrm{d}X}{\mathrm{d}t}{m}_{\star }{c}^2.\end{array}$$(7)

Inverting the above relation provides us with the duration of each timestep, given the variation of hydrogen and the luminosity at that timestep. This allows us to express the surface quantities of stars as functions of time rather than hydrogen fraction. We also obtain the lifetime of each star by equating it to the time at which the star reaches X = 0.

Finally, we consider the stars to be perfect black bodies of temperature T_eff. The rate of photon production as a function of time and stellar mass is given by the integral over the relevant range of wavelengths of the black-body spectral energy distribution, multiplied by the surface of the star:

$$\begin{array}{c}\hfill {\dot{N}}_{\gamma }(\mathrm{\Delta }\lambda ,t,m_{\star })={\displaystyle {\int }_{\mathrm{\Delta }\lambda }}\mathrm{d}\lambda F_{\lambda }(T_{\mathrm{eff}}(m_{\star },t))\frac{\lambda }{\mathit{hc}}4\pi R_{\star }^2(m_{\star },t),\end{array}$$(8)

where F_λ(T_eff(m_⋆, t)) = πB_λ(T_eff(m_⋆, t)) is the outgoing spectral flux density for the black body, Δλ is the wavelength range of the photons we consider, and R_⋆(m_⋆, t) is the stellar radius for the star of mass m_⋆ at time t.

3. Ionisation model

Once Ṅ_γ(Δλ,t,m_⋆) is known, one can calculate the reionisation history (at z ≳ 15, dominated by Pop III stars) in two steps. The first step is to calculate the (average) escaping ionising photon number Q_esc(Δλ, M_h, z) per Pop III star-forming halo for different halo masses M_h and redshifts z, in the wavelength ranges of interest. The second step is to calculate the evolution of the IGM ionisation fraction Q_ion(z) by combining the contribution of all haloes hosting Pop III stars.

3.1. Ionising photons

We consider three types of ionising photons: those capable of ionising H I (E >  13.6 eV, thereafter called ‘H II photons’), those capable of ionising He I (E >  24.6 eV), thereafter ‘He II photons’, and those capable of ionising He II (E >  54.4 eV, thereafter ‘He III photons’). Since computing the interactions of the produced photons with the surrounding medium properly with radiative transfer is beyond the scope of this analytical work, we instead consider an idealised scenario where photons ionise the species with the highest ionisation potential possible:

$$\begin{array}{c}\hfill \{\begin{array}{c}13.6\mathrm{eV}<E_{\gamma }<24.6\mathrm{eV}:\mathrm{H}\mathrm{I}\to \mathrm{H}\mathrm{II},\hfill \\ 24.6\mathrm{eV}<E_{\gamma }<54.4\mathrm{eV}:\mathrm{He}\mathrm{I}\to \mathrm{He}\mathrm{II},\hfill \\ E_{\gamma }>54.4\mathrm{eV}:\mathrm{He}\mathrm{II}\to \mathrm{He}\mathrm{III}.\hfill \end{array}\end{array}$$

We also compute the number of escaping Lyman-Werner (LW, 11.2 eV <  E_γ <  13.6 eV) photons. These photons can dissociate H₂ molecules and thus are important to produce atomic-cooling haloes in which one expects the possible formation of very and even supermassive stars (see e.g. the review by Haemmerlé et al. 2020, and references therein). We consider that those photons only escape if hydrogen is fully ionised (see e.g. Schauer et al. 2015, 2017). Since they are not ionising photons, they escape the clusters at the same rate at which they are produced (in essence we assume that the H II photons take care of the recombined atoms and the LW photons escape freely).

3.2. Escape of ionising photons from Pop III hosts

Assuming ionisation equilibrium in the ionised medium, the number of escape ionising photons equals the number of ionising photons produced by stars minus the number of recombinations, as recombinations have to be undone by absorbing incoming photons to keep the medium ionised.

3.2.1. Production: Pop III cluster model

Considering the small number of Pop III stars formed per halo (especially for low-mass haloes), instead of computing the IMF-averaged escaping photon number for a given halo mass M_h and redshift z, we do the computation for N_cl = 100 stochastic Pop III clusters and then take the average escape rate at each timestep. In each cluster, we sample stars from a power-law initial mass function (IMF) $\mathrm{d}N/\mathrm{d}{m}_{\star }\sim m_{\star }^{-\alpha }$ between m_min = 9 M_⊙ and m_max = 300 M_⊙, for a total stellar mass M_⋆ = ϵ_⋆(Ω_b/Ω_m)M_h¹. The only free parameters here are the IMF power-law index α and star formation efficiency (SFE) ϵ_⋆. In total, we explore four scenarios: the ‘Fiducial’ has α = −1, ϵ_⋆ = 0.001; the ‘HighSFE’ has α = −1, ϵ_⋆ = 0.002; the ‘Salpeter’ has α = −2.35, ϵ_⋆ = 0.001; and the ‘Salpeter & HighSFE’ has α = −2.35, ϵ_⋆ = 0.002. We also generate clusters where the maximum stellar mass is m_max = 120 M_⊙, but the effect, especially on chemically homogeneous models, is very minor. This is because the number of emitted photons per baryon over the stellar lifetime is mostly flat for m_⋆ ≳ 100 M_⊙. For this reason we do not show the results for clusters with this smaller m_max.

The sampling method is the following: as long as the total stellar mass M_tot <  M_⋆, we keep sampling stars from the IMF, bounded by m_min and m_max. The sampling is finished when we obtain M_⋆ ≤ M_tot <  1.05M_⋆. If we overshoot (M_tot >  1.05M_⋆), we remove all stars and start again. Then, for each cluster i ∈ [1, N_cl] we can compute the ionising photon production rate in the wavelength range Δλ at time t:

$$\begin{array}{c}\hfill {\dot{Q}}_{\star ,i}={\dot{Q}}_{\star ,i}(\mathrm{\Delta }\lambda ,t)={\displaystyle \sum _{j=1}^{N_{\star }(t)}}{\dot{N}}_{\gamma }(\mathrm{\Delta }\lambda ,t,m_{\star ,j}),\end{array}$$(9)

where N_⋆(t) is the total number of stars still alive at time t and Ṅγ(Δλ,t,m_⋆,j) the production rate of photons in the wavelength range Δλ at time t for star j, given by Eq. (8).

3.2.2. Absorption: The champagne flow solution

After being emitted by Pop III stars, ionising photons will be absorbed by the interstellar medium (ISM) of the host halo and only a fraction of them can escape into the IGM. To model this process, we use the idealised semi-analytical model for Pop III ionisation feedback in Alvarez et al. (2006) under spherical symmetry, which is meant to capture the late-stage overall behaviour of the ISM. Below we summarise the basic ideas of this model and the readers are referred to Alvarez et al. (2006) for details. We note that the propagation of ionising photons in Pop III star-forming clouds and haloes has been studied with more complex semi-analytical models and simulations (see e.g. Fernandez & Shull 2011; Paardekooper et al. 2015; Tanaka et al. 2018; Sugimura et al. 2020; Jaura et al. 2022). We adopt the simple model in Alvarez et al. (2006) to efficiently investigate the parameter space of Pop III star formation and demonstrate the effects of CHE. In the next section, we compare our results with previous studies to justify our approach.

For a Pop III cluster i embedded in a halo of mass M_h at redshift z, the ionisation feedback is modelled as a champagne flow bounded by a (D-type)² shock running into the primordial star-forming cloud whose initial density profile follows a singular isothermal sphere (SIS)³ with a temperature T₁ (see Eq. (1) in Alvarez et al. 2006). The shock heats the gas up to a temperature T₂ ∼ 2 × 10⁴ K such that the high pressure homogenises the downstream medium (i.e. the medium between the star and the inoisation front) and drives an outflow. The initial (upstream) temperature T₁ is related to the (normalisation of) gas density profile ρ(r) assuming that in the central region dominated by gas (i.e. ignoring dark matter), gravity is balanced by the thermal pressure of gas⁴. Considering the total mass of gas in the halo $4\pi {\int }_0^{R_{\mathrm{v}\mathrm{i}\mathrm{r}}}\rho (r|T_1)r^2\mathrm{d}r=({\mathrm{\Omega }}_{\mathrm{b}}/{\mathrm{\Omega }}_{\mathrm{m}})M_h$, we can express T₁ with halo properties:

$$\begin{array}{c}\hfill \left(\frac{T_1}{300\mathrm{K}}\right)=({\mathrm{\Omega }}_{\mathrm{b}}/{\mathrm{\Omega }}_{\mathrm{m}})\left(\frac{M_h}{950M_{\odot }}\right){\left(\frac{R_{\mathrm{vir}}}{1\mathrm{pc}}\right)}^{-1},\end{array}$$(10)

where R_vir ≡ R_vir(M_h, z) is the virial radius.

Under this setup, the (number) density profile of the gas cloud n(r, t)≡n(r, t|M_h, z) is described by the champagne flow solution⁵ in Shu et al. (2002), governed by the parameter

$$\begin{array}{c}\hfill ϵ={\left(\frac{c_{s,1}}{c_{s,2}}\right)}^2=\frac{T_1}{T_2}\frac{{\mathrm{\mu }}_2}{{\mathrm{\mu }}_1}.\end{array}$$(11)

Here we have μ₂ ≃ 0.63 for the shocked ionised downstream region, while for the unshocked region, we have μ₁ ≃ 1.22 (for neutral gas). The quantities c_s, 1 and c_s, 2 are respectively the sound speed in the neutral, ionised region. The densities of the neutral and ionised medium are obtained using Table 1 of Shu et al. (2002) linking values of ϵ to these densities.

In the early stage, the downstream medium is not fully ionised (due to the high density) and no ionising photons can escape the halo. Later on, as density decreases by the outflow, the entire downstream is ionised. A part of the ionising photons emitted by the star is used to keep the ionised downstream region ionised, while the rest can ionise the outer region beyond the shock front and even escape the halo. The rate (luminosity) of escaping ionising photons in the wavelength range Δλ can be written as:

$$\begin{array}{cc}\hfill {\dot{Q}}_{\star ,\mathrm{esc},i}& \equiv {\dot{Q}}_{\star ,\mathrm{esc},i}(\mathrm{\Delta }\lambda ,t,M_h,z)\hfill \\ \hfill & =max\{{\dot{Q}}_{\star ,i}-4\pi {\alpha }_{\mathrm{k}}^{\mathrm{B}}{\displaystyle {\int }_0^{10R_{\mathrm{vir}}}}{r}^2{n}_{\mathrm{k}}^2(r,t)\mathrm{d}r,0\}.\hfill \end{array}$$(12)

Here the first term on the right hand side is the production rate of ionising photons by all the stars in the cluster, defined in the previous section. The second term is the recombination rate within 10 times the halo virial radius R_vir, which is dominated by the downstream medium within the shock front. Here we chose 10R_vir as the boundary between the halo and the IGM. n_k(r, t)≡x_kn(r, t) is the density profile and ${\alpha }_{\mathrm{k}}^{\mathrm{B}}$ is the case-B recombination coefficient for species k that corresponds to Δλ. For our three species (H II, He II, He III), the abundances are x_H = 0.93, x_He = 0.07, and we use the formulas of ${\alpha }_{\mathrm{k}}^{\mathrm{B}}$ in Appendix E of Rosdahl et al. (2013), all in units of [cm³ s⁻¹]:

$$\begin{array}{cc}\hfill {\alpha }_{\mathrm{H}\mathrm{II}}^{\mathrm{B}}& =2.753\times {10}^{-14}\frac{{\lambda }_{\mathrm{H}\mathrm{I}}^{1.5}}{{[1+{({\lambda }_{\mathrm{H}\mathrm{I}}/2.74)}^{0.407}]}^{2.242}},\hfill \\ \hfill {\alpha }_{\mathrm{He}\mathrm{II}}^{\mathrm{B}}& =1.26\times {10}^{-14}{\lambda }_{\mathrm{He}\mathrm{I}}^{0.75},\hfill \\ \hfill {\alpha }_{\mathrm{He}\mathrm{III}}^{\mathrm{B}}& =5.506\times {10}^{-14}\frac{{\lambda }_{\mathrm{He}\mathrm{II}}^{1.5}}{{[1+{({\lambda }_{\mathrm{He}\mathrm{II}}/2.74)}^{0.407}]}^{2.242}},\hfill \end{array}$$

where ${\lambda }_{\mathrm{H}\mathrm{I}}=\frac{315614\mathrm{K}}{T[\mathrm{K}]}$, ${\lambda }_{\mathrm{He}\mathrm{I}}=\frac{570670\mathrm{K}}{T[\mathrm{K}]}$, and ${\lambda }_{\mathrm{He}\mathrm{I}\mathrm{I}}=\frac{1263030\mathrm{K}}{T[\mathrm{K}]}$, and T = T₂ = 2 × 10⁴ K in our case.

The whole procedure described above is applied to N_cl clusters to obtain the average (over all clusters) luminosity of escaping photons ${\dot{Q}}_{\star ,\mathrm{esc}}(\mathrm{\Delta }\lambda ,t,M_h,z)={\displaystyle \sum _{i=1}^{N_{\mathrm{cl}}}}{\dot{Q}}_{\star ,\mathrm{esc},i}/N_{\mathrm{cl}}$. At last, the average escaping ionising photon number per halo can be derived as $Q_{\star ,\mathrm{esc}}(\mathrm{\Delta }\lambda ,M_h,z)={\displaystyle \int {\dot{Q}}_{\star ,\mathrm{esc}}}(\mathrm{\Delta }\lambda ,t,M_h,z)\mathrm{d}t$, where we integrate over the cluster lifetime (that is, the lifetime of the longest-lived star in the cluster). The average escape fraction of ionising photons as a function of time is $f_{\mathrm{esc}}(\mathrm{\Delta }\lambda ,t,M_h,z)\equiv {\dot{Q}}_{\star ,\mathrm{esc}}(\mathrm{\Delta }\lambda ,t,M_h,z)/{\dot{Q}}_{\star }(\mathrm{\Delta }\lambda ,t,M_h,z)$, where ${\dot{Q}}_{\star }(\mathrm{\Delta }\lambda ,t,M_h,z)$ is the average production rate of ionising photons, and the overall escape fraction is ${\widehat{f}}_{\mathrm{esc}}(\mathrm{\Delta }\lambda ,M_h,z)\equiv Q_{\star ,\mathrm{esc}}(\mathrm{\Delta }\lambda ,M_h,z)/Q_{\star }(\mathrm{\Delta }\lambda ,M_h,z)$ given the total number of photons produced $Q_{\star }(\mathrm{\Delta }\lambda ,M_h,z)={\displaystyle \int {\dot{Q}}_{\star }}(\mathrm{\Delta }\lambda ,t,M_h,z)\mathrm{d}t$. In other words:

$$\begin{array}{c}\hfill {\widehat{f}}_{\mathrm{esc}}(\mathrm{\Delta }\lambda ,M_h,z)=\frac{\int {\dot{Q}}_{\star ,\mathrm{esc}}(\mathrm{\Delta }\lambda ,t,M_h,z)\mathrm{d}t}{\int {\dot{Q}}_{\star }(\mathrm{\Delta }\lambda ,t,M_h,z)\mathrm{d}t}·\end{array}$$(13)

We can then obtain the cosmic escape fraction as a function of redshift:

$$\begin{array}{c}\hfill \overline{f_{\mathrm{esc}}}(\mathrm{\Delta }\lambda ,z)=\frac{{\int }_{M_1}^{M_2}\frac{d{n}_h(M_h,z)}{dln{M}_h}{f}_{\mathrm{esc}}(\mathrm{\Delta }\lambda ,M_h,z)\mathrm{d}ln{M}_h}{{\int }_{M_1}^{M_2}\frac{d{n}_h(M_h,z)}{dln{M}_h}\mathrm{d}ln{M}_h},\end{array}$$(14)

where n_h(M_h, z) is the halo mass function (Tinker et al. 2008), for which we consider haloes in the mass range defined by M₁ = 10⁶ M_⊙ and M₂ = 10⁸ M_⊙ as the typical hosts of Pop III stars⁶.

3.3. IGM ionisation fraction

The second step is to calculate the evolution of the IGM ionisation fraction Q_ion by combining the contribution of all haloes hosting Pop III stars. We define Q_ion in the following way (where we decouple singly and doubly ionised helium for simplicity):

$$\begin{array}{c}\hfill \{\begin{array}{c}{Q}_{\mathrm{ion},\mathrm{H}\mathrm{II}}=\frac{n_{\mathrm{H}\mathrm{II}}}{n_{\mathrm{H}}},\hfill \\ Q_{\mathrm{ion},\mathrm{He}\mathrm{II}}=\frac{n_{\mathrm{He}\mathrm{III}}}{n_{\mathrm{He}}},\hfill \\ Q_{\mathrm{ion},\mathrm{He}\mathrm{III}}=\frac{n_{\mathrm{He}\mathrm{III}}}{n_{\mathrm{He}}}.\hfill \end{array}\end{array}$$

For species k, the evolution of Q_ion, k follows (see e.g. Sect. 3.1 of Hartwig et al. 2015):

$$\begin{array}{c}\hfill \frac{d{Q}_{\mathrm{ion},\mathrm{k}}(z)}{\mathit{dz}}={(x_{\mathrm{k}}{n}_{\mathrm{b}})}^{-1}\frac{d{n}_{\mathrm{ion},\mathrm{k}}(z)}{\mathit{dz}}-\left[\frac{Q_{\mathrm{ion},\mathrm{k}}(z)}{t_{\mathrm{rec},\mathrm{k}}(z)}\right]\left|\frac{\mathrm{d}t}{\mathrm{d}z}\right|,\end{array}$$(15)

given the recombination timescale

$$\begin{array}{c}\hfill t_{\mathrm{rec},\mathrm{k}}(z)=\frac{1}{C_{\mathrm{IGM}}(z){\alpha }_{\mathrm{k}}^{\mathrm{B}}{x}_{\mathrm{k}}{n}_{\mathrm{b}}{(1+z)}^3},\end{array}$$(16)

with the case-B recombination coefficient ${\alpha }_{\mathrm{k}}^{\mathrm{B}}$ evaluated at a temperature of T_IGM = 10⁴ K for ionised IGM, and

$$\begin{array}{c}\hfill \frac{d{n}_{\mathrm{ion},\mathrm{k}}(z)}{\mathit{dz}}={\displaystyle {\int }_{M_1}^{M_2}}\frac{\mathrm{d}{n}_h(M_h,z)}{\mathrm{d}z}{Q}_{\star ,\mathrm{esc}}(\mathrm{\Delta }\lambda ,M_h,z)\mathrm{d}{M}_h.\end{array}$$(17)

Here n_b ≃ 2 × 10⁻⁷ cm⁻³ is the cosmic comoving number density of baryons, n_h(M_h, z) is the halo mass function, and C_IGM(z) is the redshift-dependent clumping factor, from Haardt & Madau (2012):

$$\begin{array}{c}\hfill C_{\mathrm{IGM}}(z)=1+43z^{-1.71}.\end{array}$$(18)

We solve the differential Eq. (14) from z = 30 to z = 15, taking as initial values Q_ion,HI(z = 30) = 10⁻⁴ and Q_ion,HeII(z =3 0) = Q_ion,HeIII(z = 30)=0 (see Fig. 1 in Galli & Palla 2013).

4. Results

4.1. Stellar models

We compare the evolution in the HR diagram, lifetimes of the MS, as well as photon production over time, for stars computed with our analytical Pop III CHE model and standard Pop III stars computed with GENEC by Murphy et al. (2021). Figure 1 shows the evolution in the HR diagram of the CHE (solid) and standard (dashed) Pop III models during their MS. Right at the onset of hydrogen burning the standard models have a chemically homogeneous composition, so we would expect the starting position of CHE and classical models of the same mass to be at the same position in the HR diagram. However, there are small differences between the CHE and classical models, which come from the fact that polytropes are only providing an approximation of the structure of a star. Let us remind that using a polytropic model implicitely assumes that the ratio of the gas to the total pressure β is uniform in the star, which is not strictly the case, although β does not vary much over most of the mass of the star.

Fig. 1. HR diagram showing the MS evolution of the Pop III chemically homogeneous models (solid lines) and that of the Pop III standard models computed with GENEC by Murphy et al. (2021, dashed lines). The numbers indicate the initial mass in M⊙ of each model.

The most striking feature of the CHE models is their blue and upwards evolution in the HRD. Indeed, the very strong mixing characteristic of CHE models removes the core-envelope differentiation. This means that the traditional hydrogen envelope expansion that occurs during the MS does not happen for those models. The chemically homogeneous models do expand slightly during this phase; both their luminosity and effective temperatures increase, but T_eff does not grow fast enough to keep the radius constant. That expansion however is limited to a factor of 1.7 between the initial and final radii, whereas for the standard models it can reach a factor of 100. Interestingly, this factor decreases with mass for the CHE models (the final to initial radius ratio being 1.7 for the 9 M_⊙ star and reaching 1.07 for the 300 M_⊙ star) whereas it increases with mass for the classical models (the final to initial radius ratio also being 1.7 for the 9 M_⊙ star, and reaching 100 for the 300 M_⊙ star).

Apart for the 9 M_⊙ and 12 M_⊙ stars where the CHE models have much longer MS lifetimes, both sets of models exhibit very similar MS durations. This can be explained by the competition of two effects. On the one hand, the CHE models have a much larger fuel supply since the entire initial mass of hydrogen will burn into helium: this should increase the lifetime of the MS. On the other hand, the CHE models are hotter and more luminous than the classical models, which should decrease the lifetime of the MS. The result, as stated above, is that masses above 12 M_⊙ do not exhibit significant differences in their MS duration. The difference in lifetimes between the lower mass stellar models can be seen in the late evolution of the clusters, where no more photons are emitted after ∼18 Myr in the clusters using standard stars whereas the clusters with chemically homogeneous stars last until ∼27 Myr. All these reults, as well as the numbers of photons emitted by the individual stellar models over their MS lifetimes in our different wavelength ranges of interest, are summarised in Table 1.

Figure 2 shows the number of photons emitted per baryon integrated over the MS lifetime as a function of stellar mass, for the four wavelength ranges considered. This number incorporates both the photon emission rates and the duration of the MS. Comparing photon production, the CHE models emit systematically more photons than their classical counterparts, except for the LW band. This makes sense as the CHE stars are more luminous (more photons emitted in total), and have higher effective temperatures (more energetic photons emitted). For the LW photons, their wavelengths correspond to black-body temperatures of 26 − 32 × 10⁴ K (or log₁₀(T_eff[K])∼4.4 − 4.5), which are closer to the effective temperatures of the standard models than to those of the CHE models. As such, the effect of chemically homogeneous evolution on the LW feedback should be minor. Still, the production of ionising photons (whether H ionising or He singly or doubly ionising) is higher for the CHE than the standard models by a factor of ∼1.5 − 8 (decreasing with stellar mass).

Fig. 2. Number of photons emitted per baryon over the MS lifetime, for each stellar mass and for the four wavelength ranges considered (H II photons: 13.6eV <  E <  24.6eV, He II photons: 24.6eV <  E <  54.4eV, He III photons: E >  54.4eV, LW photons: 11.2eV <  E <  13.6eV). Instead of plotting individual data points corresponding to individual masses we join them for better clarity. The solid lines represent the CHE models and the dashed ones, the classical models computed with GENEC.

For the CHE models, it is notable that the rate of emitted H II and He II photons per baryon does not depend on stellar mass for m_⋆ ≳ 100 M_⊙. The overall values are consistent in an order of magnitude with those in Fig. 1 of Heger & Woosley (2010).

4.2. Escaping photons

In the above subsection, we looked at the photons emitted by stars, without accounting for the surrounding environment. Here, we focus on the effect of the absorption due to recombinations in the halo, corresponding to the second term on the right-hand side of Eq. (12), which integrates the density profile of the champagne flow solution at each time. For illustration purposes, we show in Fig. 3 such density profiles at t = 0.1, 0.3, 1, and 3 Myr, for a typical Pop III halo with M_h = 10⁶ M_⊙ at z = 20.

Fig. 3. Density profiles obtained from the champagne flow solution (Shu et al. 2002) at t = 0.1 (dashed), 0.3 (dash-dotted), 1 (dotted), and 3 Myr (long-dashed), on top of the initial density profile (solid) for a typical primordial star-forming cloud in a DM minihalo with Mh = 106 M⊙ at z = 20.

We can now take into account the recombinations happening both downstream and upstream of the shock front, removing stars that have reached the end of their MS lifetime, at every timestep for the CHE and the classical Pop III models.

In Fig. 4, we show the number of escaping photons per second in the hydrogen ionisation (left column) and helium first ionisation (right column) wavelength ranges, for a dark matter minihalo of mass M_h = 10⁶ M_⊙ (upper panel of each subplot; CHE models in red, standard models in blue) and an atomic-cooling halo of mass M_h = 10⁸ M_⊙ (lower panel of each subplot, same colours), at redshifts z = 30, 25, 20, 15. The four cluster scenarios’ parameters are given in Table 2. We recall that the free parameters of the clusters are the IMF slope and star formation efficiency. The panels’ vertical axes for haloes of the same mass in the same scenario use the same scale, for easier comparison between H II and He II photon escape.

Fig. 4. Number of escaping H II and He II photons per second (${\dot{Q}}_{\star ,\mathrm{esc}}(\mathrm{\Delta }\lambda ,t,M_h,z)$) as a function of time, for clusters in a minihalo (Mh = 106 M⊙) and clusters in an atomic-cooling halo (Mh = 108 M⊙). At each integer redshift value between 30 and 15 we generate 100 stochastic clusters and show the average number of escaping photons from these clusters at each timestep. Each of the four rows represents a different set of model parameters (see Table 2).

The difference in the escape rate of H II photons between the CHE and standard models can be interpreted with the higher emission rates of ionising photons and shorter lifetimes of Pop III stars under CHE compared with those under the standard stellar evolution. The former effect is more significant for lower-mass stars while the latter is only important for massive stars (m_⋆ >  120 M_⊙). In low-mass haloes (M_h = 10⁶ M_⊙), where escape of ionising photons is relatively easy, the first effect (higher emission rates) dominates so that the total number of escaping H II photons is larger for the CHE models. However, for massive haloes (M_h = 10⁸ M_⊙), particularly under the top-heavy IMF (Fiducial and HighSFE), the second effect (shorter lifetimes) becomes dominant, such that the peak of escape rate is reached a few Myr earlier for the CHE models than for the standard ones. Since at later stages, the gas density in the halo (and consequently the rate of recombinations) has decreased, more photons escape in the standard case than with CHE.

The effect of the IMF slope is very striking as the bottom-heavy scenarios yield rates of escaping photons that are much flatter over time than the top-heavy IMFs. Indeed, high mass stars produce more ionising photons over shorter lifetimes than low-mass ones, so the scenarios which create them in larger quantities lead to a more pronounced peak of ${\dot{Q}}_{\star ,\mathrm{esc},\mathrm{H}\mathrm{I}\mathrm{I}}$ around 3 Myr (the MS lifetime of those high mass stars). Interestingly, CHE clusters with the bottom-heavy IMFs have more escaping photons when integrating over their whole lifetime than the top-heavy ones. This is because lower-mass stars have a longer lifetime and keep emitting photons at late cluster times when there are very few recombinations.

The star formation efficiency ϵ_⋆ also has a strong influence on the number of escaping HII photons, which is understandable since the HighSFE scenarios double the total stellar mass in the clusters but keep the number of recombinations (linked only to the halo mass and redshift) constant.

The effect of halo mass, all other parameters being equal, should be to both increase the number of recombinations for higher halo masses, while also increasing the total production of photons (because those haloes produce more stars). It turns out that the effects do not compensate each other: there are about 100 times as many ionising photons produced by the haloes with M_h = 10⁸ M_⊙ than by those with M_h = 10⁶ M_⊙. However, the escape fractions are much lower for the high mass haloes than for the low-mass ones (see Fig. 5).

Fig. 5. Escape fractions for H II photons averaged over the lifetimes of the clusters and over the Ncl stochastic clusters, ${\widehat{f}}_{\mathrm{esc}}(\mathrm{H}\mathrm{I}\mathrm{I})$, from fiducial haloes of different masses at four different redshifts. The solid lines show the escape fractions from haloes with chemically homogeneous stars, while the dotted lines show the escape fractions from haloes with classical stars.

Finally, for all sets of parameters, we find that there are more escaping HHII photons at the lower redshifts than at higher ones. This is because redshift in our model does not influence how many photons are produced, but only affects the number of recombinations in the halo. Indeed, for the same halo mass, its virial radius R_vir increases with time (or with decreasing redshift), thus the initial upstream temperature T₁ decreases with time (or with decreasing redshift) and so is the Shu et al. (2002) parameter ϵ. As ϵ decreases towards low redshifts, so do the gas density and the number of recombinations, increasing the number of escaping photons.

The results above mostly hold for the He II wavelength range, with the relative difference between CHE and classical models being larger for the escaping HeHII photons than for the HHII photons. This is logical considering the black-body temperatures that correspond to HeHII photons overlap even more with the effective temperatures of the CHE models in comparison with the standard models, than those corresponding to HHII photons. We note also that halo mass and redshift are less impactful on the escape rate of HeHII photons relative to HHII photons. That is because the number of recombining He nuclei is very small in all cases (leading to escape fractions close to 1 for all four scenarios and both types of stellar models).

We do not show the equivalent to Fig. 4 for LW and He III photons as they are relatively less important, but we present the relevant results here. For the Fiducial and HighSFE scenarios the LW emissions are similar between classical and CHE models, whereas the CHE models dominate in the bottom-heavy Salpeter and Salpeter & HighSFE scenarios. This is linked to what is shown in Fig. 2 (purple curve). Indeed, clusters containing more lower mass stars (the Salpeter and Salpeter & HighSFE scenarios) emit more LW photons when they have CHE stars than standard ones. Conversely, the standard clusters with a top-heavy IMF actually exhibit a higher maximum number of escaping LW photons.

In the He III photons wavelength range, the classical Pop III models produce very few photons (not even enough to counterbalance the recombinations in some cases, leading to an escape rate of 0). The effects of redshift, halo mass, and cluster scenario are qualitatively the same as for other wavelength ranges discussed previously. The main features of this wavelength range are that the relative difference in escape rates between the CHE and standard clusters is the largest among all wavelength ranges, and that those escape rates are in any case very low (≲10⁴⁹ s⁻¹ for 10⁶ M_⊙ haloes, ≲100 × 10⁴⁹ s⁻¹ for 10⁸ M_⊙ haloes).

We perform a check of our ionising photons escape model by using version C17.02 of the photoionisation code Cloudy (Ferland et al. 2017) to compute the output spectra of a few synthetic clusters embedded in haloes of different masses at different redshifts. The sizes and densities of the haloes are given by the Shu et al. (2002) model. For the Cloudy modelling we assume spherical geometry, constant densities, and primordial chemical abundances, then run the code to convergence. The quantity we compare is the transmitted spectrum, and we do not take into account the diffuse spectrum. We expect significant differences in escaping photon fluxes as our ionisation model is highly idealised. To be specific, the analytical estimations were calculated in strict energy bins (see Sect. 3.1), and these bins do not take into account the fact that at energies above 24.6 eV, the ionisation cross-sections of hydrogen per hydrogen nucleus are comparable (or sometimes higher) than the helium ones (see Fig. 2 of Glatzle et al. 2019), leading to an underestimation of the amount of hydrogen ionising photons when compared to the Cloudy models. Furthermore, the analytical model does consider the case in which multiple photons with individually too low energies could together ionise atoms (e.g. two LW photons of the right energies could excite then ionise a hydrogen atom).

However, the differences turn out to be modest, thus we conclude that our simple model is appropriate. In particular, for clusters with high photon production rates (typically the HighSFE scenario, or Fiducial clusters that drew very massive stars), the escape fractions computed by Cloudy are very similar (≲2% relative difference). For clusters producing relatively few photons, we underestimate by ≲10% the H II photons’ escape fraction, and overestimate by ≲20% the He II photons’ escape fraction. This leads to similar relative differences in the IGM ionisation rates if we consider a worst case scenario where we rescale all the photon escape rates from all clusters by these amounts (even those clusters where the difference between our simple model and the Cloudy output is very small). In any case, the variation in escape rates between our model and the Cloudy outputs would apply similarly to clusters with CHE stars and to those with classical stars, so the difference in ionising power would still hold.

Figure 5 shows the variation of the escape fraction of H II photons as a function of halo mass, at redshifts z=30, 25, 20, and 15, for the Fiducial scenario. The CHE clusters are represented by solid lines while dotted lines are for the clusters containing classical Pop III stars. Interestingly, while the CHE haloes see more escaping photons at most halo masses and redshifts, their escape fractions for the large halo masses are smaller than those of the standard clusters. This is caused by a rather fine effect: the standard stars reach their maximum ionising power slightly later than the CHE models, so the gas density from Shu et al. (2002) and thus the number of recombinations are smaller for the classical haloes. This effect is particularly noticeable in haloes with large numbers of high mass stars (high mass haloes) and recombinations (high mass haloes at high redshift). Overall, we find the expected result that escape fractions decrease with increasing halo mass (as also found by e.g. Wise et al. 2014), and with increasing redshift. Indeed, halo mass scales positively with the number of recombinations, so a larger part of the emitted photons in high mass haloes must be used to reionise the recombined nuclei. On the other hand, recombination rate scales positively with matter density, so that there are fewer recombinations at a later time (lower redshift), leading to increased escape fractions.

We draw attention to the fact that the escape fractions we show in Fig. 5 are averaged over the lifetimes of the clusters, as well as over the 100 stochastic clusters. Therefore we do not show the time or the stochastic variability of the escape fractions. The time variability is quite extreme, as in a given cluster the escape fraction can take values from 0 to ∼0.9. On the other hand, the stochastic variability in the values of ${\widehat{f}}_{\mathrm{esc}}$ (see Eq. (13)) is milder, with relative differences in lifetime-integrated escape fractions of ∼10% between the maximally and minimally escaping clusters.

4.3. IGM ionisation

We can now convolve the results from stochastic clusters embedded in haloes of specific masses at specific redshifts, with the distribution of halo masses at each redshift (the halo mass function). Figure 6 shows the cosmic escape fraction $\overline{f_{\mathrm{esc}}}$ (defined in Eq. (14)) as a function of redshift, for the wavelengths corresponding to H II, He II, and He III photons (LW photons are not shown as their escape fraction is always 1 in our model). The vertical axis for He II photons goes from 0.9 to 1 so as to better see the differences between scenarios. For all three wavelength ranges we find that, as expected, the CHE models yield larger cosmic escape fractions than the standard models. Furthermore, almost all He II photons produced escape their host haloes (at least 95% of them in the worst case), regardless of the scenarios and stellar models.

Fig. 6. Cosmic escape fraction $\overline{f_{\mathrm{esc}}}(\mathrm{\Delta }\lambda ,z)$ (given by Eq. (14)). Top panel: H II photons (13.6eV <  E <  24.6eV), middle panel: He II photons (24.6eV <  E <  54.4eV), bottom panel: He III photons (E >  54.4eV). The middle panel (He II photons) is zoomed-in between 0.9 and 1 to better distinguish the different curves.

A result that at first glance may seem counter-intuitive is that the H II escape fraction is higher for the clusters computed with a Salpeter IMF than for those with a fiducial IMF with α = −1 (regardless of the SFE). Indeed, from Fig. 4 we may get the impression that the Fiducial scenario leads to more escaping H II photons than the Salpeter one. But when we integrate over the lifetime of the clusters, there are more H II photons escaping from Salpeter clusters. To understand this result, one has to combine two facts here: the Salpeter IMF yields more lower-mass and thus longer-lived stars, and the number of recombinations in a halo decreases over time as the density of the champagne flow solution decreases. As a result the late phases (t ≳ 10 Myr) of the Salpeter clusters contribute a non-negligible fraction of the total number of escaping H II photons, while in the fiducial clusters most stars will have died by then. This result holds for the chemically homogeneous clusters, but not those comprised of classical Pop III stars from the GENEC models.

Figure 7 shows the ionisation fractions of H and He (singly and doubly ionised), Q_ion,HII, Q_ion,HeII, and Q_ion,HeIII. We can see that the standard Pop III scenarios (using GENEC models with fiducial star formation efficiency and either IMF) lead to Q_ion,HII ∼ 0.1−0.2 at z ∼ 15. Keeping the same cluster parameters but switching to CHE stellar models increase Q_ion,HII to ∼0.2 − 0.3. Doubling the SFE more than doubles the hydrogen ionisation fraction. This holds true for helium ionisation. As a result, this star formation efficiency parameter ϵ_⋆ is a very important one to estimate for Pop III stars since it strongly influences the impact of these stars on the reionisation. We find slightly higher ionisation fractions for He II than for H II, but a very small doubly ionised helium fraction.

Fig. 7. Fraction of ionised species as a function of redshift Qion, k(z). The bottom panel (He III ionisation fraction) is zoomed-in between 0 and 0.1 to better distinguish the different curves.

We then compare our results with previous theoretical predictions and observational constraints on the ionisation history (of hydrogen). We only focus on Pop III stars at z ∼ 15 − 30 and do not model the transition between Pop III and metal-enriched Pop II star formation, which is expected to happen in the same redshift range (z ∼ 15 − 30, see e.g. Maio et al. 2010; Johnson et al. 2013; Xu et al. 2016; Mebane et al. 2018; Jaacks et al. 2018; Sarmento et al. 2018; Liu & Bromm 2020). So we extrapolate the ionised fraction Q_ion, HII from our model to z = 0 with a conservative estimate of the ionisation history in the Pop II-dominated era (z ≲ 15). We assume that the ionised fraction is 0.1 at z = 9, appropriate for late reionisation scenarios (e.g. Greif & Bromm 2006), and increases with decreasing z following a power-law to reach unity at z = 5.5. Between z = 9 and z = 15, we connect the high-z (predicted) and low-z (assumed) pieces with a power-law⁷. The results are shown in Fig. 8, compared with the predictions from the semi-analytical models ASLOTH (Hartwig et al. 2022), 21cmFAST (Muñoz et al. 2022) and Astraeus (Hutter et al. 2021), as well as cosmological simulations ASTRID (Bird et al. 2022), THESAN (Kannan et al. 2022) and SPHINX (Rosdahl et al. 2018). We find that with the fiducial SFE and Salpeter IMF, the GENEC model is generally consistent with the semi-analytical models that cover z ≳ 15 (e.g. Hartwig et al. 2022; Muñoz et al. 2022), while Q_ion, HII is a factor of ∼2 − 4 higher with CHE. With a higher SFE and/or a more top heavy IMF, Q_ion, HII is significantly higher for both CHE and GENEC by up to a factor of ∼10.

Fig. 8. Ionised fraction of hydrogen. Extrapolation of the predictions (at z ∼ 15 − 30) from our CHE and standard models are shown with thick solid and dashed curves. For comparison, we plot the results from the semi-analytical models ASLOTH (Hartwig et al. 2022), 21cmFAST (Muñoz et al. 2022) and Astraeus (Hutter et al. 2021), as well as cosmological simulations ASTRID (Bird et al. 2022), THESAN (Kannan et al. 2022) and SPHINX (Rosdahl et al. 2018) with thin solid and dashed curves, respectively. The vertical lines show the redshift range (z ∼ 5.5 − 15) in which our predictions at high-z are extrapolated to the post-reionisation epoch.

We also calculate the optical depth to Thomson scattering, τ, from the ionisation history, as shown in Fig. 9 on top of the value measured by Planckτ = 0.0544 ± 0.0073 (Planck Collaboration I 2020). The standard GENEC model with the fiducial SFE and Salpeter IMF well reproduces the observed τ. Under the same conditions, the CHE model slightly overpredicts τ by ∼1σ = 0.0073. τ is significantly higher than the observed value with a higher SFE and/or a more top-heavy IMF, by ∼1 − 1.5σ for the Fiducial GENEC and CHE models, and more than 3σ in the HighSFE case of CHE with either IMF and GENEC with the fiducial IMF. While in the intermediate cases with fiducial SFE for CHE and enhanced SFE for GENEC, under the Salpeter IMF, τ is overpredicted by ∼1.5 − 2σ. Our results indicate that the ionisation history and contribution to τ at z ≳ 15 (when Pop III stars can play an important role) is highly sensitive to the SFE and IMF. This also implies that the maximum SFE allowed so that the Universe is not ionised too fast is lower if most Pop III stars experience CHE; this effect is more significant for a less top-heavy IMF. Finally, we would like to point out that the above analysis is rather qualitative as τ is a measurement of the entire ionisation history and our model only covers the high-z regime and relies on extrapolation at low z. Beyond τ, the radiation fields from Pop III stars (e.g. ionising, Lyman-α, and LW photons) can be better constrained with 21-cm cosmology (see e.g. Fialkov et al. 2014; Barkana 2016; Cohen et al. 2017; Tanaka et al. 2018; Mirocha et al. 2018; Schauer et al. 2019a; Mebane et al. 2020). We plan to investigate the impact of CHE on the 21-cm signal in future work.

Fig. 9. Optical depth to Thomson scattering, based on the extrapolation of the predicted ionisation histories from our CHE and GENEC models (see Fig. 8). The recent measurement by Planck (Planck Collaboration I 2020) is shown with the horizontal line, whose 1, 2, and 3σ scatters are denoted by the shaded regions. Again, the vertical lines show the redshift range (z ∼ 5.5 − 15) in which our predictions at high-z are extrapolated to the post-reionisation epoch.

5. Discussion

The idea of this work is to remain simple in providing an upper limit to the contribution of single Pop III stars to cosmic reionisation. Still, we are aware of potential refinements that would improve the precision of our results. These improvements would also increase the complexity of this paper, and we defer them to future work on the subject. For completeness, we nevertheless discuss them here.

First, we derive analytical expressions for the evolution of the chemically homogeneous models, based on the n = 3 polytrope solution. While this seems to be a reasonable approximation for massive CHE stars (Beech & Mitalas 1989), a more accurate way to proceed would be to use a numerical code to compute the stellar evolution. However, at least in this exploratory work, the polytropic approach is very well justified since it grasps the essential physics and is well adapted to stars having homogeneous chemical composition.

In the present work, we only consider the MS phase of Pop III stellar evolution, thus omitting advanced phases. For the standard models, the contributions of the post core H-burning phase to the ionising photons budget are rather modest and actually do not significantly change the results presented here. A larger effect is expected for the CHE models, however. Indeed, cores made up of a mixture of helium, carbon and oxygen are more luminous and hotter at a given mass, compared to their counterparts made of hydrogen and helium. Such stars may therefore contribute as sources of highly energetic photons. On the other hand, those stars may have luminosity near the Eddington limit (Maeder et al. 2012), and thus lose mass. This in turn would decrease their luminosity. Also, the duration of these advanced phases corresponds to approximately 10% of the time spent on the MS, thus mitigating their importance. It would nevertheless be interesting to address this question quantitatively with CHE models in a future study.

Another potential weakness of the present models is that mass loss is not included, as Pop III stars are expected to have very weak – if any – line-driven winds, due to the lack of heavy elements at their surface. This would remain true also for homogeneous evolution during the MS phase, since during this stage no heavy elements are synthesised in the core in significant quantities. Although there are other mass loss mechanisms, some of them mentioned in Liu et al. (2021), most would occur after the MS phase and thus would not impact the results presented here.

Also on the stellar models, we do not produce detailed spectra for our stars but only consider them to emit black-body radiation. While the approximation is less of a problem for stars without metals than for those at higher metallicities, the computation of synthetic spectra would certainly increase the accuracy of our emission model. For instance, Schaerer (2002) predicted the ionising power of Pop III stars and populations based on non-LTE atmosphere models; Topping & Shull (2015) estimated the efficiency of stellar reionisation at higher metallicities using realistic spectra of GENEC models, computed with the WM-basic code (Pauldrach et al. 2001; Puls et al. 2005).

Using polytropic models does not allow to link the existence of stars following a homogeneous evolution to some special physical initial conditions. For instance, if the strong internal mixing needed for such an evolution is triggered by rotation, there are certainly lower limits on the initial mass and the initial rotation (likely depending on the initial mass in turn) above which such an evolution can indeed be realised. These limits are however very sensitive to which effects of rotation are accounted for, only hydrodynamical instabilities induced by rotation like in the theory by Zahn (1992), or magneto-hydrodynamical instabilities as proposed by Spruit (1999, 2002). In the present approach, we have explored two cases that can be seen as two extremes, either no Pop III stars follow a homogeneous evolution or all follow such an evolution. Actually, we expect that a fraction of Pop III stars would undergo CHE. This fraction can be computed once we have a better view of the physics leading to such an evolution.

Next, our ionisation model of the interactions between emitted photons and the surrounding gas is very idealised. There are two orthogonal ways to proceed should we like to increase its complexity. We could keep the model analytical, but make it more ‘physical’. An example of such improvements would be to include detailed cross-sections, thus allowing photons to ionise species over an extended wavelength range, instead of strictly segregating the interactions by a simple binning procedure. Conversely, we could use a numerical code to solve the radiative transfer and compute escape fractions from each cluster at each time much more accurately. Such computations however would necessarily reduce the number of clusters and time coverage that can be considered. One can think of a third way towards improving the ionisation model: compute escape fractions numerically (e.g. with Cloudy) for a few representative clusters at a few different times in their evolution, and interpolate the results as a function of the features of the stellar masses in the cluster (such as the average or median mass, the number of stars, as well as the minimum and maximum stellar masses in the cluster).

Modelling the complex radiation-hydrodynamical (RHD) evolution of the expanding I-front around Pop III stars with the self-similar champagne flow is clearly an idealisation. However, comparisons with sophisticated RHD simulations show that the resulting density structure is surprisingly well-matched, in particular when considering more advanced evolutionary stages (see the discussion and references in Wang et al. 2012). The density structure in turn is key to capture the propagation and possible escape from the host halo of ionising photons.

An aspect which we do not take into account is feedback by supernovae (SNe) clearing out the gas in the haloes, thus increasing the rates of escaping photons (by reducing the number of recombinations) after the deaths of the most massive stars (around 3 Myr). For the top-heavy scenarios, at least 70–80% of all escaping photons (over the lifetimes of the clusters) escape before 3 Myr, so neglecting the effect of SNe is not so much of a problem. For the bottom-heavy scenarios however, that number drops to 30–35%. In these cases, of the photons produced after 3 Myr, the fraction that have to reionise recombined atoms is ∼65% (∼45%) for the Salpeter (Salpeter & HighSFE) scenarios. This means that a non negligible increase in escape rates would occur if we cleared the haloes of their gas after the first supernovae. We aim to implement such an improvement to the model in a further work on this subject.

Finally, we assume that the era of Pop III stars is terminated uniformly at z = 15, and do not consider the effect of Pop II stars between z = 15 and z = 6 (at which point the Universe is expected to be fully reionised). This leaves a degree of freedom when extrapolating our results to observables in the local Universe such as the optical depth to Thomson scattering. The inclusion of (potentially chemically homogeneous) Pop II stars in the model would make for a significant improvement in a follow-up work, but as mentioned above, a preliminary step is to identify the physics responsible for such an evolution.

6. Conclusions

In this work, we have studied the maximal impact that single Pop III stars may have had on early reionisation, by considering the extreme case where all of them underwent chemically homogeneous evolution and comparing it to the baseline case of standard Pop III evolution. To this end, we developed both analytical polytrope models and ionisation schemes, and a semi-analytical solution for the gas density inside star clusters (Shu et al. 2002). We used stochastic sampling from an initial mass function to compute the number of escaping ionising photons from a considerable amount of clusters embedded in dark matter haloes. These haloes populate a large grid in halo mass and redshift, such that we could determine the number of escaping photons to construct cosmological predictions for high redshifts. Finally, using conservative estimates for the later ionisation history of the Universe, we extrapolated our results to the current Universe and compared them to measured observables.

Our main results are summarised as follows:

• Chemically homogeneous Pop III stars produce 2–5 times as many H I and He I-ionising photons, and up to 100 times the number of He II-ionising photons, as standard Pop III stars of the same mass. However the difference in LW photon production is more moderate, and even inverted (where standard models produce more than CHE ones) for stellar masses above 100 M_⊙.

• The number of H I-ionising photons escaping from clusters, as well as the escape fractions, are also larger for clusters containing CHE stars, except for high mass clusters. The threshold mass above which standard clusters exhibit larger escape numbers and fractions increases with cosmological time (decreasing redshift). The escape rates for He I-ionising photons are systematically about twice larger for all halo masses and redshifts.

• The escape fractions obtained analytically are in good agreement with our Cloudy models. The small differences arise from the energy bins used in this work, which do not account for the hydrogen ionisation cross-sections being comparable (or higher) than the helium ones for a small range of energies above 24.6 eV.

• The ionised fractions of hydrogen (singly-ionised helium) in the IGM reach values of 0.22–0.75 (0.40–0.99) at z = 15 for different scenarios of CHE clusters, which correspond to about twice the values for standard clusters. For a given star formation scenario, the values we report for the chemically homogeneous clusters are a good estimate for the upper limit of the ionised fraction of the IGM.

• Extrapolating our results to the present-day Universe, and computing the optical depth to Thomson scattering, rules out the scenario with high star formation efficiency and all Pop III stars undergoing chemically homogeneous evolution. The chemically homogeneous Fiducial and Salpeter scenarios are within measurement uncertainties, although the best agreement with observations is obtained for the standard Fiducial and Salpeter scenarios.

This paper highlights the importance of stellar evolution parameters on cosmological evolution, and in particular shows that chemically homogeneous evolution greatly impacts the ionising budget in the early Universe. We stress that this kind of evolution cannot happen in all Pop III stars, implying that the physical conditions under which it can occur should not be too easily realised. It remains a very interesting prospect that the details of Pop III stellar structure and evolution have important cosmological implications.

Acknowledgments

The authors thank an anonymous referee for valuable comments. The authors would also like to thank A. Verhamme and O. Attia for thought-provoking scientific discussion. YS and GM have received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 833925, project STAREX).

References

All Tables

All Figures

	Fig. 1. HR diagram showing the MS evolution of the Pop III chemically homogeneous models (solid lines) and that of the Pop III standard models computed with GENEC by Murphy et al. (2021, dashed lines). The numbers indicate the initial mass in M⊙ of each model.
In the text

Fig. 2. Number of photons emitted per baryon over the MS lifetime, for each stellar mass and for the four wavelength ranges considered (H II photons: 13.6eV <  E <  24.6eV, He II photons: 24.6eV <  E <  54.4eV, He III photons: E >  54.4eV, LW photons: 11.2eV <  E <  13.6eV). Instead of plotting individual data points corresponding to individual masses we join them for better clarity. The solid lines represent the CHE models and the dashed ones, the classical models computed with GENEC.

In the text

	Fig. 3. Density profiles obtained from the champagne flow solution (Shu et al. 2002) at t = 0.1 (dashed), 0.3 (dash-dotted), 1 (dotted), and 3 Myr (long-dashed), on top of the initial density profile (solid) for a typical primordial star-forming cloud in a DM minihalo with Mh = 106 M⊙ at z = 20.
In the text

Fig. 4. Number of escaping H II and He II photons per second (${\dot{Q}}_{\star ,\mathrm{esc}}(\mathrm{\Delta }\lambda ,t,M_h,z)$) as a function of time, for clusters in a minihalo (Mh = 106 M⊙) and clusters in an atomic-cooling halo (Mh = 108 M⊙). At each integer redshift value between 30 and 15 we generate 100 stochastic clusters and show the average number of escaping photons from these clusters at each timestep. Each of the four rows represents a different set of model parameters (see Table 2).

In the text

Fig. 5. Escape fractions for H II photons averaged over the lifetimes of the clusters and over the Ncl stochastic clusters, ${\widehat{f}}_{\mathrm{esc}}(\mathrm{H}\mathrm{I}\mathrm{I})$, from fiducial haloes of different masses at four different redshifts. The solid lines show the escape fractions from haloes with chemically homogeneous stars, while the dotted lines show the escape fractions from haloes with classical stars.

In the text

	Fig. 6. Cosmic escape fraction $\overline{f_{\mathrm{esc}}}(\mathrm{\Delta }\lambda ,z)$ (given by Eq. (14)). Top panel: H II photons (13.6eV <  E <  24.6eV), middle panel: He II photons (24.6eV <  E <  54.4eV), bottom panel: He III photons (E >  54.4eV). The middle panel (He II photons) is zoomed-in between 0.9 and 1 to better distinguish the different curves.
In the text

	Fig. 7. Fraction of ionised species as a function of redshift Qion, k(z). The bottom panel (He III ionisation fraction) is zoomed-in between 0 and 0.1 to better distinguish the different curves.
In the text

Fig. 8. Ionised fraction of hydrogen. Extrapolation of the predictions (at z ∼ 15 − 30) from our CHE and standard models are shown with thick solid and dashed curves. For comparison, we plot the results from the semi-analytical models ASLOTH (Hartwig et al. 2022), 21cmFAST (Muñoz et al. 2022) and Astraeus (Hutter et al. 2021), as well as cosmological simulations ASTRID (Bird et al. 2022), THESAN (Kannan et al. 2022) and SPHINX (Rosdahl et al. 2018) with thin solid and dashed curves, respectively. The vertical lines show the redshift range (z ∼ 5.5 − 15) in which our predictions at high-z are extrapolated to the post-reionisation epoch.

In the text

Fig. 9. Optical depth to Thomson scattering, based on the extrapolation of the predicted ionisation histories from our CHE and GENEC models (see Fig. 8). The recent measurement by Planck (Planck Collaboration I 2020) is shown with the horizontal line, whose 1, 2, and 3σ scatters are denoted by the shaded regions. Again, the vertical lines show the redshift range (z ∼ 5.5 − 15) in which our predictions at high-z are extrapolated to the post-reionisation epoch.

In the text