X-RAY BACKGROUND AT HIGH REDSHIFTS FROM POP III REMNANTS: RESULTS FROM POP III STAR FORMATION RATES IN THE RENAISSANCE SIMULATIONS

One of the most important problems in astrophysics and cosmology in both theory and observation is to obtain the thermal and ionization history of the intergalactic medium (IGM) prior, during, and after the epoch of reionization. It is critical to understand cosmological observations, such as high-redshift 21 cm neutral hydrogen line signatures, the Thomson scattering optical depth of the CMB, and properties of the Lyα forest.

The appearance of the first luminous objects—likely the first stars—marks the end of the cosmic dark ages after recombination, and the beginning of the last cosmic phase transition of reionization, a prolonged process that ionizes and heats the IGM. The first generation of stars (Population III stars) are believed to have formed from metal-free gas in small dark matter halos and having a large (tens to a few hundreds ${M}_{\odot }$) characteristic mass (Abel et al. 2002; Bromm et al. 2002; O'Shea & Norman 2007; Turk et al. 2009; Greif et al. 2012; Hirano et al. 2015; Hosokawa et al. 2016). Pop III stars have short lifetimes (Schaerer 2002) and may explode as supernovae (SNe), enriching their surrounding IGM. Once the gas metallicity passes some critical threshold, ∼10⁻⁶ ${Z}_{\odot }$ if dust cooling is efficient (Omukai et al. 2005; Schneider et al. 2006; Clark et al. 2008) or ∼${10}^{-3.5}$ Z${}_{\odot }$ otherwise (Bromm et al. 2001; Smith et al. 2009), the gas can cool efficiently and lower its Jeans mass and form (Pop II) stars in clusters at much higher rates. This makes the Pop III period of a given galaxy short. Consequently, it is believed that Pop III stars are not a major contributor to the UV photons responsible for reionization (e.g., Fan et al. 2006). (However, see Wise & Abel 2008 and Ahn et al. 2012 for their contribution to the start of reionization.)

X-ray radiation with its much longer mean free path than UV radiation has been considered a good candidate for the pre-ionization and pre-heating of the IGM much earlier than the end of the reionization at $z\sim 6$ (e.g., Oh 2001; Venkatesan et al. 2001; Ricotti & Ostriker 2004; Ricotti et al. 2005), and for smooth heating of the IGM (e.g., Haiman 2011). Its impact on the high-redshift 21 cm signal has been explored by many authors (Pritchard & Furlanetto 2007; Mirabel et al. 2011; Pritchard & Loeb 2012; Visbal et al. 2012; Fialkov & Barkana 2014; Pacucci et al. 2014; Mirocha et al. 2016), and observational constraints on the cosmic X-ray background have been investigated by Fialkov et al. (2016). Recent cosmological simulations (Turk et al. 2009; Stacy et al. 2010; Susa et al. 2014) have found that metal-free star-forming clouds might fragment to form binaries in a non-negligible fraction of Pop III star-forming events, making Pop III stars and remnants possible strong X-ray sources at high redshifts. Pop III in the approximate mass range $40\mbox{--}140\,{{\rm{M}}}_{\odot }$ and $\gt 260\,{{\rm{M}}}_{\odot }$ may directly collapse to form black holes (BHs; Heger et al. 2003). Any strong accretion onto these massive Pop III seed BHs would lead to X-ray radiation at high redshifts (Kuhlen & Madau 2005; Alvarez et al. 2009; Tanaka et al. 2012; Visbal et al. 2012). X-rays from Pop III binaries have been suggested to produce a pre-heated IGM (e.g., Haiman 2011; Mirabel et al. 2011; Visbal et al. 2012; Mesinger et al. 2013; Fialkov et al. 2014; Madau & Fragos 2016), and to partially ionize the IGM in large volumes (e.g., Ostriker & Gnedin 1996; Pritchard & Furlanetto 2007). Xu et al. (2014) showed that Pop III binaries might be the dominant X-ray sources at high redshifts and discussed how the IGM might be heated locally by hundred eV photons and globally by the X-ray background from keV photons. Homogeneously heating of the IGM on tens to hundreds megaparsec scales by an X-ray background may be important to understand the line widths of Lyα forest spectra (Tytler et al. 2009) and high-z 21 cm neutral hydrogen line signatures (Ahn et al. 2015).

Recently, Xu et al. (2016a) showed, using numerical simulations, that Pop III stars can form continuously at significant rates from redshifts higher than 20 down to the end of the reionization. This suggests that X-rays from Pop III binaries might continue to be important X-ray sources until the end of the reionization, and the X-ray background may heat the IGM for a long enough period to significantly heat the IGM globally. In this Letter, we estimate the X-ray luminosities from Pop III binaries using the results of the Renaissance Simulations (Xu et al. 2013, 2016b; O'Shea et al. 2015) and calculate the X-ray background using simple models. This Letter extends to lower redshifts the results first presented in Xu et al. (2014) concerning the X-ray background itself; we defer to a future paper a detailed calculation of the heating and ionization of the IGM and its 21 cm signature. We first describe the simulations and X-ray source and background models in Section 2. In Section 3, we present our results for the X-ray luminosity density, and the background intensity and spectra. We discuss the implications of this strong X-ray background from Pop III stars in Section 4.

The Renaissance Simulations are a suite of adaptive mesh refinement (AMR) first galaxy formation simulations carried out using the Enzo ⁵ code (Bryan et al. 2014) on the Blue Waters supercomputer operated by NCSA/UIUC. Enzo has been used extensively for high-redshift structure formation research (e.g., Abel et al. 2002; O'Shea & Norman 2007; Xu et al. 2008, 2013; Turk et al. 2009; Wise et al. 2012a, 2012b; Wise et al. 2014). Notably, UV ionizing radiation from stellar populations is followed using the Moray adaptive ray-tracing radiation transport solver (Wise & Abel 2011). Lyman–Werner radiation transport from stellar sources is calculated without considering absorption using ray tracing as well. However, we included the attenuation of LW photons, by a gray opacity approximation, the "picket-fence" modulation factor, as a function of the comoving distance between the source and the observer (Ahn et al. 2009). The baryonic properties are calculated using a nine-species primordial non-equilibrium gas chemistry and cooling network (Abel et al. 1997), supplemented by metal-dependent cooling tables (Smith et al. 2009). Prescriptions for Population III and metal-enriched star formation and feedback are employed, which have been fully described in Xu et al. (2013) and Wise et al. (2012b). Here, we briefly summarize the Pop III prescriptions because they are the focus of this Letter. Individual Pop III star particles are created when the overdensity, cooling time, H₂ fraction, and metallicity thresholds given in Wise et al. (2012a) are obeyed. The stellar mass, which dictates its lifetime and fate is drawn from a top-heavy IMF with a characteristic mass of $40{M}_{\odot }$. The UV hydrogen ionizing and LW photon luminosities and lifetimes of the Pop III stars are determined by the mass-dependent model from Schaerer (2002). Depending on their mass, they can end their lives as Type II or pair instability SNe, or collapse directly to BHs (Heger et al. 2003). Kinetic and metal yields for the SNe are taken from Heger et al. (2003). Pop III stellar remnants' positions are tracked for the duration of the simulation as they become incorporated into more massive halos. Their population statistics were examined in Xu et al. (2013).

The simulated box is a region of the universe with 28.4 comoving Mpc/h on a side (or 40 Mpc with h = 0.71). This is more than adequate to fairly sample the large fluctuations in the space density of Pop III star-forming minihalos (Barkana & Loeb 2004). We do not include the streaming velocity effect that modulates Pop III star formation on yet larger scales (Fialkov et al. 2013) because our Lyman–Werner regulated Pop III star formation occurs in minihalos with ${M}_{h}\gt 5\times {10}^{6}{M}_{\odot }$ (Xu et al. 2013, 2016a) where the streaming effect is negligible (Greif et al. 2011; Stacy et al. 2011). The simulations use the WMAP7 ΛCDM+SZ+LENS best-fit cosmology (Komatsu et al. 2011): ${{\rm{\Omega }}}_{M}=0.266$, ${{\rm{\Omega }}}_{{\rm{\Lambda }}}=0.734$, ${{\rm{\Omega }}}_{b}=0.0449$, h = 0.71, ${\sigma }_{8}=0.81$, and n = 0.963. Initial conditions were generated at z = 99 using MUSIC (Hahn & Abel 2011), and a low-resolution (512³ root grid) simulation was run to z = 6 to find regions suitable for re-simulation. The simulation volume was then smoothed on a physical scale of 6 comoving Mpc, and regions of high ($\langle \delta \rangle \equiv \langle \rho \rangle /({{\rm{\Omega }}}_{M}{\rho }_{C})-1\simeq 0.68$), average ($\langle \delta \rangle \simeq 0.09$), and low ($\langle \delta \rangle \simeq -0.26$) mean density (at $z=15,12.5,\mathrm{and}\,8$, respectively) were chosen for re-simulation. These subvolumes, hereafter designated the "Rarepeak," "Normal," and "Void" regions, with comoving volumes of 133.6, 220.5, and 220.5 Mpc³, were resimulated with three additional static nested grids, resulting in an effective initial resolution of 4096³ grid cells and particles in the region of interest, translating to a dark matter mass resolution of $2.9\times {10}^{4}$ ${M}_{\odot }$ in the same region. We allowed further local refinement (and then, the star formation) in the Lagrangian volume of the finest nested grid based on baryon or dark matter overdensity for up to 12 total levels of refinement (i.e., a comoving resolution of 19 pc in the finest cells). For more details about the calculations and scientific results about Renaissance Simulations, see Xu et al. (2013, 2014), Chen et al. (2014), and O'Shea et al. (2015).

The X-ray source and background models have been employed in Xu et al. (2014) and Ahn et al. (2015). Here, we summarize the basic assumptions and approach. The X-ray source model assumes that a fraction ${\epsilon }_{\mathrm{bin}}=0.5$ of Pop III stars form in binaries and evolve to become high=mass X-ray binaries (XRBs) that emit X-rays for a fixed lifetime τ = 30 Myr. We use the same luminosities as in Xu et al. (2014) to calculate the X-rays emitted in the AMR-refined volumes. To calculate the X-ray background, we need to estimate the X-ray emission coming from the non-refined regions, in which star formation is not resolved. We do this in two different ways. In Method 1, we adopt the same simple prescription described in Xu et al. (2014); namely, for every root grid cell in the refinement volume containing X-ray emitting sources, we sum up their baryon density ρ and X-ray luminosity L. Dividing by the number of emitting cells we have their average values $\bar{\rho }$ and $\bar{L}$. Then, in the unrefined portion of the simulation volume, if the root grid cell's baryon density $\rho \geqslant \bar{\rho }$, we assign it an X-ray luminosity $\bar{L}$. The defect of this model is that it predicts different amounts of total emission from the unresolved region, which is mostly common to each simulation. This is due to the different mean densities of the refined regions. To deal with this problem we introduce Method 2, in which we combine the results from all three refined regions to predict the X-ray emission at every redshift. For each root grid cell in the three survey volumes at each redshift, we record $L\geqslant 0,\rho$, and a smoothed density over a ∼ 1 Mpc sphere ${\rho }_{s}$ centered on the cell. From this data we apply machine-learning algorithms from the Python toolkit scikit-learn to derive two functions: ${P}_{X}(\rho ,{\rho }_{s},z)$ and ${L}_{X}(\rho ,{\rho }_{s},z)$ . P_X is the normalized probability that a cell with values $(\rho ,{\rho }_{s},z)$ contains one or more X-ray sources. We assume the probability can be represented as a second-order polynomial in $z,\rho$, and ${\rho }_{s}$ and fit it using all data from three simulations at all redshifts. L_X is a second-order polynomial fit to the dimensional X-ray luminosity for only those cells with $L\gt 0$. Armed with these two functions, we now walk through the cells in the unrefined region and draw a random number to compare with P_X to decide whether it contains X-ray sources. If it does, we assign it an X-ray luminosity using L_X. In this way, we generate the X-ray emissivity field in the unrefined root grid for all redshift outputs in each of the three simulations. As we will see in the next section, Method 1 provides only a rough estimate, for the reason stated, and is not recommended. Method 2 is more robust because it is based on fitting far more data from different redshifts and environments. It is nonetheless worthwhile presenting both as it provides a quantitative estimate for the error.

Given the X-ray emissivity everywhere in the 40 Mpc periodic volume, we perform the same calculations as in Xu et al. (2014, Section 5.2) to estimate the X-ray background. We locate the 3D field of X-ray luminosities frozen at the most recent past in the full box periodically, and calculate the contribution at every location of the box by summing over the full contribution from all the X-ray sources but within the corresponding lookback time. We simply treat the background sources to be uniformly distributed, since inhomogeneity in the source distribution at large lookback times will be observed to be almost uniform inside the box. Here, we assume that the whole simulation box is a good representation of the average universe and take its luminosity as the mean, globally averaged luminosity, and we also consider the absorption of photons by H i and He i with cosmic mean density.

Simulation data were processed and analyzed using yt ⁶ (Turk et al. 2011) on the Blue Waters system at NCSA, the Gordon system at SDSC, and the Wrangler system at TACC.

Before we show the results of X-ray generation and background, we illustrate the Pop III star formation rate densities from the three simulated volumes in Figure 1. The Pop III star formation rate density is simply correlated with the mean matter overdensity. The Rarepeak simulation has earlier Pop III star formation and higher formation rate. The Void region is the only simulation we can run to single-digit redshifts. Its Pop III star formation starts late at z ∼ 20, and the rate density increases to ∼ 10⁻⁵ ${M}_{\odot }$ yr⁻¹ Mpc⁻³ at z ∼ 13, and then gradually decreases below z ∼ 10. The Pop III star formation remains roughly constant for more than 400 Myr. The Pop III star formation of the Normal region also appears to plateau at z ∼ 13, albeit at a higher rate than the Void region. As discussed in Xu et al. (2013, 2016a), the Pop III formation is regulated through suppression of their formation in low-mass halos by the Lyman–Werner radiation from stellar feedback. However, their formation is not suppressed in halos of mass $\sim {10}^{7}{M}_{\odot }$ and above, despite strong photodissociating radiation. The plateaus reflect the rates at which pristine halos of sufficient mass are forming in these different environments.

Zoom In Zoom Out Reset image size

Figure 2 shows the evolution of the X-ray luminosity densities of the refined regions and the total X-ray luminosities of the whole box. The X-ray luminosity densities simply follow the Pop III star formation rate densities. Consequently, at z = 15, the cosmic variance effect on the predicted X-ray emissivity density is about an order of magnitude. Using Method 1, the total X-ray luminosities are underestimated for the Rarepeak, and overestimated for the Void run. To better estimate the X-ray background, we scale down the total luminosity from the Void region to match that from the Normal simulation. Using Method 2, we obtain a very smooth X-ray luminosity evolution. It grows from 1.7 × 10⁴³ to 3 × 10⁴⁴ erg s⁻¹ from z = 20 to z = 7.6, as the luminosity density grows from 2.7 × 10³⁸ to 4.7 ×10³⁹ erg s⁻¹ Mpc⁻³. This is a factor of 30(${\epsilon }_{\mathrm{bin}}/0.5)$ higher than estimates of the X-ray luminosities from normal XRBs and AGNs (Fragos et al. 2013; Madau & Fragos 2016). Taking the result from Method 2 as our reference, we would conclude that Method 1 is accurate to ±50% at the redshifts where they overlap.

Zoom In Zoom Out Reset image size

Mean background X-ray intensities are shown in Figure 3. We assume the X-ray sources are monochromatic and calculate the background for six different photon energies. The X-ray intensity at an observer redshift is the total intensity of observed redshifted X-rays from all earlier redshifts. In the upper panel, we show the X-ray background intensities from Void and Normal simulations using X-ray emissions calculated using Method 1. In the lower panel, we do the same thing using Method 2. In both cases the X-ray background gradually builds up, and is significantly dependent on the emitted photon energy. As discussed in Xu et al. (2014), only the high-energy (keV) photons accumulate to form a global background since sub-keV photons are absorbed near the source. The X-ray background intensity of 3 keV photons from Method 2, which we consider more reliable than Method 1, is over 10⁻⁶ erg s⁻¹ cm⁻² sr⁻¹.

Zoom In Zoom Out Reset image size

Figure 4 illustrates the X-ray background spectra for different source photon energies for different models at z = 12 and z = 7.6 showing how X-ray photons with different energies from higher redshifts contribute to the X-ray background. Figure 5 shows the synthetic X-ray spectra at z = 12 and z = 7.6 by assuming the source X-ray energies are equally distributed over the six energy bands. The high-energy keV photons with low absorption produce continuous broadband observed spectra by cosmological redshifting. Except for the gap between 1 and 1.5 keV for the z = 12 spectrum due to the redshift cutoff, the X-ray spectra have similar strengths at z = 12 and z = 7.6, as the Pop III star formation rate does not evolve much. At z = 7.6, keV X-rays from early times redshift to sub-keV energies and then heat and ionize the IGM. Using a one-zone calculation, Xu et al. (2014) showed that the X-rays of 1 and 3 keV in this spectrum have similar heating and ionizing effects to the IGM. The efficient absorption of X-rays happens below ∼500–600 eV, and any X-rays with photon energies lower than this are absorbed locally and only contribute a narrowband to the spectrum.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this Letter, we present estimates of the X-ray emission and background from Pop III binaries derived from the Pop III star formation histories in Renaissance Simulations beginning at $z\sim 25$ down to z = 7.6 near the end of reionization. The key results are summarized as follows:

• 1.  
  A significant amount of X-rays may be produced by Pop III binaries down to the end of reionization, as a direct result of Pop III stars forming continuously with steady rates since they begin forming at redshifts of 20–30.
• 2.  
  The X-ray background builds up gradually once the first stars begin forming. The X-ray background intensity depends on the assumed source photon energies. Only keV photons contribute to the X-ray background and become a broadband spectrum in the background.
• 3.  
  X-rays from Pop III binaries might continuously dominate X-rays from normal stellar XRBs and AGNs until the end of reionization, although this conclusion depends on the assumed fraction of Pop III stars that become X-ray binaries.

The results of this Letter may have important implications on the understanding of the reionization process, the properties of the post-reionization IGM, and many cosmological observations, like the Thomson scattering optical depth of the cosmic microwave background, neutral hydrogen 21 cm transition line, and the Lyα forest.

The Pop III binary X-ray luminosity per unit comoving volume is estimated to be $\mathrm{log}{I}_{{\rm{X}}}\simeq ({\epsilon }_{\mathrm{bin}}/0.5)[40.87-0.12(1+z)]$ in units of erg s⁻¹ Mpc⁻³, which is likely higher than those from other X-ray sources of AGNs and normal XRBs at redshifts higher than 6 (Fragos et al. 2013). This linear fit is only valid at $z\geqslant 10$ and has a maximum error of 0.5% compared to the computed values. This is based on our assumption that 50% of all Pop III stars form in binaries and evolve through a high-mass XRB phase. Since it takes a long time for the X-ray background to build up, the contribution from Pop III binaries may continue to dominate the X-ray background down to much lower redshifts.

Recently, Fialkov et al. (2016) used the unresolved soft X-ray background to place constraints on the nature of high-redshift X-ray sources. They derived constraints on the X-ray efficiency parameter f_x for three source types and two reionization scenarios, where f_x is defined implicitly by the relation (Furlanetto et al. 2006; Fragos et al. 2013)

$$\begin{eqnarray}&&\displaystyle \frac{{L}_{X}}{\mathrm{SFR}}=3\times {10}^{40}{f}_{X}\,\mathrm{erg}\,{{\rm{s}}}^{-1}{M}_{\odot }^{-1}\,\mathrm{yr},\end{eqnarray} \tag{ 1 }$$

where here SFR refers to the normal (Pop II) star formation rate density. Using the Method 2 value of ${L}_{X}\approx 5\times {10}^{39}\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-3}$ at z = 7.6, and the observed SFR of ${10}^{-2}{M}_{\odot }\,{\mathrm{yr}}^{-1}\,{\mathrm{Mpc}}^{-3}$ (Bouwens et al. 2015), we get ${f}_{x}\approx 16$, well within observational constraints for all cases considered.

It is interesting to compare the relative importance of ionizing UV and X-ray photons during the course of the simulations. To z = 7.6, both Pop III and metal-enriched stars produce about 1.4 × 10⁷¹ erg in ionizing UV photons (or about 6.3 × 10⁶⁹ photons) inside the refined region, of which less than 1% (10⁶⁹ erg) is produced by Pop III stars. There are about 7.6 × 10⁶⁸ hydrogen atoms in the survey volume of the Void run. At z = 7.6, only ∼ 13.2% hydrogen atoms have been ionized by the UV ionizing photons in the simulation (Xu et al. 2016b). Therefore, most of the ionizing UV photons are absorbed by hydrogen atoms that later recombine and are consumed to maintain the ionization state of high-density regions. On the contrary, a large fraction of X-rays are expected to penetrate to low-density regions that have recombination times longer than the Hubble time. During the same period, the Pop III binaries of the Void simulation produced ∼6.8 × 10⁶⁹ eV of energy in X-rays inside the refined region. It is about 5% of energy of ionizing UV photons and about 9 eV per hydrogen atom. Using Method 2, we estimate that Pop III binaries produce in total 2.1 × 10⁷² eV of energy in X-rays in the entire simulation box. As there are about 3.6 × 10⁷¹ hydrogen atoms in the simulation box, this amounts to about 6 eV of X-ray energy per hydrogen atom.

This photon budget analysis shows that the X-ray component is not negligible. There is enough energy in X-rays from Pop III binaries to significantly change the thermal and ionization states of the IGM. It is very hard to know the ratio of X-ray energy used in ionizing and heating the IGM without knowing the SED of the X-ray sources and doing a detailed simulation (Shull & van Steenberg 1985). Xu et al. (2014) used a one-zone model to calculate heating and ionizing effects for fixed X-ray photon energies and fluxes and to estimate the temperature and hydrogen ionization fraction of the Rarepeak simulation volume by simply assuming the X-ray background is constant since z = 15. This predicted that the IGM may be heated to several hundred Kelvin and be ionized to about $0.5 \%$ by z = 6. We have shown that the X-ray luminosity and background actually increase until the end of the simulation at z = 7.6 and are at much higher (∼10 times) levels. The electron fraction due to X-ray ionization might therefore reach a few percent by simply extrapolating the results in Xu et al. (2014), or even more since X-rays are also efficient at helium ionization. This physics is important for understanding the optical depth of CMB. In addition, by the end of reionization, the IGM may have been homogeneously heated by the keV photons of the X-ray background for hundreds of Myr. The additional heat input above and beyond the UV photoheating may reach thousands of Kelvin (Xu et al. 2014), and that is enough to have significant impacts on high-redshift 21 cm signatures (Ahn et al. 2015) and the Lyα forest spectra (Tytler et al. 2009) at lower redshifts.

This research is part of the Blue Waters project, which a joint effort supported by the NSF (award number ACI-1238993) and the state of Illinois, using NSF PRAC OCI-0832662. This research was supported by NSF grants PHY-0941373, AST-1109243, AST-1211626, and AST-1333360, and NASA grants NNX12AC98G, NNX15AP39G, HST-AR-13261.01-A, HST-AR-13895.001, HST-AR-14315.001-A, and HST-AR-14326.001. K.A. was supported by NRF-2012K1A3A7A03049606 and NRF-2014R1A1A2059811. This work was performed using the open-source Enzo and yt codes, which are the products of collaborative efforts of many independent scientists from institutions around the world. Their commitment to open science has helped make this work possible.