# The epoch of cosmic heating by early sources of X-rays

The epoch of cosmic heating by early sources of X-rays Abstract Observations of the 21 cm line from neutral hydrogen indicate that an epoch of heating (EoH) might have preceded the later epoch of reionization. Here we study the effects on the ionization state and the thermal history of the intergalactic medium (IGM) during the EoH induced by different assumptions on ionizing sources in the high-redshift Universe: (i) stars; (ii) X-ray binaries (XRBs); (iii) thermal bremsstrahlung of the hot interstellar medium (ISM); and (iv) accreting nuclear black holes (BHs). To this aim, we post-process outputs from the (100 h−1 comoving Mpc)3 hydrodynamical simulation MassiveBlack-II with the cosmological 3D radiative transfer code crash, which follows the propagation of ultraviolet and X-ray photons, computing the thermal and ionization state of hydrogen and helium through the EoH. We find that stars determine the fully ionized morphology of the IGM, while the spectrally hard XRBs pave way for efficient subsequent heating and ionization by the spectrally softer ISM. With the seeding prescription in MassiveBlack-II, BHs do not contribute significantly to either ionization or heating. With only stars, most of the IGM remains in a cold state (with a median T = 11 K at z = 10), however, the presence of more energetic sources raises the temperature of regions around the brightest and more clustered sources above that of the cosmic microwave background, opening the possibility to observing the 21 cm signal in emission. radiative transfer, dark ages, reionization, first stars 1 INTRODUCTION A central, outstanding question in cosmology is the evolution and nature of the hydrogen and helium constituting the intergalactic medium (IGM), particularly at high redshift. Hydrogen was expected to reveal itself either as a deep trough blueward of the H i 21 cm line (Davies & Jennison 1964) or the H i Lyα line (Gunn & Peterson 1965) in the spectra of distant, bright radio galaxies (Field 1962; Penzias & Scott 1968; Allen 1969). However, this did not happen, leading to the conclusion that the IGM must presently be highly ionized and hot. Not until recently was the H i trough finally detected, but at higher redshifts (e.g. Fan et al. 2000, 2006; McGreer, Mesinger & D'Odorico 2015). The IGM must once have been neutral, and then have undergone a phase transformation, the Cosmic Dawn, into its present-day state. Persisting questions are how this transformation occurred, and what drove it. Answering them will shed light on the physical properties of the majority of the baryonic Universe, and open doors to understanding the interplay between the IGM and the large- and small-scale structures that stem from it. The IGM is thus a highly uncharted terrain (see e.g. Meiksin 2009 for a review). In this work, we investigate the properties and signatures of the sources of ionizing, and possibly heating, radiation that flared up as structures formed, marking the transition from the Dark Ages to the Cosmic Dawn (Rees 1998). Early cosmological 21 cm surveys have begun to report that heating may precede full ionization, resulting in a separate epoch of heating (EoH) preceding the epoch of reionization (EoR). The recent constraint on the global signal with SARAS 2 (Singh et al. 2017) disfavours late heating, disentangling the EoH from the EoR (Fialkov, Barkana & Visbal 2014b). These observations also align themselves with other surveys, as the early results of Bowman, Rogers & Hewitt (2008), that were followed by upper limits on the scale and magnitude of brightness temperature fluctuations [e.g. from Giant Metrewave Radio Telescope (GMRT), Paciga et al. 2011; Murchison Widefield Array (MWA), Dillon et al. 2014; Ewall-Wice et al. 2016; Precision Array for Probing the Epoch of Reionization (PAPER), Parsons et al. 2014; Ali et al. 2015; and Low-Frequency Array (LOFAR), Patil et al. 2017]. However, these surveys are impeded at very high redshift by the ionospheric contamination from Earth, which is nearly fully opaque to redshifted 21 cm signal. To mitigate for this, a seminal Chinese–Dutch mission scheduled for 2018 plans to obtain the global 21 cm signal from the far side of the Moon.1 Besides 21 cm tomography, which is still in its infancy after its revival as a probe by Madau, Meiksin & Rees (1997), there are other observational means for examining the ionization and thermal state of the IGM. The spectra of bright, single sources (as quasars or gamma-ray bursts; see e.g. review by Ciardi & Ferrara 2005) reveal the properties of the intervening IGM along the line of sight. The increasing abundance of neutral H i at earlier z has been indicated from large-scale surveys, as well (e.g. Matthee et al. 2015; Zheng et al. 2017; Ouchi et al. 2018 or see review by Dijkstra 2014). Hydrogen ionization is accompanied by production of free electrons. Satellite-based experiments (Komatsu et al. 2011; Planck Collaboration XLVII 2016) have provided tight constraints on the scattering of the cosmic microwave background (CMB) off these electrons, as the CMB becomes dampened and linearly polarized. These observations are however in slight tension, as the continuing lowering of the Thomson scattering optical depth requires a later completion of cosmic reionization. We do not expect a single source type of heating and ionization to be able to account for the diversity in the observed signatures. Kakiichi et al. (2017) instead found the interplay between different ionizing sources to be essential in establishing the thermal and ionization state of the IGM and to leave different signatures. Possible viable candidates include the following. Stars residing within galaxies. These are thought to be the primary drivers of reionization, as argued by Madau, Haardt & Rees (1999). The recent low Thomson scattering optical depth reported by the Planck Collaboration XLVII (2016), combined with constraints from Hubble (Robertson et al. 2015), has rejuvenated the interest in their effect on reionization. Structure formation is also affected by the thermal and ionizing feedback from stars, as discussed by e.g. Couchman & Rees (1986), Cen (1992), and Fukugita & Kawasaki (1994). Depending on their initial mass function, mass and metallicity, they will produce copious amounts of ionizing photons, as shown in the review by Ciardi & Ferrara (2005). Accreting nuclear black holes (or simply black holes, BHs). Hereafter this term will refer to active galactic nuclei (AGNs) and the brighter quasi-stellar objects (QSOs). While these were early candidates for IGM ionization (Rees & Setti 1969; Arons & McCray 1970; Bergeron & Salpeter 1970), assessing their importance is strongly hampered by lacking data on the evolution of the faint end of the QSO luminosity function (QLF) at higher redshifts. With the detection of 22 faint AGNs by Giallongo et al. (2015), the possibility of them producing amounts of ionizing photons sufficient to drive reionization independently from stars was once again examined. Madau & Haardt (2015) and Mitra, Choudhury & Ferrara (2018) found this to be possible, but it would also doubly ionize helium before redshift 4, in contrast to recent observations of an extended He ii reionization at z ∼ 2.8 (Worseck et al. 2016). Qin et al. (2017) report that the models of Madau & Haardt (2015) and Mitra et al. (2018) overpredict the contribution from AGN, and do not find them to dominate over stars, a conclusion also shared by Hassan et al. (2018) and Oñorbe et al. (2017). Updated QLF by Parsa, Dunlop & McLure (2018) or Onoue et al. (2017) also indicates the contribution of BHs to be insufficient to drive reionization. However, the contribution from AGNs was found by Chardin, Puchwein & Haehnelt (2017) and D'Aloisio et al. (2017) to be vital to the large-scale opacity fluctuations observed in the Lyα forest, which was strengthened by the recent observations of a through spanning 240 h−1 comoving Mpc (cMpc) by Barnett et al. (2017). Galactic X-ray binary (XRB) systems comprising a neutron star or a BH devouring a companion star. Among such systems, the majority of the ionizing luminosity at high-z originates from massive [high-mass X-ray binaries (HMXBs)] rather than low-mass [low-mass X-ray binaries (LMXBs)] binary systems (Mirabel et al. 2011; Fragos et al. 2013a,b; Madau & Fragos 2017; Sazonov & Khabibullin 2017). Mineo, Gilfanov & Sunyaev (2012a,b) found the spectra of XRBs to be too hard to account for the soft X-ray flux of galaxies, while they become dominant at higher energies. Thermal bremsstrahlung in the diffuse gas of the heated interstellar medium (ISM). The heating mechanism of the diffuse gas has been a topic under investigation for decades (see e.g. the review by Fabbiano 1989). Interactions between supernova-driven galactic superwinds and clouds (e.g. Chevalier & Clegg 1985) is the preferred explanation (Pacucci et al. 2014). Shocked gas in the galactic halo and disc (e.g. Suchkov et al. 1994), and hot galactic winds (e.g. Strickland & Stevens 2000) are also processes that could yield predominantly soft X-rays. Other candidates, which will not be examined further in this paper because of their secondary role, are e.g. low-energy cosmic rays (Ginzburg & Ozernoi 1966; Nath & Biermann 1993; Sazonov & Sunyaev 2015; Leite et al. 2017), self-annihilation or decay of dark matter (e.g. Liu, Slatyer & Zavala 2016), and plasma beam instabilities in TeV blazars (e.g. Chang, Broderick & Pfrommer 2012; Puchwein et al. 2012). There are different approaches for simulating the physics of the Cosmic Dawn. Ideally, structure formation can be coupled to radiative transfer (RT; e.g. Gnedin 2014; So et al. 2014; O'Shea et al. 2015; Ocvirk et al. 2016; Pawlik et al. 2017; Semelin et al. 2017). However, this is computationally expensive and requires simplifications. For example helium physics is often not treated, or the number of frequency bins used may be limited to a few, albeit the cross-section of hydrogen varies approximately as ν−3. The spatial extent is also often well below the required ∼100 h−1 cMpc needed to provide a consistent reionization history (Iliev et al. 2014). Seminumerical codes as 21cmfast (Mesinger, Furlanetto & Cen 2011) or SimFast21 (Santos et al. 2010; Hassan et al. 2016) can be applied on large scales using an excursion-set formalism (Furlanetto, Zaldarriaga & Hernquist 2004). This allows for fast parameter-space exploration, but may lack the treatment of some important physical processes, such as helium reionization, temperature evolution, partial ionization, and X-ray implementation. Nevertheless, this approach has provided to be rewarding in studies of possible consequences of radiative feedback on structure formation and its 21 cm signatures (see e.g. Fialkov et al. 2014a). Radiative post-processing of hydrodynamical or N-body simulations is also possible (e.g. Baek et al. 2010; Ciardi et al. 2012; Graziani et al. 2015; Ross et al. 2017) and allows for the dedication of computational power to e.g. treat heating, multifrequency RT, as well as the effects of helium. Here we examine the effects of stars, BHs, XRBs, and the ISM on both the thermal and ionization states of the IGM during the EoH. Our approach is the following: we rely on the hydrodynamical cosmological structure formation simulation MassiveBlack-II (MBII, described in Section 2.1; Khandai et al. 2015), which includes baryonic physics and feedback processes, to provide us with the physical environment of the IGM (temperature and gas density), as well as the location and properties of the sources. We assign spectra and ionizing luminosity to star forming particles based on their star formation rates (SFRs), masses, ages, and metallicities, and to BH particles based on their accretion rates. We then perform radiative post-processing of the MBII simulations with the cosmological RT code crash (see Section 2.2) that gives us the evolution of the thermal and ionization state of the IGM. The properties of the ionizing sources are based on empirical relations introduced in Section 2.3. We present our findings in Section 3 and relate them to comparable studies also summarizing our conclusions in Sections 4 and 5. 2 METHODOLOGY Here we describe how we combine the outputs of the cosmological hydrodynamical simulation MBII (Section 2.1) with population synthesis modelling of ionizing sources (Section 2.3), and finally perform multifrequency RT in the same cosmological volume with crash (Section 2.2). 2.1 Cosmological hydrodynamic simulation MBII (Khandai et al. 2015, hereafter K15) is a high-resolution cosmological smoothed particle hydrodynamics (SPH) simulation tracking stellar populations, galaxies, accreting, and dormant BHs, and their properties (as position, age, metallicity, mass, accretion rate, and SFR). The simulation has been run using p-gadget, a newer version of gadget-3 (see Springel 2005, for an earlier version). It accounts for baryonic physics and feedback effects of the sources on their environment following earlier works in its approach to feedback (Di Matteo et al. 2008, 2012; Croft et al. 2009; Degraf, Di Matteo & Springel 2010), subgrid treatment of star formation (Springel & Hernquist 2003), and seeding and evolution of BHs (Di Matteo, Springel & Hernquist 2005; Springel, Di Matteo & Hernquist 2005). The simulation has a box length of 100 h−1 cMpc and is performed in the 7-year Wilkinson Microwave Anisotropy Probe (WMAP7) Λ cold dark matter (ΛCDM) cosmology (Komatsu et al. 2011),2 using 2 × 17923 particles of mass mDM = 1.1 × 107 h−1 M⊙ and mgas = 2.2 × 106 h−1 M⊙ for dark matter and gas, respectively. The adopted gravitational softening length is 1.85 h−1 ckpc. As a reference, we have a total of 28 and 169 520 haloes hosting at least a star particle at z = 18 and 10, respectively, where the lowest dark matter halo masses are 2 × 108 and 9 × 107 h−1 M⊙. The highest dark matter halo masses are 2 × 109 and 1 × 1011 h−1 M⊙ at z = 18 and 10, respectively. The BHs form from seeds of 5 × 105 h−1 M⊙, growing with accretion rates in the range (106–107) h−1 M⊙ Gyr−1 to a maximum mass of 1.4 × 106 h−1 M⊙ at z = 10. The first seed BH is found at z = 13, while we have 17 BHs at z = 10. We refer the reader to K15 for more details on the simulation. We employ six snapshots from MBII, covering the evolution between redshifts z = 18 and 10, each describing the instantaneous state of the simulation. The particle distribution of each snapshot is mapped on to a Cartesian grid of $$N_{\rm c}^3$$ cells to create maps of gas number density and temperature, as well as location and properties of the ionizing sources (see Section 2.3 for a detailed description of how we convert MBII data to ionizing sources). The reference value is Nc = 256, corresponding to a spatial resolution of 391 h−1 ckpc. After gridding, the 45 (292 685) star particles present at z = 18 (10) are reduced to 26 (56 702). These are effectively our sources. When more than one particle ends up in a cell, their properties are summed up (see following section). The gas density is converted to hydrogen and helium number densities by assuming a number fraction X = 0.92 and Y = 0.08, respectively, and no metals. 2.2 Cosmological radiative transfer simulation The RT of ionizing photons is performed by post-processing the outputs of MBII with the multifrequency Monte Carlo ray-tracing code crash (Ciardi et al. 2001; Maselli, Ferrara & Ciardi 2003; Maselli, Ciardi & Kanekar 2009; Graziani, Maselli & Ciardi 2013), which calculates the ionization state of hydrogen and helium, as well as the gas temperature in each grid cell traversed by photons emitted by a radiation source. The version of crash employed here features a self-consistent treatment of ultraviolet (UV) and soft X-ray photons, in which X-ray ionization and heating and detailed secondary electron physics are included (Graziani et al., in preparation). We refer the reader to the original papers for more details on crash. The spectral shape of the sources is discretized into 82 frequency bins, where 72 are in the UV regime below hPν = 200 eV, and the rest are in the soft X-ray regime (0.2–2 keV). The bins are spaced more densely around the ionization thresholds of hydrogen (13.6 eV) and helium (24.6 and 54.4 eV). The radiation emitted by each source i at a redshift z is discretized into Nγ photon packets. Each packet, Nph,i(ν, z) (in units photons Hz−1), holds the total number of photons emitted by source i during the time step Δtem(z) (in s) in the different frequency bins:   $$N_{ {\rm ph}, i}(\nu ,z) = \frac{S_i(\nu ,z)}{L_i(z)}f_{\rm esc}(\nu )\varepsilon _i(z) \Delta t_{\rm em}(z),$$ (1)where fesc(ν) is the frequency-dependent escape fraction, εi(z) is the rate of ionizing photons the source emits (in photons s−1), $$S_i(\nu ,z)/L_i(z)=\hat{S}_i(\nu ,z)$$ is the normalized spectral shape (in Hz−1), where Si(ν, z) defines the spectral energy distribution (SED; in units of erg Hz−1 s−1) and Li(z) the luminosity (in units of erg s−1). We adopt an escape fraction3 of UV photons fesc = 15 per cent for all sources except BHs, for which photons at all frequencies are assumed to escape. Although the adoption of a single value for the escape fraction is an oversimplification as in reality it depends e.g. on the mass of dust and gas and redshift of the host galaxy, its density distribution, the mass and location of the stellar sources (see e.g. Ciardi & Ferrara 2005 and its updated version on arXiv), fesc is to a large degree an unconstrained parameter at low (see e.g. Vanzella et al. 2016; Matthee et al. 2017) and high redshift, where physical conditions could have promoted large values of fesc (see e.g. Kitayama et al. 2004; Yoshida et al. 2007; Safarzadeh & Scannapieco 2016 or the recent z = 4 observations by Vanzella et al. 2018). In a companion paper, focused on the EoR, we will address the effect of different choices of the escape fraction on the global emissivity and reionization process. It should be noted that cells hosting sources are treated as any other cell, i.e. their ionization and temperature evolution is dictated by the properties of the crossing rays. The rate of ionizing photons, or emissivity, can be written as   $$\varepsilon _i(z) = \int _{\rm 13.6\,eV}^{\rm 2\,keV} \frac{S_i(\nu ,z)}{h_{\rm P}^2 \nu } \,{\rm d}(h_{\rm P}\nu ).$$ (2)For simplicity, from now on, we will omit the explicit dependences on ν and z. Each of these packets is followed as it traverses the simulation volume, the photons in different frequency bins ionize, and heat the intervening hydrogen and helium, and finally it either goes extinct or reaches the cosmic volume boundary.4 Because of the large simulated volume, photon redshift is also accounted for. Our reference cases use Nγ = 105 (see Appendix A for a convergence analysis). In the following section we will describe how the spectrum, luminosity, and emissivity of the sources are modelled. 2.3 Sources of ionizing radiation In our simulations we consider four types of ionizing sources: regular stars, hereafter abbreviated only as stars; neutron star/black hole X-ray binaries, hereafter XRBs; thermal bremsstrahlung from supernova-heated ISM, hereafter ISM; accreting nuclear black holes, hereafter BHs. For each source i of any of the types s we model its SED, $$S_i^s$$, and luminosity, $$L_i^s$$, as detailed in the following subsections. When several sources of different types are present in the same cell, we sum up their contributions (except for BHs, which are treated separately) to obtain $$S_i=\sum _s S^s_i$$ and $$L_i=\sum _s L^s_i$$. Using equations (2) and (1), we can then evaluate the photon content of each packet emitted. In Fig. 1 we plot the evolution of the emissivities. In the upper panel we plot the evolution of the median, the 25th–75th percentiles, and the minimum/maximum values of the emissivity for the various source types, as evaluated from equation (2). In the lower panel, we plot the comoving volume averaged emissivities. These reflect the abundance of the sources, unlike the individual emissivities plotted in the upper panel, showing that stars provide the bulk of the ionizing photons throughout the redshift range considered here. Even though individual BHs have higher emissivities than the majority of the stars in galaxies, they are much fewer in number, and are therefore not dominating the overall emissivity budget. The median values of the emissivities for the different source types remain also fairly constant. This applies to the maximum values too: the brightest stellar-type sources at z = 14 are less than an order of magnitude brighter at z = 10. The brightest BHs, on the other hand, do not have higher emissivities compared to the brightest galaxies, which can emit an order of magnitude more ionizing photons. The slow evolution in emissivities is related to the physical properties (e.g. stellar masses, SFRs, and metallicities) that go into determining the SEDs. We note that the median stellar emissivity decreases with decreasing redshift. Although the median stellar mass increases, suggesting higher luminosities, this is counteracted by an increment of both the median stellar mass weighted metallicity and age of the stars. Nevertheless, stars are the dominant producer of ionizing photons at all redshifts. The evolution of the emissivities of the XRBs and the ISM is dictated by the slow evolution of the SFRs. The sudden drop in the lower limit of the emissivity of XRBs is due to galaxies that only hosts the fainter LMXBs. The first BH arises at z = 13, and not until z ≤ 11 do we have more than one present in our volume. Figure 1. View largeDownload slide Upper panel: redshift evolution of the emissivity $$\varepsilon _i^{s}$$ of individual sources (upper panel), grouped per source type s. Stars are plotted in black (single dot), BHs in blue (two dots), XRBs in purple (three dot), and the ISM in red (four dots). The solid, central lines are the median values, enclosed in shaded regions that show the 25th–75th percentile of the values. The span between the minimum and maximum emissivities are given from the dashed lines, with ‘+’ and ‘o’ denoting upper and lower limits, respectively. Note that these values are intrinsic and have not been scaled by an escape fraction. Lower panel: redshift evolution of the volume averaged emissivity εs/V per source type s. Here, the escape fractions are accounted for. Figure 1. View largeDownload slide Upper panel: redshift evolution of the emissivity $$\varepsilon _i^{s}$$ of individual sources (upper panel), grouped per source type s. Stars are plotted in black (single dot), BHs in blue (two dots), XRBs in purple (three dot), and the ISM in red (four dots). The solid, central lines are the median values, enclosed in shaded regions that show the 25th–75th percentile of the values. The span between the minimum and maximum emissivities are given from the dashed lines, with ‘+’ and ‘o’ denoting upper and lower limits, respectively. Note that these values are intrinsic and have not been scaled by an escape fraction. Lower panel: redshift evolution of the volume averaged emissivity εs/V per source type s. Here, the escape fractions are accounted for. While crash can handle a different spectrum for each single source, for the sake of simplicity we adopt at each redshift z an average spectral shape for all sources but BHs. More specifically, for each source i we evaluate $$S_i=S^{\rm stars}_i+S^{\rm XRB}_i+S^{\rm ISM}_i$$. We then use the average $$\bar{S}=\langle S_i\rangle$$ in equation (1) rather than Si. Whenever a BH is present in a cell, this is added as a separate source having the same spatial coordinates. Its spectrum is similarly calculated as $$\bar{S}^{\rm BH}=\langle S^{\rm BH}_i\rangle$$. In Fig. 2 we plot $$\bar{S}$$ at different redshifts z. We see that stars dominate at energies hPν ≲ 60 eV, while the ISM contribution is relevant above the He ii ionization threshold, i.e. into hard UV and the soft X-rays. The XRBs provide the harder X-rays. The weak redshift evolution is due to a combination of the averaging effect effectively preferring brighter sources, and the mentioned slow evolution in the underlying physical properties determining the spectra. Figure 2. View largeDownload slide Globally averaged spectral energy distributions (SEDs) $$\bar{S}(\nu ,z)$$ of the different source types: stars (short dashed lines), XRBs (long dashed), ISM (short–long dashed), total galactic (solid), and BHs (long/double short dashed). The faint, vertical grey lines indicate the ionization thresholds for hydrogen (13.6 eV), neutral helium (24.6 eV), and singly ionized helium (54.4 eV). Different colours indicate different redshifts. Figure 2. View largeDownload slide Globally averaged spectral energy distributions (SEDs) $$\bar{S}(\nu ,z)$$ of the different source types: stars (short dashed lines), XRBs (long dashed), ISM (short–long dashed), total galactic (solid), and BHs (long/double short dashed). The faint, vertical grey lines indicate the ionization thresholds for hydrogen (13.6 eV), neutral helium (24.6 eV), and singly ionized helium (54.4 eV). Different colours indicate different redshifts. In the following sections we will describe in more detail how we evaluate the luminosity $$L^s_i$$ and normalized spectral shape $$\hat{S}^s_i$$ for the various source types. 2.3.1 Stars We model the ionizing radiation from stars by using stellar particles identified in MBII. From the age, metallicity, and mass of each stellar particle p, we obtain $$\hat{S}^{\rm stars}_p$$ and $$L^{\rm stars}_p$$ using the stellar population synthesis code bpass (Eldridge & Stanway 2012). We adopt the instantaneous starburst prescription of star formation and the evolution model that does not account for interactions in binary systems.5 As a single halo i may comprise several particles p, in this case we sum up the contributions from all the particles to obtain $$L_i^{\rm stars} = \sum _{p \in i} L_p^{\rm stars}$$ and $$\hat{S}_i^{\rm stars} = \sum _{p \in i} \hat{S}_p^{\rm stars}$$. 2.3.2 X-ray binaries To account for ionizing radiation coming from XRB systems, we combine the galactic properties provided by MBII with scaling relations from the Fragos et al. (2013a,b) XRB population synthesis model, recently updated by Madau & Fragos (2017). We can thus capture both the metallicity evolution of the HMXBs and the age and stellar mass dependence of the LMXBs, as well as the spectral shape evolution in redshift. For each stellar particle p, we obtain its spectral shape $$\hat{S}^{\rm XRB}_p$$ from Fragos et al. (2013a,b) and its luminosity $$L^{\rm XRB}_p$$ as   $$L^{\rm XRB}_{p} = L^{\rm HMXB}_{p} + L^{\rm LMXB}_{p}.$$ (3)The contribution from HMXBs can be found using equation (3) in Fragos et al. (2013b) with updated coefficients β0–4 from Madau & Fragos (2017):   \begin{eqnarray} \log \left(L_{p}^{\rm HMXB}/{\rm SFR}_i\right) &=& \beta _0 \!+\! \beta _1 Z_p + \beta _2 Z_p^2 + \beta _3 Z_p^3 + \beta _4 Z_p^4 \nonumber \\ &&\times \,({\rm erg s^{-1}} {\rm M}_{\odot }^{-1} \rm yr) \quad {\rm for}\,\, Z \in [ 0,0.025], \end{eqnarray} (4)where SFRi (in M⊙ yr−1) is the star formation rate of the halo i hosting particle p, and Zp is the metallicity of the stellar particle. The contribution to the luminosity from the LMXBs from each stellar particle p is found using equation (4) in Fragos et al. (2013a), also with updated coefficients γ(0–4) from Madau & Fragos (2017):   \begin{eqnarray} \log \left(L_{p}^{\rm LMXB}/M_p\right) &=& \gamma _0 + \gamma _1 \log (t_p/{\rm Gyr}) + \gamma _2 \log (t_p/{\rm Gyr})^2 \nonumber \\ &=& + \,\gamma _3 \log (t_p/{\rm Gyr})^3 + \gamma _4 \log (t_p/{\rm Gyr})^4 \nonumber \\ && \times \,({\rm erg\, s^{-1}} 10^{10} {\rm M}_{\odot }) \quad {\rm for}\,\, t_p \in [0,13.7] {\rm Gyr}, \nonumber\\ \end{eqnarray} (5)where tp (in Gyr) and Mp (in 1010 M⊙) are the age and mass of the stellar particle, respectively. The particles luminosities and spectral shapes are added as described in the previous section to obtain the corresponding source characteristics. Finally note that we adopt coefficients and spectral shapes that are not attenuated from interstellar absorption. 2.3.3 Thermal X-rays from hot ISM We also include ionizing radiation from the diffuse ISM of galaxies. The spectral shape $$\hat{S}_i^{\rm ISM} (\nu )$$ is assumed to be that of thermal bremsstrahlung and constant in redshift (Pacucci et al. 2014):   $$\hat{S}_i^{\rm ISM} (\nu ) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}C & {\rm for }\,\, h_{\rm P} \nu \le kT_{\rm ISM}, \\ C \left(h_{\rm P} \nu /kT_{\rm ISM}\right)^{-3} & {\rm for }\, \,h_{\rm P} \nu > kT_{\rm ISM}, \end{array}\right.$$ (6)where C is a normalization constant that ensures the correct units of Hz−1, and kTISM is the thermal energy of the ISM in eV. From the spectral analysis of observations of the diffuse gas in galaxies by Mineo et al. (2012b), we use kTISM = 240 eV, translating into a characteristic temperature of the heated ISM of ∼106 K. Each halo i has an individual luminosity $$L^{\rm ISM}_i$$, which we evaluate using equation (3) in Mineo et al. (2012b):   \begin{eqnarray} L_{i}^{\rm ISM} ({0.3\hbox{--}10\rm \, keV}) / {\rm SFR}_i &=& \left( 7.3 \pm 1.3 \right) \times 10^{39} \nonumber \\ &&\times\, ( {\rm erg\, s^{-1}} \,{\rm M}_{\odot }^{-1}\,\rm yr), \end{eqnarray} (7)where 0.3–10 keV indicates the photon energy range this relation was obtained for. Note that $$L^{\rm ISM}_i$$ is corrected to be free of interstellar attenuation. We also rescale $$L^{\rm ISM}_i$$ to match our frequency band, which has the lower and upper limits of 13.6 eV and 2 keV, respectively. 2.3.4 Accreting nuclear black holes To account for the ionizing photons originated from accretion discs surrounding nuclear BHs, we identify BH particles in MBII and use their accretion rates to determine the production of ionizing photons. The bolometric luminosity of a BH i is (Shakura & Sunyaev 1973)   $$L_{i}^{\rm BH} = \eta \dot{M}_{i} c^2,$$ (8)where η is an efficiency parameter, $$\dot{M}_{i}$$ is the accretion rate, and c is the speed of light. Consistent with the BH evolution and feedback in MBII, we choose η = 0.1. As a spectral shape, we adopt the observationally derived mean QSO SED of Krawczyk et al. (2013), which is based on 108 104 QSOs sampled at 0.064 < z < 5.46. When no observational data are available between 13.6 and 200 eV, this is derived as interpolation between the mean SEDs for which they have sufficient observations at both higher and lower energies. For energies greater than 200 eV, the spectral shape is modelled as a power law,   $$\hat{S}^{\rm BH}_i \left( h_{\rm P} \nu > 200 \, {\rm eV} \right) \propto \nu ^{-1}.$$ (9)Finally note that no evolution of the SED with redshift is assumed. 3 RESULTS A central question of this study is: can any of our source types drive an EoH separate from an EoR? In other words, is there heating preceding substantial ionization? To examine this, we first discuss the qualitative and quantitative (in terms of average physical quantities) impact of the different source types on the IGM and its ionization and temperature histories (Section 3.1). Second, we focus on the thermal and ionization state of the IGM at z = 12 and 10, to discuss the details of these processes by using global maps, difference maps, and phase diagrams (Section 3.2). For the sake of clarity, throughout this paper we will use the following nomenclature to define the various ionization states of the IGM: neutral when $$x_{\rm H\,\small {II}} < 10^{-5}$$, partially ionized when $$10^{-5} \le x_{\rm H\,\small {II}} < 10^{-1}$$, highly ionized when $$10^{-1} \le x_{\rm H\,\small {II}} < 0.9$$, and fully ionized when $$x_{\rm H\,\small {II}} \ge 0.9$$. 3.1 Global history 3.1.1 Average ionization fractions and temperatures In this section we present both mass and volume average values of various physical quantities, which are summarized in Table 1. Table 1. Thermal and ionization state of the IGM at z = 14, 12, and 10 for different combinations of source types. Note that ionization fractions below 10−5 are denoted with ‘<’. Source type  Ta (K)  $$x_{\rm H\,\small {II}}$$  $$x_{\rm He\,\small {iii}}$$    Median  Volume  Mass  Neutral  Median  Volume  Mass  Median  Volume  Mass      avg.  avg.  avg.    avg.  avg.    avg.  avg.  z = 14  Stars  10  16  23  16  <  7.144 × 10−5  1.514 × 10−4  <  <  <  Stars, BHs  10  16  23  16  <  7.146 × 10−5  1.514 × 10−4  <  <  <  Stars, XRBs  10  16  23  16  <  7.165 × 10−5  1.516 × 10−4  <  <  <  Stars, ISM  10  16  23  16  <  7.201 × 10−5  1.522 × 10−4  <  <  <  Stars, BHs, XRBs, ISM  10  16  23  16  <  7.217 × 10−5  1.523 × 10−4  <  <  <  Stars, BHs, XRBs, ISM (X = 1, Y = 0)  10  17  23  16  <  7.273 × 10−5  1.533 × 10−4  –  –  –  Stars, BHs, XRBs, ISM (X = 0.92, Y = 0)  10  17  24  16  <  8.304 × 10−5  1.720 × 10−4  –  –  –  z = 12  Stars  10  73  142  51  <  1.627 × 10−3  3.456 × 10−3  <  <  <  Stars, BHs  10  73  142  51  <  1.630 × 10−3  3.459 × 10−3  <  <  <  Stars, XRBs  10  73  142  51  <  1.633 × 10−3  3.462 × 10−3  <  <  <  Stars, ISM  11  74  143  52  <  1.644 × 10−3  3.478 × 10−3  <  <  <  Stars, BHs, XRBs, ISM  11  74  143  52  1.004 × 10−5  1.652 × 10−3  3.486 × 10−3  <  <  <  Stars, BHs, XRBs, ISM (X = 1, Y = 0)  11  75  142  54  1.381 × 10−5  1.688 × 10−3  3.545 × 10−3  –  –  –  Stars, BHs, XRBs, ISM (X = 0.92, Y = 0)  11  79  149  56  1.524 × 10−5  1.923 × 10−3  3.960 × 10−3  –  –  –  z = 10  Stars  11  496  886  214  <  2.009 × 10−2  3.581 × 10−2  <  <  1.073 × 10−5  Stars, BHs  11  497  887  214  <  2.012 × 10−2  3.584 × 10−2  <  <  1.238 × 10−5  Stars, XRBs  12  498  888  215  4.442 × 10−5  2.014 × 10−2  3.586 × 10−2  <  <  1.311 × 10−5  Stars, ISM  16  504  897  219  1.008 × 10−4  2.028 × 10−2  3.604 × 10−2  <  2.915 × 10−5  8.005 × 10−5  Stars, BHs, XRBs, ISM  18  506  900  221  1.464 × 10−4  2.036 × 10−2  3.612 × 10−2  <  3.092 × 10−5  8.408 × 10−5  Stars, BHs, XRBs, ISM (X = 1, Y = 0)  19  510  876  244  1.927 × 10−4  2.073 × 10−2  3.664 × 10−2  –  –  –  Stars, BHs, XRBs, ISM (X = 0.92, Y = 0)  19  558  941  252  2.111 × 10−4  2.374 × 10−2  4.100 × 10−2  –  –  –  Source type  Ta (K)  $$x_{\rm H\,\small {II}}$$  $$x_{\rm He\,\small {iii}}$$    Median  Volume  Mass  Neutral  Median  Volume  Mass  Median  Volume  Mass      avg.  avg.  avg.    avg.  avg.    avg.  avg.  z = 14  Stars  10  16  23  16  <  7.144 × 10−5  1.514 × 10−4  <  <  <  Stars, BHs  10  16  23  16  <  7.146 × 10−5  1.514 × 10−4  <  <  <  Stars, XRBs  10  16  23  16  <  7.165 × 10−5  1.516 × 10−4  <  <  <  Stars, ISM  10  16  23  16  <  7.201 × 10−5  1.522 × 10−4  <  <  <  Stars, BHs, XRBs, ISM  10  16  23  16  <  7.217 × 10−5  1.523 × 10−4  <  <  <  Stars, BHs, XRBs, ISM (X = 1, Y = 0)  10  17  23  16  <  7.273 × 10−5  1.533 × 10−4  –  –  –  Stars, BHs, XRBs, ISM (X = 0.92, Y = 0)  10  17  24  16  <  8.304 × 10−5  1.720 × 10−4  –  –  –  z = 12  Stars  10  73  142  51  <  1.627 × 10−3  3.456 × 10−3  <  <  <  Stars, BHs  10  73  142  51  <  1.630 × 10−3  3.459 × 10−3  <  <  <  Stars, XRBs  10  73  142  51  <  1.633 × 10−3  3.462 × 10−3  <  <  <  Stars, ISM  11  74  143  52  <  1.644 × 10−3  3.478 × 10−3  <  <  <  Stars, BHs, XRBs, ISM  11  74  143  52  1.004 × 10−5  1.652 × 10−3  3.486 × 10−3  <  <  <  Stars, BHs, XRBs, ISM (X = 1, Y = 0)  11  75  142  54  1.381 × 10−5  1.688 × 10−3  3.545 × 10−3  –  –  –  Stars, BHs, XRBs, ISM (X = 0.92, Y = 0)  11  79  149  56  1.524 × 10−5  1.923 × 10−3  3.960 × 10−3  –  –  –  z = 10  Stars  11  496  886  214  <  2.009 × 10−2  3.581 × 10−2  <  <  1.073 × 10−5  Stars, BHs  11  497  887  214  <  2.012 × 10−2  3.584 × 10−2  <  <  1.238 × 10−5  Stars, XRBs  12  498  888  215  4.442 × 10−5  2.014 × 10−2  3.586 × 10−2  <  <  1.311 × 10−5  Stars, ISM  16  504  897  219  1.008 × 10−4  2.028 × 10−2  3.604 × 10−2  <  2.915 × 10−5  8.005 × 10−5  Stars, BHs, XRBs, ISM  18  506  900  221  1.464 × 10−4  2.036 × 10−2  3.612 × 10−2  <  3.092 × 10−5  8.408 × 10−5  Stars, BHs, XRBs, ISM (X = 1, Y = 0)  19  510  876  244  1.927 × 10−4  2.073 × 10−2  3.664 × 10−2  –  –  –  Stars, BHs, XRBs, ISM (X = 0.92, Y = 0)  19  558  941  252  2.111 × 10−4  2.374 × 10−2  4.100 × 10−2  –  –  –  aThe neutral average is calculated by weighting the values by the neutral fraction $$x_{\rm H\,\small {I}}$$ of the cell. More neutral cells then have larger impact. This quantity is relevant for studies of the 21 cm signal. View Large First, we turn to the ionization fractions of hydrogen and helium. We limit our discussion to $$x_{\rm H\,\small {II}}$$ and $$x_{\rm He\,\small {iii}}$$, as the behaviour of $$x_{\rm He\,\small {II}}$$ is very similar to the one of $$x_{\rm H\,\small {II}}$$. From z = 18 to 10, they increase logarithmically with redshift, and the main driver is the stars, which dominate the emissivity budget. The scarcity of BHs makes them irrelevant on cosmic scales at z ≥ 10. Similarly, the other source types, albeit being widespread, do not contribute significantly to the statistics as they are far less luminous than stars. At z = 18, the volume and mass averaged ionization fractions are well below our convergence limits of 10−5, independently from the combination of source types. This is true also for the median value (which provides insight into the state of most of the IGM, as also discussed by Ross et al. 2017), as the vast majority of the IGM is fully neutral. At z = 14, the first redshift shown in Table 1, the median value is still below 10−5 for all combinations of source types, while the volume and mass averaged ionization fractions now have increased to ∼10−5 and 10−4, respectively. At z = 12 (i.e. 70 Myr later) the median $$x_{\rm H\,\small {II}}$$ is still below 10−5, except for the case when all sources are present, for which it is 1.0 × 10−5. The volume and mass averaged $$x_{\rm H\,\small {II}}$$ are higher (∼10−3), with slight variations depending on the combination of source types. At z = 12, when all source types are present, we thus observe that the IGM is in the beginning stages of a cosmic, partial ionization phase. Another 100 Myr later, at z = 10, the IGM has been under the influence of a factor of 5 more stellar sources, as well as one additional BH. The ionization fractions have increased by an order of magnitude, reaching a volume average $$x_{\rm H\,\small {II}} \sim 2 \times 10^{-2}$$ independent of source types, and a mass average that is about twice as large. The median hydrogen ionization fraction ranges from below 10−5 (with stars or BHs), to 4 × 10−5 with XRBs, to >10−4 with the ISM. We can understand the small differences in ionization fractions to come from their sensitivity to highly or fully ionized cells, whose prominence is totally dominated by the contribution from stars. The median, on the other hand, is not as prone to biasing by the highly/fully ionized cells. When the ISM is present in addition to the stars, helium is also doubly ionized on cosmic scales, albeit only weakly. At z = 10 the IGM is predominantly fully neutral when only stars and BHs are present, while the majority of the gas is partially ionized in the presence of more numerous and diffuse energetic sources. We now investigate the temperature evolution. At z = 18 the volume and mass averaged temperatures are 11 and 12 K, respectively, irrespective of combination of source types. The median temperature is also comparable, being 10 K. The above statistics are identical to those we obtain directly from the hydrodynamical simulation MBII, i.e. they are unaffected by photoionization. At z = 12, the temperature statistics differ. The volume (mass) average is ∼70 (140) K, independent of the combination of source types, compared to a gas temperature of 38 (78) K in MBII. Weighing the temperature by the hydrogen neutral fraction we find a value ∼50 K. The small differences between combinations of source types indicate that stars are the main driver of the evolution of the averaged temperatures, as they dominate the global emissivity budget. The median temperatures, however, are much lower, around 10 K. Having both XRBs and the ISM present in addition to stars and BHs does not raise the median temperature by z = 12. At z = 10, we have larger differences in the temperature statistics. This coincides with the widespread partial (but low) ionization we found. The volume, mass, and neutral fraction averaged temperatures are now in the range (496–506), (886–900), and (214–221) K, respectively, depending on the increasing combination of source types. The median temperature is however the least biased indicator of the state of the majority of the IGM. With only stars (and BHs), it is 11 K, but with all source types present, it is 18 K. This is an order of magnitude higher than what we obtain from the MBII, where the median temperature at z = 10 is 6 K. To summarize, we find that the IGM at z = 12 is showing its first signs of transitioning into one of two possible states, depending on source populations. With the most conservative assumption, i.e. in the presence of only stars (and BHs), the IGM is mostly fully neutral and cold (10 K), and it remains such down to z = 10. On the other hand, if XRBs and the radiation from the ISM in galaxies are also accounted for, the IGM becomes mainly partially ionized by z = 10. This very low partial ionization is accompanied by a further temperature increase of ∼10 K in all statistics. However, they also show that different regions of the IGM will have different temperatures. This in turn can have observable consequences. Thus, there is a modest EoH, and it largely coincides with large-scale low partial ionization from energetic sources. We now turn to examine the details of this epoch. 3.1.2 Evolution on light-cones To better understand the morphological impact of the source types, in Fig. 3 we illustrate how the hydrogen ionization state, $$x_{\rm H\,\small {II}}$$, evolves in light-cones from z = 16 to 10. Figure 3. View largeDownload slide Light-cones showing the evolution of the ionized hydrogen fraction $$x_{\rm H\,\small {II}}$$ in the simulations with different combinations of source types, as indicated in the labels. The vertical size is 100 h−1 cMpc and the aspect ratio between the axes is given at z = 10, hence the evolution in the angular diameter distance is not taken into account. The combined effect of all source types leaves the IGM without fully neutral regions at z ≲ 11.5. Figure 3. View largeDownload slide Light-cones showing the evolution of the ionized hydrogen fraction $$x_{\rm H\,\small {II}}$$ in the simulations with different combinations of source types, as indicated in the labels. The vertical size is 100 h−1 cMpc and the aspect ratio between the axes is given at z = 10, hence the evolution in the angular diameter distance is not taken into account. The combined effect of all source types leaves the IGM without fully neutral regions at z ≲ 11.5. We find that stars are determining the existence and shape of the fully ionized regions: clusters of galaxies produce ionization bubbles with sharp ionization fronts. These bubbles grow mostly separately in our redshift range. At z = 10, the IGM is still predominantly neutral, but ionized bubbles are present within distances of tens of cMpc between each other. With stars we thus have a clear dichotomy with fully ionized bubbles residing in an otherwise fully neutral IGM. As with the seeding procedure adopted in the MBII simulation only a handful of BHs is present at z > 10, the global history is unaffected by them. The first BH appears at z = 13, while the brightest is at z = 11 with MAB = −17.6. Its ionizing photon production rate is εBH = 2.05 × 1053 photons s−1, its mass is MBH = 1.17 × 106 h−1 M⊙, and it accretes with a rate of $$\dot{M}_{\rm BH} = 10^{-2} \,{\rm M}_{\odot }$$ yr−1. The stars within its host galaxy are even brighter, with MAB = −18.7, and produce εstars = 3.22 × 1053 photons s−1. In a companion paper focused on the EoR we will investigate the BHs impact on the later stages of cosmic reionization, as well as the effect of a different seeding procedure. Including more energetic galactic sources (XRBs and the ISM), we note a significant change in the ionization history of the IGM, which becomes partially ionized ($$10^{-5} < x_{\rm H\,\small {II}} < 10^{-1}$$) several cMpc outside the fully ionized regions. On the other hand, these sources do not significantly contribute to extending the fully ionized regions, which are mainly governed by stars. The spectra of the XRBs peak at ∼2 keV (Madau & Fragos 2017), leaving fewer photons nearer the H i ionization threshold. Their long mean free path means they interact after being redshifted significantly, leaving smooth partial signatures of the order $$x_{\rm H\,\small {II}} \sim 10^{-5}$$ tens of cMpc outside their sources. The spectrum of the ISM is a broken power law, with a high-energy tail fainter than that of the XRBs, but more photons in the hard UV regime, which are more readily absorbed closer to the source. With the ISM, we thus see stronger partial ionization gradients, $$x_{\rm H\,\small {II}} \sim 10^{-3}$$–10−5, and a more patchy ionization morphology compared to XRBs. The combined effect of having all source types present in galaxies is shown in the lowermost panel of Fig. 3. The stars determine the extent of the fully ionized regions that, at z = 10, still make up a very small part of the IGM. The extent of the partially ionized regions goes well beyond those provided by XRBs and the ISM alone. Combined, they leave a much larger fraction of the IGM in a partially ionized state. In Fig. 4 we show the temperature evolution light-cones as well. The stars effectively determine the temperature of 104 K in the regions we recognize from the previous plot to be fully ionized. With additional source types, the temperature of these regions does not markedly change, but we do however have heating of the otherwise cold, T < 10 K, IGM when it is subject to partial and low ionization. Combined, the XRBs and the ISM provide heating that extends further than either can provide alone, indicating that their combined effect is not simply dominated by the largest of each. At z = 10, there is thus a significant difference between the thermal states the IGM can be in, depending on the source types. Figure 4. View largeDownload slide Light-cones showing the evolution of the IGM gas temperature under the presence of different source types as explained in Fig. 3. Figure 4. View largeDownload slide Light-cones showing the evolution of the IGM gas temperature under the presence of different source types as explained in Fig. 3. 3.2 The IGM at $$\boldsymbol {z=12}$$ and $$\boldsymbol {z=10}$$ As we have found that the IGM may transition into a predominantly partially ionized, lightly heated state by z = 10 under the influence of XRBs and the ISM, here we examine its state at z = 12, when this transition has begun. This also coincides with the redshift at which we can investigate the effects of the first BH in our simulation volume, and allows us to quantify its impact, if any. The reader can refer to Table 1 for some numbers. 3.2.1 Morphology In Fig. 5 we plot maps of (from top to bottom) $$x_{\rm H\,\small {II}}$$, $$x_{\rm He\,\small {II}}$$, $$x_{\rm He\,\small {iii}}$$, and temperature at z = 12 for different combinations of source types. These are slices through the volume, chosen to contain the only BH in the simulation box, residing in a bright galaxy in the upper right-hand corner. This BH accretes with a rate of 10−2 M⊙ yr−1 and has a mass of 5.1 h−1 × 105 M⊙. Its absolute magnitude is MAB = −16.2, and it has an ionizing emissivity of εBH = 5.6 × 1052 photons s−1. It resides in a galaxy where the stars have a higher emissivity, εS = 4.6 × 1054 photons s−1, while the XRBs contribute with εXRBs = 1.4 × 1050 photons s−1 and the ISM with εISM = 4.1 × 1051 photons s−1. Figure 5. View largeDownload slide Maps (slices through the volume) in the plane of the only BH at z = 12. The columns indicate different combinations of source types. The upper rows show the quantity in question, and the lower rows show the difference with respect to the first panel, which refer to simulations with stars only. In the maps for $$x_{\rm H\,\small {II}}$$, we also show the location of the sources as circles and their emissivities (denoted with the colour). The only BH is shown as a triangle in the upper right-hand corner. The maps are 100 h−1 cMpc wide. Figure 5. View largeDownload slide Maps (slices through the volume) in the plane of the only BH at z = 12. The columns indicate different combinations of source types. The upper rows show the quantity in question, and the lower rows show the difference with respect to the first panel, which refer to simulations with stars only. In the maps for $$x_{\rm H\,\small {II}}$$, we also show the location of the sources as circles and their emissivities (denoted with the colour). The only BH is shown as a triangle in the upper right-hand corner. The maps are 100 h−1 cMpc wide. We also indicate the position and emissivity of the sources with coloured circles in the maps of $$x_{\rm H\,\small {II}}$$, while the BH is represented by a filled triangle. To highlight the impact of different source types, we also plot absolute differences of the quantities with respect to the simulations with stars only. Note that we create the difference maps by truncating values below our numerical convergence limit of xi = 10−5 for $$i={{\rm H}\,\small {II}}, \, {\rm {He\,\small {II}}}$$, and $${\rm He\,\small {iii}}$$ (see Appendix A). The stars are responsible for the fully ionized H ii bubbles and their morphology, as we also saw in the light-cones. He ii follows H ii due to the similar ionization potentials, but as He i has an ionization cross-section progressively larger than that of H i at photon energies hPν ≥ 24.6 eV, it is then more readily ionized, and its spatial distribution is generally more extended. Also note that sources having in their spectra higher energy photons (see the cases stars+XRBs and stars+ISM), which are redshifted at the simulated scales, provide a contribution first to He i ionization and then to H i. These combined effects give us a larger extent of partially ionized $$x_{\rm He\,\small {II}}$$ than $$x_{\rm H\,\small {II}}$$, while there is still little overlap between the fully ionized regions, leaving a very patchy H ii and He ii ionization morphology. Where sources are more grouped together, they leave imprints on their surroundings similar to what one would expect of separate, but much brighter sources. The extent of ionization must thus be considered an effect of the ionizing photons arising within a volume, rather than from a single source. This helps the partially ionized regions to grow considerably in the cases with XRBs or the ISM. The only BH at z = 12 accounts for partial ionization in its vicinity, spanning several cMpc, beginning with a smoothing of the otherwise sharp ionization front already created by its host and surrounding galaxies. However, a partially ionized tail is not an unique feature of BHs, as it is also seen outside galaxies that have XRBs and the ISM, which both contribute with high-energy photons. As in this epoch, galaxies are more common than BHs, their XRBs and ISM completely dominate the partial ionization at cosmic scales, while the contribution by the ISM component is certainly more ubiquitous than that from the XRBs. The last panels of both H ii and He ii cases finally show the concerted impact of all sources, resulting in a larger overlap of red areas (He ii, in particular, is found in the process of merging into a single one) and in a smoother transition from the fully ionized to the neutral regions. He ii ionization, on the other hand, is much more sensitive to high-energy UV photons (54.4 ≤ hPν ≤ 200 eV), while its cross-section at soft X-ray energies decreases by two orders of magnitude. This is the reason why partially ionized He iii regions with ionizations up to $$x_{\rm He\,\small {iii}} \sim 10^{-3}$$ are found only within the ionized bubble of the BH or correlating with the brightest galactic sources: the strongest XRBs, and the strongest/most clustered galaxies with a hot ISM contribution. Maps of the IGM temperature resulting from the various source combinations are finally found in the bottom panels. The thermal state of the IGM shows signatures both from the hydrodynamic simulation and the heating from the sources. The ideal gas of the hydrodynamic simulations is heated in overdense regions, whereas shocks occur at scales below the resolution of our grid (∼400 h−1 ckpc). As fully ionized H ii and He ii regions are mainly driven by stars, gas at photoionization temperatures of ∼104 K strictly correlates with them, as shown in the bottom left-hand panel by white areas. Extended blue areas, corresponding to temperatures of ∼500 K, are found along filaments connecting the stellar sources. These regions, on the other hand, do not have a clear counterpart in the hydrogen ionization pattern because their ionization fraction is below the convergence threshold xi < 10−5. Difference maps show that the addition of other source types does not significantly change the global temperature pattern as set-up by the stars. A local increase in T (from 100 up to 1000 K) is found only around the BH and around galaxies having a hot ISM. Interestingly, we note that source clustering still plays a relevant role in changing T but its pattern is less extended compared to the one of the corresponding H ii region. Also note that the temperature boost expected in the He iii regions from the release of He ii electrons, as found in Kakiichi et al. (2017), cannot be appreciated around the first quasars because of the lower emissivity and softer spectral shape. In fact, their BH had an emissivity εBH = 1.4 × 1056 photons s−1, and a power-law spectrum with spectral index −1.5, shifting more ionizing photons closer to the ionization potential of He ii. By comparing the relative contribution of different sources we immediately see that XRBs do not heat appreciably the IGM that they partially ionize. We find instead that the ISM emission, mainly providing softer photons, contributes to a heating of some tens of degrees in the regions that it partially ionizes. These do not extend more than a few cMpc outside the fully ionized regions. The different impact of XRBs and ISM can be easily explained in terms of their spectral distributions (see Fig. 2). The energy released by an absorbed X-ray photon is in fact shared between direct photoionization, secondary ionization, and gas heating. The fraction of the energy that goes into heating increases with the ionization fraction of the gas. Heating is more effective for less energetic photons (Dalgarno, Yan & Liu 1999), and the ISM, which has a softer spectrum, is therefore more effective in heating the gas. Combining the effect of XRBs and the ISM, we expect a more efficient heating: XRBs should in fact drive pre-ionization in regions that can be influenced by successive contribution of ISM photons. This concerted effect is clearly visible as an increase of the extent of the heated regions created by the presence of all source types. Also check the difference map to have a better feeling of the difference introduced by the concerted, self-consistent, impact of all the source types. To summarize this section: our models predict an IGM at z = 12 that is dichotomous, with many, patchy fully ionized and hot regions, which extent and abundance are determined by the stars. These regions reside in an IGM that is either neutral and cold, or partially ionized when XRBs and the hot ISM of the galaxies irradiate harder UV and X-ray photons. They can heat the IGM only up to a few tens of degrees, and although this heating is certainly boosted when the XRBs pre-ionize the IGM, it is the ISM that sustains the heating. 3.2.2 IGM phase statistics In Fig. 6 we show the phase state of the IGM at z = 12 under the influence of different combinations of source types. Each panel refers to the distributions of the total 1.68 × 107 cells with a given temperature as a function of overdensity δ.6 Figure 6. View largeDownload slide Phase diagrams showing the distribution (in fraction of the total 2563 = 1.68 × 107 cells) of the IGM at z = 12 in different thermal (T, y-axes), overdensity ($$\delta \equiv n/\bar{n}$$, x-axis), and ionization ($$x_{\rm H\,\small {II}}$$, rows) states for various combinations of source types, as indicated by the labels. The horizontal solid (dashed) lines indicate the mean (median) temperature in the volume, and the vertical solid (dashed) lines indicate the mean (median) overdensity δ of the part of the volume with the given ionization state. The horizontal dotted lines indicate the CMB temperature. The percentages indicate the fraction of the volume that is in the given ionization state. Figure 6. View largeDownload slide Phase diagrams showing the distribution (in fraction of the total 2563 = 1.68 × 107 cells) of the IGM at z = 12 in different thermal (T, y-axes), overdensity ($$\delta \equiv n/\bar{n}$$, x-axis), and ionization ($$x_{\rm H\,\small {II}}$$, rows) states for various combinations of source types, as indicated by the labels. The horizontal solid (dashed) lines indicate the mean (median) temperature in the volume, and the vertical solid (dashed) lines indicate the mean (median) overdensity δ of the part of the volume with the given ionization state. The horizontal dotted lines indicate the CMB temperature. The percentages indicate the fraction of the volume that is in the given ionization state. The upper row shows the highly and fully ionized IGM to be independent from the type of sources. The IGM in this state is hot, with mean and median temperatures of 1.2 × 104 K at a median overdensity δ = 2. The stars that drive full ionization primarily affect the overdense regions. As reflected in the values shown in Table 1, a negligible fraction of the gas is in this high ionization state at z = 12. The middle row shows the phase state of cells with partially ionized IGM. We see that a very small fraction is in this state with stars alone, and adding a single BH does not change this. Including XRBs (ISM) increases the percentage of volume in this partially ionized state to 7 per cent (28 per cent), while having all sources present, 50 per cent of the volume becomes partially ionized. A comparison with the lower row, which refers to the neutral IGM, clearly shows the effect of the additional energetic photons. The IGM under the presence of only stars is mostly cold, but with a fraction of cells that are both overdense (log δ > 0) and hotter than the majority (T > 10 K). A larger fraction of these cells transition to a partially ionized state when including the XRBs. An even larger fraction, both from the cold and underdense and the overdense, hotter gas, becomes partially ionized under the influence of the ISM photons. There is no apparent preference in terms of δ or T on the IGM that becomes partially ionized. This results from a combination of factors such as the location (overdense regions), emissivity (this is lower for the XRBs than for the ISM), and spectral shape (the XRBs have spectra harder than the ISM) of the sources. More than 99 per cent of the IGM is fully neutral with stars only, while this fraction is strongly decreased with XRBs (93 per cent), the ISM (72 per cent), or all source types (50 per cent). It should be noted that, as already pointed out by other authors (see e.g. Ross et al. 2017), the partially ionized and warm cells found in the presence of stellar type sources might only arise from a lack of spatial resolution. More specifically, 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, while in reality part of the gas in the cell should be neutral and cold, and part fully ionized and hot. We will discuss this issue in more detail in a companion paper. In Fig. 7 we present histograms of the temperature and ionization fractions at z = 12 and 10, showing a quantitative evolution of the physical state of the IGM. Figure 7. View largeDownload slide Fraction of cells in different ionization and thermal states for different combinations of source types (stars: solid black, no ticks; stars and BHs: blue, one tick; stars and XRBs: purple, two ticks; stars and ISM: red, three ticks; all: yellow, four ticks). We show the states at z = 12 (upper row) and at z = 10 (lower row). The percentages indicate the fraction of the IGM that is shown in the histogram. The remaining IGM has an ionization fraction xi lower than 10−5 (for $$i={{\rm H}\,\small {II}}, \, {\rm {He\,\small {II}}}$$, and $${{\rm He}\,\small {iii}}$$). Figure 7. View largeDownload slide Fraction of cells in different ionization and thermal states for different combinations of source types (stars: solid black, no ticks; stars and BHs: blue, one tick; stars and XRBs: purple, two ticks; stars and ISM: red, three ticks; all: yellow, four ticks). We show the states at z = 12 (upper row) and at z = 10 (lower row). The percentages indicate the fraction of the IGM that is shown in the histogram. The remaining IGM has an ionization fraction xi lower than 10−5 (for $$i={{\rm H}\,\small {II}}, \, {\rm {He\,\small {II}}}$$, and $${{\rm He}\,\small {iii}}$$). As a baseline, we compare the distributions to the case with stars alone, where the curves for $$x_{\rm H\,\small {II}}$$ and $$x_{\rm He\,\small {II}}$$ are almost flat, with only a surplus of cells in a highly or fully ionized and hot (T ∼ 104 K) state, while the large majority is neutral and cold (T ∼ 10 K), both at z = 12 and 10. The fraction that is highly/fully ionized and hot has increased an order of magnitude by z = 10. With the addition of other source types, partial ionization is strongly increased, especially at $$x_{\rm H\,\small {II}} < 10^{-3}$$. Only at z = 10 a discernible difference is visible in the temperature distributions, with a shift of the low temperature peak by several tens of degrees when all source types are present. This is also the redshift where there is largest difference in the fraction of the IGM that has $$x_{\rm H\,\small {II}} > 10^{-5}$$, denoted as percentages in the figures. At z = 12, 1 per cent of the IGM has $$x_{\rm H\,\small {II}} > 10^{-5}$$ when only stars are present, while this fraction has increased to 50 per cent with all sources. At z = 10, this difference is even more significant, as the percentage goes from 6 to 100 per cent when including either XRBs, the ISM, or both. The behaviour of $$x_{\rm He\,\small {II}}$$ closely follows the one of $$x_{\rm H\,\small {II}}$$. The distribution and presence of doubly ionized helium changes significantly depending on the source types. While the presence of XRBs has an impact only at $$x_{\rm He\,\small {iii}} < 10^{-4}$$, the ISM ionizes more and to higher values of $$x_{\rm He\,\small {iii}}$$. At z = 12, less than 1 per cent of the IGM has $$x_{\rm He\,\small {iii}} > 10^{-5}$$. This fraction, as well as the distribution across all $$x_{\rm He\,\small {iii}}$$, increases by an order of magnitude for all source types at z = 10, leaving 4 per cent of the helium in the IGM in a partially ionized state. 3.3 Simulations without helium The presence of helium in our simulations affects the ionization state of hydrogen and the temperature evolution. The ionization cross-section of hydrogen decays roughly as ν−3, and at the ionization threshold of 24.6 eV for He i, the cross-section of helium is a factor of 6 larger than that of hydrogen. Similarly, at 54.4 eV, the ionization potential of He ii, the helium cross-section is 16 (neutral) and 13 (singly ionized) times larger than that of hydrogen. Sources that have spectra that allow for helium to get ionized twice will result in the release of free electrons that can further heat the IGM, as found by Ciardi et al. (2012). We have run simulations without helium, i.e. Y = 0, including all source types. We have either replaced the now absent helium with hydrogen, effectively increasing its number fraction from X = 0.92 to X = 1, or we have simply not solved for helium chemistry, but kept using X = 0.92. Ionization fractions and temperature statistics are presented in Table 1. When X = 0.92, we effectively have more photons to ionize hydrogen atoms and heat the IGM, resulting in ionization fractions and temperatures higher at any given redshift compared to the case in which helium is included. This effect can be mitigated to some extent by increasing the number fraction of hydrogen to X = 1. In Fig. 8 we have histograms showing the distribution of $$x_{\rm H\,\small {II}}$$ and T at z = 12 and 10. Dashed and dashed–dot lines indicate simulations without helium. The increase in ionization fractions mentioned above is accompanied by a uniform shift at $$x_{\rm H\,\small {II}} < 10^{-1}$$ towards higher $$x_{\rm H\,\small {II}}$$ at both z, while the distributions of highly and fully ionized cells do not change. The absence of helium also contributes to another interesting effect: the maximum temperature reached in the IGM is lowered by approximately 0.2 dex. Without helium, at z = 12 (z = 10) the maximum temperature is T = 22 986 (36 141) K, while with helium it is T = 32 313 (38 709) K. We also have some additional heating of the cold IGM at z = 10, as seen from the increase in the median temperatures in Table 1. Figure 8. View largeDownload slide Distributions of volume averaged temperature and ionization fraction of hydrogen for simulations with (solid; black line) or without helium (dashed; purple lines are for X = 0.92 and Y = 0, while dash–dotted yellow lines are for X = 1 and Y = 0) at z = 12 and 10. The percentages, with the same colours and reference as the lines, indicate what fraction of the IGM has $$x_{\rm H\,\small {II}} \ge 10^{-5}$$, and hence is partially or strongly ionized. Figure 8. View largeDownload slide Distributions of volume averaged temperature and ionization fraction of hydrogen for simulations with (solid; black line) or without helium (dashed; purple lines are for X = 0.92 and Y = 0, while dash–dotted yellow lines are for X = 1 and Y = 0) at z = 12 and 10. The percentages, with the same colours and reference as the lines, indicate what fraction of the IGM has $$x_{\rm H\,\small {II}} \ge 10^{-5}$$, and hence is partially or strongly ionized. 4 DISCUSSION Our simulations show that stars are the principal source of IGM ionization throughout the EoH: in fact they fully dominate the UV part of our spectra (see Fig. 2) and have a large spatial coverage. They are responsible for the volume averaged ionization fractions of both hydrogen and helium, which at z = 10 are $$x_{\rm H\,\small {II}} \sim x_{\rm He\,\small {II}} \sim 2.0 \times 10^{-2}$$. The contribution of stars also drives the shape and extent of local H ii and He ii bubbles, having sharp ionization fronts due to the short mean free path of their UV photons. As a rapid phase transformation is not expected until $$x_{\rm H\,\small {II}} \sim 0.1$$ (Furlanetto & Oh 2016), at z = 10 the IGM is still patchy, with its vast majority being fully neutral, and only a small fraction fully ionized. Only in the presence of more energetic sources, does the morphology and state of the ionized regions change. The seeding procedure adopted in MBII, plants seeds of mass 105 M⊙ in haloes with M > 1010 M⊙. This results in a scarcity of BHs at z > 10, which makes their contribution to IGM ionization and heating negligible. This conclusion, though, is strongly model dependent. We will discuss in more detail the effect of the seeding procedure and the possibility of populating with BHs also in smaller haloes in a companion paper focused on the EoR. It is also important to note that the signature of BHs is degenerate with that of other, more abundant, energetic sources. The XRBs, accounting for both the subdominant LMXBs and the SFR-tracing HMXBs, have spectra that peak at keV scales. They provide a smooth, long-range (several cMpc) partial ionization of hydrogen and helium ($$x_{\rm H\,\small {II}} = x_{\rm He\,\small {II}} \gtrsim 10^{-5}$$). Otherwise, they contribute negligibly, in line with the findings of the semi-analytical works of Madau & Fragos (2017) and Sazonov & Khabibullin (2017). Our model of XRBs thus ionizes less than what found by Ross et al. (2017), who employ a comparable 3D RT post-processing approach, albeit on a dark matter N-body simulation. It should be noted though that their emissivities and source locations were obtained through halo mass scaling relations, which differ from our approach relying on a consistent evolution of baryonic sources (and their properties) from a hydrodynamical simulation. In particular, by z = 13 Ross et al. (2017) find a partial ionization of the cold IGM that is orders of magnitudes higher than what indicated by our computation. They also find a higher volume averaged $$x_{\rm H\,\small {II}} \sim 10^{-2}$$ independent of whether XRBs were included or not. However, while our XRB spectra peak at keV-scales, theirs is a power law with index −1.5. The luminosity of our XRBs scales with either the SFR and metallicity (HMXBs) or the stellar age and mass (LMXBs) in the halo where they reside, while Ross et al. (2017) assume a constant X-ray production efficiency of each halo. Our results are more in line with those of Meiksin et al. (2017), who found the HMXBs to heat the IGM to T = 22 K by z = 10. In their work, though, the increment in temperature from a model with stars only is ∼20 K, compared to our ∼1 K. This could be due to our RT approach capturing only the effect of X-rays up to distances of the order of the box length. On the other hand, Meiksin et al. (2017) do not consider photons with energies below 200 eV, which are more relevant on smaller scales. Note also that we are still likely to overestimate the already small contribution from XRBs. Das et al. (2017), for example, found the XRB spectra to be attenuated by the ISM, albeit not as strongly as advocated by Fragos et al. (2013b), but still more than what we have assumed with our constant escape fraction of fesc = 15 per cent for photons with energies lower than 200 eV. Our multifrequency RT highlighted the importance of the ISM contribution as its spectrum provides large abundance of photons in the hard UV and soft X-rays. We found that this spectral difference is of great importance, a conclusion also shared by Pacucci et al. (2014), or by Fialkov et al. (2014b), who similarly showed that the hardness of the spectral shape of the XRBs would make them inefficient in heating the IGM. There is however some uncertainty on the extent of the contribution from the ISM, which could be larger, as the supernova-driven heating process of the ISM gas will likely result in emission outside the galactic planes (see e.g. Chevalier & Clegg 1985 or Strickland & Stevens 2000). In this case, it is less likely that the ionizing UV photons are strongly absorbed by the host galaxy, which is what we assume with our escape fraction. However, although the production rate of UV photons is uncertain, the ISM is known to supply galaxies ubiquitously with observable amounts of soft X-rays (Mineo et al. 2012b). The Meiksin et al. (2017) work mentioned above includes also the contribution of the ISM, which raises the IGM temperature at z = 10 to T = 6 or 34 K (compared to our median T = 16 K with stars and the ISM), depending on the wind model adopted. One should be careful to attribute long range partial ionization and heating signatures, if observed, solely to the ISM. These could possibly be degenerate with those from cosmic rays. While Leite et al. (2017) found that they should heat the IGM in the immediate vicinity of star-forming haloes, this happens because of the cosmic rays confinement. If this confinement is weaker than what estimated by Leite et al. (2017), their contribution could be very similar to that of the ISM. Another RT study examining the effect of X-rays on heating and ionization is that of Baek et al. (2010). They include X-rays from all of their sources, labelling this contribution as QSOs, but noting that XRBs and supernova remnants also fall into this category. Their scaling of the X-ray luminosity with the SFR allows for some comparison to our case where all source types are present. Their mean ionization fraction at z = 10 is $$x_{\rm H\,\small {II}} \sim 10^{-2}$$ for simulations with helium, but without X-rays. They note that including QSOs that account for up to 10 per cent of the total flux does not change this result, similarly to the findings of Ross et al. (2017), and now also ours. We have run simulations without helium, finding that the ionization fraction, contrary to Baek et al. (2010, cf. their models S1 and S3), increases to $$x_{\rm H\,\small {II}} = 2.4 \times 10^{-2}$$ (keeping X = 0.92), or $$x_{\rm H\,\small {II}} = 2.1 \times 10^{-2}$$ (with X = 1). This is caused by the additional photons that are available for H i ionization when helium is absent, aligning our results with those of Ciardi et al. (2012). Including helium thus provides a slower progression of the cosmic heating and reionization. As for He iii, we find that only faint, diffuse regions having $$x_{\rm He\,\small {iii}} \sim 10^{-3}$$ start to appear near bright energetic sources. We do not find the widespread existence of highly/fully ionized He iii that Ross et al. (2017) found at higher redshifts. Consistently with Madau & Fragos (2017), we find that XRBs are not able to heat significantly the IGM due to their hard spectra. Our results show that also the ISM, which has a much softer spectrum, is not able to raise the IGM temperature significantly by z = 10. Our median temperatures are in fact T = 11 and 18 K for simulations with stars only and all sources, respectively. This is somewhat higher than the value reported in Baek et al. (2010), where the temperature of the IGM with an ionization fraction $$x_{\rm H\,\small {II}} < 10^{-2}$$ was ranging between 3 and 7 K without and with QSOs, respectively. Our volume averaged temperatures are, on the other hand, much lower than those of Ross et al. (2017), who found T ∼ 600 K at z = 13 with stars only and T ∼ 1000 K with HMXBs, compared to our volume averaged T = 34 K at z = 13 with stars and XRBs. Predicting the thermal state of the IGM during the EoH is of crucial importance for 21 cm experiments because an IGM with a temperature higher than that of the CMB will be seen in emission, and if lower, it will be seen in absorption. Interestingly, at z = 10 the CMB temperature is TCMB = 30 K. The margins between a 21 cm signal seen in emission or absorption are thus close. Note, though, that while our mean temperatures at this z are much higher than TCMB, they do not represent the state of the vast majority of the IGM, which is either fully neutral (with stars only) or partially ionized (with all sources). We find partially ionized regions that are warmer than the CMB already at z = 12, as indicated from the neutral and mass averaged temperatures. We defer to a future work further investigation of the effect that different source types leave on the 21 cm signal. Observationally, very late (after z = 8.4) heating of the IGM has been ruled out by measurements of 21 cm line (Pober et al. 2015). Less constrained is a period of early and highly efficient heating, which would suppress small-scale structure formation by increasing the Jeans mass (Ostriker & Gnedin 1996), disfavouring H2 creation associated with X-ray fluxes that do not dominate the ionizing budget (Oh 2001), and increase the Thomson-scattering optical depth of the CMB (Ricotti & Ostriker 2004). We observe that both the distribution of the sources and their surrounding ionized regions trace the underlying density field. This can be seen particularly well in the phase maps of Fig. 6, where the fully ionized regions reside where log δ > 0. We also see clustering of the sources in the global maps (Fig. 5). This clustering induces (a) ionized regions that are more extended and (b) heating of the nearby neutral IGM, when comparing the surroundings of sources with similar emissivities but that exist either in a clustered environment or alone. This distinction (i.e. clustered sources versus single source) is crucial in the interpretation of observations aiming at investigating the physical properties of the high-redshift IGM. 5 CONCLUSION In this paper we have examined the sources that could end the Dark Ages, initiate the Cosmic Dawn, and drive the EoH, when the IGM starts transitioning from a cold and neutral phase to a hot and ionized one. The Cosmic Dawn concludes with the EoR at z < 10, which we will explore in a companion paper. We have investigated how stars, XRBs, and the shock heated ISM of galaxies, as well as accreting nuclear BHs, could heat and ionize the IGM at cosmic scales. Each source type has spectral characteristics determined by underlying physical processes. The luminosity of the stars is determined by their masses, ages, and metallicities, while the ISM and XRBs are sensitive to the SFR in galaxies, as well as their masses and metallicities. The BHs, on the other hand, shine with the rate they accrete matter. We have used the hydrodynamical simulations MBii (K15) to provide us with the physical properties of the sources, their location, and abundance, as well as their environment (temperature and baryonic density). Having established the sources, their spectra, and ionizing emissivities, we post-processed MBII with the RT code crash (e.g. Ciardi et al. 2001; Graziani et al. 2013). Our main findings can be summarized as follows. Stars drive the extent, shape, abundance, and temperature of the fully ionized H ii and He ii regions, but do not leave the remaining IGM in any heated or partially ionized state. Only 7 per cent of the overall IGM is in a non-neutral state at z = 10 with stars alone. With our seeding prescription, nuclear BHs are scarce and do not contribute significantly to heating or ionization at z > 10. XRBs contribute to partial, uniform ionization ($$x_{\rm H\,\small {II}} \sim 10^{-5}$$) of the dilute IGM on cMpc scales, but do not significantly heat it by z = 10. This can be attributed to their hard, keV-peaked spectra. The ISM of the galaxies contributes to a larger extent of partial ionization than the XRBs do. Their softer spectra provide heating and ionization, with $$x_{\rm H\,\small {II}} \sim 10^{-3}$$ and T ≳ 10 K up to a few cMpc around the fully ionized regions at z = 12 and lower partial ionization further out. In concert, the ISM and the XRBs induce an ionization fraction $$x_{\rm H\,\small {II}} \ge 10^{-5}$$ in 50 per cent (100 per cent) of the IGM at z = 12 (z = 10). However, at z = 10 the IGM is still predominantly cold, with median temperatures in the range (11–19) K under the influence of stars to all source types. At z = 10 the IGM will be seen in both emission and absorption in 21 cm as the neutral averaged temperatures are T > 200 K for all combinations of source types. Helium reionization has begun on cosmic scales at z = 10 in the presence of XRBs and/or the ISM. If an EoH takes place at z > 10, our findings indicate that it will be a modest one. Acknowledgements We thank Rupert Croft for enlightening comments on the bright BHs; Andrea Ferrara for rewarding discussions; Garrelt Mellema, Philipp Busch, Martin Glatzle, Max Gronke, and Aniket Bhagwat for their insights into various parts of the project. MBE thanks Piero Madau for his hospitality. We thank the anonymous referee for constructive comments. LG acknowledges the support of the DFG Priority Program 1573 and from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement n. 306476. This work made use of a number of open source software packages, in particular numpy (Van Der Walt, Colbert & Varoquaux 2011), cython (Behnel et al. 2011), matplotlib (Hunter 2007), and GNU parallel (Tange 2011). Footnotes 1 https://www.isispace.nl/dutch-radio-antenna-depart-moon-chinese-mission/; accessed 2017 July 13. 2 σ8 = 0.816, ns = 0.968, $$\Omega _\Lambda =0.725$$, Ωm = 0.275, Ωb = 0.046, h = 0.701. 3 crash allows us to specify the escape fraction in different frequency bands individually for each source. 4 In our reference simulations we do not use periodic boundary conditions, although this option is available. 5 The effect of not including the evolution of binary systems is a reduction in ionizing flux. 6 In this paper the gas overdensity is defined as $$\delta \equiv n/\bar{n}$$, where n is the particle number density with units cm−3. REFERENCES Ali Z. S. et al.  , 2015, ApJ , 809, 61 https://doi.org/10.1088/0004-637X/809/1/61 CrossRef Search ADS   Allen R. J., 1969, A&A , 3, 382 Arons J., McCray R., 1970, Astrophys. Lett. , 5, 123 Baek S., Semelin B., Di Matteo P., Revaz Y., Combes F., 2010, A&A , 523, A4 https://doi.org/10.1051/0004-6361/201014347 CrossRef Search ADS   Barnett R., Warren S. J., Becker G. D., Mortlock D. J., Hewett P. C., McMahon R. G., Simpson C., Venemans B. P., 2017, A&A , 601, A16 https://doi.org/10.1051/0004-6361/201630258 CrossRef Search ADS   Behnel S., Bradshaw R., Citro C., Dalcin L., Seljebotn D. S., Smith K., 2011, Comput. Sci. Eng. , 13, 31 https://doi.org/10.1109/MCSE.2010.118 CrossRef Search ADS   Bergeron J., Salpeter E. E., 1970, Astrophys. Lett. , 7, 115 Bowman J. D., Rogers A. E. E., Hewitt J. N., 2008, ApJ , 676, 1 https://doi.org/10.1086/528675 CrossRef Search ADS   Cen R., 1992, ApJS , 78, 341 https://doi.org/10.1086/191630 CrossRef Search ADS   Chang P., Broderick A. E., Pfrommer C., 2012, ApJ , 752, 23 https://doi.org/10.1088/0004-637X/752/1/23 CrossRef Search ADS   Chardin J., Puchwein E., Haehnelt M. G., 2017, MNRAS , 465, 3429 https://doi.org/10.1093/mnras/stw2943 CrossRef Search ADS   Chevalier R. A., Clegg A. W., 1985, Nature , 317, 44 https://doi.org/10.1038/317044a0 CrossRef Search ADS   Ciardi B., Ferrara A., 2005, Space Sci. Rev. , 116, 625 https://doi.org/10.1007/s11214-005-3592-0 CrossRef Search ADS   Ciardi B., Ferrara A., Marri S., Raimondo G., 2001, MNRAS , 324, 381 https://doi.org/10.1046/j.1365-8711.2001.04316.x CrossRef Search ADS   Ciardi B., Bolton J. S., Maselli A., Graziani L., 2012, MNRAS , 423, 558 https://doi.org/10.1111/j.1365-2966.2012.20902.x CrossRef Search ADS   Couchman H. M. P., Rees M. J., 1986, MNRAS , 221, 53 https://doi.org/10.1093/mnras/221.1.53 CrossRef Search ADS   Croft R. A. C., Di Matteo T., Springel V., Hernquist L., 2009, MNRAS , 400, 43 https://doi.org/10.1111/j.1365-2966.2009.15446.x CrossRef Search ADS   Dalgarno A., Yan M., Liu W., 1999, ApJS , 125, 237 https://doi.org/10.1086/313267 CrossRef Search ADS   D'Aloisio A., Sanderbeck P. R. U., McQuinn M., Trac H., Shapiro P. R., 2017, MNRAS , 468, 4691 https://doi.org/10.1093/mnras/stx711 CrossRef Search ADS   Das A., Mesinger A., Pallottini A., Ferrara A., Wise J. H., 2017, MNRAS , 469, 1166 https://doi.org/10.1093/mnras/stx943 CrossRef Search ADS   Datta K. K., Mellema G., Mao Y., Iliev I. T., Shapiro P. R., Ahn K., 2012, MNRAS , 424, 1877 https://doi.org/10.1111/j.1365-2966.2012.21293.x CrossRef Search ADS   Davies R. D., Jennison R. C., 1964, MNRAS , 128, 123 https://doi.org/10.1093/mnras/128.2.123 CrossRef Search ADS   Degraf C., Di Matteo T., Springel V., 2010, MNRAS , 402, 1927 https://doi.org/10.1111/j.1365-2966.2009.16018.x CrossRef Search ADS   Dijkstra M., 2014, Publ. Astron. Soc. Aust. , 31, 26 https://doi.org/10.1017/pasa.2014.33 CrossRef Search ADS   Dillon J. S. et al.  , 2014, Phys. Rev. D , 89, 023002 https://doi.org/10.1103/PhysRevD.89.023002 CrossRef Search ADS   Di Matteo T., Springel V., Hernquist L., 2005, Nature , 433, 604 https://doi.org/10.1038/nature03335 CrossRef Search ADS PubMed  Di Matteo T., Colberg J., Springel V., Hernquist L., Sijacki D., 2008, ApJ , 676, 33 https://doi.org/10.1086/524921 CrossRef Search ADS   Di Matteo T., Khandai N., DeGraf C., Feng Y., Croft R. A. C., Lopez J., Springel V., 2012, ApJ , 745, L29 https://doi.org/10.1088/2041-8205/745/2/L29 CrossRef Search ADS   Eldridge J. J., Stanway E. R., 2012, MNRAS , 419, 479 https://doi.org/10.1111/j.1365-2966.2011.19713.x CrossRef Search ADS   Ewall-Wice A. et al.  , 2016, MNRAS , 460, 4320 https://doi.org/10.1093/mnras/stw1022 CrossRef Search ADS   Fabbiano G., 1989, ARA&A , 27, 87 https://doi.org/10.1146/annurev.aa.27.090189.000511 CrossRef Search ADS   Fan X. et al.  , 2000, AJ , 120, 1167 https://doi.org/10.1086/301534 CrossRef Search ADS   Fan X. et al.  , 2006, AJ , 132, 117 https://doi.org/10.1086/504836 CrossRef Search ADS   Fialkov A., Barkana R., Pinhas A., Visbal E., 2014a, MNRAS , 437, L36 https://doi.org/10.1093/mnrasl/slt135 CrossRef Search ADS   Fialkov A., Barkana R., Visbal E., 2014b, Nature , 506, 197 https://doi.org/10.1038/nature12999 CrossRef Search ADS   Field G. B., 1962, ApJ , 135, 684 https://doi.org/10.1086/147312 CrossRef Search ADS   Fragos T. et al.  , 2013a, ApJ , 764, 41 https://doi.org/10.1088/0004-637X/764/1/41 CrossRef Search ADS   Fragos T., Lehmer B. D., Naoz S., Zezas A., Basu-Zych A., 2013b, ApJ , 776, L31 https://doi.org/10.1088/2041-8205/776/2/L31 CrossRef Search ADS   Fukugita M., Kawasaki M., 1994, MNRAS , 269, 563 https://doi.org/10.1093/mnras/269.3.563 CrossRef Search ADS   Furlanetto S. R., Oh S. P., 2016, MNRAS , 457, 1813 https://doi.org/10.1093/mnras/stw104 CrossRef Search ADS   Furlanetto S., Zaldarriaga M., Hernquist L., 2004, ApJ , 613, 1 https://doi.org/10.1086/423025 CrossRef Search ADS   Giallongo E. et al.  , 2015, A&A , 578, A83 https://doi.org/10.1051/0004-6361/201425334 CrossRef Search ADS   Ginzburg V., Ozernoi L., 1966, SvA , 9, 726 Giri S. K., Mellema G., Dixon K. L., Iliev I. T., 2018, MNRAS , 473, 2949 https://doi.org/10.1093/mnras/stx2539 CrossRef Search ADS   Gnedin N. Y., 2014, ApJ , 793, 29 https://doi.org/10.1088/0004-637X/793/1/29 CrossRef Search ADS   Graziani L., Maselli A., Ciardi B., 2013, MNRAS , 431, 722 https://doi.org/10.1093/mnras/stt206 CrossRef Search ADS   Graziani L., Salvadori S., Schneider R., Kawata D., de Bennassuti M., Maselli A., 2015, MNRAS , 449, 3137 https://doi.org/10.1093/mnras/stv494 CrossRef Search ADS   Gunn J. E., Peterson B. A., 1965, ApJ , 142, 1633 https://doi.org/10.1086/148444 CrossRef Search ADS   Hassan S., Davé R., Finlator K., Santos M. G., 2016, MNRAS , 457, 1550 https://doi.org/10.1093/mnras/stv3001 CrossRef Search ADS   Hassan S., Davé R., Mitra S., Finlator K., Ciardi B., Santos M. G., 2018, MNRAS , 473, 227 https://doi.org/10.1093/mnras/stx2194 CrossRef Search ADS   Hunter J. D., 2007, Comput. Sci. Eng. , 9, 99 https://doi.org/10.1109/MCSE.2007.55 CrossRef Search ADS   Iliev I. T., Mellema G., Ahn K., Shapiro P. R., Mao Y., Pen U. L., 2014, MNRAS , 439, 725 https://doi.org/10.1093/mnras/stt2497 CrossRef Search ADS   Kakiichi K., Graziani L., Ciardi B., Meiksin A., Compostella M., Eide M. B., Zaroubi S., 2017, MNRAS , 468, 3718 https://doi.org/10.1093/mnras/stx603 CrossRef Search ADS   Khandai N., Di Matteo T., Croft R., Wilkins S., Feng Y., Tucker E., DeGraf C., Liu M.-S., 2015, MNRAS , 450, 1349 (K15) https://doi.org/10.1093/mnras/stv627 CrossRef Search ADS   Kitayama T., Yoshida N., Susa H., Umemura M., 2004, ApJ , 613, 631 https://doi.org/10.1086/423313 CrossRef Search ADS   Komatsu E. et al.  , 2011, ApJS , 192, 18 https://doi.org/10.1088/0067-0049/192/2/18 CrossRef Search ADS   Krawczyk C. M., Richards G. T., Mehta S. S., Vogeley M. S., Gallagher S. C., Leighly K. M., Ross N. P., Schneider D. P., 2013, ApJS , 206, 4 https://doi.org/10.1088/0067-0049/206/1/4 CrossRef Search ADS   Leite N., Evoli C., D'Angelo M., Ciardi B., Sigl G., Ferrara A., 2017, MNRAS , 469, 416 https://doi.org/10.1093/mnras/stx805 CrossRef Search ADS   Liu H., Slatyer T. R., Zavala J., 2016, Phys. Rev. D , 94, 063507 https://doi.org/10.1103/PhysRevD.94.063507 CrossRef Search ADS   McGreer I. D., Mesinger A., D'Odorico V., 2015, MNRAS , 447, 499 https://doi.org/10.1093/mnras/stu2449 CrossRef Search ADS   Madau P., Fragos T., 2017, ApJ , 840, 39 https://doi.org/10.3847/1538-4357/aa6af9 CrossRef Search ADS   Madau P., Haardt F., 2015, ApJ , 813, L8 https://doi.org/10.1088/2041-8205/813/1/L8 CrossRef Search ADS   Madau P., Meiksin A., Rees M. J., 1997, ApJ , 475, 429 https://doi.org/10.1086/303549 CrossRef Search ADS   Madau P., Haardt F., Rees M. J., 1999, ApJ , 514, 648 https://doi.org/10.1086/306975 CrossRef Search ADS   Maselli A., Ferrara A., Ciardi B., 2003, MNRAS , 345, 379 https://doi.org/10.1046/j.1365-8711.2003.06979.x CrossRef Search ADS   Maselli A., Ciardi B., Kanekar A., 2009, MNRAS , 393, 171 https://doi.org/10.1111/j.1365-2966.2008.14197.x CrossRef Search ADS   Matthee J., Sobral D., Santos S., Röttgering H., Darvish B., Mobasher B., 2015, MNRAS , 451, 400 https://doi.org/10.1093/mnras/stv947 CrossRef Search ADS   Matthee J., Sobral D., Best P., Khostovan A. A., Oteo I., Bouwens R., Röttgering H., 2017, MNRAS , 465, 3637 https://doi.org/10.1093/mnras/stw2973 CrossRef Search ADS   Meiksin A. A., 2009, Rev. Modern Phys. , 81, 1405 https://doi.org/10.1103/RevModPhys.81.1405 CrossRef Search ADS   Meiksin A., Khochfar S., Paardekooper J.-P., Dalla Vecchia C., Kohn S., 2017, MNRAS , 471, 3632 https://doi.org/10.1093/mnras/stx1857 CrossRef Search ADS   Mellema G., Iliev I. T., Pen U.-L., Shapiro P. R., 2006, MNRAS , 372, 679 https://doi.org/10.1111/j.1365-2966.2006.10919.x CrossRef Search ADS   Mesinger A., Furlanetto S., Cen R., 2011, MNRAS , 411, 955 https://doi.org/10.1111/j.1365-2966.2010.17731.x CrossRef Search ADS   Mineo S., Gilfanov M., Sunyaev R., 2012a, MNRAS , 419, 2095 https://doi.org/10.1111/j.1365-2966.2011.19862.x CrossRef Search ADS   Mineo S., Gilfanov M., Sunyaev R., 2012b, MNRAS , 426, 1870 https://doi.org/10.1111/j.1365-2966.2012.21831.x CrossRef Search ADS   Mirabel I. F., Dijkstra M., Laurent P., Loeb A., Pritchard J. R., 2011, A&A , 528, A149 https://doi.org/10.1051/0004-6361/201016357 CrossRef Search ADS   Mitra S., Choudhury T. R., Ferrara A., 2018, MNRAS , 473, 1416 https://doi.org/10.1093/mnras/stx2443 CrossRef Search ADS   Nath B. B., Biermann P. L., 1993, MNRAS , 265, 241 https://doi.org/10.1093/mnras/265.1.241 CrossRef Search ADS   Ocvirk P. et al.  , 2016, MNRAS , 463, 1462 https://doi.org/10.1093/mnras/stw2036 CrossRef Search ADS   Oh S. P., 2001, ApJ , 553, 499 https://doi.org/10.1086/320957 CrossRef Search ADS   Oñorbe J., Hennawi J. F., Lukić Z., Walther M., 2017, ApJ , 847, 63 https://doi.org/10.3847/1538-4357/aa898d CrossRef Search ADS   Onoue M. et al.  , 2017, ApJ , 847, L15 https://doi.org/10.3847/2041-8213/aa8cc6 CrossRef Search ADS   O'Shea B. W., Wise J. H., Xu H., Norman M. L., 2015, ApJ , 807, L12 https://doi.org/10.1088/2041-8205/807/1/L12 CrossRef Search ADS   Ostriker J. P., Gnedin N. Y., 1996, ApJ , 472, L63 https://doi.org/10.1086/310375 CrossRef Search ADS   Ouchi M. et al.  , 2018, PASJ , 70, S13 https://doi.org/10.1093/pasj/psx074 CrossRef Search ADS   Paciga G. et al.  , 2011, MNRAS , 413, 1174 https://doi.org/10.1111/j.1365-2966.2011.18208.x CrossRef Search ADS   Pacucci F., Mesinger A., Mineo S., Ferrara A., 2014, MNRAS , 443, 678 https://doi.org/10.1093/mnras/stu1240 CrossRef Search ADS   Parsa S., Dunlop J. S., McLure R. J., 2018, MNRAS , 474, 2904 https://doi.org/10.1093/mnras/stx2887 CrossRef Search ADS   Parsons A. R. et al.  , 2014, ApJ , 788, 106 https://doi.org/10.1088/0004-637X/788/2/106 CrossRef Search ADS   Patil A. H. et al.  , 2017, ApJ , 838, 65 https://doi.org/10.3847/1538-4357/aa63e7 CrossRef Search ADS   Pawlik A. H., Rahmati A., Schaye J., Jeon M., Vecchia C. D., 2017, MNRAS , 466, 960 https://doi.org/10.1093/mnras/stw2869 CrossRef Search ADS   Penzias A. A., Scott E. H., 1968, ApJ , 153, L7 https://doi.org/10.1086/180209 CrossRef Search ADS   Planck Collaboration XLVII, 2016, A&A , 596, A108 https://doi.org/10.1051/0004-6361/201628897 CrossRef Search ADS   Pober J. C. et al.  , 2015, ApJ , 809, 62 https://doi.org/10.1088/0004-637X/809/1/62 CrossRef Search ADS   Puchwein E., Pfrommer C., Springel V., Broderick A. E., Chang P., 2012, MNRAS , 423, 149 https://doi.org/10.1111/j.1365-2966.2012.20738.x CrossRef Search ADS   Qin Y. et al.  , 2017, MNRAS , 472, 2009 https://doi.org/10.1093/mnras/stx1909 CrossRef Search ADS   Rees M. J., 1998, Proc. Natl. Acad. Sci. USA , 95, 47 CrossRef Search ADS   Rees M. J., Setti G., 1969, A&A , 8, 410 Ricotti M., Ostriker J. P., 2004, MNRAS , 352, 547 https://doi.org/10.1111/j.1365-2966.2004.07942.x CrossRef Search ADS   Robertson B. E., Ellis R. S., Furlanetto S. R., Dunlop J. S., 2015, ApJ , 802, L19 https://doi.org/10.1088/2041-8205/802/2/L19 CrossRef Search ADS   Ross H. E., Dixon K. L., Iliev I. T., Mellema G., 2017, MNRAS , 468, 3785 https://doi.org/10.1093/mnras/stx649 CrossRef Search ADS   Safarzadeh M., Scannapieco E., 2016, ApJ , 832, 4 https://doi.org/10.3847/2041-8205/832/1/L9 CrossRef Search ADS   Santos M. G., Ferramacho L., Silva M. B., Amblard A., Cooray A., 2010, MNRAS , 406, 2421 https://doi.org/10.1111/j.1365-2966.2010.16898.x CrossRef Search ADS   Sazonov S., Khabibullin I., 2017, Astron. Lett. , 43, 211 https://doi.org/10.1134/S1063773717040077 CrossRef Search ADS   Sazonov S., Sunyaev R., 2015, MNRAS , 454, 3464 https://doi.org/10.1093/mnras/stv2255 CrossRef Search ADS   Schaepman M., 2009, in Warner T., Duane Nellis M., Foody G., eds, The SAGE Handbook of Remote Sensing . Sage, London, p. 166 Google Scholar CrossRef Search ADS   Semelin B., Eames E., Bolgar F., Caillat M., 2017, MNRAS , 472, 4508 https://doi.org/10.1093/mnras/stx2274 CrossRef Search ADS   Shakura N. I., Sunyaev R. A., 1973, A&A , 24, 337 Singh S. et al.  , 2017, ApJ , 845, L12 https://doi.org/10.3847/2041-8213/aa831b CrossRef Search ADS   So G. C., Norman M. L., Reynolds D. R., Wise J. H., 2014, ApJ , 789, 149 https://doi.org/10.1088/0004-637X/789/2/149 CrossRef Search ADS   Springel V., 2005, MNRAS , 364, 1105 https://doi.org/10.1111/j.1365-2966.2005.09655.x CrossRef Search ADS   Springel V., Hernquist L., 2003, MNRAS , 339, 289 https://doi.org/10.1046/j.1365-8711.2003.06206.x CrossRef Search ADS   Springel V., Di Matteo T., Hernquist L., 2005, MNRAS , 361, 776 https://doi.org/10.1111/j.1365-2966.2005.09238.x CrossRef Search ADS   Strickland D. K., Stevens I. R., 2000, MNRAS , 314, 511 https://doi.org/10.1046/j.1365-8711.2000.03391.x CrossRef Search ADS   Suchkov A. A., Balsara D. S., Heckman T. M., Leitherner C., 1994, ApJ , 430, 511 https://doi.org/10.1086/174427 CrossRef Search ADS   Tange O., 2011, The USENIX Magazine , 36, 42 Van Der Walt S., Colbert S. C., Varoquaux G., 2011, Comput. Sci. Eng. , 13, 22 https://doi.org/10.1109/MCSE.2011.37 CrossRef Search ADS   Vanzella E. et al.  , 2016, ApJ , 825, 41 https://doi.org/10.3847/0004-637X/825/1/41 CrossRef Search ADS   Vanzella E. et al.  , 2018, MNRAS , in press Worseck G., Prochaska J. X., Hennawi J. F., McQuinn M., 2016, ApJ , 825, 144 https://doi.org/10.3847/0004-637X/825/2/144 CrossRef Search ADS   Yoshida N., Oh S. P., Kitayama T., Hernquist L., 2007, ApJ , 663, 687 https://doi.org/10.1086/518227 CrossRef Search ADS   Zheng Z.-Y. et al.  , 2017, ApJ , 842, L22 https://doi.org/10.3847/2041-8213/aa794f CrossRef Search ADS   APPENDIX A: CONVERGENCE To check the convergence of our results, we run two simulations including all source types that differ only in the number of photon packets emitted per source, i.e. Nγ = 104 and 105. In Fig. A1 we show the fraction of cells at z = 12 with a given T, $$x_{\rm H\,\small {II}}$$, $$x_{\rm He\,\small {II}}$$, and $$x_{\rm He\,\small {iii}}$$, together with the relative difference between the results of the two simulations, i.e. D/Dref − 1, where D = T, $$x_{\rm H\,\small {II}}$$, $$x_{\rm He\,\small {II}}$$, and $$x_{\rm He\,\small {iii}}$$ and Dref represents the values obtained with the highest value of Nγ. More than 50 per cent of the cells in the Nγ = 104 simulation are very well converged, with differences at the per cent level, except for $$x_{\rm H\,\small {II}} = x_{\rm He\,\small {II}} \sim 10^{-3}$$ where the results deviate up to ∼ 20 per cent. For the temperatures, the differences are less than 5 per cent for the 67 per cent of the cells that have 10 < T ≤ 104 K. Large differences are observed only in a handful of cells at high temperature or very low ionization fraction $$x_{\rm H\,\small {II}} = x_{\rm He\,\small {II}} < 10^{-5}$$, which we set as the lower limit in this work. Figure A1. View largeDownload slide Upper panels: from left to right, fraction of cells at z = 12 with a given $$x_{\rm H\,\small {II}}$$, $$x_{\rm He\,\small {II}}$$, $$x_{\rm He\,\small {iii}}$$, and T for simulations with Nγ = 104 (red dashed lines) and 105 (blue solid lines). Lower panels: relative difference in per cent with respect to the case with Nγ = 105. Figure A1. View largeDownload slide Upper panels: from left to right, fraction of cells at z = 12 with a given $$x_{\rm H\,\small {II}}$$, $$x_{\rm He\,\small {II}}$$, $$x_{\rm He\,\small {iii}}$$, and T for simulations with Nγ = 104 (red dashed lines) and 105 (blue solid lines). Lower panels: relative difference in per cent with respect to the case with Nγ = 105. APPENDIX B: LIGHT-CONES We construct light-cones inspired by the methods applied by Mellema et al. (2006), Datta et al. (2012), and Giri et al. (2018). Our approach though needs to take into account the non-periodic boundary conditions of our simulations. We create a path Δl(ti) ≡ Δli corresponding to the light travel distance cΔti = kiΔxi covered by a ray travelling through ki cells of the volume, each with sides of physical length,   $$\Delta x_i = \frac{1}{1+z_i} \frac{100\,h^{-1} }{256} \,{\rm pMpc},$$ (B1)from one crash output i to the next, i + 1, having a difference in cosmological age Δti. We hence assume the time steps to be sufficiently small to neglect the evolution between them when calculating the path Δli. We linearly interpolate between each output, i.e. the ki cells that are covered between two crash snapshots have contributions from both snapshots. The position along the ki cells that make up the path between two snapshots can be written as lp, where p = 0, …, ki − 1. For each point lp, which corresponds to a position in both physical (xp, yp) and redshift (zp) space, we also have a third spatial dimension, $$\boldsymbol {c}_p(l_p)$$, which is plotted as the y-axis in the light-cones. It has contributions from the two snapshots at times ti and ti + 1, weighed by the position (ki − 1 − p)/(ki − 1) and p/(ki − 1), respectively, along the path between them. For the ionization fraction we have   $$\boldsymbol {c}_p(l_p) = \frac{1}{k_i - 1}\left[\left(k_i - 1 -p\right) \boldsymbol {x}_{\rm H\,\small {II}}^{t_i}(l_p) + p \boldsymbol {x}_{\rm H\,\small {II}}^{t_{i+1}}(l_p) \right],$$ (B2)where we write the ionization fractions $$\boldsymbol {x}_{\rm H\,\small {II}}^{t_i}(l_p)$$ from the two contributing snapshots as vectors to indicate that we only use the cells along the third spatial dimension at this point in (zp, xp, yp) space. There are various possible approaches for deciding the paths Δli. Instead of choosing random paths, which could overlap, we instead try fixed patterns, and attempt at covering as much of the simulations volumes as possible in the light-cones. For this purpose, we choose a continuous whisk broom scanning approach (see e.g. Schaepman 2009) with decreasing horizontal scan width and a several pixels wide (slightly decreasing) vertical scan offset between each horizontal scan. Both are reset to maximum widths once we reach the vertical boundary. This way we prevent, as far as possible, the same objects at the same cosmological ages to appear in the light-cones. We only use the interior 236 cells (out of a maximum of 256) along each scanning axis to prevent boundary effects on the light-cones. The parameters of this approach (the number of pixels to decrease the scan width and offset, vertical scan offset) must however be tuned to obtain a satisfying result. The approach is thus optimal for a visualization, rather than for an objective quantification, of the temporal progression. © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Monthly Notices of the Royal Astronomical Society Oxford University Press

# The epoch of cosmic heating by early sources of X-rays

, Volume 476 (1) – May 1, 2018
17 pages

/lp/ou_press/the-epoch-of-cosmic-heating-by-early-sources-of-x-rays-DHBPLt9sGb
Publisher
Oxford University Press
ISSN
0035-8711
eISSN
1365-2966
D.O.I.
10.1093/mnras/sty272
Publisher site
See Article on Publisher Site

### Abstract

Abstract Observations of the 21 cm line from neutral hydrogen indicate that an epoch of heating (EoH) might have preceded the later epoch of reionization. Here we study the effects on the ionization state and the thermal history of the intergalactic medium (IGM) during the EoH induced by different assumptions on ionizing sources in the high-redshift Universe: (i) stars; (ii) X-ray binaries (XRBs); (iii) thermal bremsstrahlung of the hot interstellar medium (ISM); and (iv) accreting nuclear black holes (BHs). To this aim, we post-process outputs from the (100 h−1 comoving Mpc)3 hydrodynamical simulation MassiveBlack-II with the cosmological 3D radiative transfer code crash, which follows the propagation of ultraviolet and X-ray photons, computing the thermal and ionization state of hydrogen and helium through the EoH. We find that stars determine the fully ionized morphology of the IGM, while the spectrally hard XRBs pave way for efficient subsequent heating and ionization by the spectrally softer ISM. With the seeding prescription in MassiveBlack-II, BHs do not contribute significantly to either ionization or heating. With only stars, most of the IGM remains in a cold state (with a median T = 11 K at z = 10), however, the presence of more energetic sources raises the temperature of regions around the brightest and more clustered sources above that of the cosmic microwave background, opening the possibility to observing the 21 cm signal in emission. radiative transfer, dark ages, reionization, first stars 1 INTRODUCTION A central, outstanding question in cosmology is the evolution and nature of the hydrogen and helium constituting the intergalactic medium (IGM), particularly at high redshift. Hydrogen was expected to reveal itself either as a deep trough blueward of the H i 21 cm line (Davies & Jennison 1964) or the H i Lyα line (Gunn & Peterson 1965) in the spectra of distant, bright radio galaxies (Field 1962; Penzias & Scott 1968; Allen 1969). However, this did not happen, leading to the conclusion that the IGM must presently be highly ionized and hot. Not until recently was the H i trough finally detected, but at higher redshifts (e.g. Fan et al. 2000, 2006; McGreer, Mesinger & D'Odorico 2015). The IGM must once have been neutral, and then have undergone a phase transformation, the Cosmic Dawn, into its present-day state. Persisting questions are how this transformation occurred, and what drove it. Answering them will shed light on the physical properties of the majority of the baryonic Universe, and open doors to understanding the interplay between the IGM and the large- and small-scale structures that stem from it. The IGM is thus a highly uncharted terrain (see e.g. Meiksin 2009 for a review). In this work, we investigate the properties and signatures of the sources of ionizing, and possibly heating, radiation that flared up as structures formed, marking the transition from the Dark Ages to the Cosmic Dawn (Rees 1998). Early cosmological 21 cm surveys have begun to report that heating may precede full ionization, resulting in a separate epoch of heating (EoH) preceding the epoch of reionization (EoR). The recent constraint on the global signal with SARAS 2 (Singh et al. 2017) disfavours late heating, disentangling the EoH from the EoR (Fialkov, Barkana & Visbal 2014b). These observations also align themselves with other surveys, as the early results of Bowman, Rogers & Hewitt (2008), that were followed by upper limits on the scale and magnitude of brightness temperature fluctuations [e.g. from Giant Metrewave Radio Telescope (GMRT), Paciga et al. 2011; Murchison Widefield Array (MWA), Dillon et al. 2014; Ewall-Wice et al. 2016; Precision Array for Probing the Epoch of Reionization (PAPER), Parsons et al. 2014; Ali et al. 2015; and Low-Frequency Array (LOFAR), Patil et al. 2017]. However, these surveys are impeded at very high redshift by the ionospheric contamination from Earth, which is nearly fully opaque to redshifted 21 cm signal. To mitigate for this, a seminal Chinese–Dutch mission scheduled for 2018 plans to obtain the global 21 cm signal from the far side of the Moon.1 Besides 21 cm tomography, which is still in its infancy after its revival as a probe by Madau, Meiksin & Rees (1997), there are other observational means for examining the ionization and thermal state of the IGM. The spectra of bright, single sources (as quasars or gamma-ray bursts; see e.g. review by Ciardi & Ferrara 2005) reveal the properties of the intervening IGM along the line of sight. The increasing abundance of neutral H i at earlier z has been indicated from large-scale surveys, as well (e.g. Matthee et al. 2015; Zheng et al. 2017; Ouchi et al. 2018 or see review by Dijkstra 2014). Hydrogen ionization is accompanied by production of free electrons. Satellite-based experiments (Komatsu et al. 2011; Planck Collaboration XLVII 2016) have provided tight constraints on the scattering of the cosmic microwave background (CMB) off these electrons, as the CMB becomes dampened and linearly polarized. These observations are however in slight tension, as the continuing lowering of the Thomson scattering optical depth requires a later completion of cosmic reionization. We do not expect a single source type of heating and ionization to be able to account for the diversity in the observed signatures. Kakiichi et al. (2017) instead found the interplay between different ionizing sources to be essential in establishing the thermal and ionization state of the IGM and to leave different signatures. Possible viable candidates include the following. Stars residing within galaxies. These are thought to be the primary drivers of reionization, as argued by Madau, Haardt & Rees (1999). The recent low Thomson scattering optical depth reported by the Planck Collaboration XLVII (2016), combined with constraints from Hubble (Robertson et al. 2015), has rejuvenated the interest in their effect on reionization. Structure formation is also affected by the thermal and ionizing feedback from stars, as discussed by e.g. Couchman & Rees (1986), Cen (1992), and Fukugita & Kawasaki (1994). Depending on their initial mass function, mass and metallicity, they will produce copious amounts of ionizing photons, as shown in the review by Ciardi & Ferrara (2005). Accreting nuclear black holes (or simply black holes, BHs). Hereafter this term will refer to active galactic nuclei (AGNs) and the brighter quasi-stellar objects (QSOs). While these were early candidates for IGM ionization (Rees & Setti 1969; Arons & McCray 1970; Bergeron & Salpeter 1970), assessing their importance is strongly hampered by lacking data on the evolution of the faint end of the QSO luminosity function (QLF) at higher redshifts. With the detection of 22 faint AGNs by Giallongo et al. (2015), the possibility of them producing amounts of ionizing photons sufficient to drive reionization independently from stars was once again examined. Madau & Haardt (2015) and Mitra, Choudhury & Ferrara (2018) found this to be possible, but it would also doubly ionize helium before redshift 4, in contrast to recent observations of an extended He ii reionization at z ∼ 2.8 (Worseck et al. 2016). Qin et al. (2017) report that the models of Madau & Haardt (2015) and Mitra et al. (2018) overpredict the contribution from AGN, and do not find them to dominate over stars, a conclusion also shared by Hassan et al. (2018) and Oñorbe et al. (2017). Updated QLF by Parsa, Dunlop & McLure (2018) or Onoue et al. (2017) also indicates the contribution of BHs to be insufficient to drive reionization. However, the contribution from AGNs was found by Chardin, Puchwein & Haehnelt (2017) and D'Aloisio et al. (2017) to be vital to the large-scale opacity fluctuations observed in the Lyα forest, which was strengthened by the recent observations of a through spanning 240 h−1 comoving Mpc (cMpc) by Barnett et al. (2017). Galactic X-ray binary (XRB) systems comprising a neutron star or a BH devouring a companion star. Among such systems, the majority of the ionizing luminosity at high-z originates from massive [high-mass X-ray binaries (HMXBs)] rather than low-mass [low-mass X-ray binaries (LMXBs)] binary systems (Mirabel et al. 2011; Fragos et al. 2013a,b; Madau & Fragos 2017; Sazonov & Khabibullin 2017). Mineo, Gilfanov & Sunyaev (2012a,b) found the spectra of XRBs to be too hard to account for the soft X-ray flux of galaxies, while they become dominant at higher energies. Thermal bremsstrahlung in the diffuse gas of the heated interstellar medium (ISM). The heating mechanism of the diffuse gas has been a topic under investigation for decades (see e.g. the review by Fabbiano 1989). Interactions between supernova-driven galactic superwinds and clouds (e.g. Chevalier & Clegg 1985) is the preferred explanation (Pacucci et al. 2014). Shocked gas in the galactic halo and disc (e.g. Suchkov et al. 1994), and hot galactic winds (e.g. Strickland & Stevens 2000) are also processes that could yield predominantly soft X-rays. Other candidates, which will not be examined further in this paper because of their secondary role, are e.g. low-energy cosmic rays (Ginzburg & Ozernoi 1966; Nath & Biermann 1993; Sazonov & Sunyaev 2015; Leite et al. 2017), self-annihilation or decay of dark matter (e.g. Liu, Slatyer & Zavala 2016), and plasma beam instabilities in TeV blazars (e.g. Chang, Broderick & Pfrommer 2012; Puchwein et al. 2012). There are different approaches for simulating the physics of the Cosmic Dawn. Ideally, structure formation can be coupled to radiative transfer (RT; e.g. Gnedin 2014; So et al. 2014; O'Shea et al. 2015; Ocvirk et al. 2016; Pawlik et al. 2017; Semelin et al. 2017). However, this is computationally expensive and requires simplifications. For example helium physics is often not treated, or the number of frequency bins used may be limited to a few, albeit the cross-section of hydrogen varies approximately as ν−3. The spatial extent is also often well below the required ∼100 h−1 cMpc needed to provide a consistent reionization history (Iliev et al. 2014). Seminumerical codes as 21cmfast (Mesinger, Furlanetto & Cen 2011) or SimFast21 (Santos et al. 2010; Hassan et al. 2016) can be applied on large scales using an excursion-set formalism (Furlanetto, Zaldarriaga & Hernquist 2004). This allows for fast parameter-space exploration, but may lack the treatment of some important physical processes, such as helium reionization, temperature evolution, partial ionization, and X-ray implementation. Nevertheless, this approach has provided to be rewarding in studies of possible consequences of radiative feedback on structure formation and its 21 cm signatures (see e.g. Fialkov et al. 2014a). Radiative post-processing of hydrodynamical or N-body simulations is also possible (e.g. Baek et al. 2010; Ciardi et al. 2012; Graziani et al. 2015; Ross et al. 2017) and allows for the dedication of computational power to e.g. treat heating, multifrequency RT, as well as the effects of helium. Here we examine the effects of stars, BHs, XRBs, and the ISM on both the thermal and ionization states of the IGM during the EoH. Our approach is the following: we rely on the hydrodynamical cosmological structure formation simulation MassiveBlack-II (MBII, described in Section 2.1; Khandai et al. 2015), which includes baryonic physics and feedback processes, to provide us with the physical environment of the IGM (temperature and gas density), as well as the location and properties of the sources. We assign spectra and ionizing luminosity to star forming particles based on their star formation rates (SFRs), masses, ages, and metallicities, and to BH particles based on their accretion rates. We then perform radiative post-processing of the MBII simulations with the cosmological RT code crash (see Section 2.2) that gives us the evolution of the thermal and ionization state of the IGM. The properties of the ionizing sources are based on empirical relations introduced in Section 2.3. We present our findings in Section 3 and relate them to comparable studies also summarizing our conclusions in Sections 4 and 5. 2 METHODOLOGY Here we describe how we combine the outputs of the cosmological hydrodynamical simulation MBII (Section 2.1) with population synthesis modelling of ionizing sources (Section 2.3), and finally perform multifrequency RT in the same cosmological volume with crash (Section 2.2). 2.1 Cosmological hydrodynamic simulation MBII (Khandai et al. 2015, hereafter K15) is a high-resolution cosmological smoothed particle hydrodynamics (SPH) simulation tracking stellar populations, galaxies, accreting, and dormant BHs, and their properties (as position, age, metallicity, mass, accretion rate, and SFR). The simulation has been run using p-gadget, a newer version of gadget-3 (see Springel 2005, for an earlier version). It accounts for baryonic physics and feedback effects of the sources on their environment following earlier works in its approach to feedback (Di Matteo et al. 2008, 2012; Croft et al. 2009; Degraf, Di Matteo & Springel 2010), subgrid treatment of star formation (Springel & Hernquist 2003), and seeding and evolution of BHs (Di Matteo, Springel & Hernquist 2005; Springel, Di Matteo & Hernquist 2005). The simulation has a box length of 100 h−1 cMpc and is performed in the 7-year Wilkinson Microwave Anisotropy Probe (WMAP7) Λ cold dark matter (ΛCDM) cosmology (Komatsu et al. 2011),2 using 2 × 17923 particles of mass mDM = 1.1 × 107 h−1 M⊙ and mgas = 2.2 × 106 h−1 M⊙ for dark matter and gas, respectively. The adopted gravitational softening length is 1.85 h−1 ckpc. As a reference, we have a total of 28 and 169 520 haloes hosting at least a star particle at z = 18 and 10, respectively, where the lowest dark matter halo masses are 2 × 108 and 9 × 107 h−1 M⊙. The highest dark matter halo masses are 2 × 109 and 1 × 1011 h−1 M⊙ at z = 18 and 10, respectively. The BHs form from seeds of 5 × 105 h−1 M⊙, growing with accretion rates in the range (106–107) h−1 M⊙ Gyr−1 to a maximum mass of 1.4 × 106 h−1 M⊙ at z = 10. The first seed BH is found at z = 13, while we have 17 BHs at z = 10. We refer the reader to K15 for more details on the simulation. We employ six snapshots from MBII, covering the evolution between redshifts z = 18 and 10, each describing the instantaneous state of the simulation. The particle distribution of each snapshot is mapped on to a Cartesian grid of $$N_{\rm c}^3$$ cells to create maps of gas number density and temperature, as well as location and properties of the ionizing sources (see Section 2.3 for a detailed description of how we convert MBII data to ionizing sources). The reference value is Nc = 256, corresponding to a spatial resolution of 391 h−1 ckpc. After gridding, the 45 (292 685) star particles present at z = 18 (10) are reduced to 26 (56 702). These are effectively our sources. When more than one particle ends up in a cell, their properties are summed up (see following section). The gas density is converted to hydrogen and helium number densities by assuming a number fraction X = 0.92 and Y = 0.08, respectively, and no metals. 2.2 Cosmological radiative transfer simulation The RT of ionizing photons is performed by post-processing the outputs of MBII with the multifrequency Monte Carlo ray-tracing code crash (Ciardi et al. 2001; Maselli, Ferrara & Ciardi 2003; Maselli, Ciardi & Kanekar 2009; Graziani, Maselli & Ciardi 2013), which calculates the ionization state of hydrogen and helium, as well as the gas temperature in each grid cell traversed by photons emitted by a radiation source. The version of crash employed here features a self-consistent treatment of ultraviolet (UV) and soft X-ray photons, in which X-ray ionization and heating and detailed secondary electron physics are included (Graziani et al., in preparation). We refer the reader to the original papers for more details on crash. The spectral shape of the sources is discretized into 82 frequency bins, where 72 are in the UV regime below hPν = 200 eV, and the rest are in the soft X-ray regime (0.2–2 keV). The bins are spaced more densely around the ionization thresholds of hydrogen (13.6 eV) and helium (24.6 and 54.4 eV). The radiation emitted by each source i at a redshift z is discretized into Nγ photon packets. Each packet, Nph,i(ν, z) (in units photons Hz−1), holds the total number of photons emitted by source i during the time step Δtem(z) (in s) in the different frequency bins:   $$N_{ {\rm ph}, i}(\nu ,z) = \frac{S_i(\nu ,z)}{L_i(z)}f_{\rm esc}(\nu )\varepsilon _i(z) \Delta t_{\rm em}(z),$$ (1)where fesc(ν) is the frequency-dependent escape fraction, εi(z) is the rate of ionizing photons the source emits (in photons s−1), $$S_i(\nu ,z)/L_i(z)=\hat{S}_i(\nu ,z)$$ is the normalized spectral shape (in Hz−1), where Si(ν, z) defines the spectral energy distribution (SED; in units of erg Hz−1 s−1) and Li(z) the luminosity (in units of erg s−1). We adopt an escape fraction3 of UV photons fesc = 15 per cent for all sources except BHs, for which photons at all frequencies are assumed to escape. Although the adoption of a single value for the escape fraction is an oversimplification as in reality it depends e.g. on the mass of dust and gas and redshift of the host galaxy, its density distribution, the mass and location of the stellar sources (see e.g. Ciardi & Ferrara 2005 and its updated version on arXiv), fesc is to a large degree an unconstrained parameter at low (see e.g. Vanzella et al. 2016; Matthee et al. 2017) and high redshift, where physical conditions could have promoted large values of fesc (see e.g. Kitayama et al. 2004; Yoshida et al. 2007; Safarzadeh & Scannapieco 2016 or the recent z = 4 observations by Vanzella et al. 2018). In a companion paper, focused on the EoR, we will address the effect of different choices of the escape fraction on the global emissivity and reionization process. It should be noted that cells hosting sources are treated as any other cell, i.e. their ionization and temperature evolution is dictated by the properties of the crossing rays. The rate of ionizing photons, or emissivity, can be written as   $$\varepsilon _i(z) = \int _{\rm 13.6\,eV}^{\rm 2\,keV} \frac{S_i(\nu ,z)}{h_{\rm P}^2 \nu } \,{\rm d}(h_{\rm P}\nu ).$$ (2)For simplicity, from now on, we will omit the explicit dependences on ν and z. Each of these packets is followed as it traverses the simulation volume, the photons in different frequency bins ionize, and heat the intervening hydrogen and helium, and finally it either goes extinct or reaches the cosmic volume boundary.4 Because of the large simulated volume, photon redshift is also accounted for. Our reference cases use Nγ = 105 (see Appendix A for a convergence analysis). In the following section we will describe how the spectrum, luminosity, and emissivity of the sources are modelled. 2.3 Sources of ionizing radiation In our simulations we consider four types of ionizing sources: regular stars, hereafter abbreviated only as stars; neutron star/black hole X-ray binaries, hereafter XRBs; thermal bremsstrahlung from supernova-heated ISM, hereafter ISM; accreting nuclear black holes, hereafter BHs. For each source i of any of the types s we model its SED, $$S_i^s$$, and luminosity, $$L_i^s$$, as detailed in the following subsections. When several sources of different types are present in the same cell, we sum up their contributions (except for BHs, which are treated separately) to obtain $$S_i=\sum _s S^s_i$$ and $$L_i=\sum _s L^s_i$$. Using equations (2) and (1), we can then evaluate the photon content of each packet emitted. In Fig. 1 we plot the evolution of the emissivities. In the upper panel we plot the evolution of the median, the 25th–75th percentiles, and the minimum/maximum values of the emissivity for the various source types, as evaluated from equation (2). In the lower panel, we plot the comoving volume averaged emissivities. These reflect the abundance of the sources, unlike the individual emissivities plotted in the upper panel, showing that stars provide the bulk of the ionizing photons throughout the redshift range considered here. Even though individual BHs have higher emissivities than the majority of the stars in galaxies, they are much fewer in number, and are therefore not dominating the overall emissivity budget. The median values of the emissivities for the different source types remain also fairly constant. This applies to the maximum values too: the brightest stellar-type sources at z = 14 are less than an order of magnitude brighter at z = 10. The brightest BHs, on the other hand, do not have higher emissivities compared to the brightest galaxies, which can emit an order of magnitude more ionizing photons. The slow evolution in emissivities is related to the physical properties (e.g. stellar masses, SFRs, and metallicities) that go into determining the SEDs. We note that the median stellar emissivity decreases with decreasing redshift. Although the median stellar mass increases, suggesting higher luminosities, this is counteracted by an increment of both the median stellar mass weighted metallicity and age of the stars. Nevertheless, stars are the dominant producer of ionizing photons at all redshifts. The evolution of the emissivities of the XRBs and the ISM is dictated by the slow evolution of the SFRs. The sudden drop in the lower limit of the emissivity of XRBs is due to galaxies that only hosts the fainter LMXBs. The first BH arises at z = 13, and not until z ≤ 11 do we have more than one present in our volume. Figure 1. View largeDownload slide Upper panel: redshift evolution of the emissivity $$\varepsilon _i^{s}$$ of individual sources (upper panel), grouped per source type s. Stars are plotted in black (single dot), BHs in blue (two dots), XRBs in purple (three dot), and the ISM in red (four dots). The solid, central lines are the median values, enclosed in shaded regions that show the 25th–75th percentile of the values. The span between the minimum and maximum emissivities are given from the dashed lines, with ‘+’ and ‘o’ denoting upper and lower limits, respectively. Note that these values are intrinsic and have not been scaled by an escape fraction. Lower panel: redshift evolution of the volume averaged emissivity εs/V per source type s. Here, the escape fractions are accounted for. Figure 1. View largeDownload slide Upper panel: redshift evolution of the emissivity $$\varepsilon _i^{s}$$ of individual sources (upper panel), grouped per source type s. Stars are plotted in black (single dot), BHs in blue (two dots), XRBs in purple (three dot), and the ISM in red (four dots). The solid, central lines are the median values, enclosed in shaded regions that show the 25th–75th percentile of the values. The span between the minimum and maximum emissivities are given from the dashed lines, with ‘+’ and ‘o’ denoting upper and lower limits, respectively. Note that these values are intrinsic and have not been scaled by an escape fraction. Lower panel: redshift evolution of the volume averaged emissivity εs/V per source type s. Here, the escape fractions are accounted for. While crash can handle a different spectrum for each single source, for the sake of simplicity we adopt at each redshift z an average spectral shape for all sources but BHs. More specifically, for each source i we evaluate $$S_i=S^{\rm stars}_i+S^{\rm XRB}_i+S^{\rm ISM}_i$$. We then use the average $$\bar{S}=\langle S_i\rangle$$ in equation (1) rather than Si. Whenever a BH is present in a cell, this is added as a separate source having the same spatial coordinates. Its spectrum is similarly calculated as $$\bar{S}^{\rm BH}=\langle S^{\rm BH}_i\rangle$$. In Fig. 2 we plot $$\bar{S}$$ at different redshifts z. We see that stars dominate at energies hPν ≲ 60 eV, while the ISM contribution is relevant above the He ii ionization threshold, i.e. into hard UV and the soft X-rays. The XRBs provide the harder X-rays. The weak redshift evolution is due to a combination of the averaging effect effectively preferring brighter sources, and the mentioned slow evolution in the underlying physical properties determining the spectra. Figure 2. View largeDownload slide Globally averaged spectral energy distributions (SEDs) $$\bar{S}(\nu ,z)$$ of the different source types: stars (short dashed lines), XRBs (long dashed), ISM (short–long dashed), total galactic (solid), and BHs (long/double short dashed). The faint, vertical grey lines indicate the ionization thresholds for hydrogen (13.6 eV), neutral helium (24.6 eV), and singly ionized helium (54.4 eV). Different colours indicate different redshifts. Figure 2. View largeDownload slide Globally averaged spectral energy distributions (SEDs) $$\bar{S}(\nu ,z)$$ of the different source types: stars (short dashed lines), XRBs (long dashed), ISM (short–long dashed), total galactic (solid), and BHs (long/double short dashed). The faint, vertical grey lines indicate the ionization thresholds for hydrogen (13.6 eV), neutral helium (24.6 eV), and singly ionized helium (54.4 eV). Different colours indicate different redshifts. In the following sections we will describe in more detail how we evaluate the luminosity $$L^s_i$$ and normalized spectral shape $$\hat{S}^s_i$$ for the various source types. 2.3.1 Stars We model the ionizing radiation from stars by using stellar particles identified in MBII. From the age, metallicity, and mass of each stellar particle p, we obtain $$\hat{S}^{\rm stars}_p$$ and $$L^{\rm stars}_p$$ using the stellar population synthesis code bpass (Eldridge & Stanway 2012). We adopt the instantaneous starburst prescription of star formation and the evolution model that does not account for interactions in binary systems.5 As a single halo i may comprise several particles p, in this case we sum up the contributions from all the particles to obtain $$L_i^{\rm stars} = \sum _{p \in i} L_p^{\rm stars}$$ and $$\hat{S}_i^{\rm stars} = \sum _{p \in i} \hat{S}_p^{\rm stars}$$. 2.3.2 X-ray binaries To account for ionizing radiation coming from XRB systems, we combine the galactic properties provided by MBII with scaling relations from the Fragos et al. (2013a,b) XRB population synthesis model, recently updated by Madau & Fragos (2017). We can thus capture both the metallicity evolution of the HMXBs and the age and stellar mass dependence of the LMXBs, as well as the spectral shape evolution in redshift. For each stellar particle p, we obtain its spectral shape $$\hat{S}^{\rm XRB}_p$$ from Fragos et al. (2013a,b) and its luminosity $$L^{\rm XRB}_p$$ as   $$L^{\rm XRB}_{p} = L^{\rm HMXB}_{p} + L^{\rm LMXB}_{p}.$$ (3)The contribution from HMXBs can be found using equation (3) in Fragos et al. (2013b) with updated coefficients β0–4 from Madau & Fragos (2017):   \begin{eqnarray} \log \left(L_{p}^{\rm HMXB}/{\rm SFR}_i\right) &=& \beta _0 \!+\! \beta _1 Z_p + \beta _2 Z_p^2 + \beta _3 Z_p^3 + \beta _4 Z_p^4 \nonumber \\ &&\times \,({\rm erg s^{-1}} {\rm M}_{\odot }^{-1} \rm yr) \quad {\rm for}\,\, Z \in [ 0,0.025], \end{eqnarray} (4)where SFRi (in M⊙ yr−1) is the star formation rate of the halo i hosting particle p, and Zp is the metallicity of the stellar particle. The contribution to the luminosity from the LMXBs from each stellar particle p is found using equation (4) in Fragos et al. (2013a), also with updated coefficients γ(0–4) from Madau & Fragos (2017):   \begin{eqnarray} \log \left(L_{p}^{\rm LMXB}/M_p\right) &=& \gamma _0 + \gamma _1 \log (t_p/{\rm Gyr}) + \gamma _2 \log (t_p/{\rm Gyr})^2 \nonumber \\ &=& + \,\gamma _3 \log (t_p/{\rm Gyr})^3 + \gamma _4 \log (t_p/{\rm Gyr})^4 \nonumber \\ && \times \,({\rm erg\, s^{-1}} 10^{10} {\rm M}_{\odot }) \quad {\rm for}\,\, t_p \in [0,13.7] {\rm Gyr}, \nonumber\\ \end{eqnarray} (5)where tp (in Gyr) and Mp (in 1010 M⊙) are the age and mass of the stellar particle, respectively. The particles luminosities and spectral shapes are added as described in the previous section to obtain the corresponding source characteristics. Finally note that we adopt coefficients and spectral shapes that are not attenuated from interstellar absorption. 2.3.3 Thermal X-rays from hot ISM We also include ionizing radiation from the diffuse ISM of galaxies. The spectral shape $$\hat{S}_i^{\rm ISM} (\nu )$$ is assumed to be that of thermal bremsstrahlung and constant in redshift (Pacucci et al. 2014):   $$\hat{S}_i^{\rm ISM} (\nu ) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}C & {\rm for }\,\, h_{\rm P} \nu \le kT_{\rm ISM}, \\ C \left(h_{\rm P} \nu /kT_{\rm ISM}\right)^{-3} & {\rm for }\, \,h_{\rm P} \nu > kT_{\rm ISM}, \end{array}\right.$$ (6)where C is a normalization constant that ensures the correct units of Hz−1, and kTISM is the thermal energy of the ISM in eV. From the spectral analysis of observations of the diffuse gas in galaxies by Mineo et al. (2012b), we use kTISM = 240 eV, translating into a characteristic temperature of the heated ISM of ∼106 K. Each halo i has an individual luminosity $$L^{\rm ISM}_i$$, which we evaluate using equation (3) in Mineo et al. (2012b):   \begin{eqnarray} L_{i}^{\rm ISM} ({0.3\hbox{--}10\rm \, keV}) / {\rm SFR}_i &=& \left( 7.3 \pm 1.3 \right) \times 10^{39} \nonumber \\ &&\times\, ( {\rm erg\, s^{-1}} \,{\rm M}_{\odot }^{-1}\,\rm yr), \end{eqnarray} (7)where 0.3–10 keV indicates the photon energy range this relation was obtained for. Note that $$L^{\rm ISM}_i$$ is corrected to be free of interstellar attenuation. We also rescale $$L^{\rm ISM}_i$$ to match our frequency band, which has the lower and upper limits of 13.6 eV and 2 keV, respectively. 2.3.4 Accreting nuclear black holes To account for the ionizing photons originated from accretion discs surrounding nuclear BHs, we identify BH particles in MBII and use their accretion rates to determine the production of ionizing photons. The bolometric luminosity of a BH i is (Shakura & Sunyaev 1973)   $$L_{i}^{\rm BH} = \eta \dot{M}_{i} c^2,$$ (8)where η is an efficiency parameter, $$\dot{M}_{i}$$ is the accretion rate, and c is the speed of light. Consistent with the BH evolution and feedback in MBII, we choose η = 0.1. As a spectral shape, we adopt the observationally derived mean QSO SED of Krawczyk et al. (2013), which is based on 108 104 QSOs sampled at 0.064 < z < 5.46. When no observational data are available between 13.6 and 200 eV, this is derived as interpolation between the mean SEDs for which they have sufficient observations at both higher and lower energies. For energies greater than 200 eV, the spectral shape is modelled as a power law,   $$\hat{S}^{\rm BH}_i \left( h_{\rm P} \nu > 200 \, {\rm eV} \right) \propto \nu ^{-1}.$$ (9)Finally note that no evolution of the SED with redshift is assumed. 3 RESULTS A central question of this study is: can any of our source types drive an EoH separate from an EoR? In other words, is there heating preceding substantial ionization? To examine this, we first discuss the qualitative and quantitative (in terms of average physical quantities) impact of the different source types on the IGM and its ionization and temperature histories (Section 3.1). Second, we focus on the thermal and ionization state of the IGM at z = 12 and 10, to discuss the details of these processes by using global maps, difference maps, and phase diagrams (Section 3.2). For the sake of clarity, throughout this paper we will use the following nomenclature to define the various ionization states of the IGM: neutral when $$x_{\rm H\,\small {II}} < 10^{-5}$$, partially ionized when $$10^{-5} \le x_{\rm H\,\small {II}} < 10^{-1}$$, highly ionized when $$10^{-1} \le x_{\rm H\,\small {II}} < 0.9$$, and fully ionized when $$x_{\rm H\,\small {II}} \ge 0.9$$. 3.1 Global history 3.1.1 Average ionization fractions and temperatures In this section we present both mass and volume average values of various physical quantities, which are summarized in Table 1. Table 1. Thermal and ionization state of the IGM at z = 14, 12, and 10 for different combinations of source types. Note that ionization fractions below 10−5 are denoted with ‘<’. Source type  Ta (K)  $$x_{\rm H\,\small {II}}$$  $$x_{\rm He\,\small {iii}}$$    Median  Volume  Mass  Neutral  Median  Volume  Mass  Median  Volume  Mass      avg.  avg.  avg.    avg.  avg.    avg.  avg.  z = 14  Stars  10  16  23  16  <  7.144 × 10−5  1.514 × 10−4  <  <  <  Stars, BHs  10  16  23  16  <  7.146 × 10−5  1.514 × 10−4  <  <  <  Stars, XRBs  10  16  23  16  <  7.165 × 10−5  1.516 × 10−4  <  <  <  Stars, ISM  10  16  23  16  <  7.201 × 10−5  1.522 × 10−4  <  <  <  Stars, BHs, XRBs, ISM  10  16  23  16  <  7.217 × 10−5  1.523 × 10−4  <  <  <  Stars, BHs, XRBs, ISM (X = 1, Y = 0)  10  17  23  16  <  7.273 × 10−5  1.533 × 10−4  –  –  –  Stars, BHs, XRBs, ISM (X = 0.92, Y = 0)  10  17  24  16  <  8.304 × 10−5  1.720 × 10−4  –  –  –  z = 12  Stars  10  73  142  51  <  1.627 × 10−3  3.456 × 10−3  <  <  <  Stars, BHs  10  73  142  51  <  1.630 × 10−3  3.459 × 10−3  <  <  <  Stars, XRBs  10  73  142  51  <  1.633 × 10−3  3.462 × 10−3  <  <  <  Stars, ISM  11  74  143  52  <  1.644 × 10−3  3.478 × 10−3  <  <  <  Stars, BHs, XRBs, ISM  11  74  143  52  1.004 × 10−5  1.652 × 10−3  3.486 × 10−3  <  <  <  Stars, BHs, XRBs, ISM (X = 1, Y = 0)  11  75  142  54  1.381 × 10−5  1.688 × 10−3  3.545 × 10−3  –  –  –  Stars, BHs, XRBs, ISM (X = 0.92, Y = 0)  11  79  149  56  1.524 × 10−5  1.923 × 10−3  3.960 × 10−3  –  –  –  z = 10  Stars  11  496  886  214  <  2.009 × 10−2  3.581 × 10−2  <  <  1.073 × 10−5  Stars, BHs  11  497  887  214  <  2.012 × 10−2  3.584 × 10−2  <  <  1.238 × 10−5  Stars, XRBs  12  498  888  215  4.442 × 10−5  2.014 × 10−2  3.586 × 10−2  <  <  1.311 × 10−5  Stars, ISM  16  504  897  219  1.008 × 10−4  2.028 × 10−2  3.604 × 10−2  <  2.915 × 10−5  8.005 × 10−5  Stars, BHs, XRBs, ISM  18  506  900  221  1.464 × 10−4  2.036 × 10−2  3.612 × 10−2  <  3.092 × 10−5  8.408 × 10−5  Stars, BHs, XRBs, ISM (X = 1, Y = 0)  19  510  876  244  1.927 × 10−4  2.073 × 10−2  3.664 × 10−2  –  –  –  Stars, BHs, XRBs, ISM (X = 0.92, Y = 0)  19  558  941  252  2.111 × 10−4  2.374 × 10−2  4.100 × 10−2  –  –  –  Source type  Ta (K)  $$x_{\rm H\,\small {II}}$$  $$x_{\rm He\,\small {iii}}$$    Median  Volume  Mass  Neutral  Median  Volume  Mass  Median  Volume  Mass      avg.  avg.  avg.    avg.  avg.    avg.  avg.  z = 14  Stars  10  16  23  16  <  7.144 × 10−5  1.514 × 10−4  <  <  <  Stars, BHs  10  16  23  16  <  7.146 × 10−5  1.514 × 10−4  <  <  <  Stars, XRBs  10  16  23  16  <  7.165 × 10−5  1.516 × 10−4  <  <  <  Stars, ISM  10  16  23  16  <  7.201 × 10−5  1.522 × 10−4  <  <  <  Stars, BHs, XRBs, ISM  10  16  23  16  <  7.217 × 10−5  1.523 × 10−4  <  <  <  Stars, BHs, XRBs, ISM (X = 1, Y = 0)  10  17  23  16  <  7.273 × 10−5  1.533 × 10−4  –  –  –  Stars, BHs, XRBs, ISM (X = 0.92, Y = 0)  10  17  24  16  <  8.304 × 10−5  1.720 × 10−4  –  –  –  z = 12  Stars  10  73  142  51  <  1.627 × 10−3  3.456 × 10−3  <  <  <  Stars, BHs  10  73  142  51  <  1.630 × 10−3  3.459 × 10−3  <  <  <  Stars, XRBs  10  73  142  51  <  1.633 × 10−3  3.462 × 10−3  <  <  <  Stars, ISM  11  74  143  52  <  1.644 × 10−3  3.478 × 10−3  <  <  <  Stars, BHs, XRBs, ISM  11  74  143  52  1.004 × 10−5  1.652 × 10−3  3.486 × 10−3  <  <  <  Stars, BHs, XRBs, ISM (X = 1, Y = 0)  11  75  142  54  1.381 × 10−5  1.688 × 10−3  3.545 × 10−3  –  –  –  Stars, BHs, XRBs, ISM (X = 0.92, Y = 0)  11  79  149  56  1.524 × 10−5  1.923 × 10−3  3.960 × 10−3  –  –  –  z = 10  Stars  11  496  886  214  <  2.009 × 10−2  3.581 × 10−2  <  <  1.073 × 10−5  Stars, BHs  11  497  887  214  <  2.012 × 10−2  3.584 × 10−2  <  <  1.238 × 10−5  Stars, XRBs  12  498  888  215  4.442 × 10−5  2.014 × 10−2  3.586 × 10−2  <  <  1.311 × 10−5  Stars, ISM  16  504  897  219  1.008 × 10−4  2.028 × 10−2  3.604 × 10−2  <  2.915 × 10−5  8.005 × 10−5  Stars, BHs, XRBs, ISM  18  506  900  221  1.464 × 10−4  2.036 × 10−2  3.612 × 10−2  <  3.092 × 10−5  8.408 × 10−5  Stars, BHs, XRBs, ISM (X = 1, Y = 0)  19  510  876  244  1.927 × 10−4  2.073 × 10−2  3.664 × 10−2  –  –  –  Stars, BHs, XRBs, ISM (X = 0.92, Y = 0)  19  558  941  252  2.111 × 10−4  2.374 × 10−2  4.100 × 10−2  –  –  –  aThe neutral average is calculated by weighting the values by the neutral fraction $$x_{\rm H\,\small {I}}$$ of the cell. More neutral cells then have larger impact. This quantity is relevant for studies of the 21 cm signal. View Large First, we turn to the ionization fractions of hydrogen and helium. We limit our discussion to $$x_{\rm H\,\small {II}}$$ and $$x_{\rm He\,\small {iii}}$$, as the behaviour of $$x_{\rm He\,\small {II}}$$ is very similar to the one of $$x_{\rm H\,\small {II}}$$. From z = 18 to 10, they increase logarithmically with redshift, and the main driver is the stars, which dominate the emissivity budget. The scarcity of BHs makes them irrelevant on cosmic scales at z ≥ 10. Similarly, the other source types, albeit being widespread, do not contribute significantly to the statistics as they are far less luminous than stars. At z = 18, the volume and mass averaged ionization fractions are well below our convergence limits of 10−5, independently from the combination of source types. This is true also for the median value (which provides insight into the state of most of the IGM, as also discussed by Ross et al. 2017), as the vast majority of the IGM is fully neutral. At z = 14, the first redshift shown in Table 1, the median value is still below 10−5 for all combinations of source types, while the volume and mass averaged ionization fractions now have increased to ∼10−5 and 10−4, respectively. At z = 12 (i.e. 70 Myr later) the median $$x_{\rm H\,\small {II}}$$ is still below 10−5, except for the case when all sources are present, for which it is 1.0 × 10−5. The volume and mass averaged $$x_{\rm H\,\small {II}}$$ are higher (∼10−3), with slight variations depending on the combination of source types. At z = 12, when all source types are present, we thus observe that the IGM is in the beginning stages of a cosmic, partial ionization phase. Another 100 Myr later, at z = 10, the IGM has been under the influence of a factor of 5 more stellar sources, as well as one additional BH. The ionization fractions have increased by an order of magnitude, reaching a volume average $$x_{\rm H\,\small {II}} \sim 2 \times 10^{-2}$$ independent of source types, and a mass average that is about twice as large. The median hydrogen ionization fraction ranges from below 10−5 (with stars or BHs), to 4 × 10−5 with XRBs, to >10−4 with the ISM. We can understand the small differences in ionization fractions to come from their sensitivity to highly or fully ionized cells, whose prominence is totally dominated by the contribution from stars. The median, on the other hand, is not as prone to biasing by the highly/fully ionized cells. When the ISM is present in addition to the stars, helium is also doubly ionized on cosmic scales, albeit only weakly. At z = 10 the IGM is predominantly fully neutral when only stars and BHs are present, while the majority of the gas is partially ionized in the presence of more numerous and diffuse energetic sources. We now investigate the temperature evolution. At z = 18 the volume and mass averaged temperatures are 11 and 12 K, respectively, irrespective of combination of source types. The median temperature is also comparable, being 10 K. The above statistics are identical to those we obtain directly from the hydrodynamical simulation MBII, i.e. they are unaffected by photoionization. At z = 12, the temperature statistics differ. The volume (mass) average is ∼70 (140) K, independent of the combination of source types, compared to a gas temperature of 38 (78) K in MBII. Weighing the temperature by the hydrogen neutral fraction we find a value ∼50 K. The small differences between combinations of source types indicate that stars are the main driver of the evolution of the averaged temperatures, as they dominate the global emissivity budget. The median temperatures, however, are much lower, around 10 K. Having both XRBs and the ISM present in addition to stars and BHs does not raise the median temperature by z = 12. At z = 10, we have larger differences in the temperature statistics. This coincides with the widespread partial (but low) ionization we found. The volume, mass, and neutral fraction averaged temperatures are now in the range (496–506), (886–900), and (214–221) K, respectively, depending on the increasing combination of source types. The median temperature is however the least biased indicator of the state of the majority of the IGM. With only stars (and BHs), it is 11 K, but with all source types present, it is 18 K. This is an order of magnitude higher than what we obtain from the MBII, where the median temperature at z = 10 is 6 K. To summarize, we find that the IGM at z = 12 is showing its first signs of transitioning into one of two possible states, depending on source populations. With the most conservative assumption, i.e. in the presence of only stars (and BHs), the IGM is mostly fully neutral and cold (10 K), and it remains such down to z = 10. On the other hand, if XRBs and the radiation from the ISM in galaxies are also accounted for, the IGM becomes mainly partially ionized by z = 10. This very low partial ionization is accompanied by a further temperature increase of ∼10 K in all statistics. However, they also show that different regions of the IGM will have different temperatures. This in turn can have observable consequences. Thus, there is a modest EoH, and it largely coincides with large-scale low partial ionization from energetic sources. We now turn to examine the details of this epoch. 3.1.2 Evolution on light-cones To better understand the morphological impact of the source types, in Fig. 3 we illustrate how the hydrogen ionization state, $$x_{\rm H\,\small {II}}$$, evolves in light-cones from z = 16 to 10. Figure 3. View largeDownload slide Light-cones showing the evolution of the ionized hydrogen fraction $$x_{\rm H\,\small {II}}$$ in the simulations with different combinations of source types, as indicated in the labels. The vertical size is 100 h−1 cMpc and the aspect ratio between the axes is given at z = 10, hence the evolution in the angular diameter distance is not taken into account. The combined effect of all source types leaves the IGM without fully neutral regions at z ≲ 11.5. Figure 3. View largeDownload slide Light-cones showing the evolution of the ionized hydrogen fraction $$x_{\rm H\,\small {II}}$$ in the simulations with different combinations of source types, as indicated in the labels. The vertical size is 100 h−1 cMpc and the aspect ratio between the axes is given at z = 10, hence the evolution in the angular diameter distance is not taken into account. The combined effect of all source types leaves the IGM without fully neutral regions at z ≲ 11.5. We find that stars are determining the existence and shape of the fully ionized regions: clusters of galaxies produce ionization bubbles with sharp ionization fronts. These bubbles grow mostly separately in our redshift range. At z = 10, the IGM is still predominantly neutral, but ionized bubbles are present within distances of tens of cMpc between each other. With stars we thus have a clear dichotomy with fully ionized bubbles residing in an otherwise fully neutral IGM. As with the seeding procedure adopted in the MBII simulation only a handful of BHs is present at z > 10, the global history is unaffected by them. The first BH appears at z = 13, while the brightest is at z = 11 with MAB = −17.6. Its ionizing photon production rate is εBH = 2.05 × 1053 photons s−1, its mass is MBH = 1.17 × 106 h−1 M⊙, and it accretes with a rate of $$\dot{M}_{\rm BH} = 10^{-2} \,{\rm M}_{\odot }$$ yr−1. The stars within its host galaxy are even brighter, with MAB = −18.7, and produce εstars = 3.22 × 1053 photons s−1. In a companion paper focused on the EoR we will investigate the BHs impact on the later stages of cosmic reionization, as well as the effect of a different seeding procedure. Including more energetic galactic sources (XRBs and the ISM), we note a significant change in the ionization history of the IGM, which becomes partially ionized ($$10^{-5} < x_{\rm H\,\small {II}} < 10^{-1}$$) several cMpc outside the fully ionized regions. On the other hand, these sources do not significantly contribute to extending the fully ionized regions, which are mainly governed by stars. The spectra of the XRBs peak at ∼2 keV (Madau & Fragos 2017), leaving fewer photons nearer the H i ionization threshold. Their long mean free path means they interact after being redshifted significantly, leaving smooth partial signatures of the order $$x_{\rm H\,\small {II}} \sim 10^{-5}$$ tens of cMpc outside their sources. The spectrum of the ISM is a broken power law, with a high-energy tail fainter than that of the XRBs, but more photons in the hard UV regime, which are more readily absorbed closer to the source. With the ISM, we thus see stronger partial ionization gradients, $$x_{\rm H\,\small {II}} \sim 10^{-3}$$–10−5, and a more patchy ionization morphology compared to XRBs. The combined effect of having all source types present in galaxies is shown in the lowermost panel of Fig. 3. The stars determine the extent of the fully ionized regions that, at z = 10, still make up a very small part of the IGM. The extent of the partially ionized regions goes well beyond those provided by XRBs and the ISM alone. Combined, they leave a much larger fraction of the IGM in a partially ionized state. In Fig. 4 we show the temperature evolution light-cones as well. The stars effectively determine the temperature of 104 K in the regions we recognize from the previous plot to be fully ionized. With additional source types, the temperature of these regions does not markedly change, but we do however have heating of the otherwise cold, T < 10 K, IGM when it is subject to partial and low ionization. Combined, the XRBs and the ISM provide heating that extends further than either can provide alone, indicating that their combined effect is not simply dominated by the largest of each. At z = 10, there is thus a significant difference between the thermal states the IGM can be in, depending on the source types. Figure 4. View largeDownload slide Light-cones showing the evolution of the IGM gas temperature under the presence of different source types as explained in Fig. 3. Figure 4. View largeDownload slide Light-cones showing the evolution of the IGM gas temperature under the presence of different source types as explained in Fig. 3. 3.2 The IGM at $$\boldsymbol {z=12}$$ and $$\boldsymbol {z=10}$$ As we have found that the IGM may transition into a predominantly partially ionized, lightly heated state by z = 10 under the influence of XRBs and the ISM, here we examine its state at z = 12, when this transition has begun. This also coincides with the redshift at which we can investigate the effects of the first BH in our simulation volume, and allows us to quantify its impact, if any. The reader can refer to Table 1 for some numbers. 3.2.1 Morphology In Fig. 5 we plot maps of (from top to bottom) $$x_{\rm H\,\small {II}}$$, $$x_{\rm He\,\small {II}}$$, $$x_{\rm He\,\small {iii}}$$, and temperature at z = 12 for different combinations of source types. These are slices through the volume, chosen to contain the only BH in the simulation box, residing in a bright galaxy in the upper right-hand corner. This BH accretes with a rate of 10−2 M⊙ yr−1 and has a mass of 5.1 h−1 × 105 M⊙. Its absolute magnitude is MAB = −16.2, and it has an ionizing emissivity of εBH = 5.6 × 1052 photons s−1. It resides in a galaxy where the stars have a higher emissivity, εS = 4.6 × 1054 photons s−1, while the XRBs contribute with εXRBs = 1.4 × 1050 photons s−1 and the ISM with εISM = 4.1 × 1051 photons s−1. Figure 5. View largeDownload slide Maps (slices through the volume) in the plane of the only BH at z = 12. The columns indicate different combinations of source types. The upper rows show the quantity in question, and the lower rows show the difference with respect to the first panel, which refer to simulations with stars only. In the maps for $$x_{\rm H\,\small {II}}$$, we also show the location of the sources as circles and their emissivities (denoted with the colour). The only BH is shown as a triangle in the upper right-hand corner. The maps are 100 h−1 cMpc wide. Figure 5. View largeDownload slide Maps (slices through the volume) in the plane of the only BH at z = 12. The columns indicate different combinations of source types. The upper rows show the quantity in question, and the lower rows show the difference with respect to the first panel, which refer to simulations with stars only. In the maps for $$x_{\rm H\,\small {II}}$$, we also show the location of the sources as circles and their emissivities (denoted with the colour). The only BH is shown as a triangle in the upper right-hand corner. The maps are 100 h−1 cMpc wide. We also indicate the position and emissivity of the sources with coloured circles in the maps of $$x_{\rm H\,\small {II}}$$, while the BH is represented by a filled triangle. To highlight the impact of different source types, we also plot absolute differences of the quantities with respect to the simulations with stars only. Note that we create the difference maps by truncating values below our numerical convergence limit of xi = 10−5 for $$i={{\rm H}\,\small {II}}, \, {\rm {He\,\small {II}}}$$, and $${\rm He\,\small {iii}}$$ (see Appendix A). The stars are responsible for the fully ionized H ii bubbles and their morphology, as we also saw in the light-cones. He ii follows H ii due to the similar ionization potentials, but as He i has an ionization cross-section progressively larger than that of H i at photon energies hPν ≥ 24.6 eV, it is then more readily ionized, and its spatial distribution is generally more extended. Also note that sources having in their spectra higher energy photons (see the cases stars+XRBs and stars+ISM), which are redshifted at the simulated scales, provide a contribution first to He i ionization and then to H i. These combined effects give us a larger extent of partially ionized $$x_{\rm He\,\small {II}}$$ than $$x_{\rm H\,\small {II}}$$, while there is still little overlap between the fully ionized regions, leaving a very patchy H ii and He ii ionization morphology. Where sources are more grouped together, they leave imprints on their surroundings similar to what one would expect of separate, but much brighter sources. The extent of ionization must thus be considered an effect of the ionizing photons arising within a volume, rather than from a single source. This helps the partially ionized regions to grow considerably in the cases with XRBs or the ISM. The only BH at z = 12 accounts for partial ionization in its vicinity, spanning several cMpc, beginning with a smoothing of the otherwise sharp ionization front already created by its host and surrounding galaxies. However, a partially ionized tail is not an unique feature of BHs, as it is also seen outside galaxies that have XRBs and the ISM, which both contribute with high-energy photons. As in this epoch, galaxies are more common than BHs, their XRBs and ISM completely dominate the partial ionization at cosmic scales, while the contribution by the ISM component is certainly more ubiquitous than that from the XRBs. The last panels of both H ii and He ii cases finally show the concerted impact of all sources, resulting in a larger overlap of red areas (He ii, in particular, is found in the process of merging into a single one) and in a smoother transition from the fully ionized to the neutral regions. He ii ionization, on the other hand, is much more sensitive to high-energy UV photons (54.4 ≤ hPν ≤ 200 eV), while its cross-section at soft X-ray energies decreases by two orders of magnitude. This is the reason why partially ionized He iii regions with ionizations up to $$x_{\rm He\,\small {iii}} \sim 10^{-3}$$ are found only within the ionized bubble of the BH or correlating with the brightest galactic sources: the strongest XRBs, and the strongest/most clustered galaxies with a hot ISM contribution. Maps of the IGM temperature resulting from the various source combinations are finally found in the bottom panels. The thermal state of the IGM shows signatures both from the hydrodynamic simulation and the heating from the sources. The ideal gas of the hydrodynamic simulations is heated in overdense regions, whereas shocks occur at scales below the resolution of our grid (∼400 h−1 ckpc). As fully ionized H ii and He ii regions are mainly driven by stars, gas at photoionization temperatures of ∼104 K strictly correlates with them, as shown in the bottom left-hand panel by white areas. Extended blue areas, corresponding to temperatures of ∼500 K, are found along filaments connecting the stellar sources. These regions, on the other hand, do not have a clear counterpart in the hydrogen ionization pattern because their ionization fraction is below the convergence threshold xi < 10−5. Difference maps show that the addition of other source types does not significantly change the global temperature pattern as set-up by the stars. A local increase in T (from 100 up to 1000 K) is found only around the BH and around galaxies having a hot ISM. Interestingly, we note that source clustering still plays a relevant role in changing T but its pattern is less extended compared to the one of the corresponding H ii region. Also note that the temperature boost expected in the He iii regions from the release of He ii electrons, as found in Kakiichi et al. (2017), cannot be appreciated around the first quasars because of the lower emissivity and softer spectral shape. In fact, their BH had an emissivity εBH = 1.4 × 1056 photons s−1, and a power-law spectrum with spectral index −1.5, shifting more ionizing photons closer to the ionization potential of He ii. By comparing the relative contribution of different sources we immediately see that XRBs do not heat appreciably the IGM that they partially ionize. We find instead that the ISM emission, mainly providing softer photons, contributes to a heating of some tens of degrees in the regions that it partially ionizes. These do not extend more than a few cMpc outside the fully ionized regions. The different impact of XRBs and ISM can be easily explained in terms of their spectral distributions (see Fig. 2). The energy released by an absorbed X-ray photon is in fact shared between direct photoionization, secondary ionization, and gas heating. The fraction of the energy that goes into heating increases with the ionization fraction of the gas. Heating is more effective for less energetic photons (Dalgarno, Yan & Liu 1999), and the ISM, which has a softer spectrum, is therefore more effective in heating the gas. Combining the effect of XRBs and the ISM, we expect a more efficient heating: XRBs should in fact drive pre-ionization in regions that can be influenced by successive contribution of ISM photons. This concerted effect is clearly visible as an increase of the extent of the heated regions created by the presence of all source types. Also check the difference map to have a better feeling of the difference introduced by the concerted, self-consistent, impact of all the source types. To summarize this section: our models predict an IGM at z = 12 that is dichotomous, with many, patchy fully ionized and hot regions, which extent and abundance are determined by the stars. These regions reside in an IGM that is either neutral and cold, or partially ionized when XRBs and the hot ISM of the galaxies irradiate harder UV and X-ray photons. They can heat the IGM only up to a few tens of degrees, and although this heating is certainly boosted when the XRBs pre-ionize the IGM, it is the ISM that sustains the heating. 3.2.2 IGM phase statistics In Fig. 6 we show the phase state of the IGM at z = 12 under the influence of different combinations of source types. Each panel refers to the distributions of the total 1.68 × 107 cells with a given temperature as a function of overdensity δ.6 Figure 6. View largeDownload slide Phase diagrams showing the distribution (in fraction of the total 2563 = 1.68 × 107 cells) of the IGM at z = 12 in different thermal (T, y-axes), overdensity ($$\delta \equiv n/\bar{n}$$, x-axis), and ionization ($$x_{\rm H\,\small {II}}$$, rows) states for various combinations of source types, as indicated by the labels. The horizontal solid (dashed) lines indicate the mean (median) temperature in the volume, and the vertical solid (dashed) lines indicate the mean (median) overdensity δ of the part of the volume with the given ionization state. The horizontal dotted lines indicate the CMB temperature. The percentages indicate the fraction of the volume that is in the given ionization state. Figure 6. View largeDownload slide Phase diagrams showing the distribution (in fraction of the total 2563 = 1.68 × 107 cells) of the IGM at z = 12 in different thermal (T, y-axes), overdensity ($$\delta \equiv n/\bar{n}$$, x-axis), and ionization ($$x_{\rm H\,\small {II}}$$, rows) states for various combinations of source types, as indicated by the labels. The horizontal solid (dashed) lines indicate the mean (median) temperature in the volume, and the vertical solid (dashed) lines indicate the mean (median) overdensity δ of the part of the volume with the given ionization state. The horizontal dotted lines indicate the CMB temperature. The percentages indicate the fraction of the volume that is in the given ionization state. The upper row shows the highly and fully ionized IGM to be independent from the type of sources. The IGM in this state is hot, with mean and median temperatures of 1.2 × 104 K at a median overdensity δ = 2. The stars that drive full ionization primarily affect the overdense regions. As reflected in the values shown in Table 1, a negligible fraction of the gas is in this high ionization state at z = 12. The middle row shows the phase state of cells with partially ionized IGM. We see that a very small fraction is in this state with stars alone, and adding a single BH does not change this. Including XRBs (ISM) increases the percentage of volume in this partially ionized state to 7 per cent (28 per cent), while having all sources present, 50 per cent of the volume becomes partially ionized. A comparison with the lower row, which refers to the neutral IGM, clearly shows the effect of the additional energetic photons. The IGM under the presence of only stars is mostly cold, but with a fraction of cells that are both overdense (log δ > 0) and hotter than the majority (T > 10 K). A larger fraction of these cells transition to a partially ionized state when including the XRBs. An even larger fraction, both from the cold and underdense and the overdense, hotter gas, becomes partially ionized under the influence of the ISM photons. There is no apparent preference in terms of δ or T on the IGM that becomes partially ionized. This results from a combination of factors such as the location (overdense regions), emissivity (this is lower for the XRBs than for the ISM), and spectral shape (the XRBs have spectra harder than the ISM) of the sources. More than 99 per cent of the IGM is fully neutral with stars only, while this fraction is strongly decreased with XRBs (93 per cent), the ISM (72 per cent), or all source types (50 per cent). It should be noted that, as already pointed out by other authors (see e.g. Ross et al. 2017), the partially ionized and warm cells found in the presence of stellar type sources might only arise from a lack of spatial resolution. More specifically, 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, while in reality part of the gas in the cell should be neutral and cold, and part fully ionized and hot. We will discuss this issue in more detail in a companion paper. In Fig. 7 we present histograms of the temperature and ionization fractions at z = 12 and 10, showing a quantitative evolution of the physical state of the IGM. Figure 7. View largeDownload slide Fraction of cells in different ionization and thermal states for different combinations of source types (stars: solid black, no ticks; stars and BHs: blue, one tick; stars and XRBs: purple, two ticks; stars and ISM: red, three ticks; all: yellow, four ticks). We show the states at z = 12 (upper row) and at z = 10 (lower row). The percentages indicate the fraction of the IGM that is shown in the histogram. The remaining IGM has an ionization fraction xi lower than 10−5 (for $$i={{\rm H}\,\small {II}}, \, {\rm {He\,\small {II}}}$$, and $${{\rm He}\,\small {iii}}$$). Figure 7. View largeDownload slide Fraction of cells in different ionization and thermal states for different combinations of source types (stars: solid black, no ticks; stars and BHs: blue, one tick; stars and XRBs: purple, two ticks; stars and ISM: red, three ticks; all: yellow, four ticks). We show the states at z = 12 (upper row) and at z = 10 (lower row). The percentages indicate the fraction of the IGM that is shown in the histogram. The remaining IGM has an ionization fraction xi lower than 10−5 (for $$i={{\rm H}\,\small {II}}, \, {\rm {He\,\small {II}}}$$, and $${{\rm He}\,\small {iii}}$$). As a baseline, we compare the distributions to the case with stars alone, where the curves for $$x_{\rm H\,\small {II}}$$ and $$x_{\rm He\,\small {II}}$$ are almost flat, with only a surplus of cells in a highly or fully ionized and hot (T ∼ 104 K) state, while the large majority is neutral and cold (T ∼ 10 K), both at z = 12 and 10. The fraction that is highly/fully ionized and hot has increased an order of magnitude by z = 10. With the addition of other source types, partial ionization is strongly increased, especially at $$x_{\rm H\,\small {II}} < 10^{-3}$$. Only at z = 10 a discernible difference is visible in the temperature distributions, with a shift of the low temperature peak by several tens of degrees when all source types are present. This is also the redshift where there is largest difference in the fraction of the IGM that has $$x_{\rm H\,\small {II}} > 10^{-5}$$, denoted as percentages in the figures. At z = 12, 1 per cent of the IGM has $$x_{\rm H\,\small {II}} > 10^{-5}$$ when only stars are present, while this fraction has increased to 50 per cent with all sources. At z = 10, this difference is even more significant, as the percentage goes from 6 to 100 per cent when including either XRBs, the ISM, or both. The behaviour of $$x_{\rm He\,\small {II}}$$ closely follows the one of $$x_{\rm H\,\small {II}}$$. The distribution and presence of doubly ionized helium changes significantly depending on the source types. While the presence of XRBs has an impact only at $$x_{\rm He\,\small {iii}} < 10^{-4}$$, the ISM ionizes more and to higher values of $$x_{\rm He\,\small {iii}}$$. At z = 12, less than 1 per cent of the IGM has $$x_{\rm He\,\small {iii}} > 10^{-5}$$. This fraction, as well as the distribution across all $$x_{\rm He\,\small {iii}}$$, increases by an order of magnitude for all source types at z = 10, leaving 4 per cent of the helium in the IGM in a partially ionized state. 3.3 Simulations without helium The presence of helium in our simulations affects the ionization state of hydrogen and the temperature evolution. The ionization cross-section of hydrogen decays roughly as ν−3, and at the ionization threshold of 24.6 eV for He i, the cross-section of helium is a factor of 6 larger than that of hydrogen. Similarly, at 54.4 eV, the ionization potential of He ii, the helium cross-section is 16 (neutral) and 13 (singly ionized) times larger than that of hydrogen. Sources that have spectra that allow for helium to get ionized twice will result in the release of free electrons that can further heat the IGM, as found by Ciardi et al. (2012). We have run simulations without helium, i.e. Y = 0, including all source types. We have either replaced the now absent helium with hydrogen, effectively increasing its number fraction from X = 0.92 to X = 1, or we have simply not solved for helium chemistry, but kept using X = 0.92. Ionization fractions and temperature statistics are presented in Table 1. When X = 0.92, we effectively have more photons to ionize hydrogen atoms and heat the IGM, resulting in ionization fractions and temperatures higher at any given redshift compared to the case in which helium is included. This effect can be mitigated to some extent by increasing the number fraction of hydrogen to X = 1. In Fig. 8 we have histograms showing the distribution of $$x_{\rm H\,\small {II}}$$ and T at z = 12 and 10. Dashed and dashed–dot lines indicate simulations without helium. The increase in ionization fractions mentioned above is accompanied by a uniform shift at $$x_{\rm H\,\small {II}} < 10^{-1}$$ towards higher $$x_{\rm H\,\small {II}}$$ at both z, while the distributions of highly and fully ionized cells do not change. The absence of helium also contributes to another interesting effect: the maximum temperature reached in the IGM is lowered by approximately 0.2 dex. Without helium, at z = 12 (z = 10) the maximum temperature is T = 22 986 (36 141) K, while with helium it is T = 32 313 (38 709) K. We also have some additional heating of the cold IGM at z = 10, as seen from the increase in the median temperatures in Table 1. Figure 8. View largeDownload slide Distributions of volume averaged temperature and ionization fraction of hydrogen for simulations with (solid; black line) or without helium (dashed; purple lines are for X = 0.92 and Y = 0, while dash–dotted yellow lines are for X = 1 and Y = 0) at z = 12 and 10. The percentages, with the same colours and reference as the lines, indicate what fraction of the IGM has $$x_{\rm H\,\small {II}} \ge 10^{-5}$$, and hence is partially or strongly ionized. Figure 8. View largeDownload slide Distributions of volume averaged temperature and ionization fraction of hydrogen for simulations with (solid; black line) or without helium (dashed; purple lines are for X = 0.92 and Y = 0, while dash–dotted yellow lines are for X = 1 and Y = 0) at z = 12 and 10. The percentages, with the same colours and reference as the lines, indicate what fraction of the IGM has $$x_{\rm H\,\small {II}} \ge 10^{-5}$$, and hence is partially or strongly ionized. 4 DISCUSSION Our simulations show that stars are the principal source of IGM ionization throughout the EoH: in fact they fully dominate the UV part of our spectra (see Fig. 2) and have a large spatial coverage. They are responsible for the volume averaged ionization fractions of both hydrogen and helium, which at z = 10 are $$x_{\rm H\,\small {II}} \sim x_{\rm He\,\small {II}} \sim 2.0 \times 10^{-2}$$. The contribution of stars also drives the shape and extent of local H ii and He ii bubbles, having sharp ionization fronts due to the short mean free path of their UV photons. As a rapid phase transformation is not expected until $$x_{\rm H\,\small {II}} \sim 0.1$$ (Furlanetto & Oh 2016), at z = 10 the IGM is still patchy, with its vast majority being fully neutral, and only a small fraction fully ionized. Only in the presence of more energetic sources, does the morphology and state of the ionized regions change. The seeding procedure adopted in MBII, plants seeds of mass 105 M⊙ in haloes with M > 1010 M⊙. This results in a scarcity of BHs at z > 10, which makes their contribution to IGM ionization and heating negligible. This conclusion, though, is strongly model dependent. We will discuss in more detail the effect of the seeding procedure and the possibility of populating with BHs also in smaller haloes in a companion paper focused on the EoR. It is also important to note that the signature of BHs is degenerate with that of other, more abundant, energetic sources. The XRBs, accounting for both the subdominant LMXBs and the SFR-tracing HMXBs, have spectra that peak at keV scales. They provide a smooth, long-range (several cMpc) partial ionization of hydrogen and helium ($$x_{\rm H\,\small {II}} = x_{\rm He\,\small {II}} \gtrsim 10^{-5}$$). Otherwise, they contribute negligibly, in line with the findings of the semi-analytical works of Madau & Fragos (2017) and Sazonov & Khabibullin (2017). Our model of XRBs thus ionizes less than what found by Ross et al. (2017), who employ a comparable 3D RT post-processing approach, albeit on a dark matter N-body simulation. It should be noted though that their emissivities and source locations were obtained through halo mass scaling relations, which differ from our approach relying on a consistent evolution of baryonic sources (and their properties) from a hydrodynamical simulation. In particular, by z = 13 Ross et al. (2017) find a partial ionization of the cold IGM that is orders of magnitudes higher than what indicated by our computation. They also find a higher volume averaged $$x_{\rm H\,\small {II}} \sim 10^{-2}$$ independent of whether XRBs were included or not. However, while our XRB spectra peak at keV-scales, theirs is a power law with index −1.5. The luminosity of our XRBs scales with either the SFR and metallicity (HMXBs) or the stellar age and mass (LMXBs) in the halo where they reside, while Ross et al. (2017) assume a constant X-ray production efficiency of each halo. Our results are more in line with those of Meiksin et al. (2017), who found the HMXBs to heat the IGM to T = 22 K by z = 10. In their work, though, the increment in temperature from a model with stars only is ∼20 K, compared to our ∼1 K. This could be due to our RT approach capturing only the effect of X-rays up to distances of the order of the box length. On the other hand, Meiksin et al. (2017) do not consider photons with energies below 200 eV, which are more relevant on smaller scales. Note also that we are still likely to overestimate the already small contribution from XRBs. Das et al. (2017), for example, found the XRB spectra to be attenuated by the ISM, albeit not as strongly as advocated by Fragos et al. (2013b), but still more than what we have assumed with our constant escape fraction of fesc = 15 per cent for photons with energies lower than 200 eV. Our multifrequency RT highlighted the importance of the ISM contribution as its spectrum provides large abundance of photons in the hard UV and soft X-rays. We found that this spectral difference is of great importance, a conclusion also shared by Pacucci et al. (2014), or by Fialkov et al. (2014b), who similarly showed that the hardness of the spectral shape of the XRBs would make them inefficient in heating the IGM. There is however some uncertainty on the extent of the contribution from the ISM, which could be larger, as the supernova-driven heating process of the ISM gas will likely result in emission outside the galactic planes (see e.g. Chevalier & Clegg 1985 or Strickland & Stevens 2000). In this case, it is less likely that the ionizing UV photons are strongly absorbed by the host galaxy, which is what we assume with our escape fraction. However, although the production rate of UV photons is uncertain, the ISM is known to supply galaxies ubiquitously with observable amounts of soft X-rays (Mineo et al. 2012b). The Meiksin et al. (2017) work mentioned above includes also the contribution of the ISM, which raises the IGM temperature at z = 10 to T = 6 or 34 K (compared to our median T = 16 K with stars and the ISM), depending on the wind model adopted. One should be careful to attribute long range partial ionization and heating signatures, if observed, solely to the ISM. These could possibly be degenerate with those from cosmic rays. While Leite et al. (2017) found that they should heat the IGM in the immediate vicinity of star-forming haloes, this happens because of the cosmic rays confinement. If this confinement is weaker than what estimated by Leite et al. (2017), their contribution could be very similar to that of the ISM. Another RT study examining the effect of X-rays on heating and ionization is that of Baek et al. (2010). They include X-rays from all of their sources, labelling this contribution as QSOs, but noting that XRBs and supernova remnants also fall into this category. Their scaling of the X-ray luminosity with the SFR allows for some comparison to our case where all source types are present. Their mean ionization fraction at z = 10 is $$x_{\rm H\,\small {II}} \sim 10^{-2}$$ for simulations with helium, but without X-rays. They note that including QSOs that account for up to 10 per cent of the total flux does not change this result, similarly to the findings of Ross et al. (2017), and now also ours. We have run simulations without helium, finding that the ionization fraction, contrary to Baek et al. (2010, cf. their models S1 and S3), increases to $$x_{\rm H\,\small {II}} = 2.4 \times 10^{-2}$$ (keeping X = 0.92), or $$x_{\rm H\,\small {II}} = 2.1 \times 10^{-2}$$ (with X = 1). This is caused by the additional photons that are available for H i ionization when helium is absent, aligning our results with those of Ciardi et al. (2012). Including helium thus provides a slower progression of the cosmic heating and reionization. As for He iii, we find that only faint, diffuse regions having $$x_{\rm He\,\small {iii}} \sim 10^{-3}$$ start to appear near bright energetic sources. We do not find the widespread existence of highly/fully ionized He iii that Ross et al. (2017) found at higher redshifts. Consistently with Madau & Fragos (2017), we find that XRBs are not able to heat significantly the IGM due to their hard spectra. Our results show that also the ISM, which has a much softer spectrum, is not able to raise the IGM temperature significantly by z = 10. Our median temperatures are in fact T = 11 and 18 K for simulations with stars only and all sources, respectively. This is somewhat higher than the value reported in Baek et al. (2010), where the temperature of the IGM with an ionization fraction $$x_{\rm H\,\small {II}} < 10^{-2}$$ was ranging between 3 and 7 K without and with QSOs, respectively. Our volume averaged temperatures are, on the other hand, much lower than those of Ross et al. (2017), who found T ∼ 600 K at z = 13 with stars only and T ∼ 1000 K with HMXBs, compared to our volume averaged T = 34 K at z = 13 with stars and XRBs. Predicting the thermal state of the IGM during the EoH is of crucial importance for 21 cm experiments because an IGM with a temperature higher than that of the CMB will be seen in emission, and if lower, it will be seen in absorption. Interestingly, at z = 10 the CMB temperature is TCMB = 30 K. The margins between a 21 cm signal seen in emission or absorption are thus close. Note, though, that while our mean temperatures at this z are much higher than TCMB, they do not represent the state of the vast majority of the IGM, which is either fully neutral (with stars only) or partially ionized (with all sources). We find partially ionized regions that are warmer than the CMB already at z = 12, as indicated from the neutral and mass averaged temperatures. We defer to a future work further investigation of the effect that different source types leave on the 21 cm signal. Observationally, very late (after z = 8.4) heating of the IGM has been ruled out by measurements of 21 cm line (Pober et al. 2015). Less constrained is a period of early and highly efficient heating, which would suppress small-scale structure formation by increasing the Jeans mass (Ostriker & Gnedin 1996), disfavouring H2 creation associated with X-ray fluxes that do not dominate the ionizing budget (Oh 2001), and increase the Thomson-scattering optical depth of the CMB (Ricotti & Ostriker 2004). We observe that both the distribution of the sources and their surrounding ionized regions trace the underlying density field. This can be seen particularly well in the phase maps of Fig. 6, where the fully ionized regions reside where log δ > 0. We also see clustering of the sources in the global maps (Fig. 5). This clustering induces (a) ionized regions that are more extended and (b) heating of the nearby neutral IGM, when comparing the surroundings of sources with similar emissivities but that exist either in a clustered environment or alone. This distinction (i.e. clustered sources versus single source) is crucial in the interpretation of observations aiming at investigating the physical properties of the high-redshift IGM. 5 CONCLUSION In this paper we have examined the sources that could end the Dark Ages, initiate the Cosmic Dawn, and drive the EoH, when the IGM starts transitioning from a cold and neutral phase to a hot and ionized one. The Cosmic Dawn concludes with the EoR at z < 10, which we will explore in a companion paper. We have investigated how stars, XRBs, and the shock heated ISM of galaxies, as well as accreting nuclear BHs, could heat and ionize the IGM at cosmic scales. Each source type has spectral characteristics determined by underlying physical processes. The luminosity of the stars is determined by their masses, ages, and metallicities, while the ISM and XRBs are sensitive to the SFR in galaxies, as well as their masses and metallicities. The BHs, on the other hand, shine with the rate they accrete matter. We have used the hydrodynamical simulations MBii (K15) to provide us with the physical properties of the sources, their location, and abundance, as well as their environment (temperature and baryonic density). Having established the sources, their spectra, and ionizing emissivities, we post-processed MBII with the RT code crash (e.g. Ciardi et al. 2001; Graziani et al. 2013). Our main findings can be summarized as follows. Stars drive the extent, shape, abundance, and temperature of the fully ionized H ii and He ii regions, but do not leave the remaining IGM in any heated or partially ionized state. Only 7 per cent of the overall IGM is in a non-neutral state at z = 10 with stars alone. With our seeding prescription, nuclear BHs are scarce and do not contribute significantly to heating or ionization at z > 10. XRBs contribute to partial, uniform ionization ($$x_{\rm H\,\small {II}} \sim 10^{-5}$$) of the dilute IGM on cMpc scales, but do not significantly heat it by z = 10. This can be attributed to their hard, keV-peaked spectra. The ISM of the galaxies contributes to a larger extent of partial ionization than the XRBs do. Their softer spectra provide heating and ionization, with $$x_{\rm H\,\small {II}} \sim 10^{-3}$$ and T ≳ 10 K up to a few cMpc around the fully ionized regions at z = 12 and lower partial ionization further out. In concert, the ISM and the XRBs induce an ionization fraction $$x_{\rm H\,\small {II}} \ge 10^{-5}$$ in 50 per cent (100 per cent) of the IGM at z = 12 (z = 10). However, at z = 10 the IGM is still predominantly cold, with median temperatures in the range (11–19) K under the influence of stars to all source types. At z = 10 the IGM will be seen in both emission and absorption in 21 cm as the neutral averaged temperatures are T > 200 K for all combinations of source types. Helium reionization has begun on cosmic scales at z = 10 in the presence of XRBs and/or the ISM. If an EoH takes place at z > 10, our findings indicate that it will be a modest one. Acknowledgements We thank Rupert Croft for enlightening comments on the bright BHs; Andrea Ferrara for rewarding discussions; Garrelt Mellema, Philipp Busch, Martin Glatzle, Max Gronke, and Aniket Bhagwat for their insights into various parts of the project. MBE thanks Piero Madau for his hospitality. We thank the anonymous referee for constructive comments. LG acknowledges the support of the DFG Priority Program 1573 and from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement n. 306476. This work made use of a number of open source software packages, in particular numpy (Van Der Walt, Colbert & Varoquaux 2011), cython (Behnel et al. 2011), matplotlib (Hunter 2007), and GNU parallel (Tange 2011). Footnotes 1 https://www.isispace.nl/dutch-radio-antenna-depart-moon-chinese-mission/; accessed 2017 July 13. 2 σ8 = 0.816, ns = 0.968, $$\Omega _\Lambda =0.725$$, Ωm = 0.275, Ωb = 0.046, h = 0.701. 3 crash allows us to specify the escape fraction in different frequency bands individually for each source. 4 In our reference simulations we do not use periodic boundary conditions, although this option is available. 5 The effect of not including the evolution of binary systems is a reduction in ionizing flux. 6 In this paper the gas overdensity is defined as $$\delta \equiv n/\bar{n}$$, where n is the particle number density with units cm−3. REFERENCES Ali Z. S. et al.  , 2015, ApJ , 809, 61 https://doi.org/10.1088/0004-637X/809/1/61 CrossRef Search ADS   Allen R. J., 1969, A&A , 3, 382 Arons J., McCray R., 1970, Astrophys. Lett. , 5, 123 Baek S., Semelin B., Di Matteo P., Revaz Y., Combes F., 2010, A&A , 523, A4 https://doi.org/10.1051/0004-6361/201014347 CrossRef Search ADS   Barnett R., Warren S. J., Becker G. D., Mortlock D. J., Hewett P. C., McMahon R. G., Simpson C., Venemans B. P., 2017, A&A , 601, A16 https://doi.org/10.1051/0004-6361/201630258 CrossRef Search ADS   Behnel S., Bradshaw R., Citro C., Dalcin L., Seljebotn D. S., Smith K., 2011, Comput. Sci. Eng. , 13, 31 https://doi.org/10.1109/MCSE.2010.118 CrossRef Search ADS   Bergeron J., Salpeter E. E., 1970, Astrophys. Lett. , 7, 115 Bowman J. D., Rogers A. E. E., Hewitt J. N., 2008, ApJ , 676, 1 https://doi.org/10.1086/528675 CrossRef Search ADS   Cen R., 1992, ApJS , 78, 341 https://doi.org/10.1086/191630 CrossRef Search ADS   Chang P., Broderick A. E., Pfrommer C., 2012, ApJ , 752, 23 https://doi.org/10.1088/0004-637X/752/1/23 CrossRef Search ADS   Chardin J., Puchwein E., Haehnelt M. G., 2017, MNRAS , 465, 3429 https://doi.org/10.1093/mnras/stw2943 CrossRef Search ADS   Chevalier R. A., Clegg A. W., 1985, Nature , 317, 44 https://doi.org/10.1038/317044a0 CrossRef Search ADS   Ciardi B., Ferrara A., 2005, Space Sci. Rev. , 116, 625 https://doi.org/10.1007/s11214-005-3592-0 CrossRef Search ADS   Ciardi B., Ferrara A., Marri S., Raimondo G., 2001, MNRAS , 324, 381 https://doi.org/10.1046/j.1365-8711.2001.04316.x CrossRef Search ADS   Ciardi B., Bolton J. S., Maselli A., Graziani L., 2012, MNRAS , 423, 558 https://doi.org/10.1111/j.1365-2966.2012.20902.x CrossRef Search ADS   Couchman H. M. P., Rees M. J., 1986, MNRAS , 221, 53 https://doi.org/10.1093/mnras/221.1.53 CrossRef Search ADS   Croft R. A. C., Di Matteo T., Springel V., Hernquist L., 2009, MNRAS , 400, 43 https://doi.org/10.1111/j.1365-2966.2009.15446.x CrossRef Search ADS   Dalgarno A., Yan M., Liu W., 1999, ApJS , 125, 237 https://doi.org/10.1086/313267 CrossRef Search ADS   D'Aloisio A., Sanderbeck P. R. U., McQuinn M., Trac H., Shapiro P. R., 2017, MNRAS , 468, 4691 https://doi.org/10.1093/mnras/stx711 CrossRef Search ADS   Das A., Mesinger A., Pallottini A., Ferrara A., Wise J. H., 2017, MNRAS , 469, 1166 https://doi.org/10.1093/mnras/stx943 CrossRef Search ADS   Datta K. K., Mellema G., Mao Y., Iliev I. T., Shapiro P. R., Ahn K., 2012, MNRAS , 424, 1877 https://doi.org/10.1111/j.1365-2966.2012.21293.x CrossRef Search ADS   Davies R. D., Jennison R. C., 1964, MNRAS , 128, 123 https://doi.org/10.1093/mnras/128.2.123 CrossRef Search ADS   Degraf C., Di Matteo T., Springel V., 2010, MNRAS , 402, 1927 https://doi.org/10.1111/j.1365-2966.2009.16018.x CrossRef Search ADS   Dijkstra M., 2014, Publ. Astron. Soc. Aust. , 31, 26 https://doi.org/10.1017/pasa.2014.33 CrossRef Search ADS   Dillon J. S. et al.  , 2014, Phys. Rev. D , 89, 023002 https://doi.org/10.1103/PhysRevD.89.023002 CrossRef Search ADS   Di Matteo T., Springel V., Hernquist L., 2005, Nature , 433, 604 https://doi.org/10.1038/nature03335 CrossRef Search ADS PubMed  Di Matteo T., Colberg J., Springel V., Hernquist L., Sijacki D., 2008, ApJ , 676, 33 https://doi.org/10.1086/524921 CrossRef Search ADS   Di Matteo T., Khandai N., DeGraf C., Feng Y., Croft R. A. C., Lopez J., Springel V., 2012, ApJ , 745, L29 https://doi.org/10.1088/2041-8205/745/2/L29 CrossRef Search ADS   Eldridge J. J., Stanway E. R., 2012, MNRAS , 419, 479 https://doi.org/10.1111/j.1365-2966.2011.19713.x CrossRef Search ADS   Ewall-Wice A. et al.  , 2016, MNRAS , 460, 4320 https://doi.org/10.1093/mnras/stw1022 CrossRef Search ADS   Fabbiano G., 1989, ARA&A , 27, 87 https://doi.org/10.1146/annurev.aa.27.090189.000511 CrossRef Search ADS   Fan X. et al.  , 2000, AJ , 120, 1167 https://doi.org/10.1086/301534 CrossRef Search ADS   Fan X. et al.  , 2006, AJ , 132, 117 https://doi.org/10.1086/504836 CrossRef Search ADS   Fialkov A., Barkana R., Pinhas A., Visbal E., 2014a, MNRAS , 437, L36 https://doi.org/10.1093/mnrasl/slt135 CrossRef Search ADS   Fialkov A., Barkana R., Visbal E., 2014b, Nature , 506, 197 https://doi.org/10.1038/nature12999 CrossRef Search ADS   Field G. B., 1962, ApJ , 135, 684 https://doi.org/10.1086/147312 CrossRef Search ADS   Fragos T. et al.  , 2013a, ApJ , 764, 41 https://doi.org/10.1088/0004-637X/764/1/41 CrossRef Search ADS   Fragos T., Lehmer B. D., Naoz S., Zezas A., Basu-Zych A., 2013b, ApJ , 776, L31 https://doi.org/10.1088/2041-8205/776/2/L31 CrossRef Search ADS   Fukugita M., Kawasaki M., 1994, MNRAS , 269, 563 https://doi.org/10.1093/mnras/269.3.563 CrossRef Search ADS   Furlanetto S. R., Oh S. P., 2016, MNRAS , 457, 1813 https://doi.org/10.1093/mnras/stw104 CrossRef Search ADS   Furlanetto S., Zaldarriaga M., Hernquist L., 2004, ApJ , 613, 1 https://doi.org/10.1086/423025 CrossRef Search ADS   Giallongo E. et al.  , 2015, A&A , 578, A83 https://doi.org/10.1051/0004-6361/201425334 CrossRef Search ADS   Ginzburg V., Ozernoi L., 1966, SvA , 9, 726 Giri S. K., Mellema G., Dixon K. L., Iliev I. T., 2018, MNRAS , 473, 2949 https://doi.org/10.1093/mnras/stx2539 CrossRef Search ADS   Gnedin N. Y., 2014, ApJ , 793, 29 https://doi.org/10.1088/0004-637X/793/1/29 CrossRef Search ADS   Graziani L., Maselli A., Ciardi B., 2013, MNRAS , 431, 722 https://doi.org/10.1093/mnras/stt206 CrossRef Search ADS   Graziani L., Salvadori S., Schneider R., Kawata D., de Bennassuti M., Maselli A., 2015, MNRAS , 449, 3137 https://doi.org/10.1093/mnras/stv494 CrossRef Search ADS   Gunn J. E., Peterson B. A., 1965, ApJ , 142, 1633 https://doi.org/10.1086/148444 CrossRef Search ADS   Hassan S., Davé R., Finlator K., Santos M. G., 2016, MNRAS , 457, 1550 https://doi.org/10.1093/mnras/stv3001 CrossRef Search ADS   Hassan S., Davé R., Mitra S., Finlator K., Ciardi B., Santos M. G., 2018, MNRAS , 473, 227 https://doi.org/10.1093/mnras/stx2194 CrossRef Search ADS   Hunter J. D., 2007, Comput. Sci. Eng. , 9, 99 https://doi.org/10.1109/MCSE.2007.55 CrossRef Search ADS   Iliev I. T., Mellema G., Ahn K., Shapiro P. R., Mao Y., Pen U. L., 2014, MNRAS , 439, 725 https://doi.org/10.1093/mnras/stt2497 CrossRef Search ADS   Kakiichi K., Graziani L., Ciardi B., Meiksin A., Compostella M., Eide M. B., Zaroubi S., 2017, MNRAS , 468, 3718 https://doi.org/10.1093/mnras/stx603 CrossRef Search ADS   Khandai N., Di Matteo T., Croft R., Wilkins S., Feng Y., Tucker E., DeGraf C., Liu M.-S., 2015, MNRAS , 450, 1349 (K15) https://doi.org/10.1093/mnras/stv627 CrossRef Search ADS   Kitayama T., Yoshida N., Susa H., Umemura M., 2004, ApJ , 613, 631 https://doi.org/10.1086/423313 CrossRef Search ADS   Komatsu E. et al.  , 2011, ApJS , 192, 18 https://doi.org/10.1088/0067-0049/192/2/18 CrossRef Search ADS   Krawczyk C. M., Richards G. T., Mehta S. S., Vogeley M. S., Gallagher S. C., Leighly K. M., Ross N. P., Schneider D. P., 2013, ApJS , 206, 4 https://doi.org/10.1088/0067-0049/206/1/4 CrossRef Search ADS   Leite N., Evoli C., D'Angelo M., Ciardi B., Sigl G., Ferrara A., 2017, MNRAS , 469, 416 https://doi.org/10.1093/mnras/stx805 CrossRef Search ADS   Liu H., Slatyer T. R., Zavala J., 2016, Phys. Rev. D , 94, 063507 https://doi.org/10.1103/PhysRevD.94.063507 CrossRef Search ADS   McGreer I. D., Mesinger A., D'Odorico V., 2015, MNRAS , 447, 499 https://doi.org/10.1093/mnras/stu2449 CrossRef Search ADS   Madau P., Fragos T., 2017, ApJ , 840, 39 https://doi.org/10.3847/1538-4357/aa6af9 CrossRef Search ADS   Madau P., Haardt F., 2015, ApJ , 813, L8 https://doi.org/10.1088/2041-8205/813/1/L8 CrossRef Search ADS   Madau P., Meiksin A., Rees M. J., 1997, ApJ , 475, 429 https://doi.org/10.1086/303549 CrossRef Search ADS   Madau P., Haardt F., Rees M. J., 1999, ApJ , 514, 648 https://doi.org/10.1086/306975 CrossRef Search ADS   Maselli A., Ferrara A., Ciardi B., 2003, MNRAS , 345, 379 https://doi.org/10.1046/j.1365-8711.2003.06979.x CrossRef Search ADS   Maselli A., Ciardi B., Kanekar A., 2009, MNRAS , 393, 171 https://doi.org/10.1111/j.1365-2966.2008.14197.x CrossRef Search ADS   Matthee J., Sobral D., Santos S., Röttgering H., Darvish B., Mobasher B., 2015, MNRAS , 451, 400 https://doi.org/10.1093/mnras/stv947 CrossRef Search ADS   Matthee J., Sobral D., Best P., Khostovan A. A., Oteo I., Bouwens R., Röttgering H., 2017, MNRAS , 465, 3637 https://doi.org/10.1093/mnras/stw2973 CrossRef Search ADS   Meiksin A. A., 2009, Rev. Modern Phys. , 81, 1405 https://doi.org/10.1103/RevModPhys.81.1405 CrossRef Search ADS   Meiksin A., Khochfar S., Paardekooper J.-P., Dalla Vecchia C., Kohn S., 2017, MNRAS , 471, 3632 https://doi.org/10.1093/mnras/stx1857 CrossRef Search ADS   Mellema G., Iliev I. T., Pen U.-L., Shapiro P. R., 2006, MNRAS , 372, 679 https://doi.org/10.1111/j.1365-2966.2006.10919.x CrossRef Search ADS   Mesinger A., Furlanetto S., Cen R., 2011, MNRAS , 411, 955 https://doi.org/10.1111/j.1365-2966.2010.17731.x CrossRef Search ADS   Mineo S., Gilfanov M., Sunyaev R., 2012a, MNRAS , 419, 2095 https://doi.org/10.1111/j.1365-2966.2011.19862.x CrossRef Search ADS   Mineo S., Gilfanov M., Sunyaev R., 2012b, MNRAS , 426, 1870 https://doi.org/10.1111/j.1365-2966.2012.21831.x CrossRef Search ADS   Mirabel I. F., Dijkstra M., Laurent P., Loeb A., Pritchard J. R., 2011, A&A , 528, A149 https://doi.org/10.1051/0004-6361/201016357 CrossRef Search ADS   Mitra S., Choudhury T. R., Ferrara A., 2018, MNRAS , 473, 1416 https://doi.org/10.1093/mnras/stx2443 CrossRef Search ADS   Nath B. B., Biermann P. L., 1993, MNRAS , 265, 241 https://doi.org/10.1093/mnras/265.1.241 CrossRef Search ADS   Ocvirk P. et al.  , 2016, MNRAS , 463, 1462 https://doi.org/10.1093/mnras/stw2036 CrossRef Search ADS   Oh S. P., 2001, ApJ , 553, 499 https://doi.org/10.1086/320957 CrossRef Search ADS   Oñorbe J., Hennawi J. F., Lukić Z., Walther M., 2017, ApJ , 847, 63 https://doi.org/10.3847/1538-4357/aa898d CrossRef Search ADS   Onoue M. et al.  , 2017, ApJ , 847, L15 https://doi.org/10.3847/2041-8213/aa8cc6 CrossRef Search ADS   O'Shea B. W., Wise J. H., Xu H., Norman M. L., 2015, ApJ , 807, L12 https://doi.org/10.1088/2041-8205/807/1/L12 CrossRef Search ADS   Ostriker J. P., Gnedin N. Y., 1996, ApJ , 472, L63 https://doi.org/10.1086/310375 CrossRef Search ADS   Ouchi M. et al.  , 2018, PASJ , 70, S13 https://doi.org/10.1093/pasj/psx074 CrossRef Search ADS   Paciga G. et al.  , 2011, MNRAS , 413, 1174 https://doi.org/10.1111/j.1365-2966.2011.18208.x CrossRef Search ADS   Pacucci F., Mesinger A., Mineo S., Ferrara A., 2014, MNRAS , 443, 678 https://doi.org/10.1093/mnras/stu1240 CrossRef Search ADS   Parsa S., Dunlop J. S., McLure R. J., 2018, MNRAS , 474, 2904 https://doi.org/10.1093/mnras/stx2887 CrossRef Search ADS   Parsons A. R. et al.  , 2014, ApJ , 788, 106 https://doi.org/10.1088/0004-637X/788/2/106 CrossRef Search ADS   Patil A. H. et al.  , 2017, ApJ , 838, 65 https://doi.org/10.3847/1538-4357/aa63e7 CrossRef Search ADS   Pawlik A. H., Rahmati A., Schaye J., Jeon M., Vecchia C. D., 2017, MNRAS , 466, 960 https://doi.org/10.1093/mnras/stw2869 CrossRef Search ADS   Penzias A. A., Scott E. H., 1968, ApJ , 153, L7 https://doi.org/10.1086/180209 CrossRef Search ADS   Planck Collaboration XLVII, 2016, A&A , 596, A108 https://doi.org/10.1051/0004-6361/201628897 CrossRef Search ADS   Pober J. C. et al.  , 2015, ApJ , 809, 62 https://doi.org/10.1088/0004-637X/809/1/62 CrossRef Search ADS   Puchwein E., Pfrommer C., Springel V., Broderick A. E., Chang P., 2012, MNRAS , 423, 149 https://doi.org/10.1111/j.1365-2966.2012.20738.x CrossRef Search ADS   Qin Y. et al.  , 2017, MNRAS , 472, 2009 https://doi.org/10.1093/mnras/stx1909 CrossRef Search ADS   Rees M. J., 1998, Proc. Natl. Acad. Sci. USA , 95, 47 CrossRef Search ADS   Rees M. J., Setti G., 1969, A&A , 8, 410 Ricotti M., Ostriker J. P., 2004, MNRAS , 352, 547 https://doi.org/10.1111/j.1365-2966.2004.07942.x CrossRef Search ADS   Robertson B. E., Ellis R. S., Furlanetto S. R., Dunlop J. S., 2015, ApJ , 802, L19 https://doi.org/10.1088/2041-8205/802/2/L19 CrossRef Search ADS   Ross H. E., Dixon K. L., Iliev I. T., Mellema G., 2017, MNRAS , 468, 3785 https://doi.org/10.1093/mnras/stx649 CrossRef Search ADS   Safarzadeh M., Scannapieco E., 2016, ApJ , 832, 4 https://doi.org/10.3847/2041-8205/832/1/L9 CrossRef Search ADS   Santos M. G., Ferramacho L., Silva M. B., Amblard A., Cooray A., 2010, MNRAS , 406, 2421 https://doi.org/10.1111/j.1365-2966.2010.16898.x CrossRef Search ADS   Sazonov S., Khabibullin I., 2017, Astron. Lett. , 43, 211 https://doi.org/10.1134/S1063773717040077 CrossRef Search ADS   Sazonov S., Sunyaev R., 2015, MNRAS , 454, 3464 https://doi.org/10.1093/mnras/stv2255 CrossRef Search ADS   Schaepman M., 2009, in Warner T., Duane Nellis M., Foody G., eds, The SAGE Handbook of Remote Sensing . Sage, London, p. 166 Google Scholar CrossRef Search ADS   Semelin B., Eames E., Bolgar F., Caillat M., 2017, MNRAS , 472, 4508 https://doi.org/10.1093/mnras/stx2274 CrossRef Search ADS   Shakura N. I., Sunyaev R. A., 1973, A&A , 24, 337 Singh S. et al.  , 2017, ApJ , 845, L12 https://doi.org/10.3847/2041-8213/aa831b CrossRef Search ADS   So G. C., Norman M. L., Reynolds D. R., Wise J. H., 2014, ApJ , 789, 149 https://doi.org/10.1088/0004-637X/789/2/149 CrossRef Search ADS   Springel V., 2005, MNRAS , 364, 1105 https://doi.org/10.1111/j.1365-2966.2005.09655.x CrossRef Search ADS   Springel V., Hernquist L., 2003, MNRAS , 339, 289 https://doi.org/10.1046/j.1365-8711.2003.06206.x CrossRef Search ADS   Springel V., Di Matteo T., Hernquist L., 2005, MNRAS , 361, 776 https://doi.org/10.1111/j.1365-2966.2005.09238.x CrossRef Search ADS   Strickland D. K., Stevens I. R., 2000, MNRAS , 314, 511 https://doi.org/10.1046/j.1365-8711.2000.03391.x CrossRef Search ADS   Suchkov A. A., Balsara D. S., Heckman T. M., Leitherner C., 1994, ApJ , 430, 511 https://doi.org/10.1086/174427 CrossRef Search ADS   Tange O., 2011, The USENIX Magazine , 36, 42 Van Der Walt S., Colbert S. C., Varoquaux G., 2011, Comput. Sci. Eng. , 13, 22 https://doi.org/10.1109/MCSE.2011.37 CrossRef Search ADS   Vanzella E. et al.  , 2016, ApJ , 825, 41 https://doi.org/10.3847/0004-637X/825/1/41 CrossRef Search ADS   Vanzella E. et al.  , 2018, MNRAS , in press Worseck G., Prochaska J. X., Hennawi J. F., McQuinn M., 2016, ApJ , 825, 144 https://doi.org/10.3847/0004-637X/825/2/144 CrossRef Search ADS   Yoshida N., Oh S. P., Kitayama T., Hernquist L., 2007, ApJ , 663, 687 https://doi.org/10.1086/518227 CrossRef Search ADS   Zheng Z.-Y. et al.  , 2017, ApJ , 842, L22 https://doi.org/10.3847/2041-8213/aa794f CrossRef Search ADS   APPENDIX A: CONVERGENCE To check the convergence of our results, we run two simulations including all source types that differ only in the number of photon packets emitted per source, i.e. Nγ = 104 and 105. In Fig. A1 we show the fraction of cells at z = 12 with a given T, $$x_{\rm H\,\small {II}}$$, $$x_{\rm He\,\small {II}}$$, and $$x_{\rm He\,\small {iii}}$$, together with the relative difference between the results of the two simulations, i.e. D/Dref − 1, where D = T, $$x_{\rm H\,\small {II}}$$, $$x_{\rm He\,\small {II}}$$, and $$x_{\rm He\,\small {iii}}$$ and Dref represents the values obtained with the highest value of Nγ. More than 50 per cent of the cells in the Nγ = 104 simulation are very well converged, with differences at the per cent level, except for $$x_{\rm H\,\small {II}} = x_{\rm He\,\small {II}} \sim 10^{-3}$$ where the results deviate up to ∼ 20 per cent. For the temperatures, the differences are less than 5 per cent for the 67 per cent of the cells that have 10 < T ≤ 104 K. Large differences are observed only in a handful of cells at high temperature or very low ionization fraction $$x_{\rm H\,\small {II}} = x_{\rm He\,\small {II}} < 10^{-5}$$, which we set as the lower limit in this work. Figure A1. View largeDownload slide Upper panels: from left to right, fraction of cells at z = 12 with a given $$x_{\rm H\,\small {II}}$$, $$x_{\rm He\,\small {II}}$$, $$x_{\rm He\,\small {iii}}$$, and T for simulations with Nγ = 104 (red dashed lines) and 105 (blue solid lines). Lower panels: relative difference in per cent with respect to the case with Nγ = 105. Figure A1. View largeDownload slide Upper panels: from left to right, fraction of cells at z = 12 with a given $$x_{\rm H\,\small {II}}$$, $$x_{\rm He\,\small {II}}$$, $$x_{\rm He\,\small {iii}}$$, and T for simulations with Nγ = 104 (red dashed lines) and 105 (blue solid lines). Lower panels: relative difference in per cent with respect to the case with Nγ = 105. APPENDIX B: LIGHT-CONES We construct light-cones inspired by the methods applied by Mellema et al. (2006), Datta et al. (2012), and Giri et al. (2018). Our approach though needs to take into account the non-periodic boundary conditions of our simulations. We create a path Δl(ti) ≡ Δli corresponding to the light travel distance cΔti = kiΔxi covered by a ray travelling through ki cells of the volume, each with sides of physical length,   $$\Delta x_i = \frac{1}{1+z_i} \frac{100\,h^{-1} }{256} \,{\rm pMpc},$$ (B1)from one crash output i to the next, i + 1, having a difference in cosmological age Δti. We hence assume the time steps to be sufficiently small to neglect the evolution between them when calculating the path Δli. We linearly interpolate between each output, i.e. the ki cells that are covered between two crash snapshots have contributions from both snapshots. The position along the ki cells that make up the path between two snapshots can be written as lp, where p = 0, …, ki − 1. For each point lp, which corresponds to a position in both physical (xp, yp) and redshift (zp) space, we also have a third spatial dimension, $$\boldsymbol {c}_p(l_p)$$, which is plotted as the y-axis in the light-cones. It has contributions from the two snapshots at times ti and ti + 1, weighed by the position (ki − 1 − p)/(ki − 1) and p/(ki − 1), respectively, along the path between them. For the ionization fraction we have   $$\boldsymbol {c}_p(l_p) = \frac{1}{k_i - 1}\left[\left(k_i - 1 -p\right) \boldsymbol {x}_{\rm H\,\small {II}}^{t_i}(l_p) + p \boldsymbol {x}_{\rm H\,\small {II}}^{t_{i+1}}(l_p) \right],$$ (B2)where we write the ionization fractions $$\boldsymbol {x}_{\rm H\,\small {II}}^{t_i}(l_p)$$ from the two contributing snapshots as vectors to indicate that we only use the cells along the third spatial dimension at this point in (zp, xp, yp) space. There are various possible approaches for deciding the paths Δli. Instead of choosing random paths, which could overlap, we instead try fixed patterns, and attempt at covering as much of the simulations volumes as possible in the light-cones. For this purpose, we choose a continuous whisk broom scanning approach (see e.g. Schaepman 2009) with decreasing horizontal scan width and a several pixels wide (slightly decreasing) vertical scan offset between each horizontal scan. Both are reset to maximum widths once we reach the vertical boundary. This way we prevent, as far as possible, the same objects at the same cosmological ages to appear in the light-cones. We only use the interior 236 cells (out of a maximum of 256) along each scanning axis to prevent boundary effects on the light-cones. The parameters of this approach (the number of pixels to decrease the scan width and offset, vertical scan offset) must however be tuned to obtain a satisfying result. The approach is thus optimal for a visualization, rather than for an objective quantification, of the temporal progression. © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society

### Journal

Monthly Notices of the Royal Astronomical SocietyOxford University Press

Published: May 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations