Asking Fast Radio Bursts (FRBs) for More than Reionization History
Abinash Kumar Shaw Raghunath Ghara Paz Beniamini Saleem Zaroubi Pawan Kumar






Abinash Kumar Shaw
abinashkumarshaw@gmail.com

0000-0002-6123-4383]Abinash Kumar Shaw
Astrophysics Research Center of the Open University (ARCO), The Open University of Israel, 1 University Road, Ra'anana 4353701, Israel
Department of Computer Science, University of Nevada Las Vegas, 4505 S. Maryland Pkwy., Las Vegas, NV 89154, USA

0000-0001-9816-5070]Raghunath Ghara
Astrophysics Research Center of the Open University (ARCO), The Open University of Israel, 1 University Road, Ra'anana 4353701, Israel
Department of Physical Sciences, Indian Institute of Science Education and Research Kolkata, Mohanpur, WB 741 246, India

0000-0001-7833-1043]Paz Beniamini
Astrophysics Research Center of the Open University (ARCO), The Open University of Israel, 1 University Road, Ra'anana 4353701, Israel
Department of Natural Sciences, The Open University of Israel, 1 University Road, Ra'anana 4353701, Israel
Department of Physics, The George Washington University, 725 21st Street NW, Washington, DC 20052, USA


0000-0001-9121-8467]Saleem Zaroubi
Astrophysics Research Center of the Open University (ARCO), The Open University of Israel, 1 University Road, Ra'anana 4353701, Israel
Department of Natural Sciences, The Open University of Israel, 1 University Road, Ra'anana 4353701, Israel
Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700AV Gröningen, The Netherlands


Department of Astronomy, University of Texas at Austin, Austin, TX 78712, USA



§ ABSTRACT 
We propose different estimators to probe the epoch of reionization (EoR) intergalactic medium (IGM) using the dispersion measure (DM) of the FRBs. We consider three different reionization histories which we can distinguish with a total of ≲ 1000 DM measurements during EoR if their redshifts are known. We note that the redshift derivatives of DM are also directly sensitive to the reionization history. The major point of this work is exploring the variance in the DM measurements and the information encoded in them. We find that the all-sky average DM(z) gets biased from the LoS fluctuations in the DM measurements introduced by the ionization of IGM during EoR. We find that the ratio σ_ DM/ DM depends directly on the ionization bubble sizes as well as the reionization history. On the other hand, we also find that angular variance (coined as structure function) of DM encodes the information about the duration of reionization and the typical bubble sizes as well. We establish the usefulness of variances in DM using toy models of reionization and later verify it with the realistic reionization simulations.










§ INTRODUCTION

 
According to the current understanding of cosmology, our universe transitioned from being a highly cold-neutral phase in the past to almost hot-ionized phase at present. This is supposed to be a result of UV radiations from the very first objects that formed in the universe photoionizing the intergalactic medium (IGM) <cit.>. This window of transition is termed as the Epoch of Reionization (EoR). The study of the EoR is crucial for answering several questions regarding the emergence of the first sources, their properties and impact on the IGM, evolution to present day structures, the exact timeline of this epoch, etc. Despite our efforts in the last few decades, our understanding of the EoR remains limited <cit.>.

Our present understanding of the timing and duration of EoR is guided by a few indirect observations such as the measurements of the Thomson scattering optical depth from the cosmic microwave background radiation observations <cit.> and the Gunn-Peterson troughs in the high-z quasar spectra <cit.>. Additional constraints on the timeline of EoR come from the recent observations of the high-z Ly-α emitters <cit.> and their clustering measurements <cit.>, Lyman break galaxies <cit.>, and the Ly-α damping wings in the high-z quasar spectra <cit.>. These experiments attempt to constrain the reionization history by putting bounds on the global ionization fraction of the IGM during EoR. On the other hand, the measurements of the effective optical depth of Ly-α forests (using dark-gap/pixel statistics) <cit.> suggests that the end of reionization has a longer tail extending to somewhere between z = 5.5 and 5.0 instead of z ≈ 6. However, all these analyses are either model-dependent or suffer from statistical variance, thus providing only loose bounds on the EoR timeline.

Probing EoR directly using the redshifted 21-cm signal with the current instruments is also challenging because of several hindrances such as large (∼ 10^4 times) foregrounds <cit.>, thermal noise, radio frequency interference, ionospheric turbulence and other systematics. While no undisputed detection of the EoR 21-cm signal has been achieved so far, the current data from the radio-interferometric experiments have been able to provide a few upper limits on the EoR 21-cm power spectra (e.g., LOFAR: , MWA: , HERA: ), and the upper-limits are improving gradually.

A few earlier works have demonstrated the potential of highly energetic astrophysical events such as gamma-ray bursts (GRBs) <cit.> and fast radio bursts (FRBs) <cit.> during EoR as a probe to measure the reionization history. In this work we focus on the FRBs, which are luminous short-duration (∼ few ms) astrophysical radio pulses which have been detected within a frequency band of 0.1-7  GHz <cit.>. Empirical studies on the all-sky event rates of FRBs based on observations find it to be ∼ 10^3/day above the fluence limit of 5 Jy ms and at a central frequency of 500 MHz <cit.>, and the rate is expected to increase significantly at lower fluence thresholds. 

FRB signals disperse while travelling through the ionized medium. The amount of dispersion, quantified as Dispersion Measure (DM), directly depends on the free electron content along its path. The DM of a cosmological FRB is expected to have dominant contribution coming from the electrons in the IGM. During post-EoR (z ≲ 5.5), where the IGM is almost ionized, the IGM DM roughly scales directly with the distances. Therefore the DM measurements can be turned to infer the redshift-distance of the FRB <cit.>. Conversely, knowing the redshift of the FRBs accurately can be potentially used to estimate baryonic content of IGM during the post-EoR <cit.>, probe the epoch of second helium reionization <cit.> and constrain several cosmological parameters <cit.>. In this work we explore how useful they can be as detailed probes of the EoR.


Despite their enigmatic origin, a recent discovery of a galactic FRB <cit.> clearly associates at least some FRBs to magnetars. Other, less direct, evidence linking FRBs to magnetars comes from the statistical properties of the bursts, from host galaxies and offsets relative to them and from the energetics and temporal properties of the bursts <cit.>. Hence, we can expect a sufficiently large number of FRBs during EoR (z>6) which spans a much larger time compared to the life-time of the massive Pop III stars (∼  5 - 30  Myr) that leave behind neutron star (NS) remnants with large angular momentum and strong magnetic fields. There are indirect evidences which supports a large abundance of high-z FRBs <cit.>.

With the possibility of detecting high-z FRBs <cit.>, one can turn their precise DM measurements to probe the sources and IGM during reionization. Recently, a few theoretical studies <cit.> have demonstrated the feasibility of using the DM measurements to extract the reionization history and Thomson scattering optical depth (τ_ Th). <cit.> and <cit.> have used the mean DM from their synthetic FRB population to constrain the parameters of their reionization simulation. The results in most of these works are based on the assumption of knowing the precise redshift (spectroscopic or empirically) of FRBs. However, detecting the precise spectroscopic redshifts from the host galaxies is a challenging task with the current instruments. Conversely, <cit.> suggest that the maximum value of DM for bursts spanning the EoR can provide an independent estimate of the Thomson optical depth of the universe without requiring direct redshift information. They have also shown that ∼ 40 FRBs are sufficient to estimate average electron fraction in 4 z-bins (within z = 6 to 10) with 4% accuracy, if their redshifts could be determined with 10% uncertainty. Similar results have been reported in <cit.> where they also estimated τ_ Th and the mid-redshift of reionization using the average DM. <cit.> also suggested that the reionization history can be constrained from the determination of the number of FRBs during the EoR per unit DM, i.e., dN_ FRB/d DM. Whereas, <cit.> have shown that the redshift derivatives of DM have potential to directly constrain the reionization history.

In this work, we forward model the FRB signal measurements using estimators like globally-averaged DM, its redshift derivative, global dispersion and angular dispersion (along different LoS). For the first time, we demonstrate that fluctuations in the LoS in DM significantly impact measurements during the EoR, potentially biasing mean DM estimates compared to those derived using the mean ionization fraction of the IGM. We also demonstrate that d DM/dz and the scatter in DM along different LoS (as a function of z) has the potential to discern different reionization histories and morphologies of the IGM. We explore a novel aspect of angular dispersion in DM (defined as structure function in <ref>) at different redshifts. This structure function encodes information regarding the typical bubble sizes in the IGM as it probes the angular fluctuations. We validate our estimators on simple toy models of EoR light-cone signals. Later, we also apply our estimators to more realistic light-cones obtained from simulations. We also perform a comparative study between different reionization histories. For this study, we primarily assume a scenario where the redshifts of the FRBs are known to within 10% uncertainty. Later, we ignore the redshift information and compute the marginalized average DM over the EoR window.

This manuscript is arranged as follows. We define DM and structure function estimators in <ref>. In <ref>, we validate our estimators with the toy ionization field model. This is followed by <ref> where we briefly describe the details of the actual EoR simulations, and the present the corresponding results under its subsection where we discuss the impact of post-EoR IGM on the DM estimates and presented results. Finally, we summarize and conclude this exercise in <ref>. The cosmological simulation here uses the cosmological parameter values Ω_ m =0.27, Ω_Λ = 0.73, Ω_ b = 0.044 and h=0.7 adapted from <cit.>.
 




§ METHODOLOGY

We revisit the basic theory of the IGM DM in the context of FRBs and define DM estimators which will be used in this work.



 §.§ Dispersion Measure


The multifrequency FRB pulses disperse while traveling through an ionized medium due to their interaction with the free electrons along the way. The time delay in the arrival of the signal at frequency ν, Δ t ∝ν^-2 DM; where DM is the line-integral of the free electron density. The total observed time-delay/DM will have contributions from the host galaxy (DM_host), the Milky Way galaxy including the halo (DM_MW), and the IGM (). For DM_MW we have reasonably good Galactic maps, and this can be largely removed from the data. Further, DM_host is reduced by a factor of (1+z)^-1 in the observer frame and so is suppressed when considering high-z events, whereas DM_ IGM increases with z. Hence, in this work we only focus on studying  and ignore any further discussion on other DM components, unless stated otherwise.

The  for an FRB, located at an angular position θ and the redshift z, is <cit.>

    (θ, z) = c ∫_0^z dz^'n_e(θ, z^')/(1+z^')^2 H(z^') ,

where c is speed of light in vacuum and n_e(θ, z) is the proper number density of free electrons. H(z)=H_0√(Ω_m 0 (1+z)^3 + Ω_Λ 0) denotes the Hubble parameter with H_0, Ω_m 0 and Ω_Λ 0 respectively being the present day values of Hubble constant, matter density parameter and dark energy parameter. Baryons being the primary source of free electrons during and after EoR, we can write n_e in terms of the baryon density parameter Ω_ b0 and recast eq. (<ref>) as

    (θ, z) = 13/16   3 c H_0 Ω_ b0/8 π G m_ H×
       ∫_0^z dz^'(1+z^') Δ(θ, z^') (θ, z^')/√(Ω_m0(1+z^')^3 + Ω_Λ0) ,

where G is the gravitational constant, m_ H is the mass of hydrogen atom, Δ is the matter overdensity and  denotes the ionization fraction of the IGM. Here, Δ includes the effects of evolution of the underlying matter density in the IGM whereas  is controlled by the photon field responsible for ionizing the IGM. On large scales, =0 before the reionization starts and it eventually approaches unity towards the end of EoR when the IGM is almost completely ionized. We obtain both Δ and  from our simulations which we discuss in a later section. The treatment of eq. (<ref>) assumes the hydrogen and helium constitute almost the entire baryons in the universe, and the helium being 25% of it by mass. Our model of IGM also assumes that the ionization of  to  occurs concurrently with  reionization.

The LoS path which the FRB signal traverses during post-reionization is effectively very large (∼ 6000   Mpc), adding a larger contribution to total . Post-EoR contribution acts as a nuisance since we are only interested in the impact of the reionization to the DM measurements. Hence, we restrict the lower limit of the integral (eq. <ref>) to the end of reionization z_ end and define : 


    (θ, z) = 13/16   3 c H_0 Ω_ b0/8 π G m_ H×
       ∫_z_ end^z dz^'(1+z^') Δ(θ, z^') (θ, z^')/√(Ω_m0(1+z^')^3 + Ω_Λ0) ,

where z ≥ z_ end. We estimate the sky-averaged mean (z) = ⟨(θ, z) ⟩_θ and the sample variance (z) = ⟨(θ,z) ⟩_θ^2 - ^2(z) numerically for a given z during EoR using simulations.



 §.§ Structure Function




We aim to outline how FRBs can be used as a probe for the characteristic size of the ionization () bubbles in the IGM. One can certainly expect the DM measurements of the two nearby FRBs to be correlated <cit.> as they traces the underlying structures. To this end, we define the DM structure function, for a given LoS at θ and redshift z, as

    Ξ(δθ, z) = ⟨ [(θ+ δθ, z) - (θ, z)]^2 ⟩_θ, δθ ,

where δθ is a small angular separation for all nearby LoS and ⟨⋯⟩_θ, δθ denotes double average – first, over different rotations by an amount δθ around every θ, and next, average over different LoS directions θ. This definition utilizes the assumption that the sky is statistically homogeneous and isotropic at any particular z, which leaves Ξ a function of the magnitude δθ and z only. The dependence of Ξ on , which itself is an integrated quantity, makes it unsuitable to directly conclude anything about the ionized bubble sizes and their growth. Hence, we use the derivatives of Ξ which probe local IGM properties. ∂Ξ/∂δθ encompasses the information about how fast the structures decorrelates on the sky plane. However it still has integrated effects along the LoS, and we therefore compute the second order derivative . This has both the instantaneous and local information about the structures and their scale information. We will demonstrate how the landscape of  corresponds to the different reionization histories and morphologies in the (δθ, z) plane. Later, we marginalize  over δθ and z one at a time. Marginalization makes it easier to understand the behaviour of  as a function of either z or δθ irrespective of the other variable and requires fewer observed bursts to be determined observationally. We finally compute the average of the derivatives and their marginalized values over various LoS θ. This provides us with the information of the mean sizes of the ionized regions in the IGM.







§ TOY MODEL SIMULATIONS


We demonstrate the impact of the ionized bubble sizes and the rate of reionization on the estimators mentioned in <ref>, allowing us to gain intuition and test the general validity of the approach. We use an approximate and simplistic toy model of reionization to simulate the ionization () field light-cones (LCs). The toy model assumes all the ionized bubbles have the same radii and that everything inside the bubbles is completely ionized and anything outside is completely neutral. We divide the whole LC boxes into reasonably thin slices along the LoS axis and fill them with a number of spherical bubbles matching the average input ionization fraction for every z slice. We place the bubble centers in the slices uniform randomly and allow overlap between them. A LC thus created will be a binary ionization field where Δ(θ, z)=1 and (θ, z) is either 0 or 1. This field has basic differences from the LC obtained from the real simulations (see <ref>) where additional fluctuations in the free electron density arise due to perturbations in the underlying matter density field. Furthermore, as opposed to the toy model, the realistic reionization model has inherent temporal growth in the bubble sizes apart from their percolation. 

We simulate the toy model LCs (see Figure <ref>) within a comoving volume (500 × 500 × 1500) [h^-1  Mpc]^3 that is divided into (600 × 600 × 1800) cubic voxels, accordingly. This particular choice of LC volume is to roughly match with our reionization simulations as described below (<ref>). We use the asymmetric tanh(z) form of reionization history <cit.> for our toy models. This form of the history closely mimics the histories found in simulations. We fix the origin of the toy model LCs at z=6 assuming reionization ends by then. The other end of the toy model LC boxes extends up to z≈ 15.




 §.§ Dependence on Reionization History


We study the effect of different reionization histories on the  and the other derived estimators, as defined above. We generate three toy models with `Faster', `Fiducial' and `Slower' reionization histories as shown in Figure <ref> with their corresponding DMs. We mimic the ionization histories by varying the number of bubbles in each slice of the LCs with bubble radius being fixed at 10 grid units (≈ 12  Mpc). We choose the reionization mid point at z_ mid = 7.5, 8.5, 9.0 and the corresponding reionization window to be δ z = 1.0,  2.0,  2.5 (i.e. the reionization to end at the same time for these three toy models) corresponding to the `Faster', `Fiducial' and `Slower' reionization histories, respectively. We fix the asymmetry parameter α = 3.

As shown in Figure <ref>, there is a slight offset between , the mean estimated over all the 600^2 grids on the transverse plane of the box and the average DM estimated using the x̅_ into eq. <ref> and ignoring the fluctuations along different LoS. As DM is a cumulative estimator, its mean value rises rapidly with z where most of the reionization is happening. At higher z it saturates, as there are no more free electrons to contribute to it. The asymptotic difference between the two mean DM estimators (solid and dashed lines) increases monotonically from faster to slower reionization at any z. This happens because the fluctuating structures exists for a larger LoS distance in the case of slower reionization history. We also show the corresponding 1σ sample variance around . Fluctuations due to the binary ionization field are not able to cause any significant (>1σ) deviations between (z) and DM(x̅_(z)) for all three histories considered here. The deviation should be enhanced for more realistic reionization LC where both  and Δ contribute to the LoS fluctuations in . Also, contribution to σ_ DM accumulated from z<6 will increase the spread in DM_ EoR.




Figure <ref> shows the ratio σ_ DM/ that qualitatively follows  for z ≳ 7. It first increases rapidly and then saturates towards large z as the ionized regions disappear. However, we find that σ_ DM/ increases sharply towards the end of reionization (z<7). This is because  has a value around zero at z≈ 6 since we ignore the contribution from the low-redshift IGM. Starting at z≈ 7,  increases more rapidly while moving towards higher z than the fluctuations do and finally the ratio saturates. The saturation redshift varies depending on the reionization history. It is at low redshift for the faster reionization and vice-versa. σ_ DM/ is larger for the faster reionization history and vice-versa, which indicates a direct mapping between the relative (to the mean) fluctuations in the DM_ EoR with the emergence and sizes of the structures in the IGM.

Figure <ref> also shows the derivative  d/dz which directly traces the local electron density. After Gaussian smoothing, the derivatives are roughly the same towards the end of reionization (z < 6.5) where the IGM has roughly indistinguishable properties. However, the derivatives beyond z_ mid apparently encodes the information of the reionization history in its slope when plotted against z in a log-log plane. The slower history has a shallower slope and it increases gradually towards fiducial and faster histories. This is simply because there are more ionized bubbles for the slower history causing a larger change in  at higher z.





Figure <ref> shows contours of the derivatives of the structure function  on the δθ-z plane. We consider only 400 LoS at every equidistant comoving slice to compute  instead of all available 360,000 LoS per comoving slice. Our choice of 400[Our mock simulations have access to 400 FRBs/slice even at larger redshifts. However, in reality, the number of FRBs might drop significantly down at larger z slices (see Figure <ref>).] FRBs per comoving slice is closer to real observations and computationally tractable. On the other hand, we believe that it is a good number for statistical convergence of . We finally bin-average our estimates of all slices within equispaced z bins of width 0.5 for further use. We consider  as a 2D surface and plot the contours which include top 10%, 30%, 50% and 67% of the total area under the  surface. For each reionization history we depict z_ mid and the angular size corresponding to the radius of the individual bubble in our toy model. The peak (10% contour) gradually shifts to a smaller z while moving from the slower to the faster history, approximately tracking the change in z_ mid, and the contours gets more squeezed along the z axis. This clearly indicates that the peakedness of the  landscape along z is directly connected to the reionization window. There is no considerable change in the extent of the contours along the δθ axis, which is expected since we are using bubbles of fixed radii here.

We next marginalize  over δθ. The marginalized result, ∫ ()  dδθ, is shown in Figure <ref> as a function of z. ∫ ()  dδθ peaks roughly around z_ mid. The decrease towards the end of EoR is sharper and roughly independent of reionization history. However, the drop towards higher z is related to the reionization history. The drop is slower for the slower reionization history and vice-versa. That could be simply related to the rate of emergence of the ionized bubbles in the IGM, and once the reionization is roughly around its midway, the percolation of bubbles makes it insensitive to the history.




 §.§ Dependence on Bubble Sizes




We assess the impact of ionized bubble size on our estimators by fixing the reionization histories to the fiducial case and varying the bubble radius R. We consider R=5, 10 and 20 grid units which corresponds to 6, 12 and 24  Mpc, respectively. We have repeated the same analysis as in the <ref>. The three simulations perfectly agree with their common input reionization history, leading to the similar  values. However, variations in bubble size influence the fluctuations in the DM_ EoR estimates, and consequently σ_ DM. 




Figure <ref> shows the contours of . Similar to Figure <ref>, we have used 400 LoS per comoving slices and bin them within a redshift window Δ z=0.5. The contours are nearly unchanged along z-axis for the different R values. The peak of the  surface (depicted by 10% contour) shifts towards larger δθ values with increasing R, approximately tracking the change in the bubble size. The derivatives decrease for δθ greater than the angular bubble size θ_R, as the correlation between the structures decays out. Whereas for scales less than θ_R (see the right panel), the derivatives decrease as the points are tightly correlated and Ξ itself is consistently small.



We can understand this with a simple argument. We distribute the bubbles uniformly in each z slice, following a Poisson distribution while preventing any overlap. Considering a redshift slice, the mean DM would approximately be DM_1× N_ col, where DM_1 is the contribution from a single bubble and N_ col is the average number of bubbles appearing along a LoS. Hence, the corresponding spread in the DM along different LoS would be σ_ DM≈ DM_1 × N_ col^1/2. Since we keep the ionization fraction of slices (and hence mean DM) constant while changing the bubble radius, the total number of the bubbles within the slice varies as N ∝ R^-3. Considering a cubical slice N_ col∝ N^1/3∝ R^-1. Since DM_1∝ R, we see that σ_ DM∝ R^1/2. This scaling is consistent with the results plotted in the left panel of Figure <ref>.

We plot σ_ DM/ and its z derivative as a function of z in Figure <ref>. The trend is qualitatively similar to that shown in Figure <ref>. We find a sharp turnover and rapid rise in σ_ DM/ by the end of reionization (z<7) due to near zero value of . However for z>7, σ_ DM/ increases almost linearly towards higher z and gradually starts saturating beyond z ≈ 10. As the saturation point strongly depends on the reionization history, hence the saturation knee is almost at the same redshift for all the three R values. However, the variation in σ_ DM/ is solely due to variation in σ_ DM with R. Overall, σ_ DM increases with increasing R as shown in the Figure. The reason is clear as filling the IGM with the smaller bubbles would make the distribution of the ionized regions roughly uniform and homogeneous, and therefore the fluctuations between the different LoS would be less and vice-versa. 

The derivative d(σ_ DM/)/dz is a more local quantity as shown in Figure <ref>. The derivatives drop very sharply approaching -∞ for z<7 as there is a rapid increase in σ_ DM/ with decreasing z. However, we see a power-law decrement in the derivative with increasing z for z ≥ 8. The slope of this power-law drop is roughly the same for the three R values, although the magnitude scales with R in a similar way as for σ_ DM/. 

We again marginalize  here, but this time along the z-axis, to obtain ∫ ()  dz as shown in Figure <ref> for the three R values. We observe a similar qualitative trend as seen in the peaks of the contour plots. The peak of ∫ ()  dz shifts to a larger δθ for large bubble sizes. We also note that it decreases for δθ larger than θ_R. ⟨θ_R⟩ is computed by marginalizing the θ_R(z) corresponding to the respective characteristic bubble size R in the z range of interest.





§ EOR SIMULATION USING GRIZZLY





We now use the simulated EoR ionization field LCs to estimate the DM and other related quantities to examine how the different reionization models affect them. We simulated three EoR  LCs corresponding to the three different reionization histories x̅_ as shown in Figure <ref>. In the following sections, we briefly describe the method used to simulate the EoR scenarios followed by all the corresponding estimates.






 §.§ EoR Simulations


The EoR  LCs used here are constructed by stitching thin slices from the coeval  cubes simulated at several different redshifts in chronological order. We simulate  coeval boxes using a CD-EoR code, grizzly <cit.>, which is based on a 1D radiative transfer technique. This algorithm takes the dark-matter density field and the corresponding halo catalogue from a N-body simulation to produce a  map at a redshift for a particular astrophysical source model. 

We use the dark-matter density fields and the associated halo catalogues obtained from the PRACE[Partnership for Advanced Computing in Europe: <http://www.prace-ri.eu/>] project, PRACE4LOFAR. The input dark-matter distributions are generated using dark-matter only N-body simulation code cubep^3m [<http://wiki.cita.utoronto.ca/mediawiki/index.php/CubePM>] <cit.>. The 3D density cubes have comoving volume [500 h^-1  Mpc]^3 <cit.> which are gridded into [600]^3 voxels. The dark-matter particle distributions have been used to find the collapsed halos using spherical overdensity halo finder <cit.>. The minimum halo mass in the PRACE4LOFAR simulation is ∼ 10^9, which corresponds to ≈ 25 dark-matter particles. We simulated 63 coeval dark-matter cubes between a redshift range 6.1 ≲ z ≲ 15.6 with an equal time gap of 11.4  Myr.

We consider an EoR source model where the dark-matter halos with masses larger than 10^9 host UV photon-emitting galaxies.
We assume that the stellar mass of a galaxy (M_⋆) is related to the host dark matter halo mass M_ halo as M_⋆∝ M^α_s_ halo.  We tune the ionization efficiency (ζ) so that the reionization ends at z∼ 6. Note that all reionization models considered in this study are inside-out in nature. Our fiducial grizzly model corresponds to a choice of  α_s=1 and spans from redshift 6 to 15.6. We consider a rapidly evolving Faster and a slowly evolving Slower reionization scenario which correspond to α_s=2 and  α_s=0.1 respectively.

Next, we produce coeval cubes of  at 63 redshifts between 6.1 and 15.6 for the aforementioned source model. We refer the readers to <cit.> and <cit.> for more details about these calculations. Finally, we used these coeval cubes of  to create the LC which accounts for the evolution of  with redshift. The detailed method to implement the LC effect can be found in <cit.>. The reionization histories for the three different EoR scenarios are shown in Figure <ref> while we present the corresponding LCs in Figure <ref>. The LCs clearly show the difference in the patchiness of the reionization scenarios.




 §.§ grizzly Simulation Results


We present the results for the LCs (Figure <ref>) simulated using grizzly. We exclude the contributions coming from post-reionization IGM and the host galaxy, assuming these would be perfectly measured and subtracted in the future as we expect to detect a larger population of FRBs at lower redshifts (Figure <ref>). This optimistic assumption primarily allows us to investigate the impact of  bubble sizes and reionization histories on our estimators by evading low-z contributions.

In Figure <ref>, we show the  estimated from the LC simulations for `Slower', `Fiducial' and `Faster' reionization histories. The solid lines show the mean estimated using all the LoS (here 360 000 grid points) in the simulations and the shaded regions around them are the respective 1σ deviations due to cosmic variance. We consider small z-bins (Δ z = 0.5) while computing the mean and the sample variance. We overplot the  estimates calculated using the average ionization fraction (Figure <ref>) for the three reionization scenarios.

(z) begins with a near zero value (with low-z IGM contributions subtracted) around the end of reionization (z ≈ 6) and increases almost linearly with z until when ionization is sufficiently small (x̅_≈ 0.1), where it starts to plateau. This is qualitatively consistent across both the realistic and toy model reionization histories (see Figure <ref>). The rate of linear rise at lower z and the saturation value is highest for the slower reionization history, and these decrease monotonically for faster histories. The `knee' in the (z) curves occurs at a larger redshift for slower reionization histories. The saturation values of DM differ by more than 1σ across all three histories.



The (z) and DM(x̅_(z)) are significantly distinct for all three reionization histories due to IGM electron density fluctuations driven by the formation and growth of ionized bubbles, as well as the underlying density perturbations. The bias due to the fluctuations are highest at larger z (z ≳ 12) where the DM saturates, and the differences between the two estimates decrease rapidly towards the end of reionization. This is because the merging and percolation of the ionized bubbles typically wash out the fluctuations towards the mid and advanced stages of reionization. The difference is smallest (≲ 1σ) for the Faster reionization scenario, where large ionization bubbles appear suddenly at lower redshifts (see bottom panel of Figure <ref>) and quickly percolate and fill up the entire IGM. In this case, the IGM  is more patchy and hence small-scale fluctuations do not prevail. Relatively small ionized bubbles form for slower histories and takes a longer time to fill the entire IGM. This leads to relatively more small-scale fluctuations (see Figure <ref>), thereby resulting in a larger difference (≳ 2σ) between the (z) and DM(x̅_(z)) (solid and dashed lines) for the slower reionization history.



As the average  is an integrated quantity, we compute its derivative d/dz to extract instantaneous information at a given stage of reionization. We show the derivatives in Figure <ref> for the three histories. d/dz starts with a high value which is similar for the three histories. This simply indicates that the electron distribution in the IGM is roughly the same for all three scenarios nearing the end of reionization. The large values of the derivative at z ≲ 8 are an implication of the fact that x̅_ changes rapidly during the mid and end stages of the reionization. However, the difference in how fast x̅_(z) evolves causes the derivatives to distinctly vary for the three histories considered here. The derivative for the faster scenario has small values but is the steepest among the three. The value of the derivatives increases for the fiducial and further to slower reionization histories, while the slopes of the curves decreases.

Figure <ref> depicts σ_ DM/ as a function of z. This shows how the sample variance error in the measurement is related with the mean DM. The qualitative behavior of the ratio is similar for the three reionization histories. The `knee' of saturation in σ_ DM/ depends on the z_ mid as seen in toy models (top panel of Figure <ref>). However, contrary to the toy models (Figures <ref> and <ref>), the ratio σ_ DM/ decreases with z. The 1σ cosmic variance for all the three histories are roughly similar, hence the ratio scales inversely with . This could be because of the fact that in the real simulations the electron distribution in the IGM has an additional contribution coming from the underlying matter density field which is absent in our toy models. This needs further detailed investigation with a more sophisticated toy model and we defer it to a future work.



Figure <ref> shows contours of the derivatives of the structure function  on the δθ-z plane, similar to Figures <ref> and <ref>. We estimate  from 400 LoS per comoving slice and bin them within bins of Δ z = 0.5. The behaviour is qualitatively similar to the toy models. The contours are tightly squeezed towards lower z (10% contour at z≲ 7.2) for the faster reionization history and they gradually expand for the fiducial and slower models (10% contour at z≲ 8.1). The effects of the different IGM topologies are also clearly evident from the δθ extent of the contours. The faster history has larger bubbles (hence patchy IGM) resulting in extended angular correlations. In contrast, the angular correlations are smaller (depicted by the squeezed contours along δθ) for slower histories where the bubbles are smaller in sizes. This trend closely resembles what is observed in Figure <ref>.






 §.§ Impact of Non-uniform Distribution of FRBs and Post-EoR




All the analyses presented above are independent of the distribution of FRBs in a LC box. There are various factors that make the distribution of FRBs non-uniform in the sky-plane and across redshift. We repeat the analysis for the grizzly simulated LCs considering a realistic redshift distribution of the sky-averaged number of FRBs. Following the estimates in <cit.>, we used the cumulative observable redshift distribution as shown in Figure <ref>. To make our estimates more realistic, we also consider the FRBs to be correlated with the overdensities in our matter density LC. This is well motivated because the overdense regions are the places which can host radiation sources and therefore also FRBs. Ideally, one should correlate the FRBs with the halo field which is precisely the location of sources. However, here we lack the exact halo locations and instead use the matter overdensity peaks as a proxy. We denote by N_ FRB^ tot the total number of FRBs observed in the range 6≤ z≤ 15. We divide the whole EoR redshift range (z=6 to 15) into bins of Δ z = 0.5 and estimate the number of FRBs per bin according to the distribution (Figure <ref>) for a given N_ FRB^ tot. For each bin, we randomly pick up the grid points which are biased by the high-dense regions. We finally use those grid points to estimate the average DM and corresponding dispersion for all the redshift bins.



Our main findings in Figure <ref> depict the minimum number of total FRBs required to distinguish between the three different reionization histories. In every panel, the solid lines denote the ensemble mean of the sky-averaged DM varies as a function of N_ FRB^ tot. We compute the DM from a random sample of N_ FRB^ tot distributed across the EoR redshift range. Later, we consider 1000 such independent realizations of the distribution of N_ FRB^ tot to compute the ensemble mean ⟨ DM⟩ and the 3σ cosmic variance as shown by the shaded regions around the lines. Figure <ref> shows ⟨ DM(z_ mid)⟩ corresponding to z_ mid for the individual histories. This inherently assumes that we have redshift information about the FRBs (within an uncertainty of Δ z = 0.5). However, if we drop this assumption and just assume that the FRB observed is somewhere within the duration of EoR, we can compute the marginalized sky-averaged DM defined as 


    DM_ margin = ∫_z_ min^z_ max DM_ EoR (z)  dP/dz  dz/∫_z_ min^z_ max dP/dz dz ,

where dP/dz is the probability of detecting an FRB per unit redshift (as per Figure <ref>), z_ min=6.0 and z_ max=14.5 are respectively the lower and upper bounds of the EoR redshift window. This too is shown in Figure <ref>.  


We introduce the contribution from post-EoR (0 ≤ z ≤ 6) and consider three different scenarios – (i) Optimistic, where the contribution from the post-EoR is perfectly modelled and completely removed from the data, (ii) Moderate, where the post-EoR contribution is imperfectly removed from the individual DM(z) such that the residual low-z contribution remains within 0.5 σ_ DM at z=6, and (iii) Pessimistic, where the post-EoR contribution is completely present for every FRB measured. Note that these three scenarios are differentiated based on the how accurately we know the DM(θ, z<6) for every observed LoS θ above z=6. The assumption that we either know the redshifts of every FRB or at least can separate EoR FRBs from the post-EoR FRBs confidently is implicit within all the three scenarios and is itself non-trivial. Apart from the fact that the boundary between EoR and post-EoR is still fuzzy, the erroneous/indeterminate FRB redshifts can cause a spill over the boundary. This may bias our estimators and we plan to investigate this in greater detail in future studies.

We simulate a matter density LC within the redshift range 0 ≤ z ≤ 6 (comoving distance 8542  Mpc), and assume that the electrons linearly follow the underlying density contrast. We use the publicly available gadget4[<https://wwwmpa.mpa-garching.mpg.de/gadget4/>] <cit.> to simulate the dark matter density fields within comoving cubes of volume [500 h^-1  Mpc]^3 at the lower redshifts. We choose the same cosmological parameters as mentioned in the end of <ref>. We also set the resolution (for Particle-Mesh) of the box to match those in the grizzly simulations. We simulated the particle distribution at 66 redshifts within the range 0 ≤ z ≤ 7 roughly equidistant by 200  Myr. We generated the density coeval boxes on 600^3 grids and then used them to make the density LC. We use weighted linear interpolation scheme to interpolate the density fields at the desired redshift slices from the redshifts at which the coeval boxes are generated. We cannot simulate a large box spanning the whole range up to z≤ 6. Hence, we need to repeat the boxes along the LoS (z-axis) while creating the final LC. 

The contribution from low-redshifts increases both  and σ_ DM(z), making the separation between the two different histories more challenging. In Figure <ref>, ⟨ DM(z_ mid)⟩ for the three different histories starts at different N_ FRB^ tot. This is because, for the slower history, z_ mid is larger than that of the other histories, requiring a higher N_ FRB^ tot to populate the corresponding z-bin with at least two FRBs. ⟨ DM(z_ mid)⟩ converges very quickly for N_ FRB^ tot≈ 10; however the corresponding dispersion is large and it decreases with increasing N^ tot_ FRB, as expected. N_ FRB^ tot≥ 20 is sufficient to distinguish ⟨ DM(z_ mid)⟩ at 3σ for the optimistic scenario. This lower bound remains the same to distinguish the slower reionization history with the faster one for the moderate scenario. However, it takes nearly N_ FRB^ tot = 150 and 40 to discern ⟨ DM(z_ mid)⟩ of slower-fiducial and fiducial-faster pairs of histories, respectively. Considering the pessimistic case, it takes N_ FRB^ tot = 80, 150, 600 to distinguish between faster-slower, slower-fiducial and fiducial-faster pairs of histories, respectively.

Figure <ref> also shows the variation of the ⟨ DM_ margin⟩ with N_ FRB^ tot. Considering the optimistic case, we need roughly 40 FRBs identified during the EoR to discern slower and faster reionization histories at 3σ with ⟨ DM_ margin⟩. Whereas the  N_ FRB^ tot slightly increases to 60 and 90 for faster-fiducial and fiducial-slower pairs of histories, respectively. These numbers respectively increase to 90, 200 and 300 for the moderate case as shown in the middle panel. Finally for pessimistic scenario, we need N_ FRB^ tot≈ 220, 600, 1000 to discern (at 3σ) slower-faster, faster-fiducial and fiducial-slower pairs of reionization histories, respectively. The numbers we found are realistic, and it would be possible to detect many FRBs during EoR using the upcoming SKA-mid. 





§ DISCUSSION

Understanding the epoch of reionization (EoR) is a crucial step in learning about one of the most important eras in the cosmic history, when it transitioned from being devoid of any stars, consisting of cold and neutral gas, to hot, ionized, and teeming with the objects we see today. The first sources are supposed to drive the whole process of reionization, therefore studying the IGM during EoR can be connected with the properties and emergence of the first structures. There are many direct and indirect contemporary probes such as the redshifted 21-cm signal from  in the IGM during EoR and, Thomson scattering optical depth of the CMB photons, high-z quasar spectra, Ly-α systems at high-z. However these probes are limited by their own challenges to date. In this work, we propose to use the dispersion measure (DM) of the fast radio bursts (FRBs) from the high-z to probe the IGM during EoR.

The dispersion introduced in the FRB pulse, while it travels through the ionized medium, can be used in probing the electron distribution along the line-of-sight (LoS) during EoR. We demonstrated the use of the sky-averaged DM and its derivatives to discern between the different reionization histories. Beyond this, we primarily aim to make use of the sky-averaged and the angular dispersion in the DM estimates to extract information on the ionization bubbles during the EoR. Using a toy model (Figure <ref>) of the binary ionization field (within the range 6 ≤ z ≤ 15), we see that the (z) first increases (starting from the end of EoR) roughly up to the mid-point of the reionization, z_ mid, and then tend to saturate as the ionized regions decreases towards the initial stages of the reionization (Figure <ref>). The derivative d/dz directly traces how fast the reionization progresses (Figure <ref>). The all-sky variance σ_ DM/ is larger for the reionization scenario where the bubble sizes are larger and vice-versa (Figure <ref>). We also compute  where Ξ is the angular dispersion (eq. <ref>) at any redshift. We demonstrated that the contours are elongated or squeezed along z and δθ depending respectively on the varying reionization history (Figure <ref>) and bubble sizes (Figure <ref>). The impacts are clearly prominent for the marginalized structure function derivatives (Figures <ref> and <ref>).

We also analyzed a more realistic reionization light-cone (LC) (Figure <ref>) generated using grizzly 1D radiative transfer code for three different reionization histories ending at the same z. The behaviour of  is the same as in the toy model. We find that the different LoS variance plays a significant role for realistic IGM (Figure <ref>) and biases the average (z) by more than 1σ as compared to the case when the DM is computed using the averaged ionization fraction x̅_. The slope of the d/dz is directly sensitive to the reionization history (Figure <ref>). σ_ DM/ is also dependent on the reionization history as well as indirectly on the IGM morphology (bubble sizes) in an intermingled way. Figure <ref> clearly shows that  is sensitive to the reionization window as well as the typical bubble sizes. The contours are squeezed along z axis and elongated along δθ axis for faster reionization history which has relatively larger ionized bubbles and vice-versa. Our initial analyses are rather optimistic where we have considered the FRBs to be located uniformly at every grid points in our reionization LC with their redshifts known within an uncertainty of Δ z =0.5. We also assumed that the contribution from the low-redshift (z<6) has been perfectly removed from the DM measurements. 

We next consider observationally more realistic situation where the FRBs abundances vary with redshift (Figure <ref>) and they are more clustered at the highly dense regions on the sky plane. This biases  relative to our initial results and also introduces more LoS dispersion, particularly at high redshifts where the FRB abundance drops rapidly. Taking realistic estimates of the FRB rate evolution with z, one requires ≲ 100 FRBs to be distributed across the whole EoR window in order to discern the reionization histories at 3σ (see the left column of Figure <ref>) using the mean DM only. This is assuming we have removed contribution from lower redshifts (z<6) from the DM measurements which is an `optimistic' case. Considering a `pessimistic' case where the low redshift contribution is present, we find that the numbers could shoot as high as 200-600 if we focus on the mid-reionization redshift bin (within uncertainty of Δ z = 0.5. Using the marginalized DM to discern between the reionization histories at 3σ might require roughly ∼ 1000 FRBs in total during EoR (see right column of Figure <ref>).

The numbers presented above corresponds to a particular choice of telescope sensitivity and FRB population models and are expected to vary if we change them, however we expect them to stay within the order of magnitude. We plan to include these effects gradually in our future work along the line. We have successfully demonstrated the potential of the derivative of the structure function  and of σ_ DM as probes of the ionization bubble sizes along with the reionization history. The , being a derivative, suffers from large variance, and we need a large number of FRBs (N_ FRB^ tot≲ 100 000) within range 6 ≤ z ≤ 15 to suppress the variance significantly down. The computation of  and their marginalization here considers a uniform distribution of the FRBs on the regular comoving grid. This is only done for convenience, and the derivatives of the structure function are well defined also when the sample is very uneven with z and/or δθ. In reality, FRBs should be associated with the galaxies which are generally clustered around the high-density peaks that gets ionized first. Therefore, it is highly probable to find more FRBs within the ionized regions, and that might help to compute the structure function in the vicinity of the ionized bubbles. We also plan to investigate more deeply into this estimator which will be very useful in probing the ionized regions around the sources.





§ ACKNOWLEDGEMENTS

AKS, RG and SZ acknowledge support by the Israel Science Foundation (grant no. 255/18). AKS is also supported by National Science Foundation (grant no. 2206602). RG also acknowledges support from SERB, DST Ramanujan Fellowship no. RJF/2022/000141. PB is supported by a grant (no. 2020747) from the United States-Israel Binational Science Foundation (BSF), Jerusalem, Israel, by a grant (no. 1649/23) from the Israel Science Foundation and by a grant (no. 80NSSC 24K0770) from the NASA astrophysics theory program. PK’s work is  funded in part by an NSF grant AST-2009619 and a NASA grant 80NSSC24K0770. 




