Prediction of the 21-cm signal from reionization: comparison between 3D and 1D radiative transfer schemes

Prediction of the 21-cm signal from reionization: comparison between 3D and 1D radiative transfer... Abstract Three-dimensional radiative transfer simulations of the epoch of reionization can produce realistic results, but are computationally expensive. On the other hand, simulations relying on one-dimensional radiative transfer solutions are faster but limited in accuracy due to their more approximate nature. Here, we compare the performance of the reionization simulation codes grizzly and C2-ray which use 1D and 3D radiative transfer schemes, respectively. The comparison is performed using the same cosmological density fields, halo catalogues, and source properties. We find that the ionization maps, as well as the 21-cm signal maps from these two simulations are very similar even for complex scenarios which include thermal feedback on low-mass haloes. The comparison between the schemes in terms of the statistical quantities such as the power spectrum of the brightness temperature fluctuation agrees with each other within 10 per cent error throughout the entire reionization history. grizzly seems to perform slightly better than the seminumerical approaches considered in Majumdar et al. which are based on the excursion set principle. We argue that grizzly can be efficiently used for exploring parameter space, establishing observations strategies, and estimating parameters from 21-cm observations. radiative transfer, galaxies: formation, intergalactic medium, dark ages, reionization, first stars, cosmology: theory 1 INTRODUCTION The epoch of reionization (EoR), when the primordial neutral hydrogen (H i) in the intergalactic medium (IGM) was ionized by the radiation from the first sources, is one of the milestone events in the evolution of the Universe. However, many important facts such as the exact timing of this event, the sources responsible for reionization, the impact of the EoR on the later stages of the structure formation, etc. are currently not or poorly known. Various probes of the EoR such as observations of high-redshift quasars (Becker et al. 2001; Fan et al. 2006; Goto et al. 2011; Mortlock et al. 2011; McGreer, Mesinger & D'Odorico 2015), the cosmic microwave background radiation (CMB) (Zaldarriaga & Seljak 1997; Komatsu et al. 2011; Planck Collaboration XV 2016b; Planck Collaboration XIII 2016a), high-redshift galaxies (Kashikawa et al. 2011; Bouwens et al. 2015; Mitra, Choudhury & Ferrara 2015), and Ly α emitters (McQuinn et al. 2007b; Dijkstra, Mesinger & Wyithe 2011; Choudhury et al. 2015) have been used to constrain some of these unknown facts. However, these probes can only provide limited information. For example, combining the observations of high-redshift quasars and the CMB, one can constrain the probable period of the event to lie between redshifts 15 and 6 (Choudhury & Ferrara 2006; Fan et al. 2006; Malhotra & Rhoads 2006; Mitra, Choudhury & Ferrara 2011, 2012). The evolving redshifted 21-cm signal from the H i in the IGM has the potential to revolutionize our understanding of reionization (Furlanetto, Oh & Briggs 2006; Morales & Wyithe 2010; Pritchard & Loeb 2012). Therefore a range of low frequency radio telescopes such as the Giant Metrewave Radio Telescope1 (Ghosh et al. 2012; Paciga et al. 2013), the Precision Array for Probing the Epoch of Reionization2 (Parsons et al. 2014), the Murchison Widefield Array3 (Bowman et al. 2013; Tingay et al. 2013), and the Low Frequency Array4 (van Haarlem et al. 2013; Patil et al. 2017) are attempting to detect this faint signal from the EoR. The primary aim of these first generation experiments is to detect the signal in terms of statistical quantities such as the variance, skewness, power spectrum, etc. However, future radio telescopes such as the Square Kilometre Array (SKA)5 should have enough sensitivity to produce image cubes of the signal (Mellema et al. 2015; Ghara et al. 2017). A large number of theoretical models based on either analytic calculations (e.g. Furlanetto, Zaldarriaga & Hernquist 2004; Paranjape & Choudhury 2014), seminumerical simulations (Mesinger & Furlanetto 2007; Zahn et al. 2007; Santos et al. 2008; Choudhury, Haehnelt & Regan 2009), or numerical simulations with radiative transfer (Iliev et al. 2006b; Mellema et al. 2006; McQuinn et al. 2007a; Shin, Trac & Cen 2008; Baek et al. 2009; Thomas et al. 2009; Ghara, Choudhury & Datta 2015) have been used to study the redshifted 21-cm signal from the EoR. Through these works a fair amount of theoretical understanding of the impact of various physical processes on the signal has been achieved. Aspects which have been considered include the impact of different kind of sources to the redshifted 21-cm signal, thermal feedback and peculiar velocities of the gas in the IGM. The simulations are not only important to extract information from the observations, but also to design observation strategies. One thing to keep in mind is that the 21-cm observations do not provide direct information about the properties of the sources responsible for ionizing the IGM, minimum mass of the dark matter haloes forming stars, thermal feedback, the effect of peculiar velocities of the gas in the IGM, etc. Due to the complexity of the astrophysical effects, one needs to use simulations to extract this information from the observations. In principle, this can be done using parameter estimation techniques such as the Markov chain Monte Carlo (Greig & Mesinger 2015) or artificial neural networks (Shimabukuro & Semelin 2017) in which the observations expressed in quantities such as power spectrum, variance, or skewness will be compared to the outputs of the simulations. The various simulation methods each have their own advantages and drawbacks and when using their results to design observation strategies or to interpret the observational data, one needs to keep in mind the various assumptions made by the models. While simulations of the 21-cm signal from the EoR using 3D radiative transfer methods such as Iliev et al. (2006b) and Mellema et al. (2006) are assumed to be the most accurate ones, these require supercomputer facilities to run. On the other hand, simulations performed with seminumerical methods are orders of magnitudes faster, the simplifying assumptions needed for this may lead to accuracy issues. Thus, it is absolutely necessary to compare the performance of different simulation techniques with each other so as to identify the advantages and drawbacks of these methods. This will be useful for improving the methods and for correctly interpreting the observations. Studies such as Zahn et al. (2011) compared a reionization simulation obtained using a full 3D radiative transfer technique with the result from a seminumerical one. It showed excellent agreement between these two. Similarly, Majumdar et al. (2014) compared seminumerical simulations and conditional Press-Schechter models with a full 3D radiative transfer simulation performed with C2-ray (Mellema et al. 2006). That study showed that the seminumerical simulations such as Zahn et al. (2007), Mesinger & Furlanetto (2007), Choudhury et al. (2009), and Santos et al. (2010), which use an excursion-set formalism to create the ionized bubbles, produce an ionization history very similar to the C2-ray provided that the halo lists from an N-body simulation are used. However, if the haloes are found from a conditional Press–Schechter approach, as is for example done in Alvarez et al. (2009b), Mesinger, Furlanetto & Cen (2011), the fit with the C2-ray results is considerably worse. The codes such as bears and grizzly use one-dimensional radiative transfer schemes instead of solving radiative transfer equations in 3D. Thus, these simulations assumes that the effect from each sources (i.e., the ionized regions around the sources) are perfectly spherical. When these ionized regions overlap, these codes account the overlap using special recipes. Due to these simplifying assumptions they are fast but can also be expected to be less accurate than full 3D radiative transfer schemes. Previously, Thomas et al. (2009) performed a limited comparison between the 1D radiative transfer code bears and the 3D radiative transfer code crash (Ciardi et al. 2001) for a simulation volume of size 12.5 h−1 comoving megaparsec. However, before considering these 1D schemes for parameter estimations (e.g. Patil et al. 2014), their performance for larger simulation volumes and larger numbers of sources is required. In this paper, we present a detailed comparison between the performance of the 1D radiative transfer code grizzly (Ghara et al. 2015) with the 3D radiative transfer code C2-ray. We use the same initial conditions such as the density fields, halo catalogues, source models, and compare the reionization history produced by these two schemes. We consider bubble size distributions (BSDs) as well as cross-correlations between ionization maps, while for the redshifted 21-cm signal we focus on statistical quantities such as the evolution of the power spectra. We consider two different models: (i) reionization with only large mass haloes which do not suffer from thermal feedback, (ii) reionization with both high- and low-mass haloes where the latter are assumed to be sensitive to thermal feedback. Case (i) is identical to the one used in Majumdar et al. (2014) to compare seminumerical and full radiative transfer results which therefore allows the comparison between seminumerical, 1D RT, and 3D RT results. Our paper is structured in the following way. In Section 2, we describe the simulations used in this paper. In particular different subsections describe the N-body simulation and the source model used in this study as well as a brief description of the underlying algorithms of C2-ray and grizzly. In Section 3, we present our results before we conclude in Section 4. Throughout the paper, we have chosen the cosmological parameters Ωm = 0.27, $$\Omega _\Lambda =0.73$$, ΩB = 0.044, h = 0.7 consistent with the Wilkinson Microwave Anisotropy Probe (WMAP) results (Hinshaw et al. 2013) and within the error bars consistent with Planck (Planck Collaboration XVI 2014). 2 SIMULATIONS Below we discuss the radiative transfer methods along with the source model and N-body simulation used in this study. 2.1 N-body simulation Both radiative transfer methods considered in this study use the results of the same N-body simulation. This dark matter only N-body simulation was carried out using the publicly available code cubep3m6 (Harnois-Déraps et al. 2013). The details of the simulations are (i) the size of the simulation cube is 163 comoving megaparsec (cMpc), (ii) the number of dark matter particles is 30723, (iii) the fine grid used by the particle–particle–mesh method is 61443, (iv) the mass of the dark matter particles is 5.47 × 106 M⊙. The details of the simulation can also be found in Iliev et al. (2012). Although this N-body simulation is by now rather old, we use it to enable direct comparison to the results in Majumdar et al. (2014) who compared the performance of a seminumerical method to C2-ray. Snapshots of the density field are generated from redshift 20.134 to redshift 6 using equal time intervals of 11.5 megayears (Myr). The resolution of these gridded density fields is 2563. cubep3m also produces dark matter halo lists using an on the fly halo finder based on the spherical overdensity method. The lowest mass haloes produced have a mass of ∼108 M⊙ which corresponds approximately to the mass needed to initiate star formation driven by atomic cooling. 2.2 Source model Our knowledge about the sources of reionization is very limited. Different efforts, both observational and modelling, have been undertaken to improve our understanding of the impact of various sources such as primordial galaxies (Wise et al. 2014; Xu et al. 2016), mini-quasars (Bromm & Loeb 2003; Zaroubi et al. 2007; Alvarez, Wise & Abel 2009a; Tanaka, Perna & Haiman 2012), and high-mass X-ray binaries (Fialkov, Barkana & Visbal 2014; Kaaret 2014; Knevitt et al. 2014) on the expected 21-cm signal from the EoR. However, the individual contributions of these sources towards reionization remain quite uncertain. In this study, we assume the sources of ionizing radiation formed in the dark matter haloes. We further assume their spectral energy distribution can be approximated by a blackbody of an effective temperature of 50 000 K, which is appropriate for a massive Population II star. We normalize the spectrum such that the rate of production of the ionizing photons is   \begin{equation} \dot{N_{\gamma }}=g_{\gamma }\frac{M_{\rm h} \Omega _{\rm B}}{(10 \ \rm Myr)\Omega _m m_{\rm p}}, \end{equation} (1)where gγ is the source ionization efficiency coefficient, Mh and mp are the mass of the dark matter halo and proton mass, respectively. For the fiducial model considered here, we choose gγ = 21.7 for the sources formed in dark matter haloes with mass Mh ≥ 2.2 × 109 M⊙, which will produce $$\dot{N_{\gamma }}=1.33 \times 10^{43} \times \frac{M_{\rm h}}{{\rm M}_{{\odot }}} \ \rm s^{-1}$$. This is a similar choice as the L3 model in Iliev et al. (2012). In addition, we consider a more complex reionization scenario where we use all haloes of mass Mh ≥ 108 M⊙. We divide these in two populations, low-mass atomically cooling haloes (LMACHs) in the mass range 108 M⊙ ≤ Mh < 109 M⊙ and high-mass atomically cooling haloes (HMACHs) in the mass range Mh ≥ 109 M⊙. The LMACHs are assumed to be sensitive to thermal feedback and stop producing ionizing photons once they are located in an cell which is more than 10 per cent ionized. The HMACHs are taken to be insensitive to thermal feedback. While still active the LMACHs produce ionizing photons with an efficiency gγ = 130, whereas the HMACHs have gγ = 8.7. The higher efficiency of the LMACHs is motivated by either a larger contribution from Pop III stars or a higher fraction of ionizing photons escaping. In the C2-ray results these parameters result in a rather early end of reionization around z = 8.5 and a Thomson optical depth value of 0.08, consistent with the value determined by WMAP and within the 2σ range of Planck. This scenario is identical to the L1 model of Iliev et al. (2012). The details of these scenarios are given in Table 1. Table 1. Details of the scenarios considered in this paper. LMACH and HMACH represent the low-mass and high-mass atomically cooling halo, respectively. Model  Min. halo (M⊙)  Thermal feedback  gγ  gγ  Model used in Iliev et al. (2012)        LMACH  HMACH    MODEL I (fiducial)  2.2 × 109  No  0  21.7  L3  MODEL II  108  Yes  130  8.7  L1  Model  Min. halo (M⊙)  Thermal feedback  gγ  gγ  Model used in Iliev et al. (2012)        LMACH  HMACH    MODEL I (fiducial)  2.2 × 109  No  0  21.7  L3  MODEL II  108  Yes  130  8.7  L1  View Large 2.3 3D radiative transfer simulation C2-ray The 3D radiative transfer method considered in this work uses ‘Conservative Causal Ray-tracing method’ (C2-ray). The details of the method can be found in Mellema et al. (2006). The method is briefly described in the following steps. It starts with preparing the source list in a random order at each redshift. Given the SED of the sources, the total photoionization rate (Γ) is calculated at time t at each cell in the simulation box including contributions from all the sources. This step takes into account the time evolution of the neutral fraction of the cells during the time-step which affects the optical depth from the sources. The history of the ionization fraction of hydrogen ($$x_{\rm H\, \small {II}}$$) is estimated by solving the non-equilibrium equation,   \begin{equation} \frac{{\rm d}x_{\rm H\, \small {II}}}{{\rm d}t}=\left(1-x_{\rm H\, \small {II}}\right)\left(\Gamma + n_{\rm e}C_{\rm H}\right)-xn_{\rm e}C\alpha _{\rm H}, \end{equation} (2)where ne is the electron density at the cell, CH and αH represent the collisional ionization and recombination coefficients for hydrogen, respectively. The quantity C is the clumping factor (which can be defined as 〈n2〉/〈n〉2, where n is the baryon number density) which accounts the clumpiness of the IGM. Here we have taken C = 1 but in general it can be larger than 1 or even position dependent. C2-ray was compared to other fully numerical radiative transfer codes in Iliev et al. (2006b) and Iliev et al. (2009) and was found to produce accurate results. In this study we do not consider the temperature evolution of the intergalactic gas, only its ionized state. 2.4 1D radiative transfer simulation grizzly The details of the one-dimensional radiative transfer method can be found in Ghara et al. (2015). The method mainly follows the bears algorithm (Thomas et al. 2009) with a few changes compared to the original method. To express its relation to the original bears method, we have chosen to call it grizzly. The basic idea of both methods is not to solve the one-dimensional radiative transfer equations on the fly while generating the ionization maps. Rather, they use pre-generated ionization profiles to construct the ionization maps. Here, we briefly summarize the most important steps of the algorithm. We start by choosing the range for the key source parameters: redshift, density contrast of the uniform background IGM, luminosity of the source, age of the source. We then generate a large number of ionization and kinetic temperature profiles around the sources for different combinations of these parameter values. Next we create a list of radii of the H ii bubbles as a function of different parameters. We assume the radius of the H ii region to be the radius at which the ionization fraction drops to 0.5 from the location of the source. This library of 1D profiles only needs to be compiled once for a given cosmology and spectral energy distribution of sources. For a given redshift and a certain halo position, we determine the luminosity and the age of the halo according to a source recipe, see Section 2.4.1. From this we determine the size of the H ii bubble around it using the density field and the pre-generated list of the H ii bubble radii. This step is done iteratively. We start with a small radius around the centre of the source, calculate the spherically averaged overdensity of the IGM within this radius, and then estimate the radius of H ii bubble from the list. We change the initial choice of radius such that the estimated radius of the H ii region matches with the chosen radius for the same overdensity. The parameter values corresponding to that radius are then associated with the source and the corresponding profiles will be used to generate the ionization and kinetic temperature maps. The procedure described in the previous step is performed for all haloes at the redshift being considered. Although the number of overlaps between the individual H ii regions around the haloes is negligible during the initial stages of reionization, they become important at the later stages of reionization. To incorporate the effect of overlap, we estimate the unused ionizing photons for each overlapping H ii region and distribute them around the overlapping regions such that each photon ionizes one hydrogen atom. Next, we use the pre-generated 1D ionization profiles corresponding to each halo to calculate the ionization fraction at the partially ionized region beyond the H ii regions. The ionization fraction in the overlapping partially ionized regions is then expressed as   \begin{equation} x_{\rm H\, \small {II}}({\boldsymbol x}) =\frac{\sum _i x^{i}_{\rm H\, \small {II}-1D}({\boldsymbol x}) \times n^{i}_{\rm H-1D} \times \left(1- x^{i-1}_{\rm H\, \small {II}}({\boldsymbol x})\right)}{n_{\rm H}({\boldsymbol x})}, \end{equation} (3)where the summation symbol represents a sum over all overlapping sources. Here, $$x^{i-1}_{\rm H\, \small {II}}({\boldsymbol x})$$ denotes the ionization fraction obtained after considering overlapping H ii regions from i − 1 number of sources. We take $$x^{0}_{\rm H\, \small {II}}({\boldsymbol x})=0$$. The term $$\left(1- x^{i-1}_{\rm H\, \small {II}}({\boldsymbol x})\right)$$ takes care of the fact that the ith source encounters an already partially ionized IGM, due to overlap between previous i − 1 number of sources, before contributing to ionization. Next, we generate the kinetic temperature maps using a correlation of the neutral fraction and the gas temperature in the partially ionized region (for details, see Ghara et al. 2015). 2.4.1 Implementation of source model In the model scenarios considered in this paper, we assume that dark matter haloes always have ongoing star formation, unless they are suppressed through thermal feedback, see Section 2.4.2. To select a one-dimensional profile from the library we need to specify a luminosity and age for the halo. The original bears code and Ghara et al. (2015) used a fixed age for the haloes, typically ∼ 10 Myr and also used a luminosity calculated from the instant mass of a halo at the redshift under consideration. Those models assume that the stars inside the galaxies efficiently emit hydrogen ionizing photons until an age of ∼10 Myr. However, we found this source model to be inconsistent with the one used by C2-ray in terms of the cumulative number of ionizing photons being produced. We therefore implemented another source model. In this model we follow the growth of haloes and find both the age of the halo as well as the time-averaged mass (to be called the effective mass Meff). We then calculate the luminosity based on Meff using equation (1) and use this together with the age to find the appropriate profile from the library. The grizzly code thus does not track the history of the state of the IGM gas but rather uses the history of the haloes to calculate the state of the IGM at a given redshift. It is not entirely trivial to calculate the age of the haloes accurately from the halo lists as haloes can move to new cells or multiple haloes can merge. In principle, one needs to track the haloes using a merger-tree algorithm. However, as no merger tree is supplied with the halo lists, we have used a simple method to account for their growth. For every pair of redshifts z1 and z2 we first search for a halo at the same location and it is found, we associate the change in mass with the halo whose history we are tracking. If a halo disappears from its previous grid position in some time-step, we inspect bordering grid cells. If we find a newly formed heavier halo in a bordering grid cells, we assume that the identified halo must be the present state of the halo which disappeared. If we are not able to find any newly formed halo in the neighbourhood, then we assume that a merger happened. In this case, we find previously existing heavier halo at the neighbour and set its age as the maximum of the ages of the identified and disappeared halo. Although this a posteriori reconstruction of the halo histories is rather crude, it does allow us to obtain the correct cumulative number of ionizing photons produced during reionization. 2.4.2 Implementation of thermal feedback In the second C2-ray simulation we include LMACHs which are fully suppressed once their environment has been sufficiently ionized. The criterion used is that when the ionization fraction in the cell containing the LMACH is larger than 0.1, thermal feedback is assumed to stop the production of ionizing photons.7 To implement this recipe in grizzly, it is not trivial as the code only places ionized regions around active sources at the redshift of interest. We implemented the following simple prescription to take into account thermal feedback effects in grizzly. We use the ionization map from the previous time-step to select the active sources at a certain time-step using the condition on the ionization fraction. We add to these the previously active sources which became inactive at the present time-step but assign to them their effective mass and age corresponding to the last time they were active. In principle one needs to consider the recombination in the already ionized regions to correctly account for the suppression. However, it is not straightforward to implement recombination in an already ionized region in grizzly. This is also true for other seminumerical simulations such as models considered in Majumdar et al. (2014) and Sobacchi & Mesinger (2014). However, recombinations play a minor role, especially since previously suppressed LMACH haloes become active again after their haloes have grown into HMACHs and the amount of time that recombinations act upon the relic H ii regions will thus be limited. Since the criterion for suppression is $$x_\mathrm{H\, \small {II}} > 0.1$$, the treatment of partially ionized cells is important. We therefore implemented an improvement over the previous version of grizzly when estimating the ionization fraction in cells containing an ionization front. We now assign the volume averaged ionized ionization fraction to these cells calculated from the higher resolution 1D profiles. 2.5 Differences between C2-ray and grizzly It is important to realize that C2-ray and grizzly treat the calculation of the reionization history very differently. C2-ray follows the time evolution of the ionization fractions by taking the state at the end of the previously calculated time-step as initial condition for the new time-step. This means that any simulation has to be followed from the emergence of the first sources to a chosen end point, for example when the IGM reaches an average ionization fraction of 99 per cent. The advantage is that any time-dependent effects such as recombinations and changes in the source luminosity or position are treated correctly. The grizzly as well as the related bears code do not use the previous state of the IGM as an initial condition for the next time-step but instead construct the state of the IGM from scratch for every output. The time evolution is captured by considering the growth history of dark matter haloes when calculating their luminosity and in the case of suppressed haloes by using the recipe described above. This approach has as an advantage that it is not required to calculate the entire reionization history if this is not needed. The disadvantage is that time-dependent effects such as recombinations, radiative cooling and changes in the source properties have to be implemented in an approximate manner. These 1D radiative transfer reionization codes share this property with the seminumerical methods based upon the excursion set approach such as 21cmfast (Mesinger, Furlanetto & Cen 2011). 2.6 21-cm signal The differential brightness temperature δTb of the 21-cm signal can be expressed as,   \begin{eqnarray} \delta T_{\rm b}(\boldsymbol {x}, z) & = & 27 \ x_{\rm H\, \small {I}} (\boldsymbol {x}, z) [1+\delta _{\rm B}(\boldsymbol {x}, z)] \left(\frac{\Omega _{\rm B}h^2}{0.023}\right) \nonumber \\ &&\times \, \left(\frac{0.15}{\Omega _{\rm m}h^2}\frac{1+z}{10}\right)^{1/2} \left(1-\frac{T_{\gamma }}{T_{\rm S}} \right)\, \rm {mK}, \end{eqnarray} (4)where the quantities $$x_{\rm H\, \small {I}}$$, δB, and Tγ(z) = 2.73 × (1 + z) K denote the neutral hydrogen fraction, baryonic density contrast, and the CMB brightness temperature, respectively, each at position $$\boldsymbol {x}$$ and redshift z. TS represents the spin temperature of hydrogen in the IGM. In this study we will adopt the high spin temperature approximation: TS ≫ Tγ. In this case δTb becomes insensitive to the actual value of the spin temperature. This approximation is valid when the spin temperature is coupled to the IGM gas temperature through Ly α coupling and the IGM has been sufficiently heated (TK ≫ Tγ) by the X-ray sources. Most reionization scenarios reach this state well before substantial ionized regions form (Cohen et al. 2017). After generating the ionization fraction cubes as described in Sections 2.3 and 2.4, it is straightforward to generate the differential brightness temperature cubes using equation (4). Below we will compare $$x_{\rm H\, \small {II}}$$ and δTb results at fixed redshifts. We will not consider the effect of peculiar velocities in the IGM (leading to the so-called redshift-space distortions). The comparison of fully numerical and seminumerical models in Majumdar et al. (2014) suggest that redshift-space distortions will not change the main conclusions of the comparison. 3 RESULTS In this section, we compare the different outputs from the two simulation codes considered in this work, namely C2-ray and grizzly. Before evaluating the different scenarios listed in Table 1, we first consider two single source scenarios to test the impact of the fundamental difference between the two codes as described in Section 2.5. 3.1 Test scenarios: single sources To gauge the basic differences between C2-ray and grizzly we follow the growth of an ionized region for two single source scenarios. For both we take the IGM around the source to have a uniform density at the mean value of the Universe. We choose the size of the simulation box to be 20 h−1 cMpc and use a grid of size 1003. We assume that the source formed at redshift 15 and we follow the impact of the source on the IGM up to redshift 6 using time-steps of 10 Myr. TEST A: The source is hosted by a dark matter halo with constant mass which we take to be 109 M⊙. We choose an efficiency factor gγ = 21.7, as defined in equation (1), which sets the rate of production of ionizing photons to $$\dot{N_\gamma } = 1.33\times 10^{52} \ \rm s^{-1}$$. TEST B: The source is hosted by a dark matter halo with a mass which changes exponentially with redshift as (see e.g. Wechsler et al. 2002),   \begin{equation} M_h(z)=m_i \times \exp \left[\beta (z_i -z)\right], \end{equation} (5)where mi = 109 M⊙, zi = 15, β = 0.5. This scenario mimicks the typical growth of dark matter haloes and thus will show the impact of capturing source evolution through an effective mass Meff in grizzly. We choose the same efficiency gγ as in Test A. The performance of the methods for these test scenarios is shown in Fig. 1. The left-hand and middle panels show the evolution of the neutral fraction in the simulation volume for TEST A and B, respectively, whereas the right-hand panel shows relative difference between the two codes. We see that they produce very similar results for both tests with differences less than 10 per cent. However, for TEST A the ionization front grows faster in the grizzly results. The probable reason is that grizzly underestimates the impact of recombinations. While placing the H ii bubble around the sources, we assume that the source has been emitting ionizing photons at a constant rate determined by its effective mass (Meff) during its entire age and we include the effect of recombinations using the densities at the redshift which we are considering. Since the expansion of the Universe causes densities to be higher at higher redshifts, this will underestimate the impact of recombinations over the age of the source. C2-ray calculates the recombination rates at every time-step using the densities at that redshift and thus captures this evolution of the recombination rate. Figure 1. View largeDownload slide Left-hand panel: Redshift evolution of the mass averaged neutral fraction due to a single source in a simulation volume of size 20 h−1 cMpc. The solid and dashed curves represent C2-ray and grizzly respectively both including recombination processes, while the dotted curve corresponds to grizzly without any recombination processes. The emission rate of the ionizing photons is proportional to the mass of the halo which is chosen to be constant (109 M⊙) in this case. Middle-panel: Same as the left-hand panel, but for a dark matter halo whose mass changes exponentially with redshift. Right-hand panel: Relative difference between the ionization fractions calculated from C2-ray and grizzly as a function of redshift for the test scenarios shown in the left-hand and middle panels of the figure. Figure 1. View largeDownload slide Left-hand panel: Redshift evolution of the mass averaged neutral fraction due to a single source in a simulation volume of size 20 h−1 cMpc. The solid and dashed curves represent C2-ray and grizzly respectively both including recombination processes, while the dotted curve corresponds to grizzly without any recombination processes. The emission rate of the ionizing photons is proportional to the mass of the halo which is chosen to be constant (109 M⊙) in this case. Middle-panel: Same as the left-hand panel, but for a dark matter halo whose mass changes exponentially with redshift. Right-hand panel: Relative difference between the ionization fractions calculated from C2-ray and grizzly as a function of redshift for the test scenarios shown in the left-hand and middle panels of the figure. Interestingly, the match between the two codes is better for TEST B where the halo mass is increasing in time (see the middle panel of Fig. 1). Since in this case most photons are emitted at later times and the growth of the H ii region is faster at later times, the effect of underestimating recombinations is less severe in this case. Below when considering more complex simulations with multiple sources we will see further demonstrations of the impact of recombinations. We conclude that grizzly is capable of approximating the expansion of an ionization front around a growing halo to within 10 per cent accuracy, giving confidence that the effective mass approach works well. When considering the multiple sources case below, larger differences than the ones found here will have to be due to how grizzly deals with overlapping H ii regions and source suppression. 3.2 Model I : Large haloes only Now, we will consider a model where the reionization is driven by haloes with mass Mh ≥ 2.2 × 109 M⊙ in a volume 163 comoving Mpc across. These haloes are never suppressed. This model is same as the L3 model in Iliev et al. (2012) and was also used in Majumdar et al. (2014) to compare seminumerical simulation techniques with a fully numerical one. Not taking into account the time needed to construct the library of 1D ionization profiles, the simulation run time for the whole reionization history is ∼5.5 CPU hours8 for grizzly. This number is without generating the TK maps and Ly α flux maps which are not being used here. For comparison, the entire C2-ray simulation took ∼105 CPU hours. So there is a factor ∼105 gain in computing time when using grizzly. 3.2.1 Reionization history Using the source model described in Section 2.2 and following the methods described in Sections 2.3 and 2.4, we generate the ionization maps in the simulation box. The ionization histories from these two simulations are shown in the bottom panel of Fig. 2. In this scenario, the first sources appear at redshift 19 and ends around z ∼ 8.4. The Thomson scattering optical depth from the grizzly and C2-ray simulations are 0.0696 and 0.0697, respectively. As can be seen from the figure, the evolutions of the mass averaged neutral fraction as a function of redshift from these two simulations are very similar. Figure 2. View largeDownload slide Bottom panel: Redshift evolution of the mass averaged neutral fraction for Model I from C2-ray (solid curve) and grizzly (dashed curve). Top panel: Same but showing the redshift evolution of the ratio of the volume and mass averaged neutral fractions. Figure 2. View largeDownload slide Bottom panel: Redshift evolution of the mass averaged neutral fraction for Model I from C2-ray (solid curve) and grizzly (dashed curve). Top panel: Same but showing the redshift evolution of the ratio of the volume and mass averaged neutral fractions. The top panel of Fig. 2 shows the redshift evolution of the ratio of the volume and mass averaged neutral fraction for this reionization history. As the reionization scenario considered here is ‘inside-out’, high-density regions are ionized first. This leads to a higher volume averaged neutral fraction than the mass averaged neutral fraction. Also for this ratio we find very good agreement between the two codes. This suggest that the ionization maps from these two schemes will be very similar. The small differences may originate from the approximations that go into the implementation of the ionization and differences in the treatment of recombinations in the two schemes. This agreement is better than for any of the seminumerical results from Majumdar et al. (2014), see fig. 1 in that paper. To illustrate the similarity between the reionization histories from these two methods even further, we show the map of redshifts of reionization of each pixel in a chosen slice from the simulation boxes in Fig. 3. We define the redshift of reionization of a region to be the redshift when the ionization fraction first reaches 0.5. The maps in Fig. 3 from grizzly and C2-ray are visually quite similar. The same slice for the seminumerical methods is shown in fig. 2 of Majumdar et al. (2014) and shows larger differences. To quantify the similarity of these two maps, we use the Pearson cross-correlation coefficient defined as,   \begin{equation} \chi _{ab}=\frac{\sum _{i}{(a_i-\bar{a})(b_i-\bar{b})}}{\sqrt{\sum _{i}(a_i-\bar{a})^2}\sqrt{\sum _{i}(b_i-\bar{b})^2}}, \end{equation} (6)where $$\bar{a}$$ and $$\bar{b}$$ represent the mean of the maps a and b respectively and i is the pixel index. The Pearson-cross-correlation coefficient calculated for the maps at the top and bottom panels of Fig. 3 is 0.97. Figure 3. View largeDownload slide Top panel: Image showing the redshifts of reionization for pixels in a slice from the simulation volume for Model I generated using the results from C2-ray. Bottom panel: Same as the top panel but generated from the grizzly results. Figure 3. View largeDownload slide Top panel: Image showing the redshifts of reionization for pixels in a slice from the simulation volume for Model I generated using the results from C2-ray. Bottom panel: Same as the top panel but generated from the grizzly results. To compare two different maps, we define the cross-correlation coefficient of a quantity Q ($$R^{AB}_Q$$) as given below:   \begin{equation} R^{AB}_{Q} (k) = \frac{P_{AB}(k)}{\sqrt{P_{AA}(k)P_{BB}(k)}}, \end{equation} (7)where PAA and PBB are the auto power spectrum of the fields A and B respectively and the PAB is the cross power spectrum between them. Here, we will consider A and B as grizzly and C2-ray, respectively. The quantity Q could be the redshift of reionization, neutral fraction, brightness temperatures, etc. The cross-correlation coefficient for the redshift of reionization cubes Rz is found to be larger than 0.95 for scales with k ≲ 1 Mpc− 1 (see top panel of Fig. 4). This suggests that the reionization histories from these two simulations are very similar at large scales. Although the correlation drops at scales correspond to k ≳ 1 Mpc− 1, it is clearly still better than the seminumerical schemes considered in Majumdar et al. (2014) (see fig. 3 of their paper). Figure 4. View largeDownload slide Top panel: Scale dependence of the cross-correlation coefficient of the redshift of reionization maps from C2-ray and grizzly simulations for Model I. Bottom panel: Same as the top panel but for the neutral fraction maps at different stages of reionization. Figure 4. View largeDownload slide Top panel: Scale dependence of the cross-correlation coefficient of the redshift of reionization maps from C2-ray and grizzly simulations for Model I. Bottom panel: Same as the top panel but for the neutral fraction maps at different stages of reionization. 3.2.2 Morphology of H ii regions Next we compare the morphology of the ionized regions from the two simulations using maps of the neutral fraction. We show the same slice from the neutral fraction cubes from C2-ray and grizzly in the two panels of Fig. 5. Visually, the maps are quite similar. The Pearson cross-correlation coefficient calculated for these two slices is 0.93. One can notice that the ionized bubbles in the slice from grizzly simulation appear to be more spherical than those from C2-ray. This is due to the application of 1D radiative transfer which erases the directional dependence of the H ii bubbles. However, the shapes of the bubbles become irregular once there is a significant amount of overlap between the H ii regions. Figure 5. View largeDownload slide Images of the neutral fraction for the same slice through the simulation box as in Fig. 3. Black is ionized, white is neutral. The top and bottom panels show the results from C2-ray and grizzly, respectively. The results are for Model I, z = 9.026 when the mass averaged neutral fraction is ∼0.5. Figure 5. View largeDownload slide Images of the neutral fraction for the same slice through the simulation box as in Fig. 3. Black is ionized, white is neutral. The top and bottom panels show the results from C2-ray and grizzly, respectively. The results are for Model I, z = 9.026 when the mass averaged neutral fraction is ∼0.5. The scale-dependent cross-correlation coefficients for maps of the neutral fraction at different stages of reionization from grizzly and C2-ray are shown in the bottom panel of Fig. 4. One can see very high correlation at large length scales, which are the ones of prime interest to the first generations of radio telescopes. We find that the cross-correlation coefficient of the neutral maps remains very close to unity (within ∼5 per cent) at scales k ≤ 1 Mpc− 1 until the universe is 50 per cent ionized. There is a consistent shift to smaller k values in the scale at which the correlation coefficient reaches 0.9 as reionization progresses. Although the correlation of the ionization fields drops for small length scales, the values actually remain higher than those for the seminumerical models considered in Majumdar et al. (2014). For example, the cross-correlation coefficients of the neutral fractions between C2-ray and different seminumerical schemes were found to be smaller than ∼0.4 at a scale k ∼ 2 Mpc− 1 at a stage with neutral fraction $$x_{\rm H\, \small {I}} \sim 0.5$$, while the corresponding value for grizzly is ∼0.6. 3.2.3 Size distribution of H ii regions We next compare the H ii BSD from these two simulations. BSDs can provide valuable information about the growth of the ionized regions (Friedrich et al. 2011). However, due to the complex morphologies of the ionized regions, there is no unique way to define the sizes of the irregularly shaped H ii regions in simulation volume. Therefore different size determination methods have been developed and each address different aspects of the size distributions (Giri et al. 2018). In this study, we have use two such methods, namely the Friends-of-friends (FOF; Iliev et al. 2006a) and mean-free-path (MFP; Mesinger & Furlanetto 2007). First, we discuss our bubble size statistics from the FOF method. This algorithm considers two neighbouring cells with ionization fraction larger than 0.5 to be the part of the same ionized bubble and in the process finds the volume of each ionized region. This BSD algorithm focuses on connectivity. As shown in Furlanetto & Oh (2016) the percolative nature of reionization leads to the development of one dominant connected region which contains most of the ionized volume. This large connected region is known as the percolation cluster. The single percolation cluster is invisible in the normalized histogram of the volumes (VdP/dV) as this quantity gives the number of ionized bubbles of volume V. Instead we plot (V/Vion)VdP/dV which shows the fraction of the ionized volume (Vion) contained in regions of volume V. Once the percolation cluster has formed it will have a value close to 1 in this quantity. The top panel of Fig. 6 shows the FOF-BSDs from the two simulations at three different stages of reionization. The volume of the percolation cluster is identical in both the simulations. One can observe that the number and size distribution of smaller ionized regions also match quite well during the early and late stages of reionization. However, we find that grizzly produces a few larger bubbles during the intermediate stages. However, compared to the different seminumerical models considered in Majumdar et al. (2014) (as shown in fig. 5 of that paper), the FOF-BSD from grizzly matches the C2-ray results better. Figure 6. View largeDownload slide Top panel: The volume distribution of the H ii bubbles calculated by the FOF method at three different stages of reionization in Model I. The left to right curves (with colours green, red and blue) correspond to redshift 9.5, 9, and 8.7 when the neutral fractions are 0.7, 0.5, and 0.3, respectively. The solid and dashed curves represent C2-ray and grizzly simulations respectively. Bottom panel: The PDF of the H ii bubble size calculated by the MFP method. The line styles and colours have the same meaning as in the upper panel. Figure 6. View largeDownload slide Top panel: The volume distribution of the H ii bubbles calculated by the FOF method at three different stages of reionization in Model I. The left to right curves (with colours green, red and blue) correspond to redshift 9.5, 9, and 8.7 when the neutral fractions are 0.7, 0.5, and 0.3, respectively. The solid and dashed curves represent C2-ray and grizzly simulations respectively. Bottom panel: The PDF of the H ii bubble size calculated by the MFP method. The line styles and colours have the same meaning as in the upper panel. The second BSD method that we use is MFP which is built on a Monte Carlo inference of the sizes. This method selects random cells in the ionized regions and traces a ray in a random direction until it hits a neutral cell. In this case, we consider a cell to be ionized if its ionization fraction is larger than 0.5 and neutral if it is less than that. We repeat the ray tracing for enough number of times (in this case 107 times) and record the lengths of the rays. This provides the estimate of the size distribution of the ionized regions, here expressed in radii instead of volumes. The MFP-BSDs for grizzly and C2-ray show good agreement as can be seen in the bottom panel of Fig. 6. The peaks of the curves at the different stages of reionization appear at the same radius which suggests that the most probable sizes of ionized regions in both simulations are identical. Majumdar et al. (2014) did not calculate MFP-BSDs and therefore we cannot compare to their results in this case. 3.2.4 Differential brightness temperature maps The next quantity that we compare for the two simulations is the expected 21-cm signal. Fig. 7 presents δTb maps of the same slice at three different redshifts 9.5, 9, and 8.7 which correspond to the neutral fractions 0.7, 0.5, and 0.3, respectively. Visually the δTb maps are quite similar. The Pearson-cross-correlation coefficients for these maps are 0.92, 0.9, and 0.9, respectively. However, the small-scale features are different due to the spherical nature of the bubbles around isolated sources in the 1D radiative transfer scheme. This was also seen in Section 3.2.2 when we considered the morphology of ionized regions. Figure 7. View largeDownload slide The brightness temperature maps of the same slice as in Fig. 5 for Model I. Top and bottom panels show the C2-ray and grizzly results, respectively. The left-hand to right-hand panels correspond to redshift 9.5, 9, and 8.7 with neutral fraction of 0.7, 0.5, and 0.3, respectively. Figure 7. View largeDownload slide The brightness temperature maps of the same slice as in Fig. 5 for Model I. Top and bottom panels show the C2-ray and grizzly results, respectively. The left-hand to right-hand panels correspond to redshift 9.5, 9, and 8.7 with neutral fraction of 0.7, 0.5, and 0.3, respectively. To check the level of similarity, we again use the scale-dependent cross-correlation between the two δTb maps, see Fig. 8. The large values of the cross-correlation coefficients $$R_{\delta T_{\rm b}}$$ suggest that the δTb maps from these two schemes are very similar. However, we do find that the value of $$R_{\delta T_{\rm b}}$$ decreases at small scales which indicates mismatches between the δTb maps at small scales. Still, $$R_{\delta T_{\rm b}}$$ is larger than 0.7 at scales such as k ≥ 1 Mpc− 1. This suggests that δTb maps from grizzly are more similar to the maps from C2-ray than the maps from the seminumerical models considered in Majumdar et al. (2014). For example, $$R_{\delta T_{\rm b}}$$ is less than 0.7 at k ≥ 1 Mpc− 1 for the seminumerical models used in Majumdar et al. (2014) while $$R_{\delta T_{\rm b}} \sim 0.85$$ for grizzly. We also find the $$R_{\delta T_{\rm b}}$$ value to be larger than the $$R_{x_{\rm H\, \small {I}}}$$ at scales k ≥ 1 Mpc− 1 (see Fig. 4). This is due to the fact that the δTb signal is dominated by the density fluctuation at small scales and these are identical between the two simulations. Figure 8. View largeDownload slide Scale dependence of the cross-correlation coefficient of the brightness temperature maps from C2-ray and grizzly for model I. Different curves represent redshifts 11.8, 9.5, 9, and 8.7 with neutral fraction 0.9, 0.7, 0.5, and 0.3, respectively. Figure 8. View largeDownload slide Scale dependence of the cross-correlation coefficient of the brightness temperature maps from C2-ray and grizzly for model I. Different curves represent redshifts 11.8, 9.5, 9, and 8.7 with neutral fraction 0.9, 0.7, 0.5, and 0.3, respectively. 3.2.5 Spherically averaged power spectrum of the signal Next we compare the outputs of these two schemes in terms of the power spectrum of the δTb fluctuations. The spherically averaged power spectrum P(k) of the δTb fluctuations is defined as   \begin{equation} \langle \hat{\delta T_{\rm b}}({\boldsymbol k}) \hat{\delta T_{\rm b}}^{\star }({\boldsymbol k^{\prime }})\rangle = (2 \pi )^3 \delta _D({\boldsymbol k - \boldsymbol k^{\prime }}) P(k), \end{equation} (8)where $$\hat{\delta T_{\rm b}}(\boldsymbol k)$$ represents the Fourier transform of $$\delta T_{\rm b}({\boldsymbol x})$$. Here, We present the results in terms of the spherically averaged dimensionless power spectrum Δ2(k) = k3P(k)/2π2. The top panel of Fig. 9 presents the scale dependence of Δ2(k) estimated from these two simulations at different stages of reionization, together with the ratio of the two results. While the power spectrum is dominated by the density fluctuations at the early stages of reionization, it starts to deviate from the background dark matter density power spectrum as reionization progress. The power spectra from these two schemes agree with each other to within 10 per cent for most stages and scales. Only at the largest scales for $$x_\mathrm{H\,{\small {I}}}=0.7$$ the difference exceeds 20 per cent. The overall features of the power spectra are very similar between the two codes. Figure 9. View largeDownload slide Top panel: Scale dependence of the spherically averaged dimensionless power spectrum of δTb fluctuations at redshifts 11.8, 9.5, 9, and 8.7 with neutral fraction 0.9, 0.7, 0.5, and 0.3, respectively. The power spectrum corresponds to the Model I of the paper which considers ionizing photons from the massive haloes only. Bottom panel: Redshift evolution of the large-scale (k = 0.1 Mpc− 1) spherically averaged power spectrum estimated from C2-ray and grizzly which use 3D and 1D radiative transfer, respectively. The ratios of the power spectrum from C2-ray and grizzly are also shown in the upper parts of the panels. $${\Delta }^2_{\rm C}$$ and $${\Delta }^2_{\rm G}$$ represent the dimensionless power spectrum from C2-ray and grizzly, respectively. Figure 9. View largeDownload slide Top panel: Scale dependence of the spherically averaged dimensionless power spectrum of δTb fluctuations at redshifts 11.8, 9.5, 9, and 8.7 with neutral fraction 0.9, 0.7, 0.5, and 0.3, respectively. The power spectrum corresponds to the Model I of the paper which considers ionizing photons from the massive haloes only. Bottom panel: Redshift evolution of the large-scale (k = 0.1 Mpc− 1) spherically averaged power spectrum estimated from C2-ray and grizzly which use 3D and 1D radiative transfer, respectively. The ratios of the power spectrum from C2-ray and grizzly are also shown in the upper parts of the panels. $${\Delta }^2_{\rm C}$$ and $${\Delta }^2_{\rm G}$$ represent the dimensionless power spectrum from C2-ray and grizzly, respectively. The bottom panel of the figure compares the redshift evolution for the k ∼ 0.1 Mpc− 1 mode of the power spectrum. The ratio between the C2-ray and grizzly results starts to deviate from unity at redshift 14 and remains between 0.9 and 1.1 during the later stages of reionization once again indicating excellent agreement. Compared to the seminumerical methods from Majumdar et al. (2014) the grizzly results match the C2-ray ones equally well, or better. 3.3 Model II : Including low-mass haloes and thermal feedback Finally we consider a more complex model for reionization and study the similarities of the outcomes of these two codes. While our fiducial model does not incorporate the contributions from LMACHs, these low-mass haloes can provide a significant contribution to the overall ionizing photon budget. As explained in Section 2.2, in Model II we take the minimum mass of haloes which contribute to reionization to be 108 M⊙ and we suppress star formation in haloes of masses (below 109 M⊙) in ionized regions. This scenario is same as the L1 model in Iliev et al. (2012). The details can be found in Table 1. Majumdar et al. (2014) did not study this model in their comparison between seminumerical and numerical methods, so unlike for the results of Model I we are here not able to compare to seminumerical methods. Although grizzly has to apply simplifying approximations in order to include thermal feedback, the reionization history it generates is very similar to that from C2-ray as shown in the bottom panel of Fig. 10. The Thomson scattering optical depths for C2-ray and grizzly are 0.08 and 0.081, respectively. The redshift evolution of the ratio of the volume and mass averaged neutral fraction from these two simulations also matches very well as shown in the top panel of the figure. This suggests that the neutral fraction maps from these two simulations for this complex model are also very similar. We do not show these here but below will show the 21-cm images from which the morphology of the ionized regions can also be seen. Figure 10. View largeDownload slide Bottom panel: Redshift evolution of the mass averaged neutral fraction from C2-ray and grizzly for the reionization scenario Model II which includes the contributions from low mass haloes and considers thermal feedback. Top panel: Redshift evolution of the ratio of the volume and mass averaged neutral fraction from both the simulations. Figure 10. View largeDownload slide Bottom panel: Redshift evolution of the mass averaged neutral fraction from C2-ray and grizzly for the reionization scenario Model II which includes the contributions from low mass haloes and considers thermal feedback. Top panel: Redshift evolution of the ratio of the volume and mass averaged neutral fraction from both the simulations. The FOF-BSDs obtained for Model II as shown in the top panel of Fig. 11 have mostly very similar shapes. At the early stages (with neutral fraction 0.7), grizzly is seen to produce more large ionized regions, although this difference is likely due to small size differences among the larger regions. However, this possible bias vanishes once the H ii regions start percolating during the middle and later stages of reionization. The distributions of the sizes of the H ii regions calculated using the MFP method for this scenario also display good agreement (as shown in the bottom panel of Fig. 11). We find that there is no prominent peak feature in the MFP-BSD curves during the early stages of reionization. This is due to the fact that the size of the ionized regions are very small and mostly fall within the resolution limit of the simulation. Figure 11. View largeDownload slide Top panel: The volume distribution of the H ii bubbles as a function of volume at different stages of reionization for the two different simulations considered in this work. The left to right curves (with colours green, red, and blue) correspond to ionization fraction 0.7, 0.5, and 0.3 at redshift 10.5, 9.5, and 8.9, respectively. The solid and dashed curves correspond to the C2-ray and grizzly simulation, respectively. Bottom panel: The PDF of the radius of the H ii bubbles as a function of the radius of the bubbles. The BSD in this panel is estimated using the MFP method. The reionization scenario considered here is Model II. Figure 11. View largeDownload slide Top panel: The volume distribution of the H ii bubbles as a function of volume at different stages of reionization for the two different simulations considered in this work. The left to right curves (with colours green, red, and blue) correspond to ionization fraction 0.7, 0.5, and 0.3 at redshift 10.5, 9.5, and 8.9, respectively. The solid and dashed curves correspond to the C2-ray and grizzly simulation, respectively. Bottom panel: The PDF of the radius of the H ii bubbles as a function of the radius of the bubbles. The BSD in this panel is estimated using the MFP method. The reionization scenario considered here is Model II. The differential brightness temperature maps from these two simulations also visually very similar at different stages of reionization as shown in Fig. 12. The Pearson-cross-correlation coefficient calculated for the slices shown in Fig. 12 are 0.7, 0.75, and 0.77 at reionization stages with neutral fraction 0.7, 0.5, and 0.3, respectively. The scale-dependent comparison of the δTb cubes from these two simulations is also shown in terms of the cross-correlation coefficient at different stages of reionization in Fig. 13. One can see that the δTb cubes from these two simulations are highly correlated at the large scales with k ≲ 0.5 Mpc− 1 (within 10 per cent error), while the correlation drops at smaller scales. Figure 12. View largeDownload slide Top and bottom panels represent the δTb maps corresponding to C2-ray and grizzly simulations respectively which are based on 3D and 1D radiative transfer schemes. The left-hand to right-hand panels represent redshift 10.5, 9.5, and 8.9, respectively, which correspond to neutral fraction 0.7, 0.5, and 0.3, respectively. The maps correspond to Model II which include the contributions from low-mass haloes and consider thermal feedback. Figure 12. View largeDownload slide Top and bottom panels represent the δTb maps corresponding to C2-ray and grizzly simulations respectively which are based on 3D and 1D radiative transfer schemes. The left-hand to right-hand panels represent redshift 10.5, 9.5, and 8.9, respectively, which correspond to neutral fraction 0.7, 0.5, and 0.3, respectively. The maps correspond to Model II which include the contributions from low-mass haloes and consider thermal feedback. Figure 13. View largeDownload slide Scale dependence of the cross-correlation coefficient of the brightness temperature maps from C2-ray and grizzly at different stages of reionization for Model II. Figure 13. View largeDownload slide Scale dependence of the cross-correlation coefficient of the brightness temperature maps from C2-ray and grizzly at different stages of reionization for Model II. We present the comparison between the spherically averaged power spectrum from these two simulations of Model II in Fig. 14. The redshift evolution of the spherically averaged power spectrum at scale k ∼ 0.1 Mpc− 1 are quite similar to each other as shown in the bottom panel of Fig. 14. Although the power spectrum estimated from these two schemes are very similar, there are small differences too. For example, we find that the power spectra deviate from each other around redshift 10.5 when the neutral fraction is 0.7. Up to this time the LMACHs contribute substantially to reionization, so the differences are most likely connected to the way their suppression is implemented. Still, the differences are small and except at the very largest scales do not exceed 10 per cent. Figure 14. View largeDownload slide Top panel: The spherically averaged power spectrum of the 21-cm signal for Model II as a function of scale at different stages of reionization. Different curves correspond to redshift 13.5, 10.5, 9.5, and 8.9 with neutral fraction 0.9, 0.7, 0.5, and 0.3, respectively. Bottom panel: δTb power spectrum for k = 0.1 Mpc− 1 as a function of redshift. The solid and dashed curves in both the panels represent C2-ray and grizzly simulations, respectively. Upper parts of the panels show the ratio of the power spectrum from C2-ray and grizzly simulations. Here, $${\Delta }^2_{\rm C}$$ and $${\Delta }^2_{\rm G}$$ represent the dimensionless power spectrum from C2-ray and grizzly, respectively. Figure 14. View largeDownload slide Top panel: The spherically averaged power spectrum of the 21-cm signal for Model II as a function of scale at different stages of reionization. Different curves correspond to redshift 13.5, 10.5, 9.5, and 8.9 with neutral fraction 0.9, 0.7, 0.5, and 0.3, respectively. Bottom panel: δTb power spectrum for k = 0.1 Mpc− 1 as a function of redshift. The solid and dashed curves in both the panels represent C2-ray and grizzly simulations, respectively. Upper parts of the panels show the ratio of the power spectrum from C2-ray and grizzly simulations. Here, $${\Delta }^2_{\rm C}$$ and $${\Delta }^2_{\rm G}$$ represent the dimensionless power spectrum from C2-ray and grizzly, respectively. When comparing these results for Model II to the ones for Model I we see that grizzly performs somewhat less well in the former case. However, the results are still very close to the numerical ones and definitely fall within the expected measurement errors from the observations. Overall, we find that the results from these two simulations agree well with each other even for this more complex scenario. This suggest that the 1D radiative transfer can be used for the complex scenarios as well without a substantial loss of accuracy. Specifically, the 3D and 1D radiative transfer schemes match quite well at scales k ≲ 1 Mpc−1 which are the prime target of the current radio interferometers. This suggests that 1D radiative transfer codes such as bears or grizzly can be efficiently used for predicting the expected signal for various source models as well as to find new observation strategies and possibly for parameter estimation. 4 SUMMARY AND DISCUSSION In this paper, we have compared the performance of the 1D radiative transfer code grizzly to the 3D radiative transfer code C2-ray in the context of a large-scale reionization simulation. While C2-ray simulates the process more accurately, by performing both the full 3D radiative transfer and following the evolution of the ionization fractions and the source in time, this accuracy comes at a considerable computational cost. The grizzly code simplifies the calculation by placing pre-generated 1D profiles of ionization fraction around the sources and using the halo growth history to approximate the evolutionary effects, which makes it approximately a factor ∼105 faster. We have used same sets of initial conditions, i.e. the same density fields, halo catalogues and source models for both simulations. We limit our comparison to the calculation of the ionization fraction. Beside comparing results around a single source, we consider two models to compare the outcomes of these simulations in detail. In Model I (also our fiducial model), we assume that the reionization is driven by massive haloes with masses larger than 2.2 × 109 M⊙. This model was also used in Majumdar et al. (2014) to evaluate the performance of seminumerical reionization codes which rely on the excursion set formalism. We also consider a more complex model (Model II) which includes low mass haloes (down to 108 M⊙) and suppression of star formation in haloes with mass smaller than 109 M⊙ due to thermal feedback. For all the models considered in this paper, we have used the same SED which in this case is a blackbody spectrum. When calculating the differential brightness temperature of the 21-cm signal we assume that TS ≫ Tγ. Our findings from this comparison are the following. When comparing the reionization history due to an isolated source, we observe that grizzly underestimates the effect of recombinations. Since this is a cumulative effect, the impact is largest at the later stages. Still, we find that the differences between the neutral fractions estimated from the two methods for the isolated sources remain less than 10 per cent. For our fiducial simulation, the reionization histories in terms of the evolution of the mass averaged neutral fraction and ionization maps are very similar between both simulations. The cross-correlation coefficient of the neutral fraction maps remains close to unity (within 5–10 per cent) at scales k ≲ 1 Mpc− 1 before the universe reaches 50 per cent ionization and somewhat worse after this. We have considered two different methods to characterize the size distribution of ionized regions, namely the FOF and MFP methods. For both methods the size distributions match very well between the C2-ray and grizzly results. However, visual inspection of the ionization fraction maps does show that grizzly produces more spherically shaped regions compared to C2-ray. Visually, the differential brightness temperature maps from the two simulations look very similar. The cross-correlation coefficients at large scales, k ≲ 0.1 Mpc− 1 are close to unity. We find that the correlation coefficients at small scales are higher for the δTb maps than for the ionization maps. This is due to the density fluctuations in the neutral medium. The spherically averaged power spectra from the two simulations agree with each other within 10 per cent. For the more complex scenario Model II, we find that the differences between the grizzly and C2-ray are generally somewhat larger than for the simpler Model I but remain small. We conclude that grizzly is capable of handling more complex source models where sources are not always ‘on’. The results from grizzly are more similar to those from C2-ray compared to the different seminumerical simulations considered in Majumdar et al. (2014). We find for example that the various cross-correlation coefficients are higher for grizzly, especially for smaller scales. Furthermore, for the seminumerical models the source efficiency was adapted at each redshift in order to reproduce the ionization history from C2-ray, whereas grizzly is able to produce an excellent match for the reionization history without varying the source efficiency parameter. These results imply that the faster 1D radiative transfer codes such as bears or grizzly can be efficiently used for exploring parameter space, determining observational strategies or developing parameter estimation pipelines. Once the 1D profiles have been generated, these simulations are fast enough to be used for parameter estimation using techniques such as the MCMCs, machine learning, etc. However, the details of parameter estimation such as convergence time, accuracy, etc. are beyond the scope of this work and will be addressed in future works. Although grizzly is quite fast compared to the full radiative transfer simulations and generally reproduces the results well, our comparison also shows where improvements could be made. (i) grizzly shows lower accuracy at small scales. We speculate that the accuracy at small scales can be improved by optimizing the technique grizzly uses to deal with overlapping ionized regions. (ii) The comparison for a single source shows that grizzly underestimates the effect of recombinations. This could be improved by taking into account the cosmological evolution of the density field over the lifetime of the source. (iii) The halo suppression model used in grizzly is currently a very simple one. It may be possible to implement the evolution of relic H ii regions better and thereby to improve the accuracy of grizzly for scenarios where source can turn ‘off’. Besides the neutral fraction maps, grizzly can also provide the kinetic temperature and Ly α flux maps. Thus, it can be used to predict signal in presence of spin-temperature fluctuations. The 21-cm signal from the Cosmic Dawn is expected to be dominated by these fluctuations. Thus, grizzly is also very relevant for observations with interferometers which also probe this earlier era, such as SKA1-low. grizzly can for example be used to constrain the heating and Ly α coupling parameters during the Cosmic Dawn from the observations. However, a detailed comparison of the results in presence of TS fluctuations with fully numerical simulations such as Ross et al. (2017) is beyond the scope of this paper and will be addressed in the future. Acknowledgements The authors would like to thank Saleem Zaroubi and Adi Nusser for providing useful suggestions on this work. KKD would like to thank the University Grant Commission (UGC), India for support through UGC-faculty recharge scheme (UGC-FRP) vide ref.no. F.4-5(137-FRP)/2014(BSR). SM acknowledges the financial support from the European Research Council under ERC grant number 638743-FIRSTDAWN. Footnotes 1 http://www.gmrt.tifr.res.in 2 http://eor.berkeley.edu/ 3 http://www.mwatelescope.org/ 4 http://www.lofar.org/ 5 http://www.skatelescope.org/ 6 http://wiki.cita.utoronto.ca/mediawiki/index.php/CubePM 7 The suppression of the low-mass haloes is determined by the kinetic temperature of the region. For example, one can assume no star formation within haloes with masses <109 M⊙ if it is formed within a region with TK larger than 104 K. Following Iliev et al. (2012) we use $$x_{\rm H\, \small {II}}$$ as proxy for the temperature and assume that if a region is more ionized than 0.1, it will be hot enough to suppress further star formation. 8 As mentioned above, grizzly does not track the ionization history to create the ionization map at a certain redshift for simple models where thermal feedback is ignored. Thus, the run time will be much shorter if only a single snapshot of the ionization fields is needed. For example, grizzly generates the neutral fraction map at redshift 9.026 with mass averaged neutral fraction ∼0.5 in ∼10 CPU minutes. In presence of thermal feedback, the estimated time is higher which is also true for other codes such as the 21 cm FAST (see e.g. Sobacchi & Mesinger 2014). The main computational bottleneck for grizzly is to estimate the unused ionizing photons in the overlapping regions and to redistribute these around these overlapping regions. This step is executed iteratively and there are no straightforward methods to compute this more efficiently. REFERENCES Alvarez M. A., Wise J. H., Abel T., 2009a, ApJ , 701, L133 CrossRef Search ADS   Alvarez M. A., Busha M., Abel T., Wechsler R. H., 2009b, ApJ , 703, L167 CrossRef Search ADS   Baek S., Di Matteo P., Semelin B., Combes F., Revaz Y., 2009, A&A , 495, 389 CrossRef Search ADS   Becker R. H. et al.  , 2001, AJ , 122, 2850 CrossRef Search ADS   Bouwens R. J. et al.  , 2015, ApJ , 803, 34 CrossRef Search ADS   Bowman J. D. et al.  , 2013, PASA , 30, e031 CrossRef Search ADS   Bromm V., Loeb A., 2003, ApJ , 596, 34 CrossRef Search ADS   Choudhury T. R., Ferrara A., 2006, MNRAS , 371, L55 CrossRef Search ADS   Choudhury T. R., Haehnelt M. G., Regan J., 2009, MNRAS , 394, 960 CrossRef Search ADS   Choudhury T. R., Puchwein E., Haehnelt M. G., Bolton J. S., 2015, MNRAS , 452, 261 CrossRef Search ADS   Ciardi B., Ferrara A., Marri S., Raimondo G., 2001, MNRAS , 324, 381 CrossRef Search ADS   Cohen A., Fialkov A., Barkana R., Lotem M., 2017, MNRAS , 472, 1915 CrossRef Search ADS   Dijkstra M., Mesinger A., Wyithe J. S. B., 2011, MNRAS , 414, 2139 CrossRef Search ADS   Fan X. et al.  , 2006, AJ , 131, 1203 CrossRef Search ADS   Fialkov A., Barkana R., Visbal E., 2014, Nature , 506, 197 CrossRef Search ADS PubMed  Friedrich M. M., Mellema G., Alvarez M. A., Shapiro P. R., Iliev I. T., 2011, MNRAS , 413, 1353 CrossRef Search ADS   Furlanetto S. R., Oh S. P., 2016, MNRAS , 457, 1813 CrossRef Search ADS   Furlanetto S. R., Zaldarriaga M., Hernquist L., 2004, ApJ , 613, 16 CrossRef Search ADS   Furlanetto S. R., Oh S. P., Briggs F. H., 2006, Phys. Rep. , 433, 181 CrossRef Search ADS   Ghara R., Choudhury T. R., Datta K. K., 2015, MNRAS , 447, 1806 CrossRef Search ADS   Ghara R., Choudhury T. R., Datta K. K., Choudhuri S., 2017, MNRAS , 464, 2234 CrossRef Search ADS   Ghosh A., Prasad J., Bharadwaj S., Ali S. S., Chengalur J. N., 2012, MNRAS , 426, 3295 CrossRef Search ADS   Giri S. K., Mellema G., Dixon K. L., Iliev I. T., 2018, MNRAS , 473, 2949 CrossRef Search ADS   Goto T., Utsumi Y., Hattori T., Miyazaki S., Yamauchi C., 2011, MNRAS , 415, L1 CrossRef Search ADS   Greig B., Mesinger A., 2015, MNRAS , 449, 4246 CrossRef Search ADS   Harnois-Déraps J., Pen U.-L., Iliev I. T., Merz H., Emberson J. D., Desjacques V., 2013, MNRAS , 436, 540 CrossRef Search ADS   Hinshaw G. et al.  , 2013, ApJS , 208, 19 CrossRef Search ADS   Iliev I. T., Mellema G., Pen U.-L., Merz H., Shapiro P. R., Alvarez M. A., 2006a, MNRAS , 369, 1625 CrossRef Search ADS   Iliev I. T. et al.  , 2006b, MNRAS , 371, 1057 CrossRef Search ADS   Iliev I. T. et al.  , 2009, MNRAS , 400, 1283 CrossRef Search ADS   Iliev I. T., Mellema G., Shapiro P. R., Pen U.-L., Mao Y., Koda J., Ahn K., 2012, MNRAS , 423, 2222 CrossRef Search ADS   Kaaret P., 2014, MNRAS , 440, L26 CrossRef Search ADS   Kashikawa N. et al.  , 2011, ApJ , 734, 119 CrossRef Search ADS   Knevitt G., Wynn G. A., Power C., Bolton J. S., 2014, MNRAS , 445, 2034 CrossRef Search ADS   Komatsu E. et al.  , 2011, ApJS , 192, 18 CrossRef Search ADS   Majumdar S., Mellema G., Datta K. K., Jensen H., Choudhury T. R., Bharadwaj S., Friedrich M. M., 2014, MNRAS , 443, 2843 CrossRef Search ADS   Malhotra S., Rhoads J. E., 2006, ApJ , 647, L95 CrossRef Search ADS   McGreer I. D., Mesinger A., D'Odorico V., 2015, MNRAS , 447, 499 CrossRef Search ADS   McQuinn M., Lidz A., Zahn O., Dutta S., Hernquist L., Zaldarriaga M., 2007a, MNRAS , 377, 1043 CrossRef Search ADS   McQuinn M., Hernquist L., Zaldarriaga M., Dutta S., 2007b, MNRAS , 381, 75 CrossRef Search ADS   Mellema G., Iliev I. T., Pen U.-L., Shapiro P. R., 2006, MNRAS , 372, 679 CrossRef Search ADS   Mellema G., Koopmans L., Shukla H., Datta K. K., Mesinger A., Majumdar S., 2015, HI tomographic imaging of the Cosmic Dawn and Epoch of Reionization with SKA . SISSA, Trieste, PoS(AASKA14)010 Mesinger A., Furlanetto S., 2007, ApJ , 669, 663 CrossRef Search ADS   Mesinger A., Furlanetto S., Cen R., 2011, MNRAS , 411, 955 CrossRef Search ADS   Mitra S., Choudhury T. R., Ferrara A., 2011, MNRAS , 413, 1569 CrossRef Search ADS   Mitra S., Choudhury T. R., Ferrara A., 2012, MNRAS , 419, 1480 CrossRef Search ADS   Mitra S., Choudhury T. R., Ferrara A., 2015, MNRAS , 454, L76 CrossRef Search ADS   Morales M. F., Wyithe J. S. B., 2010, ARA&A , 48, 127 CrossRef Search ADS   Mortlock D. J. et al.  , 2011, Nature , 474, 616 CrossRef Search ADS PubMed  Paciga G. et al.  , 2013, MNRAS , 433, 639 CrossRef Search ADS   Paranjape A., Choudhury T. R., 2014, MNRAS , 442, 1470 CrossRef Search ADS   Parsons A. R. et al.  , 2014, ApJ , 788, 106 CrossRef Search ADS   Patil A. H. et al.  , 2014, MNRAS , 443, 1113 CrossRef Search ADS   Patil A. H. et al.  , 2017, ApJ , 838, 65 CrossRef Search ADS   Planck Collaboration XVI, 2014, A&A , 571, A16 CrossRef Search ADS   Planck Collaboration XIII, 2016a, A&A , 594, A13 CrossRef Search ADS   Planck Collaboration XV, 2016b, A&A , 594, A15 CrossRef Search ADS   Pritchard J. R., Loeb A., 2012, Rep. Prog. Phys. , 75, 086901 CrossRef Search ADS PubMed  Ross H. E., Dixon K. L., Iliev I. T., Mellema G., 2017, MNRAS , 468, 3785 CrossRef Search ADS   Santos M. G., Amblard A., Pritchard J., Trac H., Cen R., Cooray A., 2008, ApJ , 689, 1 CrossRef Search ADS   Santos M. G., Ferramacho L., Silva M. B., Amblard A., Cooray A., 2010, MNRAS , 406, 2421 CrossRef Search ADS   Shimabukuro H., Semelin B., 2017, MNRAS , 468, 3869 CrossRef Search ADS   Shin M.-S., Trac H., Cen R., 2008, ApJ , 681, 756 CrossRef Search ADS   Sobacchi E., Mesinger A., 2014, MNRAS , 440, 1662 CrossRef Search ADS   Tanaka T., Perna R., Haiman Z., 2012, MNRAS , 425, 2974 CrossRef Search ADS   Thomas R. M. et al.  , 2009, MNRAS , 393, 32 CrossRef Search ADS   Tingay S. J. et al.  , 2013, PASA , 30, 7 CrossRef Search ADS   van Haarlem M. P. et al.  , 2013, A&A , 556, A2 CrossRef Search ADS   Wechsler R. H., Bullock J. S., Primack J. R., Kravtsov A. V., Dekel A., 2002, ApJ , 568, 52 CrossRef Search ADS   Wise J. H., Demchenko V. G., Halicek M. T., Norman M. L., Turk M. J., Abel T., Smith B. D., 2014, MNRAS , 442, 2560 CrossRef Search ADS   Xu H., Wise J. H., Norman M. L., Ahn K., O'Shea B. W., 2016, ApJ , 833, 84 CrossRef Search ADS   Zahn O., Lidz A., McQuinn M., Dutta S., Hernquist L., Zaldarriaga M., Furlanetto S. R., 2007, ApJ , 654, 12 CrossRef Search ADS   Zahn O., Mesinger A., McQuinn M., Trac H., Cen R., Hernquist L. E., 2011, MNRAS , 414, 727 CrossRef Search ADS   Zaldarriaga M., Seljak U., 1997, Phys. Rev. D , 55, 1830 CrossRef Search ADS   Zaroubi S., Thomas R. M., Sugiyama N., Silk J., 2007, MNRAS , 375, 1269 CrossRef Search ADS   © 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

Prediction of the 21-cm signal from reionization: comparison between 3D and 1D radiative transfer schemes

Loading next page...
 
/lp/ou_press/prediction-of-the-21-cm-signal-from-reionization-comparison-between-3d-RXV03So3RW
Publisher
The Royal Astronomical Society
Copyright
© 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society
ISSN
0035-8711
eISSN
1365-2966
D.O.I.
10.1093/mnras/sty314
Publisher site
See Article on Publisher Site

Abstract

Abstract Three-dimensional radiative transfer simulations of the epoch of reionization can produce realistic results, but are computationally expensive. On the other hand, simulations relying on one-dimensional radiative transfer solutions are faster but limited in accuracy due to their more approximate nature. Here, we compare the performance of the reionization simulation codes grizzly and C2-ray which use 1D and 3D radiative transfer schemes, respectively. The comparison is performed using the same cosmological density fields, halo catalogues, and source properties. We find that the ionization maps, as well as the 21-cm signal maps from these two simulations are very similar even for complex scenarios which include thermal feedback on low-mass haloes. The comparison between the schemes in terms of the statistical quantities such as the power spectrum of the brightness temperature fluctuation agrees with each other within 10 per cent error throughout the entire reionization history. grizzly seems to perform slightly better than the seminumerical approaches considered in Majumdar et al. which are based on the excursion set principle. We argue that grizzly can be efficiently used for exploring parameter space, establishing observations strategies, and estimating parameters from 21-cm observations. radiative transfer, galaxies: formation, intergalactic medium, dark ages, reionization, first stars, cosmology: theory 1 INTRODUCTION The epoch of reionization (EoR), when the primordial neutral hydrogen (H i) in the intergalactic medium (IGM) was ionized by the radiation from the first sources, is one of the milestone events in the evolution of the Universe. However, many important facts such as the exact timing of this event, the sources responsible for reionization, the impact of the EoR on the later stages of the structure formation, etc. are currently not or poorly known. Various probes of the EoR such as observations of high-redshift quasars (Becker et al. 2001; Fan et al. 2006; Goto et al. 2011; Mortlock et al. 2011; McGreer, Mesinger & D'Odorico 2015), the cosmic microwave background radiation (CMB) (Zaldarriaga & Seljak 1997; Komatsu et al. 2011; Planck Collaboration XV 2016b; Planck Collaboration XIII 2016a), high-redshift galaxies (Kashikawa et al. 2011; Bouwens et al. 2015; Mitra, Choudhury & Ferrara 2015), and Ly α emitters (McQuinn et al. 2007b; Dijkstra, Mesinger & Wyithe 2011; Choudhury et al. 2015) have been used to constrain some of these unknown facts. However, these probes can only provide limited information. For example, combining the observations of high-redshift quasars and the CMB, one can constrain the probable period of the event to lie between redshifts 15 and 6 (Choudhury & Ferrara 2006; Fan et al. 2006; Malhotra & Rhoads 2006; Mitra, Choudhury & Ferrara 2011, 2012). The evolving redshifted 21-cm signal from the H i in the IGM has the potential to revolutionize our understanding of reionization (Furlanetto, Oh & Briggs 2006; Morales & Wyithe 2010; Pritchard & Loeb 2012). Therefore a range of low frequency radio telescopes such as the Giant Metrewave Radio Telescope1 (Ghosh et al. 2012; Paciga et al. 2013), the Precision Array for Probing the Epoch of Reionization2 (Parsons et al. 2014), the Murchison Widefield Array3 (Bowman et al. 2013; Tingay et al. 2013), and the Low Frequency Array4 (van Haarlem et al. 2013; Patil et al. 2017) are attempting to detect this faint signal from the EoR. The primary aim of these first generation experiments is to detect the signal in terms of statistical quantities such as the variance, skewness, power spectrum, etc. However, future radio telescopes such as the Square Kilometre Array (SKA)5 should have enough sensitivity to produce image cubes of the signal (Mellema et al. 2015; Ghara et al. 2017). A large number of theoretical models based on either analytic calculations (e.g. Furlanetto, Zaldarriaga & Hernquist 2004; Paranjape & Choudhury 2014), seminumerical simulations (Mesinger & Furlanetto 2007; Zahn et al. 2007; Santos et al. 2008; Choudhury, Haehnelt & Regan 2009), or numerical simulations with radiative transfer (Iliev et al. 2006b; Mellema et al. 2006; McQuinn et al. 2007a; Shin, Trac & Cen 2008; Baek et al. 2009; Thomas et al. 2009; Ghara, Choudhury & Datta 2015) have been used to study the redshifted 21-cm signal from the EoR. Through these works a fair amount of theoretical understanding of the impact of various physical processes on the signal has been achieved. Aspects which have been considered include the impact of different kind of sources to the redshifted 21-cm signal, thermal feedback and peculiar velocities of the gas in the IGM. The simulations are not only important to extract information from the observations, but also to design observation strategies. One thing to keep in mind is that the 21-cm observations do not provide direct information about the properties of the sources responsible for ionizing the IGM, minimum mass of the dark matter haloes forming stars, thermal feedback, the effect of peculiar velocities of the gas in the IGM, etc. Due to the complexity of the astrophysical effects, one needs to use simulations to extract this information from the observations. In principle, this can be done using parameter estimation techniques such as the Markov chain Monte Carlo (Greig & Mesinger 2015) or artificial neural networks (Shimabukuro & Semelin 2017) in which the observations expressed in quantities such as power spectrum, variance, or skewness will be compared to the outputs of the simulations. The various simulation methods each have their own advantages and drawbacks and when using their results to design observation strategies or to interpret the observational data, one needs to keep in mind the various assumptions made by the models. While simulations of the 21-cm signal from the EoR using 3D radiative transfer methods such as Iliev et al. (2006b) and Mellema et al. (2006) are assumed to be the most accurate ones, these require supercomputer facilities to run. On the other hand, simulations performed with seminumerical methods are orders of magnitudes faster, the simplifying assumptions needed for this may lead to accuracy issues. Thus, it is absolutely necessary to compare the performance of different simulation techniques with each other so as to identify the advantages and drawbacks of these methods. This will be useful for improving the methods and for correctly interpreting the observations. Studies such as Zahn et al. (2011) compared a reionization simulation obtained using a full 3D radiative transfer technique with the result from a seminumerical one. It showed excellent agreement between these two. Similarly, Majumdar et al. (2014) compared seminumerical simulations and conditional Press-Schechter models with a full 3D radiative transfer simulation performed with C2-ray (Mellema et al. 2006). That study showed that the seminumerical simulations such as Zahn et al. (2007), Mesinger & Furlanetto (2007), Choudhury et al. (2009), and Santos et al. (2010), which use an excursion-set formalism to create the ionized bubbles, produce an ionization history very similar to the C2-ray provided that the halo lists from an N-body simulation are used. However, if the haloes are found from a conditional Press–Schechter approach, as is for example done in Alvarez et al. (2009b), Mesinger, Furlanetto & Cen (2011), the fit with the C2-ray results is considerably worse. The codes such as bears and grizzly use one-dimensional radiative transfer schemes instead of solving radiative transfer equations in 3D. Thus, these simulations assumes that the effect from each sources (i.e., the ionized regions around the sources) are perfectly spherical. When these ionized regions overlap, these codes account the overlap using special recipes. Due to these simplifying assumptions they are fast but can also be expected to be less accurate than full 3D radiative transfer schemes. Previously, Thomas et al. (2009) performed a limited comparison between the 1D radiative transfer code bears and the 3D radiative transfer code crash (Ciardi et al. 2001) for a simulation volume of size 12.5 h−1 comoving megaparsec. However, before considering these 1D schemes for parameter estimations (e.g. Patil et al. 2014), their performance for larger simulation volumes and larger numbers of sources is required. In this paper, we present a detailed comparison between the performance of the 1D radiative transfer code grizzly (Ghara et al. 2015) with the 3D radiative transfer code C2-ray. We use the same initial conditions such as the density fields, halo catalogues, source models, and compare the reionization history produced by these two schemes. We consider bubble size distributions (BSDs) as well as cross-correlations between ionization maps, while for the redshifted 21-cm signal we focus on statistical quantities such as the evolution of the power spectra. We consider two different models: (i) reionization with only large mass haloes which do not suffer from thermal feedback, (ii) reionization with both high- and low-mass haloes where the latter are assumed to be sensitive to thermal feedback. Case (i) is identical to the one used in Majumdar et al. (2014) to compare seminumerical and full radiative transfer results which therefore allows the comparison between seminumerical, 1D RT, and 3D RT results. Our paper is structured in the following way. In Section 2, we describe the simulations used in this paper. In particular different subsections describe the N-body simulation and the source model used in this study as well as a brief description of the underlying algorithms of C2-ray and grizzly. In Section 3, we present our results before we conclude in Section 4. Throughout the paper, we have chosen the cosmological parameters Ωm = 0.27, $$\Omega _\Lambda =0.73$$, ΩB = 0.044, h = 0.7 consistent with the Wilkinson Microwave Anisotropy Probe (WMAP) results (Hinshaw et al. 2013) and within the error bars consistent with Planck (Planck Collaboration XVI 2014). 2 SIMULATIONS Below we discuss the radiative transfer methods along with the source model and N-body simulation used in this study. 2.1 N-body simulation Both radiative transfer methods considered in this study use the results of the same N-body simulation. This dark matter only N-body simulation was carried out using the publicly available code cubep3m6 (Harnois-Déraps et al. 2013). The details of the simulations are (i) the size of the simulation cube is 163 comoving megaparsec (cMpc), (ii) the number of dark matter particles is 30723, (iii) the fine grid used by the particle–particle–mesh method is 61443, (iv) the mass of the dark matter particles is 5.47 × 106 M⊙. The details of the simulation can also be found in Iliev et al. (2012). Although this N-body simulation is by now rather old, we use it to enable direct comparison to the results in Majumdar et al. (2014) who compared the performance of a seminumerical method to C2-ray. Snapshots of the density field are generated from redshift 20.134 to redshift 6 using equal time intervals of 11.5 megayears (Myr). The resolution of these gridded density fields is 2563. cubep3m also produces dark matter halo lists using an on the fly halo finder based on the spherical overdensity method. The lowest mass haloes produced have a mass of ∼108 M⊙ which corresponds approximately to the mass needed to initiate star formation driven by atomic cooling. 2.2 Source model Our knowledge about the sources of reionization is very limited. Different efforts, both observational and modelling, have been undertaken to improve our understanding of the impact of various sources such as primordial galaxies (Wise et al. 2014; Xu et al. 2016), mini-quasars (Bromm & Loeb 2003; Zaroubi et al. 2007; Alvarez, Wise & Abel 2009a; Tanaka, Perna & Haiman 2012), and high-mass X-ray binaries (Fialkov, Barkana & Visbal 2014; Kaaret 2014; Knevitt et al. 2014) on the expected 21-cm signal from the EoR. However, the individual contributions of these sources towards reionization remain quite uncertain. In this study, we assume the sources of ionizing radiation formed in the dark matter haloes. We further assume their spectral energy distribution can be approximated by a blackbody of an effective temperature of 50 000 K, which is appropriate for a massive Population II star. We normalize the spectrum such that the rate of production of the ionizing photons is   \begin{equation} \dot{N_{\gamma }}=g_{\gamma }\frac{M_{\rm h} \Omega _{\rm B}}{(10 \ \rm Myr)\Omega _m m_{\rm p}}, \end{equation} (1)where gγ is the source ionization efficiency coefficient, Mh and mp are the mass of the dark matter halo and proton mass, respectively. For the fiducial model considered here, we choose gγ = 21.7 for the sources formed in dark matter haloes with mass Mh ≥ 2.2 × 109 M⊙, which will produce $$\dot{N_{\gamma }}=1.33 \times 10^{43} \times \frac{M_{\rm h}}{{\rm M}_{{\odot }}} \ \rm s^{-1}$$. This is a similar choice as the L3 model in Iliev et al. (2012). In addition, we consider a more complex reionization scenario where we use all haloes of mass Mh ≥ 108 M⊙. We divide these in two populations, low-mass atomically cooling haloes (LMACHs) in the mass range 108 M⊙ ≤ Mh < 109 M⊙ and high-mass atomically cooling haloes (HMACHs) in the mass range Mh ≥ 109 M⊙. The LMACHs are assumed to be sensitive to thermal feedback and stop producing ionizing photons once they are located in an cell which is more than 10 per cent ionized. The HMACHs are taken to be insensitive to thermal feedback. While still active the LMACHs produce ionizing photons with an efficiency gγ = 130, whereas the HMACHs have gγ = 8.7. The higher efficiency of the LMACHs is motivated by either a larger contribution from Pop III stars or a higher fraction of ionizing photons escaping. In the C2-ray results these parameters result in a rather early end of reionization around z = 8.5 and a Thomson optical depth value of 0.08, consistent with the value determined by WMAP and within the 2σ range of Planck. This scenario is identical to the L1 model of Iliev et al. (2012). The details of these scenarios are given in Table 1. Table 1. Details of the scenarios considered in this paper. LMACH and HMACH represent the low-mass and high-mass atomically cooling halo, respectively. Model  Min. halo (M⊙)  Thermal feedback  gγ  gγ  Model used in Iliev et al. (2012)        LMACH  HMACH    MODEL I (fiducial)  2.2 × 109  No  0  21.7  L3  MODEL II  108  Yes  130  8.7  L1  Model  Min. halo (M⊙)  Thermal feedback  gγ  gγ  Model used in Iliev et al. (2012)        LMACH  HMACH    MODEL I (fiducial)  2.2 × 109  No  0  21.7  L3  MODEL II  108  Yes  130  8.7  L1  View Large 2.3 3D radiative transfer simulation C2-ray The 3D radiative transfer method considered in this work uses ‘Conservative Causal Ray-tracing method’ (C2-ray). The details of the method can be found in Mellema et al. (2006). The method is briefly described in the following steps. It starts with preparing the source list in a random order at each redshift. Given the SED of the sources, the total photoionization rate (Γ) is calculated at time t at each cell in the simulation box including contributions from all the sources. This step takes into account the time evolution of the neutral fraction of the cells during the time-step which affects the optical depth from the sources. The history of the ionization fraction of hydrogen ($$x_{\rm H\, \small {II}}$$) is estimated by solving the non-equilibrium equation,   \begin{equation} \frac{{\rm d}x_{\rm H\, \small {II}}}{{\rm d}t}=\left(1-x_{\rm H\, \small {II}}\right)\left(\Gamma + n_{\rm e}C_{\rm H}\right)-xn_{\rm e}C\alpha _{\rm H}, \end{equation} (2)where ne is the electron density at the cell, CH and αH represent the collisional ionization and recombination coefficients for hydrogen, respectively. The quantity C is the clumping factor (which can be defined as 〈n2〉/〈n〉2, where n is the baryon number density) which accounts the clumpiness of the IGM. Here we have taken C = 1 but in general it can be larger than 1 or even position dependent. C2-ray was compared to other fully numerical radiative transfer codes in Iliev et al. (2006b) and Iliev et al. (2009) and was found to produce accurate results. In this study we do not consider the temperature evolution of the intergalactic gas, only its ionized state. 2.4 1D radiative transfer simulation grizzly The details of the one-dimensional radiative transfer method can be found in Ghara et al. (2015). The method mainly follows the bears algorithm (Thomas et al. 2009) with a few changes compared to the original method. To express its relation to the original bears method, we have chosen to call it grizzly. The basic idea of both methods is not to solve the one-dimensional radiative transfer equations on the fly while generating the ionization maps. Rather, they use pre-generated ionization profiles to construct the ionization maps. Here, we briefly summarize the most important steps of the algorithm. We start by choosing the range for the key source parameters: redshift, density contrast of the uniform background IGM, luminosity of the source, age of the source. We then generate a large number of ionization and kinetic temperature profiles around the sources for different combinations of these parameter values. Next we create a list of radii of the H ii bubbles as a function of different parameters. We assume the radius of the H ii region to be the radius at which the ionization fraction drops to 0.5 from the location of the source. This library of 1D profiles only needs to be compiled once for a given cosmology and spectral energy distribution of sources. For a given redshift and a certain halo position, we determine the luminosity and the age of the halo according to a source recipe, see Section 2.4.1. From this we determine the size of the H ii bubble around it using the density field and the pre-generated list of the H ii bubble radii. This step is done iteratively. We start with a small radius around the centre of the source, calculate the spherically averaged overdensity of the IGM within this radius, and then estimate the radius of H ii bubble from the list. We change the initial choice of radius such that the estimated radius of the H ii region matches with the chosen radius for the same overdensity. The parameter values corresponding to that radius are then associated with the source and the corresponding profiles will be used to generate the ionization and kinetic temperature maps. The procedure described in the previous step is performed for all haloes at the redshift being considered. Although the number of overlaps between the individual H ii regions around the haloes is negligible during the initial stages of reionization, they become important at the later stages of reionization. To incorporate the effect of overlap, we estimate the unused ionizing photons for each overlapping H ii region and distribute them around the overlapping regions such that each photon ionizes one hydrogen atom. Next, we use the pre-generated 1D ionization profiles corresponding to each halo to calculate the ionization fraction at the partially ionized region beyond the H ii regions. The ionization fraction in the overlapping partially ionized regions is then expressed as   \begin{equation} x_{\rm H\, \small {II}}({\boldsymbol x}) =\frac{\sum _i x^{i}_{\rm H\, \small {II}-1D}({\boldsymbol x}) \times n^{i}_{\rm H-1D} \times \left(1- x^{i-1}_{\rm H\, \small {II}}({\boldsymbol x})\right)}{n_{\rm H}({\boldsymbol x})}, \end{equation} (3)where the summation symbol represents a sum over all overlapping sources. Here, $$x^{i-1}_{\rm H\, \small {II}}({\boldsymbol x})$$ denotes the ionization fraction obtained after considering overlapping H ii regions from i − 1 number of sources. We take $$x^{0}_{\rm H\, \small {II}}({\boldsymbol x})=0$$. The term $$\left(1- x^{i-1}_{\rm H\, \small {II}}({\boldsymbol x})\right)$$ takes care of the fact that the ith source encounters an already partially ionized IGM, due to overlap between previous i − 1 number of sources, before contributing to ionization. Next, we generate the kinetic temperature maps using a correlation of the neutral fraction and the gas temperature in the partially ionized region (for details, see Ghara et al. 2015). 2.4.1 Implementation of source model In the model scenarios considered in this paper, we assume that dark matter haloes always have ongoing star formation, unless they are suppressed through thermal feedback, see Section 2.4.2. To select a one-dimensional profile from the library we need to specify a luminosity and age for the halo. The original bears code and Ghara et al. (2015) used a fixed age for the haloes, typically ∼ 10 Myr and also used a luminosity calculated from the instant mass of a halo at the redshift under consideration. Those models assume that the stars inside the galaxies efficiently emit hydrogen ionizing photons until an age of ∼10 Myr. However, we found this source model to be inconsistent with the one used by C2-ray in terms of the cumulative number of ionizing photons being produced. We therefore implemented another source model. In this model we follow the growth of haloes and find both the age of the halo as well as the time-averaged mass (to be called the effective mass Meff). We then calculate the luminosity based on Meff using equation (1) and use this together with the age to find the appropriate profile from the library. The grizzly code thus does not track the history of the state of the IGM gas but rather uses the history of the haloes to calculate the state of the IGM at a given redshift. It is not entirely trivial to calculate the age of the haloes accurately from the halo lists as haloes can move to new cells or multiple haloes can merge. In principle, one needs to track the haloes using a merger-tree algorithm. However, as no merger tree is supplied with the halo lists, we have used a simple method to account for their growth. For every pair of redshifts z1 and z2 we first search for a halo at the same location and it is found, we associate the change in mass with the halo whose history we are tracking. If a halo disappears from its previous grid position in some time-step, we inspect bordering grid cells. If we find a newly formed heavier halo in a bordering grid cells, we assume that the identified halo must be the present state of the halo which disappeared. If we are not able to find any newly formed halo in the neighbourhood, then we assume that a merger happened. In this case, we find previously existing heavier halo at the neighbour and set its age as the maximum of the ages of the identified and disappeared halo. Although this a posteriori reconstruction of the halo histories is rather crude, it does allow us to obtain the correct cumulative number of ionizing photons produced during reionization. 2.4.2 Implementation of thermal feedback In the second C2-ray simulation we include LMACHs which are fully suppressed once their environment has been sufficiently ionized. The criterion used is that when the ionization fraction in the cell containing the LMACH is larger than 0.1, thermal feedback is assumed to stop the production of ionizing photons.7 To implement this recipe in grizzly, it is not trivial as the code only places ionized regions around active sources at the redshift of interest. We implemented the following simple prescription to take into account thermal feedback effects in grizzly. We use the ionization map from the previous time-step to select the active sources at a certain time-step using the condition on the ionization fraction. We add to these the previously active sources which became inactive at the present time-step but assign to them their effective mass and age corresponding to the last time they were active. In principle one needs to consider the recombination in the already ionized regions to correctly account for the suppression. However, it is not straightforward to implement recombination in an already ionized region in grizzly. This is also true for other seminumerical simulations such as models considered in Majumdar et al. (2014) and Sobacchi & Mesinger (2014). However, recombinations play a minor role, especially since previously suppressed LMACH haloes become active again after their haloes have grown into HMACHs and the amount of time that recombinations act upon the relic H ii regions will thus be limited. Since the criterion for suppression is $$x_\mathrm{H\, \small {II}} > 0.1$$, the treatment of partially ionized cells is important. We therefore implemented an improvement over the previous version of grizzly when estimating the ionization fraction in cells containing an ionization front. We now assign the volume averaged ionized ionization fraction to these cells calculated from the higher resolution 1D profiles. 2.5 Differences between C2-ray and grizzly It is important to realize that C2-ray and grizzly treat the calculation of the reionization history very differently. C2-ray follows the time evolution of the ionization fractions by taking the state at the end of the previously calculated time-step as initial condition for the new time-step. This means that any simulation has to be followed from the emergence of the first sources to a chosen end point, for example when the IGM reaches an average ionization fraction of 99 per cent. The advantage is that any time-dependent effects such as recombinations and changes in the source luminosity or position are treated correctly. The grizzly as well as the related bears code do not use the previous state of the IGM as an initial condition for the next time-step but instead construct the state of the IGM from scratch for every output. The time evolution is captured by considering the growth history of dark matter haloes when calculating their luminosity and in the case of suppressed haloes by using the recipe described above. This approach has as an advantage that it is not required to calculate the entire reionization history if this is not needed. The disadvantage is that time-dependent effects such as recombinations, radiative cooling and changes in the source properties have to be implemented in an approximate manner. These 1D radiative transfer reionization codes share this property with the seminumerical methods based upon the excursion set approach such as 21cmfast (Mesinger, Furlanetto & Cen 2011). 2.6 21-cm signal The differential brightness temperature δTb of the 21-cm signal can be expressed as,   \begin{eqnarray} \delta T_{\rm b}(\boldsymbol {x}, z) & = & 27 \ x_{\rm H\, \small {I}} (\boldsymbol {x}, z) [1+\delta _{\rm B}(\boldsymbol {x}, z)] \left(\frac{\Omega _{\rm B}h^2}{0.023}\right) \nonumber \\ &&\times \, \left(\frac{0.15}{\Omega _{\rm m}h^2}\frac{1+z}{10}\right)^{1/2} \left(1-\frac{T_{\gamma }}{T_{\rm S}} \right)\, \rm {mK}, \end{eqnarray} (4)where the quantities $$x_{\rm H\, \small {I}}$$, δB, and Tγ(z) = 2.73 × (1 + z) K denote the neutral hydrogen fraction, baryonic density contrast, and the CMB brightness temperature, respectively, each at position $$\boldsymbol {x}$$ and redshift z. TS represents the spin temperature of hydrogen in the IGM. In this study we will adopt the high spin temperature approximation: TS ≫ Tγ. In this case δTb becomes insensitive to the actual value of the spin temperature. This approximation is valid when the spin temperature is coupled to the IGM gas temperature through Ly α coupling and the IGM has been sufficiently heated (TK ≫ Tγ) by the X-ray sources. Most reionization scenarios reach this state well before substantial ionized regions form (Cohen et al. 2017). After generating the ionization fraction cubes as described in Sections 2.3 and 2.4, it is straightforward to generate the differential brightness temperature cubes using equation (4). Below we will compare $$x_{\rm H\, \small {II}}$$ and δTb results at fixed redshifts. We will not consider the effect of peculiar velocities in the IGM (leading to the so-called redshift-space distortions). The comparison of fully numerical and seminumerical models in Majumdar et al. (2014) suggest that redshift-space distortions will not change the main conclusions of the comparison. 3 RESULTS In this section, we compare the different outputs from the two simulation codes considered in this work, namely C2-ray and grizzly. Before evaluating the different scenarios listed in Table 1, we first consider two single source scenarios to test the impact of the fundamental difference between the two codes as described in Section 2.5. 3.1 Test scenarios: single sources To gauge the basic differences between C2-ray and grizzly we follow the growth of an ionized region for two single source scenarios. For both we take the IGM around the source to have a uniform density at the mean value of the Universe. We choose the size of the simulation box to be 20 h−1 cMpc and use a grid of size 1003. We assume that the source formed at redshift 15 and we follow the impact of the source on the IGM up to redshift 6 using time-steps of 10 Myr. TEST A: The source is hosted by a dark matter halo with constant mass which we take to be 109 M⊙. We choose an efficiency factor gγ = 21.7, as defined in equation (1), which sets the rate of production of ionizing photons to $$\dot{N_\gamma } = 1.33\times 10^{52} \ \rm s^{-1}$$. TEST B: The source is hosted by a dark matter halo with a mass which changes exponentially with redshift as (see e.g. Wechsler et al. 2002),   \begin{equation} M_h(z)=m_i \times \exp \left[\beta (z_i -z)\right], \end{equation} (5)where mi = 109 M⊙, zi = 15, β = 0.5. This scenario mimicks the typical growth of dark matter haloes and thus will show the impact of capturing source evolution through an effective mass Meff in grizzly. We choose the same efficiency gγ as in Test A. The performance of the methods for these test scenarios is shown in Fig. 1. The left-hand and middle panels show the evolution of the neutral fraction in the simulation volume for TEST A and B, respectively, whereas the right-hand panel shows relative difference between the two codes. We see that they produce very similar results for both tests with differences less than 10 per cent. However, for TEST A the ionization front grows faster in the grizzly results. The probable reason is that grizzly underestimates the impact of recombinations. While placing the H ii bubble around the sources, we assume that the source has been emitting ionizing photons at a constant rate determined by its effective mass (Meff) during its entire age and we include the effect of recombinations using the densities at the redshift which we are considering. Since the expansion of the Universe causes densities to be higher at higher redshifts, this will underestimate the impact of recombinations over the age of the source. C2-ray calculates the recombination rates at every time-step using the densities at that redshift and thus captures this evolution of the recombination rate. Figure 1. View largeDownload slide Left-hand panel: Redshift evolution of the mass averaged neutral fraction due to a single source in a simulation volume of size 20 h−1 cMpc. The solid and dashed curves represent C2-ray and grizzly respectively both including recombination processes, while the dotted curve corresponds to grizzly without any recombination processes. The emission rate of the ionizing photons is proportional to the mass of the halo which is chosen to be constant (109 M⊙) in this case. Middle-panel: Same as the left-hand panel, but for a dark matter halo whose mass changes exponentially with redshift. Right-hand panel: Relative difference between the ionization fractions calculated from C2-ray and grizzly as a function of redshift for the test scenarios shown in the left-hand and middle panels of the figure. Figure 1. View largeDownload slide Left-hand panel: Redshift evolution of the mass averaged neutral fraction due to a single source in a simulation volume of size 20 h−1 cMpc. The solid and dashed curves represent C2-ray and grizzly respectively both including recombination processes, while the dotted curve corresponds to grizzly without any recombination processes. The emission rate of the ionizing photons is proportional to the mass of the halo which is chosen to be constant (109 M⊙) in this case. Middle-panel: Same as the left-hand panel, but for a dark matter halo whose mass changes exponentially with redshift. Right-hand panel: Relative difference between the ionization fractions calculated from C2-ray and grizzly as a function of redshift for the test scenarios shown in the left-hand and middle panels of the figure. Interestingly, the match between the two codes is better for TEST B where the halo mass is increasing in time (see the middle panel of Fig. 1). Since in this case most photons are emitted at later times and the growth of the H ii region is faster at later times, the effect of underestimating recombinations is less severe in this case. Below when considering more complex simulations with multiple sources we will see further demonstrations of the impact of recombinations. We conclude that grizzly is capable of approximating the expansion of an ionization front around a growing halo to within 10 per cent accuracy, giving confidence that the effective mass approach works well. When considering the multiple sources case below, larger differences than the ones found here will have to be due to how grizzly deals with overlapping H ii regions and source suppression. 3.2 Model I : Large haloes only Now, we will consider a model where the reionization is driven by haloes with mass Mh ≥ 2.2 × 109 M⊙ in a volume 163 comoving Mpc across. These haloes are never suppressed. This model is same as the L3 model in Iliev et al. (2012) and was also used in Majumdar et al. (2014) to compare seminumerical simulation techniques with a fully numerical one. Not taking into account the time needed to construct the library of 1D ionization profiles, the simulation run time for the whole reionization history is ∼5.5 CPU hours8 for grizzly. This number is without generating the TK maps and Ly α flux maps which are not being used here. For comparison, the entire C2-ray simulation took ∼105 CPU hours. So there is a factor ∼105 gain in computing time when using grizzly. 3.2.1 Reionization history Using the source model described in Section 2.2 and following the methods described in Sections 2.3 and 2.4, we generate the ionization maps in the simulation box. The ionization histories from these two simulations are shown in the bottom panel of Fig. 2. In this scenario, the first sources appear at redshift 19 and ends around z ∼ 8.4. The Thomson scattering optical depth from the grizzly and C2-ray simulations are 0.0696 and 0.0697, respectively. As can be seen from the figure, the evolutions of the mass averaged neutral fraction as a function of redshift from these two simulations are very similar. Figure 2. View largeDownload slide Bottom panel: Redshift evolution of the mass averaged neutral fraction for Model I from C2-ray (solid curve) and grizzly (dashed curve). Top panel: Same but showing the redshift evolution of the ratio of the volume and mass averaged neutral fractions. Figure 2. View largeDownload slide Bottom panel: Redshift evolution of the mass averaged neutral fraction for Model I from C2-ray (solid curve) and grizzly (dashed curve). Top panel: Same but showing the redshift evolution of the ratio of the volume and mass averaged neutral fractions. The top panel of Fig. 2 shows the redshift evolution of the ratio of the volume and mass averaged neutral fraction for this reionization history. As the reionization scenario considered here is ‘inside-out’, high-density regions are ionized first. This leads to a higher volume averaged neutral fraction than the mass averaged neutral fraction. Also for this ratio we find very good agreement between the two codes. This suggest that the ionization maps from these two schemes will be very similar. The small differences may originate from the approximations that go into the implementation of the ionization and differences in the treatment of recombinations in the two schemes. This agreement is better than for any of the seminumerical results from Majumdar et al. (2014), see fig. 1 in that paper. To illustrate the similarity between the reionization histories from these two methods even further, we show the map of redshifts of reionization of each pixel in a chosen slice from the simulation boxes in Fig. 3. We define the redshift of reionization of a region to be the redshift when the ionization fraction first reaches 0.5. The maps in Fig. 3 from grizzly and C2-ray are visually quite similar. The same slice for the seminumerical methods is shown in fig. 2 of Majumdar et al. (2014) and shows larger differences. To quantify the similarity of these two maps, we use the Pearson cross-correlation coefficient defined as,   \begin{equation} \chi _{ab}=\frac{\sum _{i}{(a_i-\bar{a})(b_i-\bar{b})}}{\sqrt{\sum _{i}(a_i-\bar{a})^2}\sqrt{\sum _{i}(b_i-\bar{b})^2}}, \end{equation} (6)where $$\bar{a}$$ and $$\bar{b}$$ represent the mean of the maps a and b respectively and i is the pixel index. The Pearson-cross-correlation coefficient calculated for the maps at the top and bottom panels of Fig. 3 is 0.97. Figure 3. View largeDownload slide Top panel: Image showing the redshifts of reionization for pixels in a slice from the simulation volume for Model I generated using the results from C2-ray. Bottom panel: Same as the top panel but generated from the grizzly results. Figure 3. View largeDownload slide Top panel: Image showing the redshifts of reionization for pixels in a slice from the simulation volume for Model I generated using the results from C2-ray. Bottom panel: Same as the top panel but generated from the grizzly results. To compare two different maps, we define the cross-correlation coefficient of a quantity Q ($$R^{AB}_Q$$) as given below:   \begin{equation} R^{AB}_{Q} (k) = \frac{P_{AB}(k)}{\sqrt{P_{AA}(k)P_{BB}(k)}}, \end{equation} (7)where PAA and PBB are the auto power spectrum of the fields A and B respectively and the PAB is the cross power spectrum between them. Here, we will consider A and B as grizzly and C2-ray, respectively. The quantity Q could be the redshift of reionization, neutral fraction, brightness temperatures, etc. The cross-correlation coefficient for the redshift of reionization cubes Rz is found to be larger than 0.95 for scales with k ≲ 1 Mpc− 1 (see top panel of Fig. 4). This suggests that the reionization histories from these two simulations are very similar at large scales. Although the correlation drops at scales correspond to k ≳ 1 Mpc− 1, it is clearly still better than the seminumerical schemes considered in Majumdar et al. (2014) (see fig. 3 of their paper). Figure 4. View largeDownload slide Top panel: Scale dependence of the cross-correlation coefficient of the redshift of reionization maps from C2-ray and grizzly simulations for Model I. Bottom panel: Same as the top panel but for the neutral fraction maps at different stages of reionization. Figure 4. View largeDownload slide Top panel: Scale dependence of the cross-correlation coefficient of the redshift of reionization maps from C2-ray and grizzly simulations for Model I. Bottom panel: Same as the top panel but for the neutral fraction maps at different stages of reionization. 3.2.2 Morphology of H ii regions Next we compare the morphology of the ionized regions from the two simulations using maps of the neutral fraction. We show the same slice from the neutral fraction cubes from C2-ray and grizzly in the two panels of Fig. 5. Visually, the maps are quite similar. The Pearson cross-correlation coefficient calculated for these two slices is 0.93. One can notice that the ionized bubbles in the slice from grizzly simulation appear to be more spherical than those from C2-ray. This is due to the application of 1D radiative transfer which erases the directional dependence of the H ii bubbles. However, the shapes of the bubbles become irregular once there is a significant amount of overlap between the H ii regions. Figure 5. View largeDownload slide Images of the neutral fraction for the same slice through the simulation box as in Fig. 3. Black is ionized, white is neutral. The top and bottom panels show the results from C2-ray and grizzly, respectively. The results are for Model I, z = 9.026 when the mass averaged neutral fraction is ∼0.5. Figure 5. View largeDownload slide Images of the neutral fraction for the same slice through the simulation box as in Fig. 3. Black is ionized, white is neutral. The top and bottom panels show the results from C2-ray and grizzly, respectively. The results are for Model I, z = 9.026 when the mass averaged neutral fraction is ∼0.5. The scale-dependent cross-correlation coefficients for maps of the neutral fraction at different stages of reionization from grizzly and C2-ray are shown in the bottom panel of Fig. 4. One can see very high correlation at large length scales, which are the ones of prime interest to the first generations of radio telescopes. We find that the cross-correlation coefficient of the neutral maps remains very close to unity (within ∼5 per cent) at scales k ≤ 1 Mpc− 1 until the universe is 50 per cent ionized. There is a consistent shift to smaller k values in the scale at which the correlation coefficient reaches 0.9 as reionization progresses. Although the correlation of the ionization fields drops for small length scales, the values actually remain higher than those for the seminumerical models considered in Majumdar et al. (2014). For example, the cross-correlation coefficients of the neutral fractions between C2-ray and different seminumerical schemes were found to be smaller than ∼0.4 at a scale k ∼ 2 Mpc− 1 at a stage with neutral fraction $$x_{\rm H\, \small {I}} \sim 0.5$$, while the corresponding value for grizzly is ∼0.6. 3.2.3 Size distribution of H ii regions We next compare the H ii BSD from these two simulations. BSDs can provide valuable information about the growth of the ionized regions (Friedrich et al. 2011). However, due to the complex morphologies of the ionized regions, there is no unique way to define the sizes of the irregularly shaped H ii regions in simulation volume. Therefore different size determination methods have been developed and each address different aspects of the size distributions (Giri et al. 2018). In this study, we have use two such methods, namely the Friends-of-friends (FOF; Iliev et al. 2006a) and mean-free-path (MFP; Mesinger & Furlanetto 2007). First, we discuss our bubble size statistics from the FOF method. This algorithm considers two neighbouring cells with ionization fraction larger than 0.5 to be the part of the same ionized bubble and in the process finds the volume of each ionized region. This BSD algorithm focuses on connectivity. As shown in Furlanetto & Oh (2016) the percolative nature of reionization leads to the development of one dominant connected region which contains most of the ionized volume. This large connected region is known as the percolation cluster. The single percolation cluster is invisible in the normalized histogram of the volumes (VdP/dV) as this quantity gives the number of ionized bubbles of volume V. Instead we plot (V/Vion)VdP/dV which shows the fraction of the ionized volume (Vion) contained in regions of volume V. Once the percolation cluster has formed it will have a value close to 1 in this quantity. The top panel of Fig. 6 shows the FOF-BSDs from the two simulations at three different stages of reionization. The volume of the percolation cluster is identical in both the simulations. One can observe that the number and size distribution of smaller ionized regions also match quite well during the early and late stages of reionization. However, we find that grizzly produces a few larger bubbles during the intermediate stages. However, compared to the different seminumerical models considered in Majumdar et al. (2014) (as shown in fig. 5 of that paper), the FOF-BSD from grizzly matches the C2-ray results better. Figure 6. View largeDownload slide Top panel: The volume distribution of the H ii bubbles calculated by the FOF method at three different stages of reionization in Model I. The left to right curves (with colours green, red and blue) correspond to redshift 9.5, 9, and 8.7 when the neutral fractions are 0.7, 0.5, and 0.3, respectively. The solid and dashed curves represent C2-ray and grizzly simulations respectively. Bottom panel: The PDF of the H ii bubble size calculated by the MFP method. The line styles and colours have the same meaning as in the upper panel. Figure 6. View largeDownload slide Top panel: The volume distribution of the H ii bubbles calculated by the FOF method at three different stages of reionization in Model I. The left to right curves (with colours green, red and blue) correspond to redshift 9.5, 9, and 8.7 when the neutral fractions are 0.7, 0.5, and 0.3, respectively. The solid and dashed curves represent C2-ray and grizzly simulations respectively. Bottom panel: The PDF of the H ii bubble size calculated by the MFP method. The line styles and colours have the same meaning as in the upper panel. The second BSD method that we use is MFP which is built on a Monte Carlo inference of the sizes. This method selects random cells in the ionized regions and traces a ray in a random direction until it hits a neutral cell. In this case, we consider a cell to be ionized if its ionization fraction is larger than 0.5 and neutral if it is less than that. We repeat the ray tracing for enough number of times (in this case 107 times) and record the lengths of the rays. This provides the estimate of the size distribution of the ionized regions, here expressed in radii instead of volumes. The MFP-BSDs for grizzly and C2-ray show good agreement as can be seen in the bottom panel of Fig. 6. The peaks of the curves at the different stages of reionization appear at the same radius which suggests that the most probable sizes of ionized regions in both simulations are identical. Majumdar et al. (2014) did not calculate MFP-BSDs and therefore we cannot compare to their results in this case. 3.2.4 Differential brightness temperature maps The next quantity that we compare for the two simulations is the expected 21-cm signal. Fig. 7 presents δTb maps of the same slice at three different redshifts 9.5, 9, and 8.7 which correspond to the neutral fractions 0.7, 0.5, and 0.3, respectively. Visually the δTb maps are quite similar. The Pearson-cross-correlation coefficients for these maps are 0.92, 0.9, and 0.9, respectively. However, the small-scale features are different due to the spherical nature of the bubbles around isolated sources in the 1D radiative transfer scheme. This was also seen in Section 3.2.2 when we considered the morphology of ionized regions. Figure 7. View largeDownload slide The brightness temperature maps of the same slice as in Fig. 5 for Model I. Top and bottom panels show the C2-ray and grizzly results, respectively. The left-hand to right-hand panels correspond to redshift 9.5, 9, and 8.7 with neutral fraction of 0.7, 0.5, and 0.3, respectively. Figure 7. View largeDownload slide The brightness temperature maps of the same slice as in Fig. 5 for Model I. Top and bottom panels show the C2-ray and grizzly results, respectively. The left-hand to right-hand panels correspond to redshift 9.5, 9, and 8.7 with neutral fraction of 0.7, 0.5, and 0.3, respectively. To check the level of similarity, we again use the scale-dependent cross-correlation between the two δTb maps, see Fig. 8. The large values of the cross-correlation coefficients $$R_{\delta T_{\rm b}}$$ suggest that the δTb maps from these two schemes are very similar. However, we do find that the value of $$R_{\delta T_{\rm b}}$$ decreases at small scales which indicates mismatches between the δTb maps at small scales. Still, $$R_{\delta T_{\rm b}}$$ is larger than 0.7 at scales such as k ≥ 1 Mpc− 1. This suggests that δTb maps from grizzly are more similar to the maps from C2-ray than the maps from the seminumerical models considered in Majumdar et al. (2014). For example, $$R_{\delta T_{\rm b}}$$ is less than 0.7 at k ≥ 1 Mpc− 1 for the seminumerical models used in Majumdar et al. (2014) while $$R_{\delta T_{\rm b}} \sim 0.85$$ for grizzly. We also find the $$R_{\delta T_{\rm b}}$$ value to be larger than the $$R_{x_{\rm H\, \small {I}}}$$ at scales k ≥ 1 Mpc− 1 (see Fig. 4). This is due to the fact that the δTb signal is dominated by the density fluctuation at small scales and these are identical between the two simulations. Figure 8. View largeDownload slide Scale dependence of the cross-correlation coefficient of the brightness temperature maps from C2-ray and grizzly for model I. Different curves represent redshifts 11.8, 9.5, 9, and 8.7 with neutral fraction 0.9, 0.7, 0.5, and 0.3, respectively. Figure 8. View largeDownload slide Scale dependence of the cross-correlation coefficient of the brightness temperature maps from C2-ray and grizzly for model I. Different curves represent redshifts 11.8, 9.5, 9, and 8.7 with neutral fraction 0.9, 0.7, 0.5, and 0.3, respectively. 3.2.5 Spherically averaged power spectrum of the signal Next we compare the outputs of these two schemes in terms of the power spectrum of the δTb fluctuations. The spherically averaged power spectrum P(k) of the δTb fluctuations is defined as   \begin{equation} \langle \hat{\delta T_{\rm b}}({\boldsymbol k}) \hat{\delta T_{\rm b}}^{\star }({\boldsymbol k^{\prime }})\rangle = (2 \pi )^3 \delta _D({\boldsymbol k - \boldsymbol k^{\prime }}) P(k), \end{equation} (8)where $$\hat{\delta T_{\rm b}}(\boldsymbol k)$$ represents the Fourier transform of $$\delta T_{\rm b}({\boldsymbol x})$$. Here, We present the results in terms of the spherically averaged dimensionless power spectrum Δ2(k) = k3P(k)/2π2. The top panel of Fig. 9 presents the scale dependence of Δ2(k) estimated from these two simulations at different stages of reionization, together with the ratio of the two results. While the power spectrum is dominated by the density fluctuations at the early stages of reionization, it starts to deviate from the background dark matter density power spectrum as reionization progress. The power spectra from these two schemes agree with each other to within 10 per cent for most stages and scales. Only at the largest scales for $$x_\mathrm{H\,{\small {I}}}=0.7$$ the difference exceeds 20 per cent. The overall features of the power spectra are very similar between the two codes. Figure 9. View largeDownload slide Top panel: Scale dependence of the spherically averaged dimensionless power spectrum of δTb fluctuations at redshifts 11.8, 9.5, 9, and 8.7 with neutral fraction 0.9, 0.7, 0.5, and 0.3, respectively. The power spectrum corresponds to the Model I of the paper which considers ionizing photons from the massive haloes only. Bottom panel: Redshift evolution of the large-scale (k = 0.1 Mpc− 1) spherically averaged power spectrum estimated from C2-ray and grizzly which use 3D and 1D radiative transfer, respectively. The ratios of the power spectrum from C2-ray and grizzly are also shown in the upper parts of the panels. $${\Delta }^2_{\rm C}$$ and $${\Delta }^2_{\rm G}$$ represent the dimensionless power spectrum from C2-ray and grizzly, respectively. Figure 9. View largeDownload slide Top panel: Scale dependence of the spherically averaged dimensionless power spectrum of δTb fluctuations at redshifts 11.8, 9.5, 9, and 8.7 with neutral fraction 0.9, 0.7, 0.5, and 0.3, respectively. The power spectrum corresponds to the Model I of the paper which considers ionizing photons from the massive haloes only. Bottom panel: Redshift evolution of the large-scale (k = 0.1 Mpc− 1) spherically averaged power spectrum estimated from C2-ray and grizzly which use 3D and 1D radiative transfer, respectively. The ratios of the power spectrum from C2-ray and grizzly are also shown in the upper parts of the panels. $${\Delta }^2_{\rm C}$$ and $${\Delta }^2_{\rm G}$$ represent the dimensionless power spectrum from C2-ray and grizzly, respectively. The bottom panel of the figure compares the redshift evolution for the k ∼ 0.1 Mpc− 1 mode of the power spectrum. The ratio between the C2-ray and grizzly results starts to deviate from unity at redshift 14 and remains between 0.9 and 1.1 during the later stages of reionization once again indicating excellent agreement. Compared to the seminumerical methods from Majumdar et al. (2014) the grizzly results match the C2-ray ones equally well, or better. 3.3 Model II : Including low-mass haloes and thermal feedback Finally we consider a more complex model for reionization and study the similarities of the outcomes of these two codes. While our fiducial model does not incorporate the contributions from LMACHs, these low-mass haloes can provide a significant contribution to the overall ionizing photon budget. As explained in Section 2.2, in Model II we take the minimum mass of haloes which contribute to reionization to be 108 M⊙ and we suppress star formation in haloes of masses (below 109 M⊙) in ionized regions. This scenario is same as the L1 model in Iliev et al. (2012). The details can be found in Table 1. Majumdar et al. (2014) did not study this model in their comparison between seminumerical and numerical methods, so unlike for the results of Model I we are here not able to compare to seminumerical methods. Although grizzly has to apply simplifying approximations in order to include thermal feedback, the reionization history it generates is very similar to that from C2-ray as shown in the bottom panel of Fig. 10. The Thomson scattering optical depths for C2-ray and grizzly are 0.08 and 0.081, respectively. The redshift evolution of the ratio of the volume and mass averaged neutral fraction from these two simulations also matches very well as shown in the top panel of the figure. This suggests that the neutral fraction maps from these two simulations for this complex model are also very similar. We do not show these here but below will show the 21-cm images from which the morphology of the ionized regions can also be seen. Figure 10. View largeDownload slide Bottom panel: Redshift evolution of the mass averaged neutral fraction from C2-ray and grizzly for the reionization scenario Model II which includes the contributions from low mass haloes and considers thermal feedback. Top panel: Redshift evolution of the ratio of the volume and mass averaged neutral fraction from both the simulations. Figure 10. View largeDownload slide Bottom panel: Redshift evolution of the mass averaged neutral fraction from C2-ray and grizzly for the reionization scenario Model II which includes the contributions from low mass haloes and considers thermal feedback. Top panel: Redshift evolution of the ratio of the volume and mass averaged neutral fraction from both the simulations. The FOF-BSDs obtained for Model II as shown in the top panel of Fig. 11 have mostly very similar shapes. At the early stages (with neutral fraction 0.7), grizzly is seen to produce more large ionized regions, although this difference is likely due to small size differences among the larger regions. However, this possible bias vanishes once the H ii regions start percolating during the middle and later stages of reionization. The distributions of the sizes of the H ii regions calculated using the MFP method for this scenario also display good agreement (as shown in the bottom panel of Fig. 11). We find that there is no prominent peak feature in the MFP-BSD curves during the early stages of reionization. This is due to the fact that the size of the ionized regions are very small and mostly fall within the resolution limit of the simulation. Figure 11. View largeDownload slide Top panel: The volume distribution of the H ii bubbles as a function of volume at different stages of reionization for the two different simulations considered in this work. The left to right curves (with colours green, red, and blue) correspond to ionization fraction 0.7, 0.5, and 0.3 at redshift 10.5, 9.5, and 8.9, respectively. The solid and dashed curves correspond to the C2-ray and grizzly simulation, respectively. Bottom panel: The PDF of the radius of the H ii bubbles as a function of the radius of the bubbles. The BSD in this panel is estimated using the MFP method. The reionization scenario considered here is Model II. Figure 11. View largeDownload slide Top panel: The volume distribution of the H ii bubbles as a function of volume at different stages of reionization for the two different simulations considered in this work. The left to right curves (with colours green, red, and blue) correspond to ionization fraction 0.7, 0.5, and 0.3 at redshift 10.5, 9.5, and 8.9, respectively. The solid and dashed curves correspond to the C2-ray and grizzly simulation, respectively. Bottom panel: The PDF of the radius of the H ii bubbles as a function of the radius of the bubbles. The BSD in this panel is estimated using the MFP method. The reionization scenario considered here is Model II. The differential brightness temperature maps from these two simulations also visually very similar at different stages of reionization as shown in Fig. 12. The Pearson-cross-correlation coefficient calculated for the slices shown in Fig. 12 are 0.7, 0.75, and 0.77 at reionization stages with neutral fraction 0.7, 0.5, and 0.3, respectively. The scale-dependent comparison of the δTb cubes from these two simulations is also shown in terms of the cross-correlation coefficient at different stages of reionization in Fig. 13. One can see that the δTb cubes from these two simulations are highly correlated at the large scales with k ≲ 0.5 Mpc− 1 (within 10 per cent error), while the correlation drops at smaller scales. Figure 12. View largeDownload slide Top and bottom panels represent the δTb maps corresponding to C2-ray and grizzly simulations respectively which are based on 3D and 1D radiative transfer schemes. The left-hand to right-hand panels represent redshift 10.5, 9.5, and 8.9, respectively, which correspond to neutral fraction 0.7, 0.5, and 0.3, respectively. The maps correspond to Model II which include the contributions from low-mass haloes and consider thermal feedback. Figure 12. View largeDownload slide Top and bottom panels represent the δTb maps corresponding to C2-ray and grizzly simulations respectively which are based on 3D and 1D radiative transfer schemes. The left-hand to right-hand panels represent redshift 10.5, 9.5, and 8.9, respectively, which correspond to neutral fraction 0.7, 0.5, and 0.3, respectively. The maps correspond to Model II which include the contributions from low-mass haloes and consider thermal feedback. Figure 13. View largeDownload slide Scale dependence of the cross-correlation coefficient of the brightness temperature maps from C2-ray and grizzly at different stages of reionization for Model II. Figure 13. View largeDownload slide Scale dependence of the cross-correlation coefficient of the brightness temperature maps from C2-ray and grizzly at different stages of reionization for Model II. We present the comparison between the spherically averaged power spectrum from these two simulations of Model II in Fig. 14. The redshift evolution of the spherically averaged power spectrum at scale k ∼ 0.1 Mpc− 1 are quite similar to each other as shown in the bottom panel of Fig. 14. Although the power spectrum estimated from these two schemes are very similar, there are small differences too. For example, we find that the power spectra deviate from each other around redshift 10.5 when the neutral fraction is 0.7. Up to this time the LMACHs contribute substantially to reionization, so the differences are most likely connected to the way their suppression is implemented. Still, the differences are small and except at the very largest scales do not exceed 10 per cent. Figure 14. View largeDownload slide Top panel: The spherically averaged power spectrum of the 21-cm signal for Model II as a function of scale at different stages of reionization. Different curves correspond to redshift 13.5, 10.5, 9.5, and 8.9 with neutral fraction 0.9, 0.7, 0.5, and 0.3, respectively. Bottom panel: δTb power spectrum for k = 0.1 Mpc− 1 as a function of redshift. The solid and dashed curves in both the panels represent C2-ray and grizzly simulations, respectively. Upper parts of the panels show the ratio of the power spectrum from C2-ray and grizzly simulations. Here, $${\Delta }^2_{\rm C}$$ and $${\Delta }^2_{\rm G}$$ represent the dimensionless power spectrum from C2-ray and grizzly, respectively. Figure 14. View largeDownload slide Top panel: The spherically averaged power spectrum of the 21-cm signal for Model II as a function of scale at different stages of reionization. Different curves correspond to redshift 13.5, 10.5, 9.5, and 8.9 with neutral fraction 0.9, 0.7, 0.5, and 0.3, respectively. Bottom panel: δTb power spectrum for k = 0.1 Mpc− 1 as a function of redshift. The solid and dashed curves in both the panels represent C2-ray and grizzly simulations, respectively. Upper parts of the panels show the ratio of the power spectrum from C2-ray and grizzly simulations. Here, $${\Delta }^2_{\rm C}$$ and $${\Delta }^2_{\rm G}$$ represent the dimensionless power spectrum from C2-ray and grizzly, respectively. When comparing these results for Model II to the ones for Model I we see that grizzly performs somewhat less well in the former case. However, the results are still very close to the numerical ones and definitely fall within the expected measurement errors from the observations. Overall, we find that the results from these two simulations agree well with each other even for this more complex scenario. This suggest that the 1D radiative transfer can be used for the complex scenarios as well without a substantial loss of accuracy. Specifically, the 3D and 1D radiative transfer schemes match quite well at scales k ≲ 1 Mpc−1 which are the prime target of the current radio interferometers. This suggests that 1D radiative transfer codes such as bears or grizzly can be efficiently used for predicting the expected signal for various source models as well as to find new observation strategies and possibly for parameter estimation. 4 SUMMARY AND DISCUSSION In this paper, we have compared the performance of the 1D radiative transfer code grizzly to the 3D radiative transfer code C2-ray in the context of a large-scale reionization simulation. While C2-ray simulates the process more accurately, by performing both the full 3D radiative transfer and following the evolution of the ionization fractions and the source in time, this accuracy comes at a considerable computational cost. The grizzly code simplifies the calculation by placing pre-generated 1D profiles of ionization fraction around the sources and using the halo growth history to approximate the evolutionary effects, which makes it approximately a factor ∼105 faster. We have used same sets of initial conditions, i.e. the same density fields, halo catalogues and source models for both simulations. We limit our comparison to the calculation of the ionization fraction. Beside comparing results around a single source, we consider two models to compare the outcomes of these simulations in detail. In Model I (also our fiducial model), we assume that the reionization is driven by massive haloes with masses larger than 2.2 × 109 M⊙. This model was also used in Majumdar et al. (2014) to evaluate the performance of seminumerical reionization codes which rely on the excursion set formalism. We also consider a more complex model (Model II) which includes low mass haloes (down to 108 M⊙) and suppression of star formation in haloes with mass smaller than 109 M⊙ due to thermal feedback. For all the models considered in this paper, we have used the same SED which in this case is a blackbody spectrum. When calculating the differential brightness temperature of the 21-cm signal we assume that TS ≫ Tγ. Our findings from this comparison are the following. When comparing the reionization history due to an isolated source, we observe that grizzly underestimates the effect of recombinations. Since this is a cumulative effect, the impact is largest at the later stages. Still, we find that the differences between the neutral fractions estimated from the two methods for the isolated sources remain less than 10 per cent. For our fiducial simulation, the reionization histories in terms of the evolution of the mass averaged neutral fraction and ionization maps are very similar between both simulations. The cross-correlation coefficient of the neutral fraction maps remains close to unity (within 5–10 per cent) at scales k ≲ 1 Mpc− 1 before the universe reaches 50 per cent ionization and somewhat worse after this. We have considered two different methods to characterize the size distribution of ionized regions, namely the FOF and MFP methods. For both methods the size distributions match very well between the C2-ray and grizzly results. However, visual inspection of the ionization fraction maps does show that grizzly produces more spherically shaped regions compared to C2-ray. Visually, the differential brightness temperature maps from the two simulations look very similar. The cross-correlation coefficients at large scales, k ≲ 0.1 Mpc− 1 are close to unity. We find that the correlation coefficients at small scales are higher for the δTb maps than for the ionization maps. This is due to the density fluctuations in the neutral medium. The spherically averaged power spectra from the two simulations agree with each other within 10 per cent. For the more complex scenario Model II, we find that the differences between the grizzly and C2-ray are generally somewhat larger than for the simpler Model I but remain small. We conclude that grizzly is capable of handling more complex source models where sources are not always ‘on’. The results from grizzly are more similar to those from C2-ray compared to the different seminumerical simulations considered in Majumdar et al. (2014). We find for example that the various cross-correlation coefficients are higher for grizzly, especially for smaller scales. Furthermore, for the seminumerical models the source efficiency was adapted at each redshift in order to reproduce the ionization history from C2-ray, whereas grizzly is able to produce an excellent match for the reionization history without varying the source efficiency parameter. These results imply that the faster 1D radiative transfer codes such as bears or grizzly can be efficiently used for exploring parameter space, determining observational strategies or developing parameter estimation pipelines. Once the 1D profiles have been generated, these simulations are fast enough to be used for parameter estimation using techniques such as the MCMCs, machine learning, etc. However, the details of parameter estimation such as convergence time, accuracy, etc. are beyond the scope of this work and will be addressed in future works. Although grizzly is quite fast compared to the full radiative transfer simulations and generally reproduces the results well, our comparison also shows where improvements could be made. (i) grizzly shows lower accuracy at small scales. We speculate that the accuracy at small scales can be improved by optimizing the technique grizzly uses to deal with overlapping ionized regions. (ii) The comparison for a single source shows that grizzly underestimates the effect of recombinations. This could be improved by taking into account the cosmological evolution of the density field over the lifetime of the source. (iii) The halo suppression model used in grizzly is currently a very simple one. It may be possible to implement the evolution of relic H ii regions better and thereby to improve the accuracy of grizzly for scenarios where source can turn ‘off’. Besides the neutral fraction maps, grizzly can also provide the kinetic temperature and Ly α flux maps. Thus, it can be used to predict signal in presence of spin-temperature fluctuations. The 21-cm signal from the Cosmic Dawn is expected to be dominated by these fluctuations. Thus, grizzly is also very relevant for observations with interferometers which also probe this earlier era, such as SKA1-low. grizzly can for example be used to constrain the heating and Ly α coupling parameters during the Cosmic Dawn from the observations. However, a detailed comparison of the results in presence of TS fluctuations with fully numerical simulations such as Ross et al. (2017) is beyond the scope of this paper and will be addressed in the future. Acknowledgements The authors would like to thank Saleem Zaroubi and Adi Nusser for providing useful suggestions on this work. KKD would like to thank the University Grant Commission (UGC), India for support through UGC-faculty recharge scheme (UGC-FRP) vide ref.no. F.4-5(137-FRP)/2014(BSR). SM acknowledges the financial support from the European Research Council under ERC grant number 638743-FIRSTDAWN. Footnotes 1 http://www.gmrt.tifr.res.in 2 http://eor.berkeley.edu/ 3 http://www.mwatelescope.org/ 4 http://www.lofar.org/ 5 http://www.skatelescope.org/ 6 http://wiki.cita.utoronto.ca/mediawiki/index.php/CubePM 7 The suppression of the low-mass haloes is determined by the kinetic temperature of the region. For example, one can assume no star formation within haloes with masses <109 M⊙ if it is formed within a region with TK larger than 104 K. Following Iliev et al. (2012) we use $$x_{\rm H\, \small {II}}$$ as proxy for the temperature and assume that if a region is more ionized than 0.1, it will be hot enough to suppress further star formation. 8 As mentioned above, grizzly does not track the ionization history to create the ionization map at a certain redshift for simple models where thermal feedback is ignored. Thus, the run time will be much shorter if only a single snapshot of the ionization fields is needed. For example, grizzly generates the neutral fraction map at redshift 9.026 with mass averaged neutral fraction ∼0.5 in ∼10 CPU minutes. In presence of thermal feedback, the estimated time is higher which is also true for other codes such as the 21 cm FAST (see e.g. Sobacchi & Mesinger 2014). The main computational bottleneck for grizzly is to estimate the unused ionizing photons in the overlapping regions and to redistribute these around these overlapping regions. This step is executed iteratively and there are no straightforward methods to compute this more efficiently. REFERENCES Alvarez M. A., Wise J. H., Abel T., 2009a, ApJ , 701, L133 CrossRef Search ADS   Alvarez M. A., Busha M., Abel T., Wechsler R. H., 2009b, ApJ , 703, L167 CrossRef Search ADS   Baek S., Di Matteo P., Semelin B., Combes F., Revaz Y., 2009, A&A , 495, 389 CrossRef Search ADS   Becker R. H. et al.  , 2001, AJ , 122, 2850 CrossRef Search ADS   Bouwens R. J. et al.  , 2015, ApJ , 803, 34 CrossRef Search ADS   Bowman J. D. et al.  , 2013, PASA , 30, e031 CrossRef Search ADS   Bromm V., Loeb A., 2003, ApJ , 596, 34 CrossRef Search ADS   Choudhury T. R., Ferrara A., 2006, MNRAS , 371, L55 CrossRef Search ADS   Choudhury T. R., Haehnelt M. G., Regan J., 2009, MNRAS , 394, 960 CrossRef Search ADS   Choudhury T. R., Puchwein E., Haehnelt M. G., Bolton J. S., 2015, MNRAS , 452, 261 CrossRef Search ADS   Ciardi B., Ferrara A., Marri S., Raimondo G., 2001, MNRAS , 324, 381 CrossRef Search ADS   Cohen A., Fialkov A., Barkana R., Lotem M., 2017, MNRAS , 472, 1915 CrossRef Search ADS   Dijkstra M., Mesinger A., Wyithe J. S. B., 2011, MNRAS , 414, 2139 CrossRef Search ADS   Fan X. et al.  , 2006, AJ , 131, 1203 CrossRef Search ADS   Fialkov A., Barkana R., Visbal E., 2014, Nature , 506, 197 CrossRef Search ADS PubMed  Friedrich M. M., Mellema G., Alvarez M. A., Shapiro P. R., Iliev I. T., 2011, MNRAS , 413, 1353 CrossRef Search ADS   Furlanetto S. R., Oh S. P., 2016, MNRAS , 457, 1813 CrossRef Search ADS   Furlanetto S. R., Zaldarriaga M., Hernquist L., 2004, ApJ , 613, 16 CrossRef Search ADS   Furlanetto S. R., Oh S. P., Briggs F. H., 2006, Phys. Rep. , 433, 181 CrossRef Search ADS   Ghara R., Choudhury T. R., Datta K. K., 2015, MNRAS , 447, 1806 CrossRef Search ADS   Ghara R., Choudhury T. R., Datta K. K., Choudhuri S., 2017, MNRAS , 464, 2234 CrossRef Search ADS   Ghosh A., Prasad J., Bharadwaj S., Ali S. S., Chengalur J. N., 2012, MNRAS , 426, 3295 CrossRef Search ADS   Giri S. K., Mellema G., Dixon K. L., Iliev I. T., 2018, MNRAS , 473, 2949 CrossRef Search ADS   Goto T., Utsumi Y., Hattori T., Miyazaki S., Yamauchi C., 2011, MNRAS , 415, L1 CrossRef Search ADS   Greig B., Mesinger A., 2015, MNRAS , 449, 4246 CrossRef Search ADS   Harnois-Déraps J., Pen U.-L., Iliev I. T., Merz H., Emberson J. D., Desjacques V., 2013, MNRAS , 436, 540 CrossRef Search ADS   Hinshaw G. et al.  , 2013, ApJS , 208, 19 CrossRef Search ADS   Iliev I. T., Mellema G., Pen U.-L., Merz H., Shapiro P. R., Alvarez M. A., 2006a, MNRAS , 369, 1625 CrossRef Search ADS   Iliev I. T. et al.  , 2006b, MNRAS , 371, 1057 CrossRef Search ADS   Iliev I. T. et al.  , 2009, MNRAS , 400, 1283 CrossRef Search ADS   Iliev I. T., Mellema G., Shapiro P. R., Pen U.-L., Mao Y., Koda J., Ahn K., 2012, MNRAS , 423, 2222 CrossRef Search ADS   Kaaret P., 2014, MNRAS , 440, L26 CrossRef Search ADS   Kashikawa N. et al.  , 2011, ApJ , 734, 119 CrossRef Search ADS   Knevitt G., Wynn G. A., Power C., Bolton J. S., 2014, MNRAS , 445, 2034 CrossRef Search ADS   Komatsu E. et al.  , 2011, ApJS , 192, 18 CrossRef Search ADS   Majumdar S., Mellema G., Datta K. K., Jensen H., Choudhury T. R., Bharadwaj S., Friedrich M. M., 2014, MNRAS , 443, 2843 CrossRef Search ADS   Malhotra S., Rhoads J. E., 2006, ApJ , 647, L95 CrossRef Search ADS   McGreer I. D., Mesinger A., D'Odorico V., 2015, MNRAS , 447, 499 CrossRef Search ADS   McQuinn M., Lidz A., Zahn O., Dutta S., Hernquist L., Zaldarriaga M., 2007a, MNRAS , 377, 1043 CrossRef Search ADS   McQuinn M., Hernquist L., Zaldarriaga M., Dutta S., 2007b, MNRAS , 381, 75 CrossRef Search ADS   Mellema G., Iliev I. T., Pen U.-L., Shapiro P. R., 2006, MNRAS , 372, 679 CrossRef Search ADS   Mellema G., Koopmans L., Shukla H., Datta K. K., Mesinger A., Majumdar S., 2015, HI tomographic imaging of the Cosmic Dawn and Epoch of Reionization with SKA . SISSA, Trieste, PoS(AASKA14)010 Mesinger A., Furlanetto S., 2007, ApJ , 669, 663 CrossRef Search ADS   Mesinger A., Furlanetto S., Cen R., 2011, MNRAS , 411, 955 CrossRef Search ADS   Mitra S., Choudhury T. R., Ferrara A., 2011, MNRAS , 413, 1569 CrossRef Search ADS   Mitra S., Choudhury T. R., Ferrara A., 2012, MNRAS , 419, 1480 CrossRef Search ADS   Mitra S., Choudhury T. R., Ferrara A., 2015, MNRAS , 454, L76 CrossRef Search ADS   Morales M. F., Wyithe J. S. B., 2010, ARA&A , 48, 127 CrossRef Search ADS   Mortlock D. J. et al.  , 2011, Nature , 474, 616 CrossRef Search ADS PubMed  Paciga G. et al.  , 2013, MNRAS , 433, 639 CrossRef Search ADS   Paranjape A., Choudhury T. R., 2014, MNRAS , 442, 1470 CrossRef Search ADS   Parsons A. R. et al.  , 2014, ApJ , 788, 106 CrossRef Search ADS   Patil A. H. et al.  , 2014, MNRAS , 443, 1113 CrossRef Search ADS   Patil A. H. et al.  , 2017, ApJ , 838, 65 CrossRef Search ADS   Planck Collaboration XVI, 2014, A&A , 571, A16 CrossRef Search ADS   Planck Collaboration XIII, 2016a, A&A , 594, A13 CrossRef Search ADS   Planck Collaboration XV, 2016b, A&A , 594, A15 CrossRef Search ADS   Pritchard J. R., Loeb A., 2012, Rep. Prog. Phys. , 75, 086901 CrossRef Search ADS PubMed  Ross H. E., Dixon K. L., Iliev I. T., Mellema G., 2017, MNRAS , 468, 3785 CrossRef Search ADS   Santos M. G., Amblard A., Pritchard J., Trac H., Cen R., Cooray A., 2008, ApJ , 689, 1 CrossRef Search ADS   Santos M. G., Ferramacho L., Silva M. B., Amblard A., Cooray A., 2010, MNRAS , 406, 2421 CrossRef Search ADS   Shimabukuro H., Semelin B., 2017, MNRAS , 468, 3869 CrossRef Search ADS   Shin M.-S., Trac H., Cen R., 2008, ApJ , 681, 756 CrossRef Search ADS   Sobacchi E., Mesinger A., 2014, MNRAS , 440, 1662 CrossRef Search ADS   Tanaka T., Perna R., Haiman Z., 2012, MNRAS , 425, 2974 CrossRef Search ADS   Thomas R. M. et al.  , 2009, MNRAS , 393, 32 CrossRef Search ADS   Tingay S. J. et al.  , 2013, PASA , 30, 7 CrossRef Search ADS   van Haarlem M. P. et al.  , 2013, A&A , 556, A2 CrossRef Search ADS   Wechsler R. H., Bullock J. S., Primack J. R., Kravtsov A. V., Dekel A., 2002, ApJ , 568, 52 CrossRef Search ADS   Wise J. H., Demchenko V. G., Halicek M. T., Norman M. L., Turk M. J., Abel T., Smith B. D., 2014, MNRAS , 442, 2560 CrossRef Search ADS   Xu H., Wise J. H., Norman M. L., Ahn K., O'Shea B. W., 2016, ApJ , 833, 84 CrossRef Search ADS   Zahn O., Lidz A., McQuinn M., Dutta S., Hernquist L., Zaldarriaga M., Furlanetto S. R., 2007, ApJ , 654, 12 CrossRef Search ADS   Zahn O., Mesinger A., McQuinn M., Trac H., Cen R., Hernquist L. E., 2011, MNRAS , 414, 727 CrossRef Search ADS   Zaldarriaga M., Seljak U., 1997, Phys. Rev. D , 55, 1830 CrossRef Search ADS   Zaroubi S., Thomas R. M., Sugiyama N., Silk J., 2007, MNRAS , 375, 1269 CrossRef Search ADS   © 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

There are no references for this article.

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
discover and read the research
that matters to you.

Enjoy affordable access to
over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Unlimited reading

Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere.

Stay up to date

Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates.

Organize your research

It’s easy to organize your research with our built-in tools.

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.

See the journals in your area

Monthly Plan

  • Read unlimited articles
  • Personalized recommendations
  • No expiration
  • Print 20 pages per month
  • 20% off on PDF purchases
  • Organize your research
  • Get updates on your journals and topic searches

$49/month

Start Free Trial

14-day Free Trial

Best Deal — 39% off

Annual Plan

  • All the features of the Professional Plan, but for 39% off!
  • Billed annually
  • No expiration
  • For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588

$360/year

billed annually
Start Free Trial

14-day Free Trial