A Hybrid Approach Toward Simulating Reionization: Coupling Ray Tracing with Excursion Sets

As we know, the first luminous objects started to form around a redshift of 30 or so. One can perform n-body simulations to get the density distribution of dark matter at these redshifts (z ≲ 30). One can then use the excursion set formalism to identify the fraction of matter that has collapsed to form stars and galaxies from this density distribution. I have used the outputs of Illustris-2 simulations (Genel et al. 2014) to get the density distributions. The cosmological parameters used are wmap-9 (Hinshaw et al. 2013): Ω_m = 0.2726, Ω_b = 0.0456, σ₈ = 0.821, and h = 0.704.

The method used in this paper is intermediate between theoretical semianalytic methods, which are based on excursion set ideas (Furlanetto et al. 2004; Choudhury et al. 2009; Mesinger et al. 2011), and extensively numerical methods, which are based on ray-tracing algorithms (Abel et al. 1999; Razoumov & Scott 1999; Nakamoto et al. 2001). Our method is closest to the hybrid approach of Trac & Cen (2007), Lee et al. (2008), and Zheng et al. (2011). This earlier work focused more on small-scale structure and hence is computationally more intensive. In our approach, excursion set ideas are used to compute the effective number of ionizing photons produced by the collapsed objects. But, for computing the ionization status of the IGM, instead of excursion set ideas, ray-tracing is applied. The idea is to compute the ionizing photons produced in each grid cell and then propagate them in small steps. Obviously, it is not possible to propagate all the photons, but it is possible to propagate a large number of photon rays so as to get a stable brightness temperature distribution in the end. Each photon ray carries a weight that indicates the number of ionizing photons the ray carries. As a photon ray propagates through IGM it would ionize the neutral hydrogen and its weight would decrement depending upon the number of neutral hydrogen atoms it encounters. The method is based on the assumption that ionizing photons will always ionize a neutral atom almost immediately. This is alright as the cross-section for ionization is very high for the situations considered (x_{H i} ≫ 10⁻⁴).

The algorithm details are outlined in Figure 1

Zoom In Zoom Out Reset image size

After getting the overdensity distribution from N-body simulations, the first step is to compute the collapse fraction f_coll for each grid cell using excursion set formulations (Press & Schechter 1974; Bardeen et al. 1986; Bond et al. 1991; Lacey & Cole 1993; Liddle & Lyth 2000):

$$\begin{eqnarray}&&{f}_{\mathrm{coll}}=\mathrm{erfc}\left\{\displaystyle \frac{{\delta }_{c}(z)-{\delta }_{m}}{\sqrt{2\left[{\sigma }_{\min }^{2}-{\sigma }^{2}(m)\right]}}\right\}.\end{eqnarray} \tag{ 1 }$$

Here, δ_m is the mean linear overdensity of the grid cell, m is the total mass in the grid cell, σ²(m) is the variance of density fluctuations on the scale m, ${\sigma }_{\min }^{2}={\sigma }^{2}({M}_{\min })$, M_min is the minimum mass of ionizing sources and δ_c(z) is the critical density for collapse at redshift z. The grid size used is 64³ while the box size is (106.5 Mpc)³.

The second step is to compute the number of ionizing photons n_γ dumped in the IGM in a grid cell by the objects that are collapsed in that cell:

$$\begin{eqnarray}&&{n}_{\gamma }={n}_{\mathrm{ion}}\displaystyle \frac{{n}_{{\rm{H}}}}{1-{Y}_{\mathrm{He}}}{f}_{\mathrm{coll}}.\end{eqnarray} \tag{ 2 }$$

Here, n_ion is number of ionizing photons deposited in the IGM per collapsed baryon, n_H is hydrogen density and Y_He is helium fraction (Choudhury et al. 2009). Note that n_ion is closely related to the ζ parameter that is usually encountered in 21 cm literature (Furlanetto et al. 2004)

$$\begin{eqnarray}&&{n}_{\mathrm{ion}}={f}_{\mathrm{esc}}\cdot {f}_{* }\cdot {N}_{\gamma /b}.\end{eqnarray} \tag{ 3 }$$

Here, f_esc is the escape fraction of the ionizing photons, f_* is the star formation efficiency and N_γ/b is the number of ionizing photons produces per baryon in stars. Thus, n_ion is nothing but ζ without the ${\left(1+{n}_{\mathrm{rec}}\right)}^{-1}$ factor. If r is the ratio of recombination rate to ionization rate then one has $1-r\,={\left(1+{n}_{\mathrm{rec}}\right)}^{-1}$ or $r={n}_{\mathrm{rec}}/(1+{n}_{\mathrm{rec}})$. This factor is incorporated while propagation of photons. Thus, if in a certain step, n₁ number of hydrogen atoms are ionized in a grid cell then $r\cdot {n}_{1}$ are added back to the cell to account for recombination.

Next step is to create photon rays. About 10 million photon rays are created and a weight value is assigned to each ray based on the total number of photons that are created in the whole simulation box. We have weight = total photons/ray count. Each grid cell will have a number of photon rays created, given by rays = n_γ/weight. These rays are assigned initial positions and initial direction randomly inside the grid cell. The process is repeated for each grid cell.

Once photon rays are generated, they are propagated at each step by step size. The mean free path of individual photons in the IGM is very small. If the step size is assigned this value, then it would take a long time for the simulations to run. So a value of 0.05 of the grid cell size was assigned to it and stability of the simulations with respect to this step size value is considered.

After this, at each step, neutral hydrogen is ionized and photon rays are annihilated (or have their weight decremented). This is explained in detail in Figure 2. Position of each photon ray is tracked. Each photon ray carries an ID which tells you the grid cell that it has arrived into. Once this is done, photon rays with the same ID (implying that they are in the same grid cell) are considered together. Their weights are sorted and weights are decremented based on the number of neutral hydrogen atoms present in that grid cell. For example, if the cell is already completely ionized then all the photon rays that arrive in that cell will have their weights unchanged. On the other hand, if the total of all the photon ray weights that arrive in a grid cell is smaller than the total number of neutral hydrogen atoms present in the cell, then all the photon ray weights would go to zero and the neutral hydrogen count in that cell is decremented by the total of the photon ray weights. For in-between scenarios, the neutral hydrogen to be ionized is distributed equally between all the photon rays. If there are n photon rays arriving in the cell with the number of neutral hydrogen atoms N_{H i}, then the weight of each photon ray is decremented by N_{H i}/n. If for a given photon ray, its weight is less than this number then the additional load is shared equally between all the remaining photon rays. This process is well defined as we start sorting the photon ray weights and decrementing the weight from the ray with the lowest weight.

Zoom In Zoom Out Reset image size

The procedure is repeated until all the photon rays are completely annihilated. Practically, for a quick termination of the algorithm, the propagation is stopped once the total of all the photon ray weights falls below 0.01% of its initial value.

At the end of propagation one gets neutral hydrogen distribution that is left in each grid cell. This is used to compute the neutral hydrogen fraction ${x}_{{\rm{H}}{\rm{I}}}$ in that cell, and then the 21 cm brightness temperature in milli-kelvin units is computed using the well-known expression (Wouthuysen 1952; Field 1959; Madau et al. 1997; Furlanetto et al. 2006)

$$\begin{eqnarray}&&\delta {T}_{b}\approx 27{x}_{{\rm{H}}{\rm{I}}}(1+\delta )(1-{T}_{\gamma }/{T}_{{\rm{s}}}){\left(\displaystyle \frac{1+z}{10}\displaystyle \frac{0.15}{{{\rm{\Omega }}}_{M}{h}^{2}}\right)}^{1/2}\left(\displaystyle \frac{{{\rm{\Omega }}}_{b}{h}^{2}}{0.023}\right).\end{eqnarray} \tag{ 4 }$$

While computing the brightness temperature we assume T_γ ≪ T_s, that is, the cosmic microwave background temperature is much less than the spin temperature (see, e.g., Ghara et al. 2015).

Figure 3 plots 21 cm maps (δT_b) along with density field and initial photon field (${N}_{\gamma }^{0}$). As seen from Figure 3, the density field, which is proportional to ($1+\delta$), is similar due to the difference in the redshift being small. But, the collapse fraction, $\left\langle {f}_{\mathrm{coll}}\right\rangle$, is fairly different giving different ${N}_{\gamma }^{0}$ and δT_b distributions. Also, due to the presence of more neutral hydrogen, the redshift 9 maps have the maximum structure.

Zoom In Zoom Out Reset image size

After getting the brightness temperature distribution, the next step is to get the power spectrum. The power spectrum obtained for the fiducial model is depicted in Figure 4. To get the power spectrum, the 21 cm brightness temperature distribution, $\delta {T}_{b}\,\equiv {\delta }_{21}$, is Fourier transformed and then the dimensionless power spectrum is obtained using the following prescription:

$$\begin{eqnarray}&&\left\langle {\tilde{\delta }}_{21}({\boldsymbol{k}})\,{\tilde{\delta }}_{21}^{* }({\boldsymbol{k}}^{\prime} )\right\rangle ={\left(2\pi \right)}^{3}\,{\delta }_{D}({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} )\,P({\boldsymbol{k}}),\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}^{2}({\boldsymbol{k}})=\displaystyle \frac{{k}^{3}P({\boldsymbol{k}})}{2{\pi }^{2}}.\end{eqnarray} \tag{ 6 }$$

The spherically averaged power spectrum is obtained by averaging Δ²(k) over all possible angles:

$$\begin{eqnarray}&&{{\rm{\Delta }}}^{2}(k)=\displaystyle \frac{1}{2}{\int }_{-1}^{1}d\mu \,{{\rm{\Delta }}}^{2}({\boldsymbol{k}})={\int }_{0}^{1}d\mu \,{{\rm{\Delta }}}^{2}({\boldsymbol{k}}),\end{eqnarray} \tag{ 7 }$$

where $\mu ={k}_{\parallel }/k$ is the cosine of the angle that wavevector k makes with the line-of-sight direction.

Zoom In Zoom Out Reset image size

The power spectrum obtained is of expected magnitude (see, e.g., Mesinger et al. 2011) and also seems to show the expected trend with respect to the redshift.

This method of generating 21 cm maps has the following features:

• 1.  
  A simple and intuitive method where excursion set formalism is used only for getting the collapse fraction.
• 2.  
  An explicitly photon conserving method where one can compute the mean neutral hydrogen fraction once one knows, the mean collapse fraction ($\left\langle {f}_{\mathrm{coll}}\right\rangle$) and the other two reionization parameters (n_ion and n_rec).
• 3.  
  The expression for mean neutral hydrogen fraction is: ${\bar{x}}_{{\rm{H}}{\rm{I}}}=1-\tfrac{{n}_{\mathrm{ion}}\left\langle {f}_{\mathrm{coll}}\right\rangle }{(1+{n}_{\mathrm{rec}})(1-\left\langle {f}_{\mathrm{coll}}\right\rangle )(1-{Y}_{\mathrm{He}}/2)}$.
• 4.  
  Method is fairly fast. It takes about 3–4 minutes to propagate 10 million photon rays and generate a temperature cube (64 × 64 × 64) on an ordinary intel core i5 machine with 8 GB memory.
• 5.  
  It is precise even with a small number of photon rays as it automatically assigns more rays to more overdense regions or more luminous objects.

To check if the reionization procedure is giving a convergent answer, two tests were performed. In the first test, the step size was halved, twice, from the value chosen earlier, that is, it was made to be 0.025 and 0.0125 of the grid cell size. And, in the second test, the number of photon rays was doubled, again twice, to 20 million and 40 million. The effect of these two procedures on temperature maps and power spectrum is summarized in Tables 1 and 2 respectively. In the case of temperature maps the quantity listed is $\left\langle | {\rm{\Delta }}\delta {T}_{b}| \right\rangle /\left\langle \delta {T}_{b}\right\rangle$, where $| {\rm{\Delta }}\delta {T}_{b}|$ is the absolute value of the difference in brightness temperature measured with respect to the fiducial case, that is, with respect to step size of .05 and ray count of 10 million case. In this case, $\left\langle \right\rangle$ indicates the average over all the grid cells.

As seen from Table 1, the change in δT_b averaged over the whole grid cell is small (≲5%). Also the increment in going from the double to the four times case is smaller than going from the fiducial to the double case. In the case of changing step size, the percentage change is nearly the same. In the case of the power spectrum, in Table 2, the quantity listed is $\left\langle {\rm{\Delta }}{{\rm{\Delta }}}^{2}(k)/{{\rm{\Delta }}}^{2}(k)\right\rangle$, that is a fractional change in 21 cm power spectrum averaged over all k values. The smallness of values in the second table indicates that the procedure of getting the power spectrum has nearly converged. Again the increment in going from the double (half) to the four times (one-fourth) case is smaller than going from the fiducial to the double (half) case. To see the effect over different scales, the change in power spectrum for different scales is plotted in Figure 5. As seen from the figure, the change is very small (≈1 or 2%) for large scales and reasonably small (≈5%–10%) for small scales.

Zoom In Zoom Out Reset image size

Figure 6 shows the probability distribution of δT_b obtained for the three redshifts. To check whether the temperature distribution is converged or not, another number is considered. This number is half the area in between the two distributions of δT_b, one obtained in the fiducial case (step size of 0.05 and ray count of 10 million) and the other obtained in the changed parameter case. The results are depicted in Table 3. For computing the area, the total range of δT_b values was divided into about 100 bins. As seen from the table the difference is smaller than or of the order of 1% (for completely different distributions this number would be 1). Again the increment in going from the double (half) to the four times (one-fourth) case is smaller than going from the fiducial to the double (half) case. This indicates that the step size of .05 and the ray count of 10 million is reasonably adequate for the simulations.

Zoom In Zoom Out Reset image size

The parameters of this model are n_ion, n_rec, and M_min. The third parameter, M_min, directly impacts the collapse fraction and hence would have to be considered immediately after getting the overdensity from n-body simulations. The fiducial values of these parameters were 6.0, 1.0, and 10⁸M_☉ respectively. M_min ∼ 10⁸ M_☉ is a choice that is prescribed by atomic line cooling (Barkana & Loeb 2001). n_rec ∼ 1 is motivated by semianalytic models (Sobacchi & Mesinger 2014). If one chooses n_ion = 6 the mean neutral hydrogen fraction for this choice of parameters is 0.37, 0.57, and 0.7 for redshifts 7, 8, and 9 respectively. This set of values is close to (although not exactly the same as) the values that are consistent with CMBR and Lyα data (Raut & Choudhury 2018). One can ask the question, how would the power spectrum change if any one of these values is perturbed? The results are depicted in Figure 7. As seen from Figure 7, the power spectrum on small scales increases if one increases n_rec or decreases n_ion or increases M_min. This is understandable as with these changes one is effectively producing a lower number of photons and hence available neutral hydrogen is more and hence the signal is more. Such a behavior has been observed before by many authors (see, e.g., Lidz et al. 2008; Mesinger et al. 2011; Greig & Mesinger 2015; Cohen et al. 2018). On the largest scales the behavior needs to be reexamined using larger data cubes so as to suppress cosmic variance. But, it is arguable whether these modes can be measured accurately by future telescopes like SKA1-Low.¹

Zoom In Zoom Out Reset image size