Large scale simulations of H and He reionization and heating driven by stars and more energetic sources

We present simulations of cosmic reionization and reheating from $z=18$ to $z=5$, investigating the role of stars (emitting soft UV-photons), nuclear black holes (BHs, with power-law spectra), X-ray binaries (XRBs, with hard X-ray dominated spectra), and the supernova-associated thermal bremsstrahlung of the diffuse interstellar medium (ISM, with soft X-ray spectra). We post-process the hydrodynamical simulation Massive-Black II (MBII) with multifrequency ionizing radiative transfer. The source properties are directly derived from the physical environment of MBII, and our only real free parameter is the ionizing escape fraction $f_{\rm esc}$. We find that, among the models explored here, the one with an escape fraction that decreases with decreasing redshift yields results most in line with observations, such as of the neutral hydrogen fraction and the Thomson scattering optical depth. Stars are the main driver of hydrogen reionization and consequently of the thermal history of the intergalactic medium (IGM). We obtain $\langle x_{\rm HII} \rangle = 0.99998$ at $z=6$ for all source types, with volume averaged temperatures $\langle T \rangle \sim 20,000~{\rm K}$. BHs are rare and negligible to hydrogen reionization, but conversely they are the only sources which can fully ionize helium, increasing local temperatures by $\sim 10^4~{\rm K}$. The thermal and ionization state of the neutral and lowly ionized hydrogen differs significantly with different source combinations, with ISM and (to a lesser extent) XRBs, playing a significant role and, as a consequence, determining the transition from absorption to emission of the 21 cm signal from neutral hydrogen.


INTRODUCTION
The existence and characteristics of the last phase change of our Universe-the Epoch of Reionization (EoR, see e.g. Zaroubi 2013 for a review) has been a puzzle ever since the cosmological models of last century envisioned some sort of matter, likely hydrogen, to exist in between galaxies (an intergalactic medium, Hoyle 1948), which thus could appear as an absorption feature in the spectra of high redshift quasars (Davies & Fennison 1964;Gunn & Peterson 1965). Another half a century of search (e.g Penzias & Wilson 1969;Field 1972;Schneider et al. 1991) has not only ultimately uncovered the IGM's existence (e.g. Fan et al. 2006), but also coincidentally revealed the afterglow from the primordial cosmic fireball (CMB, Penzias & Wilson 1965). Various probes (e.g. the disappearance of Lyman α radiation at higher redshifts, Dijkstra 2014; CMB scattering off electrons freed in the EoR, Planck Collaboration et al. 2018; absorption in the spectra of gamma ray bursts and QSOs, Ciardi & Ferrara 2005; as well as the intrinsic 21 cm glow of the neutral IGM, Madau et al. 1997) indicate that the IGM transition from neutral to E-mail: eide@MPA-Garching.MPG.DE ionized occurred more than 12 billion years ago. The details of its occurrence, however, remain an outstanding and heavily investigated question.
Stellar type sources are believed to be the dominant driver of the EoR (Robertson et al. 2015) if a sufficiently large fraction of the ionizing stellar radiation generated in galaxies escapes (e.g. Vanzella et al. 2018;Matthee et al. 2018;Steidel et al. 2018). Nuclear black holes (BHs) are long believed to play a crucial role, particularly in recent years after the discovery of a population of faint quasars (QSOs; Giallongo et al. 2015), indications of an extended and patchy reionization of helium (Worseck et al. 2016(Worseck et al. , 2019 and enormous, opaque troughs in an otherwise partially ionized IGM (Becker et al. 2015; Barnett et al. 2017) that fail to be explained with stellar sources alone (Chardin et al. 2015(Chardin et al. , 2017. Estimating the BHs' role requires sampling of the faint end of their luminosity function at high z (e.g. Parsa et al. 2018;Wang et al. 2019;Kulkarni et al. 2019). An additional source that has been somewhat overlooked, but nevertheless of potential importance, is the thermal Bremsstrahlung from the interstellar medium (ISM), driven by the cooling of supernovashocked gas. It has been observed as a diffuse X-ray component in galaxies (e.g. Mineo et al. 2012b), and its existence may be followed To model the EoR including a variety of sources, we post-process outputs of the hydrodynamic simulation MassiveBlack-II (MBII, Khandai et al. 2015) with the 3D multifrequency radiative transfer (RT) code CRASH (Ciardi et al. 2001;Maselli et al. 2003Maselli et al. , 2009Graziani et al. 2013;Hariharan et al. 2017;Graziani et al. 2018;Glatzle et al. 2019). Our approach to modelling the EoR differs slightly from models where sources are populated solely based on the dark matter halo mass (e.g. Iliev et al. 2006;McQuinn et al. 2007;Trac & Cen 2007;Ross et al. 2017), but it is similar in its multifrequency and post-processing nature (see also Baek et al. 2010;Paardekooper et al. 2013). In comparison to monofrequency approaches in which the rays' properties (and impact) do not change as they propagate through the IGM, a multifrequency RT allows to better capture the impact of radiation on the ionization and thermal state of the IGM (e.g. Ciardi et al. 2012). It should be noted that as our RT is not fully coupled to the cosmological simulations (e.g. Gnedin 2000Gnedin , 2014Norman et al. 2015;Ocvirk et al. 2016;Pawlik et al. 2017;Katz et al. 2018;Rosdahl et al. 2018;Finlator et al. 2018), the ionization and thermal state do not provide feedback and affect subsequent galaxy formation. Our approach, though, is computationally cheaper and thus allows for a more accurate modeling of the RT (e.g. highly resolved spectra of the sources and inclusion of both H and He) and of the ionization and thermal state of the IGM. For completeness, we also mention another category of approaches employed to model cosmic reionization, namely semi-numerical methods (e.g. Mesinger & Furlanetto 2007;Choudhury et al. 2009), which are typically based on an excursion-set approach and thus particularly fast and ideal when a large parameter exploration is required at the expense of a more accurate determination of the physical properties of the IGM. In this respect a good compromise is achieved by methods based on a combination of N-body simulations and post-processing with 1D (rather than 3D) RT codes (e.g. Thomas et al. 2009;Ghara et al. 2018).
Our work is structured as follows. In § 2 we present the details of MBII, the RT and the source modeling. In § 3 we present our results in terms of discussion of the physical state of the IGM as determined by different source types, as well as of comparison to available observations. We discuss the results and give our conclusions in § 4.

METHOD
The approach followed to model cosmic reionization is the same one presented in Eide et al. (2018, in the following Eide2018). Here we outline its main characteristics, and refer the reader to the original paper for more details.
To retrieve the properties of the environment and of the sources of ionizing photons we use the cosmological SPH simulation MassiveBlack-II (MBII, Khandai et al. 2015) with sides of length 100h −1 cMpc 1 . It was run with P-GADGET (see Springel et al. 2005, for details concerning the earlier GADGET-2) using 2 × 1792 3 dark matter and gas particles of mass m DM = 1.1 × 10 7 h −1 M and m gas = 2.2 × 10 6 h −1 M , respectively. The transition of gas into stars is governed by the sub-grid physics presented in Springel & Hernquist (2003) (a comparison to observed star formation rates can be found in Fig. 6 of Khandai et al. 2012 and in Fig. 23 of Khandai et al. 2015). Black holes form following the prescriptions by Di  and Springel et al. (2005). MBII also accounts for the resulting feedback processes (Di Matteo et al. 2008;Croft et al. 2009;Degraf et al. 2010;Di Matteo et al. 2012). We grid the 15 snapshots 2 between z = 18 and 5 onto 256 3 cells, resulting in a resolution of 391h −1 ckpc. This gives us locations of potential ionizing sources through stellar particles, halos and black holes representative, respectively, of stellar populations, galaxies and galactic nuclei in various states of activity. We also grid the gas temperature and density 3 . When multiple source candidates appear within the same cell, we add together their derived spectral shapes and luminosities.
At z < 10 we further reduce the number of sources by evenly absorbing the luminosity of faint cells into all their neighbours that are at least 100 times brighter, if they themselves are not at least 100 times brighter than any neighbour. This procedure is iterated until no more cells can be absorbed into brighter neighbours. The result is a source field with a greatly reduced number of sources. The structure of this new field still closely resembles that of the original one, as only very faint extensions of larger source clusters are integrated into the brighter central cells that would outshine them in any case. Isolated faint sources are left untouched. The combination reduces the number of source by 50% at z = 10 and up to 70% at z = 7. Note that, with this procedure, at z = 7 the vast majority of sources are displaced by less than 1 h −1 cMpc and virtually all remain within 2 h −1 cMpc, i.e. much less than the typical dimension of an ionized region at the same redshift (see also Busch et al. 2020). As a reference, for the simulations with all sources included and f esc = 15% at z = 8.56 we find that clustering the sources induces a difference in the volume averaged temperature and HII and HeII fractions of 0.09%, which increases to 0.3% for the HeIII fraction. These values are within the numerical noise as they are below our convergence limits. A further comparison of cell by cell values indicates that the clustering of the sources does not affect the results presented in the paper.
Next, we detail how the spectral and luminous properties of the sources are obtained. For each source i at a redshift z, we calculate 1 The cosmological parameters adopted in the simulation are taken fom WMAP7 (Komatsu et al. 2011): σ 8 = 0.816, n s = 0.968, Ω Λ = 0.725, Ω m = 0.275, Ω b = 0.046 and h = 0.701. 2 The redshift of the snapshots is dictated by the outputs of the MBII simulations and it is z=18, 16, 14, 13, 12, 11, 10, 9, 8, 7.5, 7, 6.5, 6, 5.5, 5. 3 The gas density is translated to hydrogen and helium number densities assuming number fractions of X = 0.92 and Y = 0.08, respectively, and no metals.  the ionizing emissivity between 13.6 eV and 2 keV as where h P is the Planck's constant, ν is the frequency in Hz, and S i (in erg s −1 Hz −1 ) is the source's spectrum, which is evaluated following the procedure described in Eide2018. As in Eide2018 the emissivities of stars, XRBs and ISM sharing the same cell are summed together and assigned a spectrum which is the sum of the volume averaged spectra of each source type. BHs are instead treated separately, as not every galaxy hosts an active BH. In Fig. 1 we show the spectral energy distribution (SED) of the various sources at different redshifts. Note that this is very similar to Fig. 2 of Eide2018, but has been nevertheless included to show the spectra at all redshift for completeness. As discussed in Eide2018, the evolution in z is extremely mild. The SEDs of the various sources are discussed in the following. In Fig. 2 we show the luminosity function (LF) of the various source types at z = 6, estimated at an absolute magnitude M AB at 1450 Å 4 . We note that the LF is calculated before applying any grid shifting of the sources and estimated on a 1024 3 grid, i.e. with sources grouped together in cells with widths of 143 ckpc, slightly larger than present-day Milky Way-sizes. The LFs are discussed in the following, where we summarize the main features of the four different adopted source types.
Stars: The spectra of stars are obtained from the 2012 version of the stellar population synthesis code BPASS (Eldridge & Stanway 2012) by evaluating a spectrum for each star particle given its stellar mass, age and metallicity. We do not account for contributions from stars that have evolved into binary systems, i.e. we use the 'single   (Bouwens et al. 2015a) and BHs (Giallongo et al. 2015), respectively, while the dotted and double-dashed lines are corresponding fits.
star' version of BPASS, nor for nebular emission. As a reference, at z = 6 we have a total of 3.1 × 10 7 star particles. The volume averages of the individual galaxy spectra are strongest at HI and HeI ionizing frequencies, and fall by orders of magnitude near the HeII ionization threshold, as seen in Fig. 1. The spectra display little redshift evolution towards lower z except for a dampening of the spectra at lower energies and a hardening at energies > 54.4 eV. The LF has a sharp cut-off at both faint (M AB > −10) and bright (M AB < −18) magnitudes, with a plateau and peak near M AB ∼ −15. There is a sizeable contribution from faint stellar sources, as well as a population of sources that are brighter than most BHs, except for the very brightest BHs at M AB ∼ −23. The LF from the simulation matches observational constraints well (Bouwens et al. 2015a), showing a similar distribution of brighter stars and reproducing the sharply rising behaviour of a Schechter LF. At the fainter magnitudes, where no observational constraint is available, the LF decreases, a feature shown to be related to a lower SFR and stellar content in smaller mass halos (see e.g. O'Shea et al. 2015). A reduction in SFR in low-mass galaxies due to photoionization feedback has also been found in fully coupled reionization simulations (Dawoodbhoy et al. 2018).
Nuclear Black Holes (BHs): In MBII each BHs is modelled as a collisionless sink particle of mass 5 × 10 5 h −1 M seeded within newly formed halos with mass above 5 × 10 10 h −1 M . The BH seed grows by accreting surrounding gas or by merging with other BHs. The accretion rate, M i , depends on the BH mass, its velocity relative to the surrounding gas and the density and sound speed of the gas. M i is mildly super-Eddington. Each BH in MBII is assumed to radiate with a bolometric luminosity L BH i = η M i c 2 (in erg s −1 ) given by its accretion rate M i (in g s −1 ) scaled by an efficiency parameter η and the speed of light c (Shakura & Sunyaev 1973). The efficiency parameter η = 0.1 is chosen to be consistent with MBII. For all BHs we adopt an averaged spectrum obtained from observations of 108,104 low-z QSOs (Krawczyk et al. 2013), which is a broken power law with a spectral index of 1 for h P ν > 200 eV. We rescale it with the luminosity of each BH, which rarely exceeds the Eddington value. A notable exception is our brightest BH, found at z = 7, which has an Eddington ratio of 2.9. In general, as can be seen from Fig. 1, the volume average of the BH spectra is brighter than that of other source types. At z = 6 we have a total of 2, 745 BHs. In Fig. 2 the BH LF at z = 6 is shown, along with the observational data (and a fitting LF) from Giallongo et al. (2015). Although we are constrained by the box size and therefore do not sample the rarest, brightest BHs, we do however have a population of faint BHs. Our LF matches the range of brighter magnitudes where observational constraints are available. In a companion paper, by means of a neural network applied to the BH population of the MBII simulation, we re-estimate the faint end of the LF (Eide et al., in prep).
X-ray Binary Systems (XRBs): As with the black holes, we separately obtain a spectrum and rescale it with the individual luminosities of the sources. We look up spectra (evolving only with redshift) and luminosities (scaling with the physical properties of the sources) from the libraries of Fragos et al. (2013a,b), using the updated version presented by Madau & Fragos (2017). In particular, the luminosity scales with properties such as total mass, metallicity and age of the stellar particles, as well as the star formation rate in the halo. Our models include contributions from high-mass XRBs, which trace star formation and are dominant at z 2.5, as well as low-mass XRBs, which trace the stellar mass and dominate at lower z. Their spectra are generally hard, peaking at keV-scales, as shown in Fig. 1. As their AB magnitude is negligibly small, we calculate their z = 6 LF at a wavelength of 6.2 Å, corresponding to photon energies of 2 keV, and include it in Fig. 2. The peak of the 2 keV LF of the XRBs is at a magnitude M 2 keV ∼ 1.
Diffuse thermal bremsstrahlung from shock heated ISM (ISM): Whereas XRBs are point sources of X-rays in galaxies, a diffuse X-ray component has also been observed (Mineo et al. 2012b). We model this as a redshift independent broken power-law spectrum which is flat until the characteristic thermal energy break at k B T ISM = 240 eV (observationally determined by Mineo et al. 2012b, where k B is the Boltzmann's constant), and has a spectral index of 3 at higher energies. As the shock heating process is likely to be associated with supernovae , we scale the ISM luminosities of the galaxies with their star formation rate following Mineo et al. (2012b). As seen in Fig. 1, the spectra are fainter than the stellar ones below HeII ionizing energies, while they are considerably softer than those of XRBs. From Fig. 2 we note that the ISM is much fainter than stars and BHs, as the brightest interstellar gas corresponds to the faintest BHs at M AB ∼ −13.
As we combine scaling relations from the literature with the physical properties obtained from MBII, the only real free parameter in our simulations is the escape fraction of ionizing photons. While in Eide2018 we rescaled the emission of h P ν < 200 eV photons by the constant factor f esc = 15%, observational and theoretical investigations suggest that the averaged escape fraction of ionizing photons may have evolved with redshift along with the changing occurrence of environments that allow for their escape 5 . 5 See e.g. observations from Izotov et al. (2018), Vanzella et al. (2018) and Fletcher et al. (2019), in addition to suggestions of higher f esc due to turbulence by e.g. Safarzadeh & Scannapieco (2016), Kakiichi & Gronke (2019) and Kimm et al. (2019), or of an environment-dependent f esc as seen in simulations from e.g. Paardekooper et al. (2015), Trebitsch et al. (2017) and Katz et al. (2018). A strong dependence is expected also on the evolution . Redshift evolution of the comoving volume averaged ionizing emissivity per source type assuming either a constant UV escape fraction f esc = 15% (solid lines) or redshift-dependent f esc (z) (dashed lines). From top to bottom the sets of curves correspond to: stars (black), BHs (blue), XRBs (purple) and ISM (red). The grey area corresponds to observationally constrained 95% CI for the evolution of the cosmic ionizing emissivity (Bouwens et al. 2015b).
In this work, we examine the effect of varying the escape fraction following the models and Bayesian constraints from Planck presented by Price et al. (2016), where f esc,z=8 is the escape fraction at z = 8. We choose f esc,z=8 = 15% and β = 2.2, where the latter is slightly lower than the Price et al. (2016) lower limit. Our evolving escape fraction is essentially 100% at high-z before decreasing to single-digit values at z < 8. Our choices for the value of the escape fractions have also been guided to yield emissivities consistent with the Bouwens et al.
(2015b) measurements, as can be seen in Fig. 3, where we show the evolution of the volume averaged ionizing emissivity from the various source types. The difference between a constant and varying escape fraction is now reflected on the emissivity, which is higher at z > 8 with f esc (z). Stars dominate the budget at all z, as expected, independently from the adopted escape fraction. There is however an increasingly significant contribution from black holes at z < 10, as their volume averaged emissivity increases by three orders of magnitude from the appearance of the first BH at z = 13, and z = 6. The contribution from the ISM is similar to that of BHs between z = 13 and 11 assuming f esc (z), otherwise it is sub-dominant. The volume averaged emissivity of the XRBs is fainter by yet another order of magnitude compared to that of the ISM.
To investigate which halos contribute most to the ionization budget, in Fig. 4 we plot the stellar emissivity (i.e. produced by all stellar particles contained within a halo), ε stars , as a function of the hosting of the chemical properties of galaxies and in particular on the presence of dust (see e.g. Chisholm   halo mass, M h . Halos are identified with a friends-of-friends procedure and have a minimum mass of ∼ 9 × 10 6 h −1 M . However, we have stars that have been associated (potentially unreliably) with even lower mass halos. We do not consider this to have significant impact on the reliability on the emissivity of the stars, as their properties are derived mainly from the star particles and not the halos. As ε stars is assigned according to a number of physical properties, it has a strong degeneracy with M h . The majority of the sources (90%) are found in the narrow range of halo masses 8.5 < log M h /(M h −1 ) < 9.6 at all redshifts considered here, but their emissivities differ by eight order of magnitudes, indicating the strong contribution from e.g. stellar ages and metallicities to the stellar emissivity. The most massive halos, with log M h /(M h −1 ) > 11, are few in number, but yield exclusively high emissivities. These results are also suggestive that adopting the same escape fraction for all halos is an oversimplification, a topic which will be investigated in more detail in a separate paper.
Finally, we perform multifrequency 3D Monte Carlo radiative transfer simulations with the code CRASH (Ciardi et al. 2001;Maselli et al. 2003Maselli et al. , 2009Graziani et al. 2013Graziani et al. , 2018. We discretize the spectrum of each source into 82 frequency bins between 13.6 eV and 2 keV, spaced closer around the ionization thresholds of hydrogen (13.6 eV) and helium (24.6 and 54.4 eV), both of which CRASH tracks the ionization state of, in addition to the gas temperature 6 . The emitted packets of photons are depleted and redshifted as they propagate (see Eide2018 for a more detailed description), and they are assumed to be lost once they exit the volume. This assumption does not impact the results presented here (see Appendix A for a more extensive discussion).

Figure 5.
Volumetric rendering at z = 6.9 for a simulation with ionizing radiation from stars, BHs, XRBs and the ISM of galaxies, using a UV escape fraction f esc = 15%. Each side of the volume is 100h −1 cMpc long. The blue regions show the distribution of the neutral HI and are drawn at x HI > 0.1, whereas red contours indicate x HeIII > 0.1, which we only find around BHs. The brightest BH is located in the lower right corner.

RESULTS
We first present some key qualitative findings that manifest themselves in all our results ( § 3.1), and then turn to present detailed reionization histories, phase diagrams and thermal properties of the IGM in our simulations with various combinations of different source types and escape fractions ( § 3.2). We finally study the results in light of various observational constraints ( § 3.3).

Qualitative overview
In Fig. 5 we present a 3D rendering of the full (100h −1 cMpc) 3 volume at z = 6.9, when the hydrogen reionization process nears completion and with helium reionization well underway 7 . It shows the HI and HeII ionization state of the IGM for a simulation including all the source types and with f esc = 15%. The volume averaged ionization fraction is x HII = 0.85365 for hydrogen, and x HeIII = 0.00731 for HeII 8 The neutral IGM is perforated by ionized HII regions, not necessarily centered on the brightest sources. Instead, their extent illustrates the effect of the ionization provided by stars. Multiple galaxies contribute jointly to enlarge the HII fronts after individual HII bubbles merge. The morphology thus displays signatures of the percolation inherent to the HI reionization process (Furlanetto & Oh 2016;Bag et al. 2019), a trait which will be investigated in more detail in a companion paper (Busch et al. sub). The HII regions host HeIII bubbles of various extent, most centered on BHs, indicating their responsibility in driving HeII reionization. The figure highlights the effect of our most luminous BH, which distinguishes itself through a large HeIII region. The spikes associated to it are mainly due to recently ionized channels of lower optical HeII depth, but there is also a smaller component associated to numerical artefacts of a yet unconverged, non-equilibrium physical state (see also the physical properties of the environment of the high-z BH modelled by Kakiichi et al. 2017).
In Fig. 6 we show maps of the thermal and ionization state of the IGM at z = 7.5, centered around the brightest BH seen in Fig. 5. The BH resides in a 5.76 × 10 11 M halo, has an AB magnitude at 1450 ÃĚ M AB = −21.9, a mass M BH = 1.85 × 10 7 M and an ionizing emissivity ε = 1.13 × 10 55 phots s −1 . This is not as bright or massive as ULAS J1120+0641 at z = 7.54 (which has M AB = −26.76 ± 0.04 and a mass M BH = 7.8 +3.3 −1.9 × 10 8 M ; Bañados et al. 2018), but it nevertheless belongs to the brightest end of the BH LF that can be sampled with a 100h −1 cMpc volume. The maps illustrate a few important concepts that will be further discussed in the following sections. The extent and temperature of ∼ 10 4 K of the HII regions are determined by the stars. Additional and significant HII heating is only provided when HeII gets fully ionized by the BHs (see also Kakiichi et al. 2017 for further discussion), although partial ionization and heating around fully ionized regions is also provided by the shock heated ISM and, to a lesser extent because of the harder spectrum, the XRBs (see also Graziani et al. 2018). This will have a strong impact on 21 cm HI observations, which will be examined in a companion paper (Ma et al. in prep). Additional heating is also clearly visible within HII regions, in correspondence to a substantial presence of HeIII, i.e. in the immediate surroundings of the more energetic sources. The alternating blue/red pixels correspond to Monte Carlo noise, which is particularly evident in the temperature maps within the HII regions. This happens because it is statistically impossible for two cells with x HII = 1 to have the same exact temperature in simulations with different source types.

Reionization and reheating history
In Fig. 7 we present the reionization and thermal histories for our five combinations of source types. The onset of HI and HeI reionization occurs at z ∼ 16 for all models, and slightly earlier with an evolving escape fraction. With f esc (z) the volume averaged ionization fractions remain higher until z ∼ 7, when f esc falls below 15% and the trend is reversed.
The main driver of the evolution of x HII and x HeII are the stars, with a non negligible contribution from the ISM and, at z < 8, the BHs (see Fig. 3). The ionization fractions differ most at z ∼ 7, but they converge again to similar values towards the end of reionization, when x HII ∼ 1. We also observe that the HeII fraction below z = 7 with stars only is higher than with more energetic sources, as these are starting to produce an appreciable amount of HeIII, visibly depleting HeII.
The onset of HeII reionization is strongly dependent on the spectral hardness of the sources, occurring between z ∼ 11 (when ISM or all sources are included) and z ∼ 9 (with stars only). While XRBs increase x HeIII by a factor of about five compared to stars only, including the contribution from the ISM or the BHs increases the difference to two order of magnitudes, the latter becoming the dominant source of HeII reionization at z < 8. With the exception of the initial and final stages of reionization, the evolution of the volume averaged temperature and ionization fractions are approximately exponential.
By z ∼ 6 both HI and HeI reionization are concluded globally, with x HII ≈ x HeII ∼ 1 (see Table 1). However, the timing of reionization, z reion , varies with location, as it is clear from Fig. 8, where z reion (defined as the redshift at which x HII becomes larger than 0.9) is shown through a slice of our simulation with only stellar sources and f esc = 15% (see Fig. B2 for a comparison with the other simulations). Reionization occurs earlier closer to the sources, supporting an inside-out scenario. Hydrogen in a small fraction of the IGM gets completely ionized at z > 12, including the one surrounding the massive BH discussed earlier, but reionization occurs at z ∼ 7 for the majority of the IGM. Lowering the ionization threshold to x HII ≥ 0.1 does not alter z reion for most of the IGM, except that which is closest to the sources, where z reion could increase from z ∼ 10 to z ∼ 12, as detailed in App. B. We defer a more thorough discussion of the correlation between the timing of reionization and its sources to a separate paper.
Turning to the reheating history, its global onset, evolution and conclusion follows closely that of HI/HeI reionization. The temperature of the HII/HeII regions is determined mainly by stellar type sources, while the harder sources heat the IGM surrounding the fully ionized regions (as seen from the maps in Fig. 6). The volume averaged impact of these energetic sources is minor, but nevertheless visible, as the temperature is systematically higher in their presence (see Table 1). Similarly to what we found for HeIII, the difference with the stellar case only increases with decreasing redshift. Consistently with the trend observed in the ionization fractions, the temperatures are always higher with an evolving escape fraction. From the inset in Fig. 7, at z < 6.5 we observe a clear decline in temperature for all the simulations with f esc = 15%, except when all sources are included, in which case it rises from T = 20, 277 K at z = 6.5, whereas this cooling is less evident with f esc (z). The cooling is indicative of the transition towards an expected thermal asymptote (Hui & Gnedin 1997) which can be upset and postponed by HeII reionization, during which further heating is provided (McQuinn et al. 2009). At z < 6 we see the clear thermal signature of the latter in all our simulations except for those with stars alone, or stars and XRBs, both of which are inefficient HeII ionizers.
In Fig. 9 we plot the distributions of the various quantities characterising the physical state of the IGM at three different redshifts through reionization, z = 9, 7.5 and 6. We find pronounced bimodalities for all quantities except HeIII and while the IGM still has neutral hydrogen and helium. As previously discussed, full ionization of hydrogen and singly ionized helium, and heating of the gas to T 10 4 K, is caused by stellar type sources. The IGM that is not ionized, is coldest and most neutral with stellar type sources. Of the other source types in addition to stars, we find more ionization of HI and HeI and heating with BHs, which is further increased by the presence of XRBs. The most partial ionization and heating, though, is provided by the ISM. The distribution of HeIII is complex. With stars, log x HeIII < −1.4 at z = 9 and 7.5, and higher ionization is barely seen at z = 6. The BHs do not provide large amounts of HeIII until z = 6, when they are responsible for the vast majority of fully ionized HeIII. The XRBs produce HeIII in appreciable quantities, but at lower x HeIII values than the stars, while the ISM provides copious amounts of HeII ionization resulting in two peaks which shift towards higher values with decreasing redshift, but very little fully ionized helium.
We find gas with temperatures and ionization states in between the two peaks in the HII, HeII and T distributions. This gas may be a numerical artefact, as it is unlikely that such gas exists with stellar sources, as discussed by e.g. Ross et al. (2017) and Eide et al. (2018). Whenever the cell size is too large to resolve the sharp ionization  front expected from stellar type sources, the cell containing the front appears partially ionized and warm, as in our simulations, while in reality part of the gas in the cell should be neutral and cold, and part fully ionized and hot. We discuss this issue more in detail in Ma et al. (2020) and Ma et al. (2020b, in prep).
In Fig. 10 we show the distribution of the volume with a given temperature T and free electron fraction where n e and n H are the number densities of free electrons and of hydrogen, while X and Y are the hydrogen and helium number fractions, respectively. As a reference, x e = 1.00 when H is fully ionized and He is fully neutral, x e = 1.09 when H is fully ionized and He is singly ionized, and x e = 1.17 when both H and He are fully ionized. We plot the relations at three different redshifts, z = 9, 7.5 and 6, showing distributions in which the gas is at various stages in transforming from cold (T ∼ 10 K) and largely neutral (x e < 10 −4 ) to hot (T ∼ 10 4 K) and ionized (x e > 1). At z = 9, all the simulations have a significant amount of gas with x e < 10 −2 and temperatures ranging from 10 K to ∼ 10 3 K, but still with significant amounts below T CMB . At z = 7.5, only the simulations with either stars alone or stars and BHs have gas with T < T CMB . At z = 6, all the gas in all the simulations has transitioned into an ionized state. The highest temperatures (above ∼ 10 4.5 K) are reached in regions where both H and He are ionized. Their abundance depends on the spectral shape of the ionizing photons and it is maximum in the presence of BHs. In fact, only BHs are able to fully ionize He (see third column of Fig. 9), although as a whole their contribution to ionization and heating is marginal because of their paucity. A slight amount of gas in this physical state is also present with stars only at z = 7.5 and 6 (see also the upper panel of Fig. 6), and corresponds to cells hosting sources 9 .
The distribution of temperatures in the HII/HeII regions (1 < x e < 1.09) is very similar in all simulations, ranging between ∼8,000 K and ∼20,000 K. The similarity confirms that stars are responsible for the physical state of the gas in such regions. The spread of temperatures observed in Fig. 10 is indicative of the time elapsed since the ionization of the gas. In fact, under purely adiabatic expansion and assuming no recombination, T(z) ∝ (1 + z) 2 . Furthermore, the cooling rate increases for temperature above ∼ 10 4 K. The net effect is that gas that has recently been ionized is hotter than gas that has experienced ionization at an earlier time. This points to the gas temperature being a possible archaeological tracer of the reionization timing, supporting the work of e.g. Keating et al. (2018), who also found recently ionized regions to be hotter than earlier ionized ones.
As the temperature of the neutral and partially ionized IGM is of great importance to the cosmological HI 21 cm signal, we will 9 A negligible fraction of these cells are shock heated. discuss this further in the following. In Fig. 11 we plot the median temperature of gas with x HII < 0.9 against the volume averaged x HII for our different source types and escape fractions. The simulations reach x HII ∼ 0.1 at z = 8-10, when the CMB temperature ranges between 24.5 and 30 K. As expected, little or no heating of this largely neutral gas is observed with stars only. In spite of their harder photons' ability to penetrate deeper into the IGM, the contribution from BHs is negligible (a few degrees more than the stars) due to their paucity. The more numerous XRBs and ISM, individually fainter than BHs (see Fig. 2), are able to raise the temperature to ∼ 70 K and ∼ 220 K, respectively, as x HII approaches 0.9. An additional ∼ 50 K are gained when all sources are taken together. With stars alone, the temperature of the neutral IGM will always be below T CMB , whereas with BHs it will be patchy with temperatures above in the vicinity of BHs and below further away from these sources. As for the XRBs and ISM, we find T gas T CMB at all x HII values plotted in the figure, corresponding to redshifts z 10. This behaviour has important implications for the 21cm signal, which will be presented in Ma et al. (2020b, in prep).
Simulations with f esc = 15% have temperatures higher than those with f esc (z) at a fixed x HII as long as x HII 0.85, while at higher ionization fractions the trend is reversed, as f esc (z) at z < 8 reaches values lower than 15%, resulting in a hardening of the spectra in simulations with a varying escape fraction compared to those with f esc = 15%, i.e. higher temperatures in the neutral parts of the IGM.

Observational constraints
In this section we will discuss our results in the context of available observational constraints.

Abundance of neutral hydrogen
In Fig. 12 we plot the evolution of x HI at z < 12 for simulations including all source types and the two different escape fractions (for reference, in Table 1 the values for all simulations are reported). Although the timing of full reionization is roughly the same in both simulations, there are some differences in the residual neutral fraction towards its conclusion, with x HI = 6.5 × 10 −4 (1.8 × 10 −5 ) at z = 6 for f esc (z) ( f esc = 15%). x HI remains approximately constant from z = 12 until z = 9, when it starts to decrease visibly. Consistently with what observed in Fig. 7, x HI is lower for simulations with f esc (z) until z ∼ 7, and then the trend is reversed, with a difference at z = 6 of about two orders of magnitude, as below z = 6.5 f esc (z) has reached the single-digit percent level.
Our results are consistent with observational data at z 7, reproducing the x HI derived from absorption lines of the z = 7.54 QSO ULAS J1342+0928 by Bañados et al. (2018), the x HI statistically inferred from damping wing absorption in the spectra of ULAS J1120+0641 (Greig et al. 2017), the lower limit of Mason et al.   At z < 7 our results are well within the upper limits from the dark pixel analysis of McGreer et al. (2015), from the LAE clustering by Ouchi et al. (2010), as well as the QSO sample analysis of Gallerani et al. (2008a). Nevertheless, at z ∼ 6 we observe some tension. More specifically, the results with f esc = 15% are marginally consistent with the value of the neutral fraction derived along the line of sight to a GRB at z = 6.29 by Gallerani et al. (2008b), but both are at odds with the Fan et al. (2006) HI trough constraints at z ≤ 6 and the QSO near-zones explored by Bolton & Haehnelt (2007), suggesting e.g. that toward the end of reionization a lower escape fraction would be more appropriate. The simulation with an evolving f esc in fact seems to match better the observational constraints at z 6. We should note that our estimate of the neutral fraction at these redshifts is likely an upper limit, as we do not resolve nor model as sub-grid physics the effect of small-scale Lyman-limit systems (LLS), which are instead crucial in the final stages of reionization as well as in the post-reionization IGM, as discussed e.g. by Madau (2017). The inclusion of LLS would also help to reconcile both our models with the QSO damping wing analysis by Schroeder et al. (2013), which indicates that a substantial part of the IGM must still be neutral at z = 6. Such neutral patches may also be required to explain the long troughs observed in front of QSOs (Chardin et al. 2017).
Finally, our timing of the reionization midpoint, z 50% , when the volume average of the free electron fraction reaches x e = 0.5, is 7.5 for f esc = 15% (independent of source types) and 7.7 for f esc (z), in both cases within the constraints of Planck Collaboration

HII,HeII HeIII
All sources  et al. (2018), where the middle 68th percentile is estimated to be 6.9 < z 50% < 8.1 with a tanh-parametrisation of x e .

IGM temperature
In Table 1 we also report the values of the temperature at the mean density, T 0 . These are slightly higher than the log(T 0 /K)=4.21 ± 0.03 determined by Bolton et al. (2012) (the value is 3.9 ± 0.1 when correcting for the effect of HeII ionization within QSO near zones). It should be noted that the IGM temperature is very sensitive to the spectral shape adopted for the sources. While our approach avoids the introduction of additional degrees of freedom that comes with the assumption of a stellar spectrum (e.g. slope and shape), the results are strongly constrained by the hydrodynamic simulation 10 .
In Fig. 13 we show the cumulative energy deposited as heat per unit mass in the simulations, u 0 , calculated following Boera et al. (2019) and Nasir et al. (2016). This can be used as an independent constraint on reionization histories alongside the gas temperature. Independently from the source type, we observe a rapid increase in u 0 , which roughly corresponds to the midpoint of reionization, z 50% . This period of rapid heating is slightly shorter with a constant rather than a varying escape fraction, although in the latter case u 0 is ∼ 0.2 eV m −1 prot larger at z ∼ 8, consistently with the evolution of the gas temperature. The rate of heating declines at z < 6.6 and it eventually becomes lower with f esc (z) (see inset of the figure). At z = 6 all our simulations produce values of u 0 in the range (5.58 − 5.93) eV m −1 prot , higher than those derived by Boera et al. (2019) from measurements of the Lyα flux power spectrum, with 10 We have tested that adopting a simple power-law spectrum results in lower volume averaged temperatures T and ionization fractions x HII , whereas the neutral and partially ionized IGM is heated and ionized more.  the exception of models with stars only and stars plus BHs with an evolving f esc . On the other hand, our results at z = 6 are all below the u 0 obtained by Boera et al. (2019) from the UV background (UVB) of Puchwein et al. (2019), while at z < 6, they seem to converge towards the Boera et al. (2019) and Puchwein et al. (2019) results at z = 4.6. Our u 0 histories show a more rapid increase in heating than Boera et al. (2019) find in their example models with the underlying Haardt & Madau (2012) UVB. However, the T 0 they derive in their models and from their measurements (T 0 = 7, 600 K), as well as the T 0 = 12, 000 K of Puchwein et al. (2019), are both lower than our values. The more rapid evolution of our u 0 and the higher T 0 can be understood as a consequence of the more rapid redshift evolution of the emissivity in our simulations, compared to those embedded in the UVB underlying the above estimates. Earlier ionization also gives ionized gas time to cool to a thermal asymptote. The majority of our IGM is reionized at later times (as seen in Fig. 8) to temperatures similar to those expected by Boera et al. (2019) (who estimate the gas temperature of recently reionized gas to be 20,000 K) and has hence not undergone a similar cooling.
Our wide range of temperatures in the post-reionized regions is consistent with the thermal properties of the ionized gas in the full radiation-hydrodynamics simulation investigated by D'Aloisio et al. (2019). As mentioned earlier, we find that the source spectral properties strongly affect also the thermal distribution, which is much more uniform in the presence of a power-law spectrum. This dependence on the spectral shape is not found in D'Aloisio et al. (2019), most probably because their spectra are cropped at 4 Ryd. We note that, when employing a full RT, D'Aloisio et al. (2019) find temperatures higher than when using uniform UVB models (such as of Puchwein et al. 2019), or RT simulations with low resolution or monofrequency spectra (e.g. Keating et al. 2018). This is consistent with test cases we have investigated. However, we do also find a wide range of temperatures on all scales in the post-ionization front zones (ranging from ∼ 18, 000 K to 25, 000 K), highlighting that the simulated volume needs to be large enough for the results to be representative ). This does not rule out consistency between our findings and the lower T in works with smaller scales, such as Finlator et al. (2018). Our temperatures, though, are lower than those found by Ross et al. (2019) in their large scale RT simulations including QSOs and XRBs, where they predict T ∼ 200 K at z = 14, compared to the T ∼ 10 K of our models. However, our approaches differ in the way galaxies are populated with QSOs (as they use N-body simulations) and the choice of the spectra, so that a detailed comparison is not straightforward.
Finally, we note that the recent estimated constraints from the 21 cm observations of LOFAR at z ≈ 9.1 (Ghara et al. 2020) find that, if the gas heating remains negligible, a mean ionization fraction of the IGM 0.13 is ruled out. These results align with the temperature inferences from PAPER-64 (Greig et al. 2016) and SARAS 2 (Singh et al. 2018), which rule out a cold reionization scenario. Fig. 11 suggests that the model including all source types is in closest agreement with the LOFAR constraints.

Thomson scattering optical depth
In Fig. 14 we show the optical depth τ due to electron scattering for our simulations together with the 68% CI optical depth τ = 0.054±0.007 evaluated by the Planck Collaboration et al. (2018) and the slightly larger τ = 0.058 ± 0.012 of Planck Collaboration et al. (2016). The optical depth is calculated directly from the ionization fractions of our simulations at z > 5, while we assume x i (3 < z < 5) = x i (z = 5) (where i= HII, HeII, HeIII) and that HeII  Eide et al. (2018) found that during the Cosmic Dawn there is a subtle interplay between various sources of ionizing photons and that their collective impact on the reionization process is not simply given by a sum of the single components. This is particularly relevant for a correct determination of the IGM temperature, but also for partially ionized H and He ionization. Here we thus perform an analysis similar to Eide et al. (2018) concentrating on the Epoch of Reionization, at z 10. This is done by modelling the reionization of hydrogen and helium by post-processing the outputs of the SPH simulation MassiveBlack-II (Khandai et al. 2015) with the 3D radiative transfer code CRASH (e.g. Ciardi et al. 2001;Maselli et al. 2003;Graziani et al. 2013Graziani et al. , 2018, evaluating the physical properties of the IGM as determined by a variety of sources, namely stars, nuclear accreting black holes, X-ray binaries and shock heated ISM. The characteristics of the sources are derived directly from the properties of the stellar particles (such as stellar mass, age, metallicity for stars and XRBs), of the host galaxies (such as the star formation rate for the XRBs and the ISM) and of the BHs. Effectively, the only free parameter of the simulations is the escape fraction of UV photons.

DISCUSSION AND CONCLUSIONS
Unlike most simulations in which the radiative transfer is coupled to the hydrodynamics (e.g. Gnedin 2014;O'Shea et al. 2015;Aubert et al. 2018;Rosdahl et al. 2018), our simulations give results on cosmological scales large enough to provide representative predictions for the progression of reionization , although they are still not large enough to capture the effect of the largest and rarest QSOs and cannot match the large scales (and parameter space) explored by semi-numerical methods (e.g. Mesinger et al. 2011;Fialkov et al. 2014;Hassan et al. 2017). Being less computationally expensive, though, our post-processing approach allows for a much higher sampling of the spectra of the sources in the radiative transfer, which is essential for a correct evaluation of the helium ionization and gas temperature in the presence of sources more energetic than stars. On the other hand, it fails to capture the feedback on structure formation from reionization itself. Our approach to post-processing differs in respect to those of other authors (e.g. Ross et al. 2017;Keating et al. 2018) as we use an SPH environment where the emissivities are not derived from dark matter halo properties, but rather from the properties of the sources as given by the hydrodynamic simulations (in essence, the only free parameter in our approach is the escape fraction of UV photons), in addition to having a very fine sampling of their spectrum.
Our results can be summarized as follows.
• Consistently with previous works and with Eide2018, full hydrogen reionization is driven by stars, which create rapidly growing and overlapping bubbles. The volume average HII fraction at z ∼ 6 is in fact 0.99998 in all the scenarios investigated.
• Photons from more energetic sources (XRBs and ISM) propagate further into the IGM, inducing partial ionization and heating of the gas, in both its hydrogen and helium components, with values x HII ∼ x HeII ∼ 10 −4 − 10 −2 and temperatures as high as ∼ 10 3 K. Due to its soft spectrum, the ISM is particularly efficient. These results are consistent with what discussed by Eide2018 at higher redshift.
• While BHs do not have a strong impact on the average properties of the IGM due to their paucity, whenever active they dominate the local production of ionizing photons, increasing the extent of the HII regions. More importantly, they are the only sources capable of fully ionize helium (as it can be only partially ionized by XRBs and ISM), increasing the volume average HeIII fraction at z ∼ 6 to ∼ 0.02 compared to ∼ 0.0002 with stars only, 0.0009 with XRBs and 0.01 with ISM.
• In the vicinity of energetic sources, where values x HeIII 10 −2 are reached, the gas temperature increases as much as 10 4 K compared to the case in which only radiation from stars is considered. This confirms the importance of including the treatment of the helium component of the gas for a correct evaluation of the temperature. Similarly, because both the ionization level of helium and the value of the temperature are highly dependent on the spectral shape of the ionization field, it is equally crucial to sample the spectrum with a high accuracy.
• Among the models explored here, we find that simulations with an escape fraction that decreases with decreasing redshift reproduce more closely observational data, including the volume averaged emissivity. This also suggests that a constant value lower than the 15% adopted here could give a better fit. It should be noted, though, that the escape fraction is degenerate with properties such as the SFR and the stellar IMF. These are determined by the hydrodynamic simulations and the only free parameter thus remains f esc .
• Our models are consistent with constraints on the Thomson scattering optical depth from Planck, as well as on the volume average HI fraction down to z ∼ 6, while a sub-grid treatment of the LLSs is necessary for a better modelling of the properties of the IGM at lower redshifts.
• Finally, our thermal histories are all dominated by the HII ionization of stars, reaching mean gas temperatures of T = 19, 274 K (18, 643 K with an evolving escape fraction) with stars alone at z = 6. These temperatures only increase by at most ∼ 1, 000 K when including other source types. At z < 6 we see further heating from HeII ionization. Our temperatures are marginally larger than those deduced from observations, such as the temperature at mean density T 0 and the cumulative heating per unit mass u 0 , and we conclude that this mainly is a feature of the redshift evolution of the ionizing emissivity. Important to the HI 21 cm line, the thermal state of the neutral or lowly ionized IGM is extremely sensitive to the presence and abundance of energetic sources. With only stars or additionally the BHs, we expect a 21 cm signal in absorption even at z < 10, whereas we expect a strong (XRBs) and stronger (ISM) signal in emission with the other source types. Observations (PAPER-64, SARAS 2, and LOFAR) favour our model with all source types combined.
Although clear differences in the IGM properties emerge depending on the source combinations, the observational data used in this paper are not able to unambiguously discriminate between different source types and their relative contribution to the reionization process. Additional information is expected from the 21 cm line from neutral hydrogen (Ma et al. in prep), but, to maximize the extraction of information, it will be crucial to cross-correlate this signal with observations in different frequency bands such as of Lyα (e.g. Vrbanec et al. 2020) or [OIII] (Moriwaki et al. 2019) emitters, or of integrated quantities as X-ray (e.g. Ma et al. 2018b) and infrared (e.g. Fernandez et al. 2014) background, and the CMB (e.g. Jelić et al. 2010;Ma et al. 2018a).
Whereas redshifted [OIII] may be observed with JWST, two other hyperfine transitions holds the potential to be probes of the high-z IGM (e.g. Furlanetto et al. 2006). The 3.46 cm transition of 3 He + can be a unique signature of HeIII regions (Khullar et al. 2020), which we have found to overwhelmingly exist only in the vicinity of BHs. Furthermore, this 3 He + signal and the even more elusive DI 91.6 cm transition of deuterium may be observables less prone to the foreground contamination that occludes HI 21 cm signals.
To conclude, our large scale radiative transfer simulations of IGM reionization, anchored to detailed hydrodynamic simulations, give a new insight into the relative role of a variety of source types. The exceptional spectral resolution employed, as well as the inclusion of the helium component of the gas, also assure an accurate evaluation of the IGM temperature. port from the Amaldi Research Center funded by the MIUR program "Dipartimento di Eccellenza" (CUP:B81I18001170001).
We are grateful for the availability of open source software. In this work, we made use of numpy (van der Walt et al. 2011), cython (Behnel et al. 2011), Yt (Turk et al. 2011) and Matplotlib (Hunter 2007).

DATA AVAILABILITY
No new data were generated or analysed in support of this research. While energetic photons are important for partial ionization and heating, their mean free path is very long and thus, with our assumption of non-periodicity in the RT, they could easily escape from the simulation box and be lost 12 . In the following we show that this is not expected to have a significant impact on the results presented in this paper.
In Fig. A1 we plot at which redshift, z req , a photon of energy E needs to be emitted in order to be redshifted to 100 eV (and thus be easily absorbed) by a given redshift z. As a reference, to reach 100 eV by the end of the simulation at z ∼ 5, a photon of 300 eV should have been emitted at z ∼ 17. As the IGM has already been partially or fully ionized well before z = 5 by less energetic photons, this suggests that photons above a couple of hundred eV need times longer than those available in the simulation to possibly play a relevant role.
In Fig. A2 we plot the distribution of the mean free path (MFP) of cells in the simulation with all source types and f esc (z) at z = 7. When only highly neutral cells (i.e. with x HII ≤ 0.9) are considered (top panel), the distributions are strongly peaked and shift to larger MFP with increasing photon energy. For all photons with h P ν < 100 eV the MFP is smaller than the box dimension. When only highly ionized cells (i.e. with x HII > 0.9) are considered (bottom panel), instead, the distributions are more complex, as the presence of He (and its ionization state) becomes more relevant. We then see that the MFP of photons below the ionization threshold of HeII becomes even larger than the box size, as HI and HeI are almost fully ionized, and the corresponding distributions are much wider, reflecting the strong dependence on the detailed ionization state of the various components of the gas. Conversely, higher energy photons have distributions which still resemble those in the left panel, with a strong peak, albeit shifted towards larger values of the MFP due to the reduced absorption from HI and HeI. In Fig. A3 we show the redshift evolution of the MFPs for photons of different energies. These have been calculated by taking the median of the MFPs in cells where hydrogen is highly neutral (x HII ≤ 0.9, first panel) as well as where it is highly ionized (x HII > 0.9, second panel); by taking the mean of the MFPs in all cells (third panel) and finally by inferring the MFP from the volume averaged ionization fractions x HII , x HeII and x HeIII . The results in the highly neutral IGM reflect what was observed in the previous figure, i.e. the MFP smoothly increases with increasing photon energy. All photons with energies below a couple of hundred eV have MFPs shorter than the box size at all redshifts. Conversely, for the highest energy photons this happens only at high redshift. For example, the MFP of a 1 keV photon becomes larger than 100h −1 cMpc at z < 8. This increase in MFP with decreasing redshift observed for all photon energies is associated to the declining gas number density with cosmological expansion and is proportional to 1/z 3 .
In highly ionized regions (second panel), the MFP of photons with energies below 54.4 eV becomes much longer due to the almost lack of HI and HeI, and once reionization is well under way it becomes even larger than the box size. For these photons there is also a stronger

MFP [cMpc]
Inferred from x HII , x HeII , x HeIII and n Figure A3. Redshift evolution of the mean free path (derived from the statistics of the MFPs in each simulation cell) of photons of different energies (as indicated by the line colour) in the simulation with all source types and f esc (z). The grey dashed lines refer to the length of the simulation box. From top to bottom the panels refer to: the median of the MFPs in cells with x HII ≤ 0.9; the median in cells with x HII > 0.9; the volume average of all the MFPs; and the MFPs one would obtain by calculating them directly from the volume averaged number density n and ionization fractions x HII , x HeII and x HeIII of the simulation. dependence on redshift, with a more rapid increase observed from z ∼ 9, when x HII > 10 −1 . Photons with energies above 54.4 eV, instead, have a behaviour similar to the one observed in the previous plot, but the MFPs can be more than one order of magnitude larger. This difference means that even the highest energy photons are sensitive to the presence (or absence) of HI (and HeI), even though the ionizing cross section of hydrogen σ HI at such energies is extremely small (as σ HI ∝ (h P ν/13.6 eV) −3 ).
The volume averaged MFP (irrespective of HI ionization state, as plotted in the third panel) displays how the MFPs transition from being dominated by those found in the neutral HI at high z, to the much longer ones in HII regions at lower z. This increase is most dramatic and extended for HI and HeI ionizing photons, whose MFPs increase by more than five orders of magnitude between z ∼ 14 and the end of the simulation. For higher energy HeII ionizing photons we see the same uptick in MFPs at z < 6, when HeII ionization becomes significant.
Finally, in the last panel we show the MFPs one would obtain by calculating them directly from the volume averaged ionization fractions x HII , x HeII and x HeIII and number density. The MFPs display the same behaviour as in the third panel, albeit with a much shorter transition period at z < 7 from being dominated by MFPs in the neutral HI to the ionized HII at z = 6.
We conclude by observing that even very high energy photons have MFPs similar to our box length, as long as some HI is still present. It should be highlighted, though, that photons above a couple of hundred eV need times longer than those available in the simulation to possibly play a relevant role.

APPENDIX B: REIONIZATION TIMING
We define the reionization redshift z reion of the cell as the redshift at which the cell has reached a hydrogen ionization fraction x HII larger than a threshold value, x th HII . In our reference case we adopt x th HII = 0.9. In Fig. B1 we show how z reion changes when we lower this threshold 0.1 for the simulation with all source types and f esc = 15%. We find  Figure B2. Maps showing the differences in the timing of reionization, z reion , for simulations with different source types (as indicated by the labels), f esc = 15% and an ionization threshold x t h HII = 0.9, with respect to the simulation with only stars. The maps are 100h −1 cMpc wide. The black square highlights a region in which the features resembling Monte Carlo noise are actually due to physical reasons (see text for more details). that z reion increases by as much as ∆z reion ∼ 2 in the vicinity of the sources, while it remains unaltered for the majority of the IGM, with ∆z reion < 0.1 for 50% of the gas, and ∆z reion < 0.4 for 90% of the gas.
In Fig. B2 we show how z reion changes with different source types for a reference value of x th HII = 0.9. We find that reionization in the vicinity of BHs happens earlier, by more than a factor of 0.2 in redshift, whereas XRBs and ISM have a very small impact. Although some of the features observed in all these panels closely resemble Monte Carlo noise, we have verified that this is not the case for all cells and they are instead related to the complexity of the multi-frequency radiative transfer. We highlight one such region with a square in the 'Stars,XRBs' panel of the figure. Indeed some gas pockets experience earlier reionization due to the longer mean free path of the high energy photons emitted by these sources, but at the same time, because such photons ionize less efficiently, for some cells we observe a delayed reionization. The behaviour of the simulation including all source types reflexes that observed in the simulation with stars and BHs, with reionization occurring slightly earlier (∆z reion ≈ 0.1) in the vicinity of BHs. This paper has been typeset from a T E X/L A T E X file prepared by the author.