We propose different estimators to probe the epoch of reionization EoR intergalactic medium IGM using the dispersion measure rm DM of the FRBs. We consider three different reionization histories which we can distinguish with a total of lesssim 1000 rm DM measurements during EoR if their redshifts are known. We note that the redshift derivatives of rm DM are also directly sensitive to the reionization history. The major point of this work is exploring the variance in the rm DM measurements and the information encoded in them. We find that the allsky average overlinerm DMz gets biased from the LoS fluctuations in the rm DM measurements introduced by the ionization of IGM during EoR. We find that the ratio sigmarm DMoverlinerm 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 rm DM encodes the information about the duration of reionization and the typical bubble sizes as well. We establish the usefulness of variances in rm DM using toy models of reionization and later verify it with the realistic reionization simulations.
['Astrophysics - Cosmology and Nongalactic Astrophysics', 'Astrophysics - High Energy Astrophysical Phenomena']
Asking Fast Radio Bursts FRBs for More than Reionization History
{'Asking Fast Radio Bursts (FRBs) for More than Reionization History': "Abinash Kumar Shaw , 1, 2 Raghunath Ghara , 1, 3 Paz Beniamini , 1, 4, 5 Saleem Zaroubi , 1, 4, 6 and Pawan Kumar 7  \n- 1 Astrophysics Research Center of the Open University (ARCO), The Open University of Israel, 1 University Road, Ra'anana 4353701, Israel  \n2 Department of Computer Science, University of Nevada Las Vegas, 4505 S. Maryland Pkwy., Las Vegas, NV 89154, USA 3 Department of Physical Sciences, Indian Institute of Science Education and Research Kolkata, Mohanpur, WB 741 246, India 4 Department of Natural Sciences, The Open University of Israel, 1 University Road, Ra'anana 4353701, Israel 5 Department of Physics, The George Washington University, 725 21st Street NW, Washington, DC 20052, USA 6 Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700AV Groningen, The Netherlands 7 Department of Astronomy, University of Texas at Austin, Austin, TX 78712, USA  \n(Dated: September 6, 2024)  \nSubmitted to ApJ", '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.  \nKeywords: Reionization (1383), Radio transient sources (2008), Intergalactic medium (813), Largescale structure of the universe (902), Radio bursts (1339)', '1. INTRODUCTION': "According to the current understanding of cosmology, our universe transitioned from being a highly coldneutral 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) (see e.g., Loeb & Barkana 2001; Zaroubi 2013; Pritchard & Loeb 2012). This window of transition is termed as the Epoch of Reionization (EoR). The study of the EoR  \nCorresponding author: Abinash Kumar Shaw  \nabinashkumarshaw@gmail.com  \nis 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 (see e.g., Shaw et al. 2023).  \nOur 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 (e.g., Planck Collaboration et al. 2020) and the Gunn-Peterson troughs in the highz quasar spectra (e.g., Becker et al. 2001; Fan et al. 2006; Gallerani et al. 2006; Becker et al. 2015; Bosman et al. 2022;  \nD'Odorico et al. 2023; Gaikwad et al. 2023; Greig et al. 2024; Asthana et al. 2024; Spina et al. 2024). Additional constraints on the timeline of EoR come from the recent observations of the highz Lyα emitters (e.g., Hu et al. 2010; Kashikawa et al. 2011; Ota et al. 2017; Ishigaki et al. 2018; Morales et al. 2021; Bruton et al. 2023; Nakane et al. 2023) and their clustering measurements (e.g., Faisst et al. 2014; Santos et al. 2016; Wold et al. 2022), Lyman break galaxies (Mason et al. 2018; Hoag et al. 2019; Naidu et al. 2020; Bolan et al. 2022), and the Lyα damping wings in the highz quasar spectra (e.g., Ba˜nados et al. 2018; Davies et al. 2018; ˇ Durovˇc'ıkov'a et al. 2020; Wang et al. 2020; Yang et al. 2020; Umeda et al. 2023; ˇ Durovˇc'ıkov'a et al. 2024). 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 darkgap/pixel statistics) (e.g., McGreer et al. 2014; Keating et al. 2019; Kulkarni et al. 2019; Zhu et al. 2021, 2022; Bosman et al. 2022) 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.  \nAfew earlier works have demonstrated the potential of highly energetic astrophysical events such as gamma-ray bursts (GRBs) (e.g., Ciardi & Loeb 2000; Ioka 2003; Inoue 2004; Lidz et al. 2021) and fast radio bursts (FRBs) (e.g., Beniamini et al. 2021; Hashimoto et al. 2021; Heimersheim et al. 2022) 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 (see Petroff et al. 2022, for a review). Empirical studies on the all-  \ny 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 (The CHIME/FRB Collaboration et al. 2021), and the rate is expected to increase significantly at lower fluence thresholds.  \nDespite their enigmatic origin, a recent discovery of a galactic FRB (Bochenek et al. 2020; CHIME/FRB Collaboration et al. 2020; Bera et al. 2024) 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 (e.g., Popov & Postnov 2010; Falcke, Heino & Rezzolla, Luciano 2014; Ioka 2020; Wadiasingh et al. 2020; Zhang 2022; Beniamini & Kumar 2023; Chen et al. 2023; Totani & Tsuzuki 2023). 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 highz FRBs (see the introduction of Beniamini et al. 2021, for more details).  \nFRB 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 (Zhang 2018; Kumar & Linder 2019). Conversely, knowing the redshift of the FRBs accurately can be potentially used to estimate baryonic content of IGM during the post-EoR (e.g., McQuinn 2013; Macquart et al. 2020; Lee et al. 2022; Khrykin et al. 2024), probe the epoch of second helium reionization (e.g., Caleb et al. 2019; Linder 2020; Bhattacharya et al. 2021; Lau et al. 2021) and constrain several cosmological parameters (e.g., Deng & Zhang 2014; Zhou et al. 2014; Yang & Zhang 2017; Walters et al. 2018; Pol et al. 2019; Jaroszynski 2019; Wu et al. 2020; Wucknitz et al. 2021; James et al. 2022; Hagstotz et al. 2022). In this work we explore how useful they can be as detailed probes of the EoR.  \nWith the possibility of detecting highz FRBs (e.g., Hashimoto et al. 2020), one can turn their precise DM measurements to probe the sources and IGM during reionization. Recently, a few theoretical studies (Be-  \nProbing 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 (e.g., Ali et al. 2008; Bernardi, G. et al. 2009, 2010; Ghosh et al. 2012), 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: Patil et al. 2017; Mertens et al. 2020, MWA: Barry et al. 2019; Li et al. 2019; Trott et al. 2020, HERA: Abdurashidova et al. 2022; The HERA Collaboration et al. 2022), and the upperlimits are improving gradually.  \nniamini et al. 2021; Hashimoto et al. 2021; Heimersheim et al. 2022) have demonstrated the feasibility of using the DM measurements to extract the reionization history and Thomson scattering optical depth ( τ Th ). Pagano & Fronenberg (2021) and Maity (2024) 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, Beniamini et al. (2021) 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 Heimersheim et al. (2022) where they also estimated τ Th and the mid-redshift of reionization using the average DM. Beniamini et al. (2021) 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, Hashimoto et al. (2021) have shown that the redshift derivatives of DM have potential to directly constrain the reionization history.  \nredshift information and compute the marginalized average DM over the EoR window.  \nThis manuscript is arranged as follows. We define DM and structure function estimators in §2. In §3, we validate our estimators with the toy ionization field model. This is followed by §4 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 §5. The cosmological simulation here uses the cosmological parameter values Ω m = 0 . 27, Ω Λ = 0 . 73, Ω b = 0 . 044 and h = 0 . 7 adapted from Hinshaw et al. (2013).", '2. 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.', '2.1. 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 (DM 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 highz events, whereas DM IGM increases with z . Hence, in this work we only focus on studying DM IGM and ignore any further discussion on other DM components, unless stated otherwise.  \nThe DM IGM for an FRB, located at an angular position θ and the redshift z , is (e.g., Beniamini et al. 2021)  \nDM IGM ( θ , z ) = c z ∫ 0 dz ' n e ( θ , z ' ) (1 + z ' ) 2 H ( z ' ) , (1)  \nwhere 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  \nIn 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 §2.2) 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  \nΩ b0 and recast eq. (1) as  \nfor a given LoS at θ and redshift z , as  \nDM IGM ( θ , z ) = 13 16 3 cH 0 Ω b0 8 πGm H × z ∫ 0 dz ' (1 + z ' )∆( θ , z ' ) x H ii ( θ , z ' ) √ Ω m 0 (1 + z ' ) 3 +Ω Λ0 , (2)  \nwhere G is the gravitational constant, m H is the mass of hydrogen atom, ∆ is the matter overdensity and x H ii denotes the ionization fraction of the IGM. Here, ∆ includes the effects of evolution of the underlying matter density in the IGM whereas x H ii is controlled by the photon field responsible for ionizing the IGM. On large scales, x H ii = 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 x H ii from our simulations which we discuss in a later section. The treatment of eq. (2) 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 He i to He ii occurs concurrently with H i reionization.  \nThe LoS path which the FRB signal traverses during post-reionization is effectively very large ( ∼ 6000 Mpc), adding a larger contribution to total DM IGM . 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. 2) to the end of reionization z end and define DM EoR :  \nDM EoR ( θ , z ) = 13 16 3 cH 0 Ω b0 8 πGm H × z ∫ z end dz ' (1 + z ' )∆( θ , z ' ) x H ii ( θ , z ' ) √ Ω m 0 (1 + z ' ) 3 +Ω Λ0 , (3)  \nwhere z ≥ z end . We estimate the sky-averaged mean DM EoR ( z ) = 〈 DM EoR ( θ , z ) 〉 θ and the sample variance σ 2 DM ( z ) = 〈 DM EoR ( θ , z ) 〉 2 θ -DM EoR 2 ( z ) numerically for a given z during EoR using simulations.", '2.2. Structure Function': 'We aim to outline how FRBs can be used as a probe for the characteristic size of the ionization ( x H ii ) bubbles in the IGM. One can certainly expect the DM measurements of the two nearby FRBs to be correlated (Reischke & Hagstotz 2023) as they traces the underlying structures. To this end, we define the DM structure function ,  \nΞ( δθ, z ) = 〈 [DM EoR ( θ + δθ , z ) -DM EoR ( θ , z )] 2 〉 θ , ̂ δθ , (4)  \nwhere δθ 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 DM EoR , 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 ∂ 2 Ξ /∂z ∂δθ . This has both the instantaneous and local information about the structures and their scale information. We will demonstrate how the landscape of ∂ 2 Ξ /∂z ∂δθ corresponds to the different reionization histories and morphologies in the ( δθ, z ) plane. Later, we marginalize ∂ 2 Ξ /∂z ∂δθ over δθ and z one at a time. Marginalization makes it easier to understand the behaviour of ∂ 2 Ξ /∂z ∂δθ 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.', '3. TOY MODEL SIMULATIONS': 'We demonstrate the impact of the ionized bubble sizes and the rate of reionization on the estimators mentioned in §2, 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 ( x H ii ) 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 x H ii ( θ , z ) is either 0 or 1. This field has basic differences from the LC obtained from the real sim-  \nFigure 1. Image slice from the LC box of the fiducial EoR toy model. Both axes are in comoving units with the axis labelled x being the LoS axis. The bubble radius is 10 grid units ( ≈ 12 Mpc). Black and white regions are respectively completely ionized and completely neutral. Reionization progresses from right to left in this box.  \n<!-- image -->  \nFigure 2. Left: Asymmetrical tanh( z ) ionization histories. Different colors represent three ionization histories corresponding to the simulated toy models. Right: Mean DM estimates and 1 σ fluctuations. Colors represent different reionization histories as shown in the left panel. Solid lines and the shaded regions are the average DM and corresponding 1 σ errors estimated over all grid points available in the simulation volume. Dashed lines correspond to the DM estimated using the mean ionization fraction directly in eq. 3. Three vertical lines mark the redshift of mid reionization.  \n<!-- image -->  \nulations (see §4.1) 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.  \net al. 2024) 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.', '3.1. Dependence on Reionization History': "We study the effect of different reionization histories on the DM EoR and the other derived estimators, as defined above. We generate three toy models with 'Faster', 'Fiducial' and 'Slower' reionization histories as shown in Figure 2 with their corresponding DMs. We mimic the ionization histories by varying the number of bubbles in  \nWe simulate the toy model LCs (see Figure 1) 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 (§4.1). We use the asymmetric tanh( z ) form of reionization history (e.g., Heinrich et al. 2017; Ghara  \neach 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.  \nFigure 3 shows the ratio σ DM / DM EoR that qualitatively follows DM EoR for z ≳ 7. It first increases rapidly and then saturates towards large z as the ionized regions disappear. However, we find that σ DM / DM EoR increases sharply towards the end of reionization ( z < 7). This is because DM EoR has a value around zero at z ≈ 6 since we ignore the contribution from the low-redshift IGM. Starting at z ≈ 7, DM EoR 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 / DM EoR 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.  \nAs shown in Figure 2, there is a slight offset between DM EoR , the mean estimated over all the 600 2 grids on the transverse plane of the box and the average DM estimated using the ¯ x H ii into eq. 3 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 DM EoR . Fluctuations due to the binary ionization field are not able to cause any significant ( > 1 σ ) deviations between DM EoR ( z ) and DM(¯ x H ii ( z )) for all three histories considered here. The deviation should be enhanced for more realistic reionization LC where both x H ii and ∆ contribute to the LoS fluctuations in DM EoR . Also, contribution to σ DM accumulated from z < 6 will increase the spread in DM EoR .  \nFigure 3 also shows the derivative d DM EoR /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  \nFigure 4 shows contours of the derivatives of the structure function ∂ 2 Ξ /∂z ∂δθ on the δθ -z plane. We consider only 400 LoS at every equidistant comoving slice to compute ∂ 2 Ξ /∂z ∂δθ instead of all available 360 , 000 LoS per comoving slice. Our choice of 400 1 FRBs per comoving slice is closer to real observations and com-  \n<!-- image -->  \nz  \nFigure 3. Top: The ratio σ DM / DM EoR as a function of redshift z . The three different line styles corresponds to the three reionization histories. Bottom: The redshift derivative of DM EoR ( z ). The green lines show the fluctuating estimates whereas the orange lines are after a Gaussian smoothing.  \nderivatives 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 DM EoR at higher z .  \nFigure 4. Density distribution of the double derivative of the structure function ∂ 2 Ξ /∂z ∂δθ . The three different panels correspond to the different reionization histories considered. The different colored filled-contours here represent various percentages of the area encompassed by the ∂ 2 Ξ /∂z ∂δθ surface as shown in the labels. The horizontal dashed lines mark z mid . The dashed lines represent the angular size of bubbles varying with z .  \n<!-- image -->  \nFigure 5. The distribution shown in Figure 4 marginalized along δθ . The three solid lines represent ∫ ( ∂ 2 Ξ /∂z ∂δθ ) dδθ for the three reionization histories. The dashed vertical lines correspond to the respective z mid values.  \n<!-- image -->  \nputationally tractable. On the other hand, we believe that it is a good number for statistical convergence of ∂ 2 Ξ /∂z ∂δθ . We finally bin-average our estimates of all slices within equispaced z bins of width 0 . 5 for further use. We consider ∂ 2 Ξ /∂z ∂δθ as a 2D surface and plot the contours which include top 10%, 30%, 50% and 67% of the total area under the ∂ 2 Ξ /∂z ∂δθ 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  \nthe contours gets more squeezed along the z axis. This clearly indicates that the peakedness of the ∂ 2 Ξ /∂z ∂δθ 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.  \nWe next marginalize ∂ 2 Ξ /∂z ∂δθ over δθ . The marginalized result, ∫ ( ∂ 2 Ξ /∂z ∂δθ ) dδθ , is shown in Figure 5 as a function of z . ∫ ( ∂ 2 Ξ /∂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.", '3.2. 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 §3.1. The three simulations perfectly agree with their common input reionization history, leading to the similar DM EoR values. However, variations in bubble size influence the fluctuations in the DM EoR estimates, and consequently σ DM .  \nFigure 7 shows the contours of ∂ 2 Ξ /∂z ∂δθ . Similar to Figure 4, 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  \nFigure 6. Left: σ DM / DM EoR as a function of z . The three different linestyles correspond to the different bubble sizes chosen here with radii R = 5 , 10 and 20 (in grid units). Right: Derivative of the ratio shown in left panel w.r.t. z . The three different green lines corresponds to the fluctuating estimates from the toy model simulations, whereas the orange lines are after a Gaussian smoothing. The vertical gray line shows the redshift corresponding to z mid .  \n<!-- image -->  \nFigure 7. A similar plot of ∂ 2 Ξ /∂z ∂δθ as in Figure 4. Three different panels correspond to the models with different bubble sizes but the same reionization history. The horizontal lines mark z mid , and the dashed line represents how the bubble size varies with redshift in terms of the angular size, θ R .  \n<!-- image -->  \nR values. The peak of the ∂ 2 Ξ /∂z ∂δθ 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.  \nWe 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 ap-  \nproximately 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 1 / 2 col . 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 6.  \nWe plot σ DM / DM EoR and its z derivative as a function of z in Figure 6. The trend is qualitatively similar  \nFigure 8. The distribution shown in Figure 7, but marginalized along z . Vertical lines denote, 〈 θ R 〉 which is θ R ( z ) marginalized over the z range of interest.  \n<!-- image -->  \nto that shown in Figure 3. We find a sharp turnover and rapid rise in σ DM / DM EoR by the end of reionization ( z < 7) due to near zero value of DM EoR . However for z > 7, σ DM / DM EoR 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 / DM EoR 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.  \nWe again marginalize ∂ 2 Ξ /∂z ∂δθ here, but this time along the z -axis, to obtain ∫ ( ∂ 2 Ξ /∂z ∂δθ ) dz as shown in Figure 8 for the three R values. We observe a similar qualitative trend as seen in the peaks of the contour plots. The peak of ∫ ( ∂ 2 Ξ /∂z ∂δθ ) dz shifts to a larger δθ for large bubble sizes. We also note that it decreases for δθ larger than θ R . 〈 θ R 〉 is computed by marginalizing  \nThe derivative d ( σ DM / DM EoR ) /dz is a more local quantity as shown in Figure 6. The derivatives drop very sharply approaching -∞ for z < 7 as there is a rapid increase in σ DM / DM EoR with decreasing z . However, we see a power-law decrement in the derivative with increasing z for z ≥ 8. The slope of this powerlaw drop is roughly the same for the three R values, although the magnitude scales with R in a similar way as for σ DM / DM EoR .  \nFigure 9. Global averaged ionization fraction of the IGM during reionization. The three different lines represent different reionization histories obtained from realistic grizzly simulations of EoR.  \n<!-- image -->  \nthe θ R ( z ) corresponding to the respective characteristic bubble size R in the z range of interest.', '4. 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 x H ii LCs corresponding to the three different reionization histories ¯ x H ii as shown in Figure 9. In the following sections, we briefly describe the method used to simulate the EoR scenarios followed by all the corresponding estimates.', '4.1. EoR Simulations': "The EoR x H ii LCs used here are constructed by stitching thin slices from the coeval x H ii cubes simulated at several different redshifts in chronological order. We simulate x H ii coeval boxes using a CD-EoR code, grizzly (Ghara et al. 2015a, 2018), 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 x H ii map at a redshift for a particular astrophysical source model.  \nWe use the dark-matter density fields and the associated halo catalogues obtained from the PRACE 2 project, PRACE4LOFAR. The input dark-matter distributions are generated using dark-matter only N -  \nz  \nFigure 10. Light-cone slices of the ionization field. The three different panels correspond to the LC slices obtained from the realistic grizzly simulations for Slower , Fiducial and Faster reionization histories as labelled.  \n<!-- image -->  \nbody simulation code cubep 3 m 3 (Harnois-D'eraps et al. 2013). The 3D density cubes have comoving volume [500 h -1 Mpc] 3 (see e.g, Dixon et al. 2015; Giri et al. 2019) which are gridded into [600] 3 voxels. The darkmatter particle distributions have been used to find the collapsed halos using spherical overdensity halo finder (Watson et al. 2013). The minimum halo mass in the PRACE4LOFAR simulation is ∼ 10 9 M /circledot , 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.  \nWe consider an EoR source model where the darkmatter halos with masses larger than 10 9 M /circledot host UV photon-emitting galaxies. We assume that the stellar  \nmass of a galaxy ( M /star ) is related to the host dark matter halo mass M halo as M /star ∝ 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.  \nNext, we produce coeval cubes of x H ii at 63 redshifts between 6 . 1 and 15 . 6 for the aforementioned source model. We refer the readers to Ghara et al. (2015a) and Islam et al. (2019) for more details about these calculations. Finally, we used these coeval cubes of x H ii to create the LC which accounts for the evolution of x H ii with redshift. The detailed method to implement the LC effect can be found in Ghara et al. (2015b). The reion-  \nization histories for the three different EoR scenarios are shown in Figure 9 while we present the corresponding LCs in Figure 10. The LCs clearly show the difference in the patchiness of the reionization scenarios.", '4.2. grizzly Simulation Results': "We present the results for the LCs (Figure 10) 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 15). This optimistic assumption primarily allows us to investigate the impact of H ii bubble sizes and reionization histories on our estimators by evading lowz contributions.  \nDM EoR ( z ) begins with a near zero value (with lowz IGM contributions subtracted) around the end of reionization ( z ≈ 6) and increases almost linearly with z until when ionization is sufficiently small (¯ x H ii ≈ 0 . 1), where it starts to plateau. This is qualitatively consistent across both the realistic and toy model reionization histories (see Figure 2). 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 DM EoR ( 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.  \nIn Figure 11, we show the DM EoR 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 DM EoR estimates calculated using the average ionization fraction (Figure 9) for the three reionization scenarios.  \nThe DM EoR ( z ) and DM(¯ x H ii ( 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 10) and quickly percolate  \nFigure 11. Sky-averaged DM and the corresponding cosmic variance during EoR. The solid lines represent the DM EoR estimates computed from realistic grizzly LC simulations. Shaded regions around them are the corresponding 1 σ cosmic variance. The cyan, green and orange colors respectively represent Slower , Fiducial and Faster reionization histories. The three dashed lines represent the DM EoR computed using the globally-averaged ionization fractions ¯ x H ii shown in Figure 9.  \n<!-- image -->  \nand 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 10), thereby resulting in a larger difference ( ≳ 2 σ ) between the DM EoR ( z ) and DM(¯ x H ii ( z )) (solid and dashed lines) for the slower reionization history.  \nAs the average DM EoR is an integrated quantity, we compute its derivative d DM EoR /dz to extract instantaneous information at a given stage of reionization. We show the derivatives in Figure 12 for the three histories. d DM EoR /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 H ii changes rapidly during the mid and end stages of the reionization. However, the difference in how fast ¯ x H ii ( 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.  \nFigure 12. Derivative of DM EoR ( z ) w.r.t. redshift z . The three lines show the derivatives obtained from realistic simulations shown in Figure 10.  \n<!-- image -->  \nFigure 13 depicts σ DM / DM EoR 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 / DM EoR depends on the z mid as seen in toy models (top panel of Figure 3). However, contrary to the toy models (Figures 3 and 6), the ratio σ DM / DM EoR decreases with z . The 1 σ cosmic variance for all the three histories are roughly similar, hence the ratio scales inversely with DM EoR . 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.  \nFigure 14 shows contours of the derivatives of the structure function ∂ 2 Ξ /∂z ∂δθ on the δθ -z plane, similar to Figures 4 and 7. We estimate ∂ 2 Ξ /∂z ∂δθ 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  \nFigure 13. The ratio σ DM / DM EoR as a function of z . The ratios are computed using the estimates from the three realistic grizzly simulations.  \n<!-- image -->  \nare smaller in sizes. This trend closely resembles what is observed in Figure 7.", '4.3. 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 nonuniform 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 Beniamini et al. (2021), we used the cumulative observable redshift distribution as shown in Figure 15. 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 tot FRB 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 15) for a given N tot FRB . 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.  \nFigure 14. Density distribution of ∂ 2 Ξ /∂z ∂δθ . This plot is similar to Figure 4, but for the realistic EoR simulations. The three different panels correspond to the three different reionization histories as labelled. However, this is obtained from a realistic simulation of EoR LCs which does not have a fixed bubble size.  \n<!-- image -->  \nFigure 15. Cumulative redshift distribution of observable FRBs. The orange solid line shows the cumulative distribution of the FRB abundance computed following the model in Beniamini et al. (2021). The green shaded region here demarcates the redshift range (6 ≤ z ≤ 15) used in this work.  \n<!-- image -->  \nOur main findings in Figure 16 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 skyaveraged DM varies as a function of N tot FRB . We compute the DM from a random sample of N tot FRB distributed across the EoR redshift range. Later, we consider 1000 such independent realizations of the distribution of N tot FRB to compute the ensemble mean 〈 DM 〉 and the 3 σ cosmic variance as shown by the shaded regions around the lines. Figure 16 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 ob-  \nomewhere within the duration of EoR, we can compute the marginalized sky-averaged DM defined as  \nDM margin = ∫ z max z min DM EoR ( z ) dP/dz dz ∫ z max z min dP/dz dz , (5)  \nwhere dP/dz is the probability of detecting an FRB per unit redshift (as per Figure 15), 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 16.  \nWe 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  \nWe 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 lowz 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.  \nFigure 16. Variation of 〈 DM( z mid ) 〉 and the marginalized 〈 DM margin 〉 as a function of total number of FRBs N tot FRB observed during the EoR. Top row: This shows the ensemble averaged DM( z mid ) and its standard deviation as a function of N tot FRB for the three reionization histories. Bottom row: This shows ensemble mean of DM margin , marginalized DM over the EoR redshift range, as a function of N tot FRB for the three different reionization histories. Solid lines here are the ensemble mean and shaded regions correspond to 3 σ cosmic variances which are estimated, for every choice of N tot FRB , using 1000 independent ensembles of FRB distributions in the IGM during EoR. The left, middle and right columns correspond to the optimistic, moderate and pessimistic scenarios, respectively, as described in the text.  \n<!-- image -->  \ngadget4 4 (Springel et al. 2021) 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 §1. 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.  \nThe contribution from low-redshifts increases both DM EoR and σ DM ( z ), making the separation between the two different histories more challenging. In Figure 16, 〈 DM( z mid ) 〉 for the three different histories starts at different N tot FRB . This is because, for the slower history, z mid is larger than that of the other histories, requiring a higher N tot FRB to populate the corresponding z -bin with at least two FRBs. 〈 DM( z mid ) 〉 converges very quickly for N tot FRB ≈ 10; however the corresponding dispersion is large and it decreases with increasing N tot FRB , as expected. N tot FRB ≥ 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 tot FRB = 150 and 40 to discern 〈 DM( z mid ) 〉 of slower-fiducial and fiducialfaster pairs of histories, respectively. Considering the pessimistic case, it takes N tot FRB = 80 , 150 , 600 to distin-  \nlower, slower-fiducial and fiducialfaster pairs of histories, respectively.  \nFigure 16 also shows the variation of the 〈 DM margin 〉 with N tot FRB . 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 tot FRB 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 tot FRB ≈ 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.', '5. 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 H i in the IGM during EoR and, Thomson scattering optical depth of the CMB photons, highz quasar spectra, Lyα systems at highz . 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 highz to probe the IGM during EoR.  \nThe 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 skyaveraged 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 1) of the binary ionization field (within the range 6 ≤ z ≤ 15), we see that the DM EoR ( 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 2). The derivative d DM EoR /dz directly traces how fast the reionization progresses (Figure 3). The all-sky variance σ DM / DM EoR  \nis larger for the reionization scenario where the bubble sizes are larger and vice-versa (Figure 6). We also compute ∂ 2 Ξ /∂z ∂δθ where Ξ is the angular dispersion (eq. 4) at any redshift. We demonstrated that the contours are elongated or squeezed along z and δθ depending respectively on the varying reionization history (Figure 4) and bubble sizes (Figure 7). The impacts are clearly prominent for the marginalized structure function derivatives (Figures 5 and 8).  \nWe next consider observationally more realistic situation where the FRBs abundances vary with redshift (Figure 15) and they are more clustered at the highly dense regions on the sky plane. This biases DM EoR 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 16) 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  \nWe also analyzed a more realistic reionization lightcone (LC) (Figure 10) generated using grizzly 1D radiative transfer code for three different reionization histories ending at the same z . The behaviour of DM EoR is the same as in the toy model. We find that the different LoS variance plays a significant role for realistic IGM (Figure 11) and biases the average DM EoR ( z ) by more than 1 σ as compared to the case when the DM is computed using the averaged ionization fraction ¯ x H ii . The slope of the d DM EoR /dz is directly sensitive to the reionization history (Figure 12). σ DM / DM EoR is also dependent on the reionization history as well as indirectly on the IGM morphology (bubble sizes) in an intermingled way. Figure 14 clearly shows that ∂ 2 Ξ /∂z ∂δθ 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 viceversa. 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.  \nbetween the reionization histories at 3 σ might require roughly ∼ 1000 FRBs in total during EoR (see right column of Figure 16).  \nThe 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 ∂ 2 Ξ /∂z ∂δθ and of σ DM as probes of the ionization bubble sizes along with the reionization history. The ∂ 2 Ξ /∂z ∂δθ , being a derivative, suffers from large variance, and we need a large number of FRBs ( N tot FRB ≲ 100 000) within range 6 ≤ z ≤ 15 to suppress the variance significantly down. The computation of ∂ 2 Ξ /∂z ∂δθ 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  \nwhich 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.", 'REFERENCES': 'Abdurashidova, Z., Aguirre, J. E., Alexander, P., et al. 2022, ApJ, 925, 221, doi: 10.3847/1538-4357/ac1c78 Ali, S. S., Bharadwaj, S., & Chengalur, J. N. 2008, MNRAS, 385, 2166, doi: 10.1111/j.1365-2966.2008.12984.x Asthana, S., Haehnelt, M. G., Kulkarni, G., et al. 2024, arXiv e-prints, arXiv:2404.06548, doi: 10.48550/arXiv.2404.06548 Ba˜nados, E., Venemans, B. P., Mazzucchelli, C., et al. 2018, Nature, 553, 473, doi: 10.1038/nature25180 Barry, N., Wilensky, M., Trott, C. M., et al. 2019, ApJ, 884, 1, doi: 10.3847/1538-4357/ab40a8 Becker, G. D., Bolton, J. S., Madau, P., et al. 2015, MNRAS, 447, 3402, doi: 10.1093/mnras/stu2646 Becker, R. H., Fan, X., White, R. L., et al. 2001, AJ, 122, 2850. http://stacks.iop.org/1538-3881/122/i=6/a=2850 Beniamini, P., & Kumar, P. 2023, MNRAS, 519, 5345, doi: 10.1093/mnras/stad028 Beniamini, P., Kumar, P., Ma, X., & Quataert, E. 2021, MNRAS, 502, 5134, doi: 10.1093/mnras/stab309 Bera, A., James, C. W., Deller, A. T., et al. 2024, arXiv e-prints, arXiv:2406.13704, doi: 10.48550/arXiv.2406.13704  \nBernardi, G., de Bruyn, A. G., Brentjens, M. A., et al. 2009, A&A, 500, 965, doi: 10.1051/0004-6361/200911627'}