AMBER: A Semi-numerical Abundance Matching Box for the Epoch of Reionization

The reionization-redshift fields are quantitatively very similar when we account for the different reionization histories. To show this, we take the ${z}_{\mathrm{re}}({\boldsymbol{x}})$ fields from Sims 0 and 2 and individually scale the redshift values for each grid cell without changing their spatial rank order such that they have the same ionization fraction ${\bar{x}}_{\mathrm{0,2}}(z)={\bar{x}}_{1}(z)$ as Sim 1. We are matching the abundance of ionized mass or volume at all redshifts, which we call abundance matching the reionization history (Section 3.4.2 for more details). Alvarez & Abel (2012) use a similar technique to rescale the reionization-redshift field in their semi-numerical method. Here we consider both mass-weighted ${\bar{x}}_{{\rm{i}},{\rm{M}}}(z)$ and volume-weighted ${\bar{x}}_{{\rm{i}},{\rm{V}}}(z)$ ionization fractions.

To compare the ${z}_{\mathrm{re}}({\boldsymbol{x}})$ fields cell by cell, we define the fractional difference field,

$$\begin{eqnarray}&&{\epsilon }_{{\rm{z}}}({\boldsymbol{x}})\equiv \displaystyle \frac{[1+{z}_{\mathrm{re}}({\boldsymbol{x}})]-[1+{z}_{1}({\boldsymbol{x}})]}{1+{z}_{1}({\boldsymbol{x}})},\end{eqnarray} \tag{ 3 }$$

for Sim 0 or 2 relative to Sim 1 and then calculate the probability distribution function (PDF) of _z. Figure 2 shows the distribution of fractional differences for the ${z}_{\mathrm{re}}({\boldsymbol{x}})$ fields at the RT resolution l_rt = 98 h⁻¹kpc, which is approximately ten times smaller than the typical mesh cell size used in semi-numerical models. The difference distributions for the original Sims 0 and 2 are skewed because the reionization histories differ from that of Sim 1. After abundance matching, the distributions are approximately Gaussian with zero mean and small root-mean-square differences ${\sigma }_{{\rm{z}}}=\langle {\left({\epsilon }_{{\rm{z}}}({\boldsymbol{x}})-{\bar{\epsilon }}_{{\rm{z}}}\right)}^{2}{\rangle }^{1/2}\lesssim 0.01$. We find very similar results when matching either the mass-weighted or the volume-weighted ionization fractions.

Zoom In Zoom Out Reset image size

To correlate the fields, we define the normalized redshift field,

$$\begin{eqnarray}&&{\delta }_{{\rm{z}}}({\boldsymbol{x}})\equiv \displaystyle \frac{[1+{z}_{\mathrm{re}}({\boldsymbol{x}})]-[1+{\bar{z}}_{\mathrm{re}}]}{1+{\bar{z}}_{\mathrm{re}}},\end{eqnarray} \tag{ 4 }$$

where ${\bar{z}}_{\mathrm{re}}$ is the volume-weighted average value, and Fourier transform to get δ_z( k ) for all three simulations. Generally, the two-point power spectrum and its dimensionless version are defined as

$$\begin{eqnarray}&&{P}_{{ij}}(k)\equiv \langle {\delta }_{i}({\boldsymbol{k}}){\delta }_{j}({\boldsymbol{k}})\rangle ,\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{ij}}^{2}(k)\equiv \displaystyle \frac{{k}^{3}}{2{\pi }^{2}}{P}_{{ij}}(k),\end{eqnarray} \tag{ 6 }$$

where i = j for auto-correlations and i ≠ j for cross correlations. For comparing two different fields, we quantify with the bias and cross correlation,

$$\begin{eqnarray}&&{b}_{{ij}}(k)\equiv \sqrt{\displaystyle \frac{{P}_{{ii}}(k)}{{P}_{{jj}}(k)}},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{r}_{{ij}}(k)\equiv \displaystyle \frac{{P}_{{ij}}(k)}{\sqrt{{P}_{{ii}}(k){P}_{{jj}}(k)}}.\end{eqnarray} \tag{ 8 }$$

The field δ_i is said to be biased, unbiased, or underbiased relative to δ_j for b_ij > 1, = 1, and < 1, respectively. Similarly, the fields and their fluctuations are said to be correlated, uncorrelated, and anticorrelated for r_ij > 0, = 0, and <0, with perfect correlation or anticorrelation given by the limits r_ij = 1 and = −1, respectively.

Figure 3 shows the auto-power spectra for the ${z}_{\mathrm{re}}({\boldsymbol{x}})$ fields prior to abundance matching. On larger scales, the dimensionless spectra similarly peak at k ≈ 0.3–0.4 h Mpc⁻¹. As this is only a factor of 3 from the RadHydro simulation box limit ${k}_{\min }=2\pi /L$, the peak scale may possibly be underestimated. Nonetheless, this suggests that reionization simulations and semi-numerical models must have box lengths >20 h⁻¹Mpc in order to capture the relevant large-scale fluctuations and correlations in the reionization-redshift field. The spectra similarly decline in power-law form on intermediate scales until k ≈ 10 h Mpc⁻¹, followed by steeper declines on smaller scales down to the RT Nyquist frequency.

Zoom In Zoom Out Reset image size

We also cross-correlate Sims 0 and 2 relative to Sim 1 and plot the bias b_zz and normalized cross correlation r_zz in Figure 3. Using the normalized fields δ_z( x ) partially adjusts for the different redshift midpoints, but the shorter and longer durations result in lower and higher bias for the original Sims 0 and 2, respectively. After abundance matching, the b_zz(k) curves are unity at the box scale and deviate by only ≲10% at the RT-resolution scale. The cross correlations show that the original fields are already very highly correlated, allowing abundance matching to be applied correctly. The r_zz(k) curves are exactly unity at the largest scale and remain within a few percent on scales k ≲ 10 h Mpc⁻¹.

All of these results lead to new ideas for AMBER: there is a spatial order to the reionization process, and abundance matching can be applied to a correlated field to accurately predict the reionization-redshift field.

The redshift midpoint z_mid, duration Δ_z, and asymmetry A_z are input parameters, which in turn give three ionization points. For the early, middle, and late stages of reionization, let the redshifts z_ear > z_mid > z_lat correspond to the ionization fractions x_ear < x_mid < x_lat. We take z_mid as the redshift midpoint and define the duration,

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{\rm{z}}}\equiv {z}_{\mathrm{ear}}-{z}_{\mathrm{lat}},\end{eqnarray} \tag{ 10 }$$

and asymmetry,

$$\begin{eqnarray}&&{A}_{{\rm{z}}}\equiv \displaystyle \frac{{z}_{\mathrm{ear}}-{z}_{\mathrm{mid}}}{{z}_{\mathrm{mid}}-{z}_{\mathrm{lat}}}.\end{eqnarray} \tag{ 11 }$$

The other two redshifts are then given by

$$\begin{eqnarray}&&{z}_{\mathrm{ear}}={z}_{\mathrm{mid}}+\displaystyle \frac{{{\rm{\Delta }}}_{{\rm{z}}}{A}_{{\rm{z}}}}{1+{A}_{{\rm{z}}}},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{z}_{\mathrm{lat}}={z}_{\mathrm{mid}}-\displaystyle \frac{{{\rm{\Delta }}}_{{\rm{z}}}{A}_{{\rm{z}}}}{1+{A}_{{\rm{z}}}}={z}_{\mathrm{ear}}-{{\rm{\Delta }}}_{{\rm{z}}}.\end{eqnarray} \tag{ 13 }$$

Instantaneous reionization models would correspond to Δ_z = 0, but reionization is an extended process with finite Δ_z > 0. Symmetric reionization histories would correspond to A_z = 1, but reionization simulations typically find that the early stage of reionization takes longer than the later stage such that A_z > 1.

Table 1 lists the midpoint, duration, and asymmetry parameters for the three RadHydro simulations. In Doussot et al. (2019), we take ${\bar{x}}_{\mathrm{mid}}=0.50$ and present two practical choices for the other ionization points. In the first case, we choose quartile ionization fractions $({\bar{x}}_{\mathrm{ear}},{\bar{x}}_{\mathrm{lat}})=(0.25,0.75)$, and Δ_z is analogous to a full width half max. In the second case, we take $({\bar{x}}_{\mathrm{ear}},{\bar{x}}_{\mathrm{lat}})=(0.05,0.95)$, and Δ_z more effectively quantifies the whole EoR. We will adopt the latter as the default convention throughout this paper. AMBER users can specify other values, but extreme choices (e.g., ${\bar{x}}_{\mathrm{ear}}\lesssim 0.01$, ${\bar{x}}_{\mathrm{lat}}\gtrsim 0.99$) are not recommended because the start and end of the EoR are not well defined and are difficult to determine precisely.

Once the midpoint, duration, and asymmetry parameters are specified and used to calculate three ionization points, we use an analytical interpolating function to calculate the reionization history. In Trac (2018), we use Lagrange interpolating functions to construct an analytical function for the ionization fraction x_i(z) that exactly fits the ionization points. For three points, the interpolating polynomial is a quadratic, which has the advantage of being invertible but has the disadvantage of being nonmonotonic. A log transformation of the variables improves monotonicity over the relevant redshift range, but it is not an ideal solution.

In AMBER, we interpolate the three ionization points with a modified Weibull function (Weibull 1951),

$$\begin{eqnarray}&&{\bar{x}}_{{\rm{i}}}(z)=\exp \left[-\max {\left(\displaystyle \frac{z-{a}_{{\rm{w}}}}{{b}_{{\rm{w}}}},0\right)}^{{c}_{{\rm{w}}}}\right],\end{eqnarray} \tag{ 14 }$$

where the coefficients a_w, b_w, c_w are all positive values. Note that a_w corresponds to z_end when x_end = 1 and the max function ensures full ionization at lower redshifts. In Appendix A, we show that the coefficients can be easily determined by first solving a nonlinear equation for c_w and then substituting its value into algebraic equations for the other two coefficients. We find that solutions exist for the asymmetry range A_z ≲ 15, which is more than sufficient for parameter-space studies. We also experimented with a generalized logistic function (Richards 1959), but find that it is limited to A_z ≲ 5. The Weibull function can be analytically inverted,

$$\begin{eqnarray}&&z({\bar{x}}_{{\rm{i}}})={a}_{{\rm{w}}}+{b}_{{\rm{w}}}{\left|\mathrm{ln}{\bar{x}}_{{\rm{i}}}\right|}^{1/{c}_{{\rm{w}}}},\end{eqnarray} \tag{ 15 }$$

which is a useful feature for our abundance matching technique.

Table 1 lists the Weibull coefficients for the three RadHydro simulations. The offset a_w determines when the function reaches unity, and it decreases when reionization ends later. The divisor b_w controls the stretch in redshift, and it increases with the duration Δ_z. The power c_w controls the steepness of the function, and it is inversely related to the asymmetry A_z. In general, changing one of the reionization shape parameters while holding the other two fixed affects all three Weibull coefficients.

Figure 4 compares the evolution of the mass-weighted ionization fraction x_i(z) from the RadHydro simulations with the Weibull functions. The analytical parameterizations are excellent matches to the simulated reionization histories. The typical differences are only $| {\rm{\Delta }}{\bar{x}}_{{\rm{i}}}| \lesssim 0.01$, while the maximum differences of $| {\rm{\Delta }}{\bar{x}}_{{\rm{i}}}| \lesssim 0.02$ are found near the end of the EoR, which is not accurately captured in reionization simulations and semi-analytical models. We find that the Weibull function gives slightly better fits than Lagrange interpolation and is valid for a larger range of parameter space.

Zoom In Zoom Out Reset image size

Density fields are constructed from the LPT particles using techniques from particle-mesh (PM) methods. In N-body simulations, particles are added to a Cartesian mesh to create a density field for solving Poisson's equation with Fast Fourier Transforms (FFT). In large-scale structure analyses, galaxies or halos are added to a grid to create a density field for calculating the power spectrum. Similar techniques are also used in other semi-numerical models of reionization, but we will describe an improved implementation for AMBER.

There is a hierarchy of particle assignment schemes, and the first three are called the nearest grid point (NGP), cloud-in-cell (CIC), and triangular-shaped cells (TSC). The mass assignment process interpolates the discrete distribution of particles and can be interpreted as convolving and sampling the underlying density field ρ( x ) to produce a discretized density field,

$$\begin{eqnarray}&&\rho ({{\boldsymbol{x}}}_{{\rm{n}}})=\int \rho ({\boldsymbol{x}}){W}_{\mathrm{pm}}({{\boldsymbol{x}}}_{{\rm{n}}}-{\boldsymbol{x}}){d}^{3}x,\end{eqnarray} \tag{ 20 }$$

where W_pm( x ) is the PM assignment window function. See Hockney & Eastwood (1988) for a seminal review.

The particle assignments introduce effects such as aliasing, smoothing, and shot noise that have to be accounted for (e.g., Jing 2005). While the latter two effects can have analytical corrections, the first is especially problematic. Hockney & Eastwood (1988) have suggested that interlacing two or more staggered meshes can reduce aliasing. Interlacing has been used to improve force calculations in N-body simulations (Couchman 1991), produce more accurate overdensity fields for power spectra estimation (Sefusatti et al. 2016), and recently implemented in the nbodykit package (Hand et al. 2018).

We use the interlacing technique to produce more accurate and robust density fields in AMBER. We construct the first density field ρ₁( x ) as normal; while for the second density mesh ρ₂( x ), the particle positions are shifted by half of a grid cell spacing in each direction by adding the vector Δ x = (l_c/2, l_c/2, l_c/2), which is equivalent to shifting the mesh by the opposite vector. In Fourier space, the two transformed fields are combined into a single, effective field as

$$\begin{eqnarray}&&\rho ({\boldsymbol{k}})=\displaystyle \frac{\left[{\rho }_{1}({\boldsymbol{k}})+{e}^{i{\boldsymbol{k}}\cdot {\boldsymbol{\Delta }}{\boldsymbol{x}}}{\rho }_{2}({\boldsymbol{k}})\right]/2}{{W}_{\mathrm{pm}}({\boldsymbol{k}})}.\end{eqnarray} \tag{ 21 }$$

The twiddle factor e^{i k ·Δ x} multiplied to the second transformed field accounts for the shifted particle positions. We also deconvolve the effects of smoothing by dividing with the Fourier transform of the assignment window function,

$$\begin{eqnarray}&&{W}_{\mathrm{pm}}({\boldsymbol{k}})={\left[\displaystyle \frac{\sin ({k}_{{\rm{x}}}{l}_{{\rm{c}}}/2)\sin ({k}_{{\rm{y}}}{l}_{{\rm{c}}}/2)\sin ({k}_{{\rm{z}}}{l}_{{\rm{c}}}/2)}{({k}_{{\rm{x}}}{l}_{{\rm{c}}}/2)({k}_{{\rm{y}}}{l}_{{\rm{c}}}/2)({k}_{{\rm{z}}}{l}_{{\rm{c}}}/2)}\right]}^{p},\end{eqnarray} \tag{ 22 }$$

where p = 1, 2, and 3 for the NGP, CIC, and TSC schemes, respectively.

For comparison, we construct the AMBER and RadHydro density fields at two resolutions, using a 64³ mesh with cell size l_c = 0.8 h⁻¹Mpc and a 512³ mesh with cell size l_c = 0.1 h⁻¹Mpc. The fiducial lower resolution is more typical of semi-numerical models, while the higher resolution allows us to test the limits of LPT. The AMBER density fields are efficiently made using equal numbers of LPT particles as mesh cells. The RadHydro density fields are made by simply binning the 2048³ dark matter particles and the 2048³ gas cells down to the appropriate size meshes.

Figure 5 is a visualization of the matter density fields ρ( x ) at redshift z = 8 from RadHydro and AMBER. The LPT particles are assigned using the TSC scheme with interlacing and deconvolution. Both images are shown with the higher-resolution pixels in order to see the large-scale structure more clearly, but have been averaged along the line of sight at the lower resolution. The RadHydro matter and gas (not shown) density fields are visually indistinguishable, while the AMBER density field also has close resemblance, but minor smoothing is visible in the small-scale, high-density regions. At the lower resolution, the AMBER and RadHydro fields are indistinguishable, though the images appear quite pixelated.

Zoom In Zoom Out Reset image size

Figure 6 compares the one-point PDF of the relative density ρ/〈ρ〉 = 1 + δ at redshift z = 8. The RadHydro matter and gas distributions are in very good agreement even at the higher resolution. For AMBER, the TSC assignment scheme with interlacing and deconvolution gives the best agreement, and we only show this case to avoid overcrowding the plot. The results are also in very good agreement at the fiducial resolution, but not at the higher resolution. LPT calculated at the second order or an even higher order cannot capture the nonlinear shell-crossing on smaller scales and therefore underresolves high-density, collapsed regions. The disagreement in underdense regions is not concerning, because it is due to the much smaller number of LPT particles and the different assignment and deconvolution process for AMBER. We would find the same effects for RadHydro if we use a lower-resolution simulation and had not done the simple binning.

Zoom In Zoom Out Reset image size

Figure 7 shows the density auto-power spectra ${{\rm{\Delta }}}_{\mathrm{dd}}^{2}(k)$ at redshift z = 8 and cross correlations relative to the RadHydro matter overdensity δ_m. All of the dimensionless power spectra continue to increase in amplitude toward smaller scales, typical for cold dark-matter-dominated models. For the fiducial lower resolution, the power reaches a maximum value of only ${{\rm{\Delta }}}_{\mathrm{dd}}^{2}\approx 0.2$. As this is still in the quasi-linear regime, 2LPT is an accurate and efficient approach for modeling the mass distribution and density field.

Zoom In Zoom Out Reset image size

The RadHydro matter and gas perfectly trace each other with bias b_dd = 1 and cross correlation r_dd = 1 for k ≲ 10 h Mpc⁻¹. The suppression in the gas bias on smaller scales is due to the Jeans smoothing from the gas pressure that opposes gravitational collapse. For AMBER, the LPT bias starts out at unity on the largest scales, is still within a few percent for k ≳ 1 h Mpc⁻¹, but drops to ∼0.9 by k ∼ 10 h Mpc⁻¹. However, the normalized cross correlation is still closer to unity down to smaller scales, showing that the density perturbations are highly in-phased even though the amplitudes may differ. We find that using the TSC assignment scheme with interlacing and deconvolution gives more accurate auto-power and cross-power spectra. The NGP and CIC results tend to be overbiased and undercorrelated near the Nyquist frequency k_Ny = π/Δx because of aliasing. Our results are similar to the previous findings by Sefusatti et al. (2016), who present a thorough analysis of how interlacing can reduce aliasing and improve power spectrum estimation.

Velocity fields are constructed similarly to the density fields by assigning the LPT particles to interlaced meshes. For a given mesh cell, the velocity v is a weighted average of the individual velocities v _i from all overlapping particles,

$$\begin{eqnarray}&&{\boldsymbol{v}}=\sum _{i}\displaystyle \frac{{w}_{i}}{{w}_{\mathrm{tot}}}{{\boldsymbol{v}}}_{i}=\displaystyle \frac{\sum {w}_{i}{{\boldsymbol{v}}}_{i}}{\sum {w}_{i}},\end{eqnarray} \tag{ 23 }$$

where w_i and w_tot are the individual and total mass assignment weights, respectively. With the NGP assignment scheme, we find that some cells can have no overlapping particles and therefore undefined velocities. Both the CIC and TSC schemes work better in practice with no empty cells and missing velocities for our moderately clustered distribution of LPT particles.

Figure 8 is a visualization of the line-of-sight component of the velocity fields v( x ) at redshift z = 8 from RadHydro and AMBER. The velocity images are constructed using the same approach as for the density images (Figure 5). The shown images with the higher-resolution pixels are visually similar and have an even closer resemblance with lower-resolution pixels. The RadHydro and AMBER velocity fields show large coherence lengths that are actually underestimated because of the moderate simulation box length L = 50 h⁻¹Mpc. We will discuss this interesting feature in more detail below.

Zoom In Zoom Out Reset image size

Figure 9 compares the one-point PDF of the velocity components v ⊂ (v_x, v_y, v_z), at redshift z = 8. The RadHydro matter and gas velocity distributions are in excellent agreement, even more so than their density distributions. All of the distributions are Gaussian, and the widths decrease slightly by a few percent for the fiducial lower resolution. For AMBER, the CIC and TSC assignment schemes with interlacing and deconvolution give similarly accurate results, and we only show the latter to avoid overcrowding the plot. The differences in the tails of the PDFs are due to the RadHydro simulations having slightly faster velocities because of virialization in rare, collapsed regions.

Zoom In Zoom Out Reset image size

Figure 10 shows the velocity auto-power spectra Δ_vv(k) and the cross correlations relative to the RadHydro matter velocity v_m. We individually compare and average over each velocity component v ⊂ (v_x, v_y, v_z). The velocity power spectrum declines in amplitude toward smaller scales, unlike the density power spectrum (Figure 7). In linear theory, the transformed velocity field is related to the matter overdensity field as v ( k ) ∝ ( k /k²)δ( k ), and therefore the power is predominantly coming from larger scales. This explains the large coherence length seen in the velocity images (Figure 8). While radiation-hydrodynamic simulations are still computationally too expensive, semi-numerical methods are sufficiently efficient to afford large box sizes of ≳ 500 h⁻¹Mpc, which is necessary to capture the large-scale power and accurately model the velocity fields.

Zoom In Zoom Out Reset image size

The RadHydro matter and gas velocities perfectly trace each other with bias b_vv = 1 and cross correlation r_vv = 1 for k ≲ 20 h Mpc⁻¹, even to smaller scales compared to their densities (Figure 7). For AMBER, the NGP assignment scheme gives incorrect bias and poor stochasticity when we simply set the velocities to zero where ever it is undefined. The CIC and TSC assignment schemes have similarly accurate bias, but the latter has preferably better cross correlation. For the velocity bias, the lower-resolution results diverge from the higher-resolution ones at larger scales compared to the density bias. The velocity field is a weighted average of the particle velocities (Equation (23)), and a simple deconvolution using the mass assignment window function (Equation (22)) is only approximately correct.

In AMBER, we adopt the TSC assignment scheme with interlacing and deconvolution as the default method for constructing density and velocity fields from the LPT particles. This approach reduces aliasing and smoothing, and produces the most accurate one-point distributions and two-point correlations. Accurate density and velocity fields are essential components for modeling EoR observables. The 21 cm signal and Lyα forest both depend on the neutral hydrogen density n_HI( x ) and line-of-sight velocity field v_los( x ). For the CMB, the patchy Thomson optical depth and patchy KSZ effect depend on the electron number density n_e( x ) and line-of-sight momentum n_e( x )v_los( x ), respectively.

We construct halo mass density fields with ESF and compare them against results from a high-resolution N-body simulation. For reference, we use the halo catalog from SCORCH II (Section 2.1) and select all halos above a minimum mass ${M}_{\min }={10}^{8}\ {h}^{-1}{M}_{\odot }$, which corresponds to the atomic-cooling limit for galaxy formation at redshift z ≈ 8. To assign the halo mass to interlaced meshes, we choose the CIC scheme as it gives both reduced aliasing and strong correlation with the matter density field. The NGP scheme has the most aliasing, while the TSC scheme has the most stochasticity as it spreads the small halos over too many larger grid cells.

For ESF-L, we calculate the collapsed fraction and mass in Lagrangian space,

$$\begin{eqnarray}&&{m}_{\mathrm{coll}}({\boldsymbol{q}})={f}_{\mathrm{coll}}({\boldsymbol{q}}){m}_{{\rm{p}}},\end{eqnarray} \tag{ 33 }$$

and for LPT particles of mass m_p, move the particles to their comoving positions x (Equation (18)), and then assign the collapsed mass to the interlaced meshes (Equation (21)) to obtain the halo mass density field ρ_h( x ). For ESF-E, we calculate the collapsed fraction field in Eulerian space, and then the halo mass density field is calculated as

$$\begin{eqnarray}&&{\rho }_{{\rm{h}}}({\boldsymbol{x}})={f}_{\mathrm{coll}}({\boldsymbol{x}}){\bar{\rho }}_{{\rm{m}}}[1+D(z){\delta }_{{\rm{f}}}({\boldsymbol{x}})],\end{eqnarray} \tag{ 34 }$$

where ${\bar{\rho }}_{{\rm{m}}}$ is the average matter density, and the growth factor D(z) reverses the linear extrapolation done earlier.

ESF requires that the mass smoothing scale M_s be greater than the minimum halo mass M_min in Equation (30). The smoothing scale also cannot be lower than the particle mass resolution. For the fiducial lower resolution with 64³ particles on a 64³ mesh in the 50 h⁻¹Mpc box, the particle mass is m_p = 4.0 × 10¹⁰ h⁻¹ M_⊙, and the cubical cell size is l_c = 0.8 h⁻¹Mpc. Therefore, we will choose and vary the smoothing scales such that M_s ≥ m_p and ${R}_{{\rm{s}}}={\left[{M}_{{\rm{s}}}/(4\pi {\bar{\rho }}_{0}/3)\right]}^{1/3}\gtrsim 0.5\ {h}^{-1}\mathrm{Mpc}$.

Figure 11 is a visualization of the halo mass density fields ρ_h( x )/〈ρ_h〉 for ${M}_{\min }={10}^{8}\ {h}^{-1}{M}_{\odot }$ at z = 8. Here, we normalize each density field using the same mean 〈ρ_h〉 taken from the N-body result to preserve the relative differences in ρ_h( x ). We also use the fiducial smoothing scale M_s = m_p. The halo images appear more pixelated than the corresponding density and velocity images (Figures 5, 8) because of the lower resolution shown here. The ESF-L images with the sharp k-space and tophat filters are remarkably similar to the N-body result, but the ESF-E images are noticeably different. ESF-E produces relatively higher halo densities in prominent collapsed regions and lower values elsewhere, and this larger range of density contrast can be traced back to using the evolved density field rather than the initial density field for the filtering process.

Zoom In Zoom Out Reset image size

Figure 12 shows the halo auto-power spectra ${{\rm{\Delta }}}_{\mathrm{hh}}^{2}(k)$ and cross correlations. The halo power spectrum is approximately power law with slope $n=d\mathrm{ln}P/d\mathrm{ln}k\approx -1.5$ at z = 8. In comparison with the matter power spectrum (Figure 7), we find the scale-dependent halo-matter bias b_hm(k) with a large-scale value of ≈3.7. Similarly, the halo-matter cross correlation r_hm(k) is close to unity on the largest scales, but drops below 0.9 for k ≳ 1 h Mpc⁻¹. Therefore, we cannot accurately model the halo mass density field assuming linear bias, even if we allow for scale dependence.

Zoom In Zoom Out Reset image size

In cross correlation with the N-body halo mass density field, the ESF-L with the sharp k-space filter gives the best agreement, closely followed by with the tophat filter. Both the bias b_hh(k) and cross correlation r_hh(k) differ from unity by ≲0.05 down to the smallest scale shown. The strong agreement is remarkable given that we have adopted nominal choices for the filter radius R_sk and collapse barrier δ_c without any fine-tuning. ESF-E with either the sharp k-space filter or the tophat filter gives poorer agreement with the N-body results. The larger amplitude in the power spectrum and bias is consistent with the larger range in the density contrast seen in the halo density images (Figure 11).

We now vary the minimum halo mass M_min and smoothing mass M_s and quantify the mass-weighted average collapse fraction 〈f_coll〉 in Figure 13. As M_min is increased above the atomic-cooling limit while fixing M_s = 4 × 10¹⁰ h⁻¹ M_⊙ and R_s = 0.5 h⁻¹Mpc, we again find that ESF-L with the sharp k-space filter gives the best agreement, followed by with the tophat filter. However, ESF-E overpredicts the average collapse fraction by a factor of ≳3, with increasingly larger differences toward higher M_min. Note that we should expect small differences in the average collapse fraction. In the N-body halo finder (Trac et al. 2015), the halo mass M₂₀₀ is defined such that the average halo density ${\bar{\rho }}_{{\rm{h}}}$ is 200 times the average matter density ${\bar{\rho }}_{{\rm{m}}}(z)$, which is also equal the critical density ${\bar{\rho }}_{\mathrm{crit}}(z)$ at high redshifts. In ESF, the halo mass is not clearly defined; although in spherical collapse theory, the average halo density is usually 18π² ≈ 180 times the critical density at high redshifts or in Einstein-de Sitter cosmologies.

Zoom In Zoom Out Reset image size

As M_s and R_s are increased while fixing ${M}_{\min }={10}^{8}\ {h}^{-1}{M}_{\odot }$, we find that ESF-L gives reliably accurate and stable values for the average collapsed fraction that are within 10% of the N-body halo value. However, ESF-E gives highly variable results that are only in agreement at large smoothing radius R_s ≳ 10 h⁻¹Mpc. Similarly, we also find that the power and bias change as the smoothing is varied, which is not found for ESF-L.

The sensitivity of ESF-E to the smoothing and filtering scales is especially troublesome for the semi-numerical models of reionization that are based on this approach. For a given ionized region, the filtering scale is adaptively varied until the number of ionizing photons equals the number of hydrogen atoms. For the entire modeled volume, there is a distribution of filtering scales rather than just one value, and as a consequence, there is no one value for the average collapse fraction and large-scale halo bias. Furthermore, as one changes the reionization parameters, the distribution of filtering scales will change, and so too will the halo mass density field and power spectrum, both of which should be independent of the reionization history.

In AMBER, we adopt the ESF-L with the sharp k-space filter as the default choice for modeling the halo mass density field. In the next section, we will discuss that our abundance matching technique does not depend on the overall average density 〈ρ_h〉 nor the overall normalization of the halo bias b_hm. In future work, we can calibrate the EPS collapsed fraction relation (Equation (30)) using the N-body simulations and halo catalogs from SCORCH I and II. Accurate halo abundance and density fields are required to properly model high-redshift galaxies and quasars as radiation sources.

The radiation field is specified by the hydrogen photoionization rate,

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{HI}}({\boldsymbol{x}})=4\pi {\int }_{{\nu }_{\mathrm{HI}}}^{\infty }\displaystyle \frac{{J}_{\nu }(\nu ,{\boldsymbol{x}})}{h\nu }{\sigma }_{\mathrm{HI}}(\nu )d\nu ,\end{eqnarray} \tag{ 35 }$$

where J_ν is the specific intensity, σ_HI is the hydrogen photoionization cross-section, and ν_HI is the Lyman limit frequency. In our RadHydro simulations, the RT is calculated using an accurate but expensive raytracing algorithm (Trac & Cen 2007). However, this is an unnecessarily costly approach for calculating the radiation field just for abundance matching, which does not depend on the overall normalization. In AMBER, we model the specific intensity field as

$$\begin{eqnarray}&&{J}_{\nu }(\nu ,{\boldsymbol{x}})=\int \displaystyle \frac{{S}_{\nu }(\nu ,{{\boldsymbol{x}}}^{{\prime} })}{4\pi | {\boldsymbol{x}}-{{\boldsymbol{x}}}^{{\prime} }{| }^{2}}\exp \left(-\displaystyle \frac{| {\boldsymbol{x}}-{{\boldsymbol{x}}}^{{\prime} }| }{{l}_{\mathrm{mfp}}}\right){d}^{3}{x}^{{\prime} },\end{eqnarray} \tag{ 36 }$$

where S_ν is the source function, the 1/(4π r²) term is the inverse-square law for flux, and the ${e}^{-r/{l}_{\mathrm{mfp}}}$ term is the transmitted fraction over the radiation mean free path l_mfp.

The source field is related to the specific luminosity density and can be expressed as

$$\begin{eqnarray}&&{S}_{\nu }(\nu ,{\boldsymbol{x}})={f}_{\mathrm{esc}}{\epsilon }_{\nu }(\nu ){\rho }_{\mathrm{sfr}}({\boldsymbol{x}}),\end{eqnarray} \tag{ 37 }$$

where ρ_sfr is the star formation rate density, _ν is the radiative energy per unit star formation rate per frequency, and f_esc is the radiation escape fraction (e.g., Madau et al. 1999; Robertson et al. 2010). The star formation rate density field is often modeled using the halo mass density field,

$$\begin{eqnarray}&&{\rho }_{\mathrm{sfr}}({\boldsymbol{x}})=\displaystyle \frac{{f}_{\mathrm{star}}{\rho }_{\mathrm{halo}}({\boldsymbol{x}})}{{\tau }_{\mathrm{star}}},\end{eqnarray} \tag{ 38 }$$

where f_star and τ_star are the star formation efficiency and timescale, respectively (e.g., Cen & Ostriker 1992; Trac & Cen 2007). The source field is proportional to the halo mass density field, while the normalization depends on the product of several astrophysical parameters and functions (f_esc, f_star, τ_star, _ν ), which may have redshift dependence that further alters the reionization history. In semi-numerical methods based on the ESF, this combination is often represented by a single efficiency parameter ζ, which is generally assumed to be constant for simplicity. With our abundance matching technique, we do not need to specify the astrophysical quantities to calculate the normalization of the the source field, but instead get to cast them and their effects in terms of the the redshift midpoint, duration, and asymmetry parameters that directly set the reionization history (Section 3.1).

In AMBER, we adopt the common assumption S_ν ( x ) ∝ ρ_halo( x ) and ignore the overall normalization for the source field as it is not necessary for abundance matching. The radiation field is computed quickly using FFTs to perform the convolution in Equation (36). In practice, we can abundance match the specific intensity field instead of the photoionization rate field, as the integration over frequency only changes the overall normalization for the latter.

We use an effective mean free path to account for the attenuation of the radiation field. Our use of a mean free path parameter is more physical than the maximum bubble size imposed in some semi-numerical models based on the ESF. It is more self-consistent with RT theory (e.g., Madau et al. 1999) and in better agreement with observational measurements (e.g., Worseck et al. 2014) extrapolated to higher redshifts for the EoR. It also helps to mitigate uncertainties on radiation sinks. Recently, Davies & Furlanetto (2021) have also introduced an ionizing photon mean free path as an improved alternative to the maximum filtering scale in ESF methods like 21cmFAST. Their implementation is different in that l_mfp ultimately shows up in their scale-dependent ionization criterion, whereas we use it to directly compute an attenuated radiation field.

The reionization-redshift field ${z}_{\mathrm{re}}({\boldsymbol{x}})$ is assumed to be correlated with a given field, either the matter density field, the halo mass density field, or the radiation field. A region with higher matter density, halo density, or radiation intensity is considered to be photoionized earlier and has a higher reionization redshift. The abundance matching technique assigns redshift values such that the reionization history follows a given mass-weighted ionization fraction ${\bar{x}}_{{\rm{i}}}(z)$, specified with the redshift midpoint, duration, and asymmetry parameters and interpolated with a Weibull function (Section 3.1).

We perform the abundance matching on a correlated field at a single redshift for computational efficiency, but it can also be done tomographically using multiple redshift intervals. The correlated field is first constructed at the redshift midpoint z_mid, and the data array is ranked in descending order using a parallel quicksort. For the nth rank order cell out of a total of N, all cells with indices m ≤ n are considered ionized. The reionization redshift z_n is calculated by equating the cumulative mass fraction with the mass-weighted ionization fraction:

$$\begin{eqnarray}&&\displaystyle \frac{\sum _{m=1}^{n}1+{\delta }_{m}({z}_{n})}{\sum _{m=1}^{N}1+{\delta }_{m}({z}_{n})}={\bar{x}}_{{\rm{i}}}({z}_{n}),\end{eqnarray} \tag{ 39 }$$

where the matter overdensity for the mth cell is linearly extrapolated from the redshift midpoint,

$$\begin{eqnarray}&&{\delta }_{m}({z}_{n})=\displaystyle \frac{D({z}_{n})}{D({z}_{\mathrm{mid}})}{\delta }_{m}({z}_{\mathrm{mid}}).\end{eqnarray} \tag{ 40 }$$

The linear growth is appropriate for the modest overdensities in lower-resolution semi-numerical models, and because we have chosen the redshift midpoint as the pivot.

We construct the abundance matching models using the matter density field (AM-M), halo mass density field (AM-H), and radiation intensity fields (AM-R), and compare them with the RadHydro reionization-redshift fields (Section 2.3). For the AM-M models, we use the TSC particle assignment scheme with interlacing and deconvolution to construct the matter density field (Section 3.2.1). For the AM-H models, we use the ESF-L with sharp k − space filter for constructing the halo mass density field and set the minimum halo mass to ${M}_{\min }={10}^{8}\ {h}^{-1}{M}_{\odot }$ (Section 3.3.1). While this corresponds to the threshold for atomic-cooling halos, the RadHydro simulations are more complex because of the nonmonotonic star formation efficiency and episodic star formation (Trac et al. 2015; Doussot et al. 2019). For the AM-R models, we vary the radiation mean free path to find the best agreement between the reionization-redshift fields.

Figure 14 is a visualization of the reionization-redshift fields ${z}_{\mathrm{re}}({\boldsymbol{x}})$ from the RadHydro simulations and abundance matching models. While all of the models have the same reionization history for a given simulation, their reionization-redshift fields are visibly different. The AM-R images most closely resemble the simulation results and correctly show that large-scale regions near radiation sources are generally reionized earlier than large-scale regions far from sources. However, the AM-M and AM-H images appear too grainy. The low-density regions just outside of collapsed regions are assigned low redshifts and considered to be reionized late despite their close proximity to the radiation sources. We can possibly improve these results by smoothing the matter density and halo density fields prior to abundance matching.

Zoom In Zoom Out Reset image size

Figure 15 shows the reionization-redshift auto-power spectra ${{\rm{\Delta }}}_{\mathrm{zz}}^{2}(k)$ and cross correlations relative to the RadHydro normalized redshift δ_z (Equation (4)). We only show Sims 0 and 2 for clarity. The AM-R power spectra are most similar to the simulation results, rising in power until a characteristic peak scale and then declining in amplitude with k. However, the AM-M and AM-H power spectra generally always increase with k, similar to the mass density (Figure 7) and halo density power spectra (Figure 12). The AM-R models are the least biased with b_zz ≈ 1 for k ≲ 1 h Mpc⁻¹, while the AM-M and AM-H biases significantly deviate from unity. Furthermore, the AM-R models have improved cross correlation r_zz that approach unity on large scales. The bias and stochasticity on small scales can be attributed to three main factors. In the RadHydro simulations, the episodic star formation and spatially varying mean free path are not accounted for. Furthermore, the simulations and models have differences in smoothing near the grid scale.

Zoom In Zoom Out Reset image size

In AMBER, we choose to abundance match using the radiation field evaluated at the redshift midpoint. In future work, we can improve the spatial accuracy by incorporating varying mean free paths (e.g., Cain et al. 2021; Davies & Furlanetto 2021). We can also perform the abundance matching tomographically using multiple redshifts.

As the redshift midpoint decreases over the range 7.0 ≤ z_mid ≤ 9.0, we see, overall, later reionization in Figure 19. The one-point PDF(δ_z) becomes wider, the two-point power spectrum ${{\rm{\Delta }}}_{\mathrm{zz}}^{2}(k)$ increases overall in amplitude, but the redshift-matter bias b_zm(k) remains almost constant in Figure 20. In Equation (4) for the normalized redshift δ_z, the numerator (${z}_{\mathrm{re}}-{\bar{z}}_{\mathrm{re}}$) remains approximately constant, but the denominator ($1+{\bar{z}}_{\mathrm{re}}$) decreases. Therefore, δ_z grows by a factor of ${\bar{a}}_{\mathrm{re}}=1/(1+{\bar{z}}_{\mathrm{re}})$, whereas the matter overdensity δ grows by a factor D(z_mid) ≈ 1/(1 + z_mid). We find that the normalized redshift grows similarly to the matter overdensity as the bias is only ≈5% lower for z_mid = 7.0 than for z_mid = 9.0. Note that, while the radiation field also grows like the matter density, the abundance matching is insensitive to the overall amplitude.

As the duration increases over the range 2.0 ≤ Δ_z ≤ 6.0, we see a larger range of values in the reionization-redshift images. The PDF(δ_z) has larger variance, and the power spectrum increases overall in amplitude. Note that the variance is given by ${\sigma }_{{\rm{z}}}^{2}=\langle | {\delta }_{{\rm{z}}}({\boldsymbol{x}}){| }^{2}\rangle =\int {{\rm{\Delta }}}_{\mathrm{zz}}^{2}(k)d\mathrm{ln}k$. Therefore, we expect the large-scale bias to scale with σ_z, which in turn scales with Δ_z. We find that the large-scale bias changes by a factor ≈ 3.1 when the duration increases from Δ_z = 2.0 to Δ_z = 6.0.

The large-scale bias can differ significantly from the single value of b_B13 = 1/δ_c = 0.593 adopted in the RLS method (Battaglia et al. 2013b). Barkana & Loeb (2004) derived this value assuming that fluctuations in the reionization-redshift field are solely determined by fluctuations in the collapsed mass density field from EPS. Relating the growth of ionized regions to the growth of dark matter halos would correspond to a particular duration value. However, the reionization history can depend on other nonconstant factors, and the large-scale bias can have a range of values in general.

As the asymmetry increases over the range 1.0 ≤ A_z ≤ 5.0, the images show a smaller range of redshift values for the regions reionized after the midpoint and a larger range of redshift values for the regions reionized earlier than the midpoint. The PDF(δ_z) becomes skewed, and the power spectrum changes shape. There is less power on large scales and more power on small scales as the characteristic peak shifts to smaller scales. For a larger A_z at fixed duration, the first portion (z_mid < z < z_ear) of the reionization history becomes longer, while the second portion (z_lat < z < z_mid) becomes shorter. We have just seen that increasing (decreasing) the duration increases (decreases) the overall power. We find that lengthening the first half of the portion slightly increases the power on small scales while shortening the second portion more significantly decreases the power on large scales. Correspondingly, the large-scale bias decreases with more asymmetry in the reionization history.

Over the minimum halo mass range ${10}^{7}\leqslant {M}_{\min }/[{h}^{-1}{M}_{\odot }]\,\leqslant {10}^{10}$, there are only minor changes to the reionization-redshift field. The PDF(δ_z) only has minor changes since the redshift midpoint, duration, and asymmetry are kept fixed. There is minor suppression in ${{\rm{\Delta }}}_{\mathrm{zz}}^{2}(k)$ and b_zm(k) on small scales. For larger M_min, the collapsed fraction decreases and the halo bias increases. Correspondingly, this decreases the amplitude and increases the bias of the radiation intensity field (Equation (36)). However, the abundance matching technique is insensitive to the overall normalization of J_ν and to the overall normalization of the fluctuating component δ_J. The abundance matching only depends on the relative ranking of the field values.

On one hand, varying the minimum halo mass and halo density field will generally change the reionization history and process if astrophysical parameters such as the star formation efficiency, photon production efficiency, and radiation escape fraction are kept fixed. On the other hand, varying these astrophysical terms, which can be functions of mass and redshift, can offset the changes from varying M_min and produce the same reionization history. With AMBER, we find that when the redshift midpoint, duration, asymmetry, and radiation mean free path are kept fixed, then the reionization-redshift field has only minor dependence on the minimum halo mass. For certain applications, it may then be possible to ignore M_min and reduce the number of free parameters.

As the radiation mean free path increases over the range 1.0 ≤ l_mfp/[h⁻¹Mpc] ≤ 5.0, we see more spatial coherence in the reionization-redshift images. For a larger l_mfp, there is a larger coherence length and less small-scale structure in the radiation intensity field used for abundance matching. The PDF(δ_z) only has minor changes since the redshift midpoint, duration, and asymmetry are kept fixed. The reionization-redshift power spectrum changes shape though, having more power on large scales and less on small scales. Increasing l_mfp produces relatively larger ionized regions, resulting in the characteristic peak shifting to larger scales.

In the RLS method, Battaglia et al. (2013b) change the duration by varying a_RLS and k_RLS. But without varying b_RLS, they also change the correlation length of ionized regions at the same time. More specifically, they obtain smaller correlation lengths for longer durations. With AMBER, we find longer Δ_z increases the bias, while shorter l_mfp decreases the bias. This explains the inverse relationship between the duration and correlation length when b_RLS is fixed in the RLS method. AMBER has the flexibility to change Δ_z and l_mfp independently.