Measuring Hubble constant using localized and unlocalized fast radio bursts

The former two components, ${\rm DM_{ISM}}\$and ${\rm DM_{halo}}\$, are contributed by medium within Milky Way and also referred to as a whole, i.e. $\rm DM_{MW}=DM_{ISM}+DM_{halo}$. ${\rm DM_{ISM}}\$could be well described by galactic electron distribution models such as YMW16 by Yao et al. (2017) and NE2001 by Cordes & Lazio (2002). According to Ocker et al. (2021), YMW16 might overestimate ${\rm DM_{ISM}}\$for FRBs in the direction of anticenter and possibly other low-latitude bursts. The difference between NE2001 and YMW16 may have little influence on constraint of $H_{0}$ (Wu et al., 2022), still we apply NE2001 to estimate ${\rm DM_{ISM}}\$as previous researchers do. Prochaska & Zheng (2019) gave the constraint of ${\rm DM_{halo}}\$= 50$\sim$80 ${\rm pc\ cm^{-3}}$based on multiple observation data. It is reasonable to assume that ${\rm DM_{halo}}\$follows a Gaussian distribution with $\langle\rm DM_{halo}\rangle$ = 65 ${\rm pc\ cm^{-3}}$and $\sigma=15$ ${\rm pc\ cm^{-3}}$. To reduce calculation and avoid another level of integration, we simplify ${\rm DM_{halo}}\$to be its mean value, i.e. $\rm DM_{halo}$ = 65 ${\rm pc\ cm^{-3}}$, which should have minor influence on our estimation since the halo component appears as linear form in $\rm DM_{exc}=DM_{obs}-DM_{ISM}-DM_{halo}$.

As discussed in Section 2.1, contributions from within our Galaxy can be well estimated with galactic electron models, we subtract $\rm DM_{MW}$ from total DM to obtain the extragalactic component:

${\rm DM_{host}}\$is the dispersion measure contributed by host galaxy of FRB source, and its probability distribution can be fitted by a lognormal distribution according to Macquart et al. (2020b). The intergalactic component ${\rm DM_{IGM}}\$could be described with a Gaussian-like distribution around its mean value. In a standard $\rm\Lambda CDM$ universe model, the mean value of ${\rm DM_{IGM}}\$is given by Deng & Zhang (2014) as:

where $m_{p}$ is proton mass, $f_{\rm IGM}$ is the fraction of baryon in IGM. Previous researches prefer a value of $f_{\rm IGM}\simeq 0.84$ according to Shull et al. (2012), yet Connor et al. (2024) gives a more accurate constraint of $f_{\rm IGM}\simeq 0.93$ with data from the Deep Synoptic Array (DSA-110) using both FRB and non-FRB methods (notably what we called $f_{\rm IGM}$ corresponds to $f_{d}$ instead of $f_{\rm IGM}$ in their research). We adopt the value of $f_{\rm IGM}=0.93$ in our simulation. The integral $f_{e}(z)$ is defined as:

The cosmological parameters $\Omega_{m}$ and $\Omega_{\Lambda}$ are given by Planck 2018 results (Planck Collaboration et al., 2020), and we have $\Omega_{\Lambda}=1-\Omega_{m}$ assuming a flat universe. $\Omega_{b}$ is decided as an assumption in Planck results based on big bang nucleosynthesis (BBN) constraints and primordial deuterium abundance measurements by Cooke et al. (2018), which is always given in the form of $\Omega_{b}h^{2}$ where $h=H_{0}/100$ ${\rm km\ s^{-1}Mpc^{-1}}$, thus we modify the equation to keep $\Omega_{b}$ in the form of $\Omega_{b}{H_{0}}^{2}$. $y_{1}$ and $y_{2}$ in Equation (4) are hydrogen and helium fractions normalized to 0.75 and 0.25 respectively, which can be neglected as $y_{1}\simeq y_{2}\simeq 1$. $\chi_{e,{\rm H}}(z)$ and $\chi_{e,{\rm He}}(z)$ are ionization fraction of hydrogen and helium, which could also be considered to be $\chi_{e,{\rm H}}(z)=\chi_{e,{\rm He}}(z)=1$ at $z<3$. Equations (3) and (4) can now be further rewritten as:

To run a Markov Chain Monte Carlo (MCMC) simulation, the probability distribution of extragalactic DM components must be calculated. As described in Section 2.2, ${\rm DM_{host}}\$follows a lognormal distribution while ${\rm DM_{IGM}}\$can be fitted with a Gaussian-like distribution, which often writes as (McQuinn, 2013; Macquart et al., 2020b):

where $\alpha$ and $\beta$ are parameters indicating inner density profile of halo gas, and Macquart et al. (2020b) proposed the best-fitting result of $\alpha=\beta=3$. $A$ and $C_{0}$ are normalization parameters given by $\int p_{\rm IGM}=1$ and $\langle\Delta\rangle=1$, while $\sigma_{\rm host}\$, $\mu$ and $\sigma_{\rm IGM}\$are distribution parameters respectively. $e^{\mu}\$is generally chosen as distribution parameter instead of $\mu$, since $e^{\mu}\$indicates mean value of ${\rm DM_{host}}\$.

One possible way of further simulation to determine $H_{0}$ is to fit all undetermined parameters at the same time, and the probability function goes like:

Still, since the number of localized FRBs is limited ($n_{\rm local}<50$), the confidence may be weakened to fit a four-dimensional parameter space $(e^{\mu},\sigma_{\rm host},\sigma_{\rm IGM},H_{0})$ with current data. Another concern is that distribution parameters are not fixed constant for different FRBs, and some have shown redshift-dependent. Zhang et al. (2020) and Zhang et al. (2021) proposed the best-fitted distribution parameters of ${\rm DM_{host}}\$and ${\rm DM_{IGM}}\$derived from IllustrisTNG Simulation (Pillepich et al., 2017), which takes redshift-dependent evolution into consideration as well. Zhang et al. (2021) also provided the evolution of $A$ and $C_{0}$, and it is shown that the best-fitted value slightly deviates from normalization. To apply results above, all localized FRBs are divided roughly into three categories based on host galaxy: (a) non-repeating (one-off) bursts; (b) FRB121102-like repeating bursts: host galaxy stellar mass $\simeq 4\sim 7\times 10^{7}\ M_{\odot}$, star formation rate (SFR) $\simeq 0.4$; (c) FRB180916-like repeating bursts: host galaxy stellar mass $\simeq 10^{10}\ M_{\odot}$, SFR $\simeq 0.016$. Note that this classification only indicates qualitative properties of the host galaxy and is not absolute. Zhang et al. (2020) and Zhang et al. (2021) provided best-fitted values of distribution parameters at several typical redshifts, and we perform a monotone cubic spline interpolation for each parameter to obtain values at any given redshift.

By obtaining values for $(e^{\mu},\sigma_{\rm host},\sigma_{\rm IGM},A,C_{0})$ from IllustrisTNG, the only free parameter left is $H_{0}$, which is what we are interested in. Taking all parameters into Equation (6) and  (7), we have the likelihood function for any single FRB:

Though not wrote explicitly, the distribution parameters in Equation (9) are distinct for each burst. For the $i^{\rm th}$ FRB, the complete likelihood function is:

and the total log-likelihood function of all FRBs is:

According to Bayesian theory, we still need a prior of $H_{0}$ to perform parameter estimation, and we use a uniform distribution $H_{0}\in\mathcal{U}(0,100)$ ${\rm km\ s^{-1}Mpc^{-1}}$as its prior.

With all data of pseudo redshifts for unlocalized FRBs prepared, difference still exists between simulation of localized and unlocalized FRBs. Generally there are two methods of simulating: integrating with entire probability density, and using randomly scattered pseudo redshifts, both of which are necessary to give a comprehensive understanding of constraint on Hubble constant.

Statistical error mainly refers to the error of cosmological parameters $\Omega_{b}h^{2}$ and $\Omega_{m}$ appeared in Equation (5). The error of $f_{\rm IGM}$ will be discussed in Section 5.2, and other constant in Equation (5) such as $G$ and $m_{p}$ are already measured with extremely high precision. Planck Collaboration et al. (2020) gave $\Omega_{b}h^{2}=0.02242\pm 0.00014$ and $\Omega_{m}=0.3111\pm 0.0056$.

For $\Omega_{b}h^{2}$, it appears in Equation (5) as a linear term. Assuming a Gaussian distribution $p(x)$ around $\mu=0.02242$, we can do a rough estimation:

which shows a symmetric distribution of $\Omega_{b}h^{2}$ has little influence on $H_{0}$. Furthermore, the precision of $\Omega_{b}h^{2}$ ($\sim 0.6\%$) is slightly smaller than that of $H_{0}$ of our constraint ($\sim 0.94\%$).

For $\Omega_{m}$, which appears inside the integral in Equation (5), it has an uncertainty of $\sim 1.8\%$ and could not be ignored. To marginalize $\Omega_{m}$, the best way is to add another level of integral, which would significantly increase computation. Alternatively, we consider replacing integral with expansion. Assume a Gaussian distribution $p(x)$ with $\mu=0.3111$, $\sigma=0.0056$, the integral can be wrote as:

where $I_{0}=\int_{\mu-\sigma}^{\mu+\sigma}p(x)dx=0.683$ is the normalization factor and $f(x,z)=p(x)/\left[x(1+z)^{3}+1-x\right]^{1/2}$ is our target function. Note that $Z$ is redshift of FRB, and $z$ is our integration variable. Since both $\int f(x,z)dx$ and $\int f(x,z)dz$ could not be expressed by elementary functions, we perform a series expansion on $f(x,z)dz$ around $x=\mu$ and get $g(x,z)=\sum_{i=0}^{5}a_{i}(z)(x-\mu)^{i}$. We plot the figure of $g(x,z)$ at different $z$ and compare it with original function to ensure the expansion is acceptable. Now the inner integral in Equation (17) can be wrote as an explicit function of $z$, i.e. $h(z)=\int_{\mu-\sigma}^{\mu+\sigma}g(x,z)dx$. However, the form of $h(z)$ is still too complicated to integrate, thus we need to do another series expansion around $z=z_{0}$. To determine the best value for $z_{0}$, we plot the expanded function at different $z_{0}$ and degrees. It is found that expanding $h(z)$ to the term of $(z-z_{0})^{4}$ around $z_{0}=2.5$ provides best fit for both $z\rightarrow 0$ and $z\rightarrow 4$. Denote it as $j(z)=\sum_{i=0}^{4}b_{i}(z)(z-z_{0})^{i}$, and we can complete the whole integral: $I\simeq 1.02\ Z+0.19\ Z^{2}-0.14\ Z^{3}+0.043\ Z^{4}-0.0066\ Z^{5}+0.00042\ Z^{6}$.

Taking the new expression into Equation (5), we could run MCMC simulation with $\Omega_{m}$ marginalized. The result is shown in Fig. 6 (we use the method in Section 4.3.2). The $1\sigma$ confidence interval of $H_{0}$ is $69.00_{-0.69}^{+0.68}$ ${\rm km\ s^{-1}Mpc^{-1}}$and has little difference with previous result in Fig. 5, which indicates that the statistical error of $\Omega_{m}$ does not influence constraint on $H_{0}$ significantly.

There are several systematic errors that need to be discussed. In Equation (5), compared with $\Omega_{b}h^{2}$ and $\Omega_{m}$, the fraction of baryon in IGM $f_{\rm IGM}$ may introduce more uncertainty. However, we still know little about $f_{\rm IGM}$. Shull et al. (2012) gave the value of $\sim 0.84$ yet does not propose the errorbar. Connor et al. (2024) provided a more accurate constraint, yet it depends on several other models and theories. Furthermore, it is difficult for any constraint using Equation (5) to separate the error of $f_{\rm IGM}$ from $H_{0}$ as they appear in a coupling term $f_{\rm IGM}/H_{0}$ in Equation (5). To be precise, only constraint on $f_{\rm IGM}/H_{0}$ can be made instead of constraint on $H_{0}$. So, the value of $H_{0}$ must be determined by other approaches when constraining $f_{\rm IGM}$ from FRBs. By fixing $H_{0}$, it has been found that the uncertainty of $f_{\rm IGM}$ is about 8% (Yang et al., 2022; Connor et al., 2024). Obviously, the systematic error from $f_{\rm IGM}$ dominates the error of measured $H_{0}$ at present. On the other hand, it may varies with redshift. Researches with other methods are needed to further investigate the error of $f_{\rm IGM}$.

${\rm DM_{host}}\$includes all contributions to DM that come from FRB host galaxy. Using IllustrisTNG simulations, the probabilty distribution of ${\rm DM_{host}}\$including the local cosmic structure (e.g. filament) halo and interstellar medium of the host have been derived (Zhang et al., 2020). However, the vicinity environment of FRB progenitors have not been considered. The most promising progenitors of FRBs are young magnetars, which can be formed by core-collapse of massive stars or mergers of two compact objects (Wang et al., 2020). So they may be embedded in a magnetar wind nebula and supernova remnant (Yang & Zhang, 2017; Piro & Gaensler, 2018; Zhao & Wang, 2021). Meanwhile, the largest ${\rm DM_{host}}\$with rotation measure reversal for FRB 20190520B indicates it may reside in a binary system (Wang et al., 2022; Anna-Thomas et al., 2023). So similar FRBs should be removed when measuring $H_{0}$. For FRBs with little DM contribution from vicinity environment, precise modelling should be performed. It is important to use optical observations of the FRB host galaxy environment, combined with the rotation measure and scattering times of FRBs to constrain ${\rm DM_{host}}\$(Cordes et al., 2022).

Several selection biases should be discussed, such as SNR effect, ISM effect, selection effect of unlocalized FRBs and gridding effect (James et al., 2022). For the latter two effects, our constraint has no such error since we make use of most unlocalized FRBs in CHIME database and use continuous value for redshift and dispersion measures instead of discrete variables. For SNR effect, the log-log figure of number of events observed above SNR threshold should follow power law of -1.5 (in log-log plot, i.e. $N\propto{\rm SNR}^{-1.5}$), and James et al. (2022) found that events from CRAFT/ICS deviates from the -1.5 power law. We plot the same figure with unlocalized FRBs from CHIME database in Fig. 7, and the histogram followed the power law well, thus our constraint is not much influenced by SNR effect. Finally, for the ISM effect, James et al. (2022) claimed that ${\rm DM_{ISM}}\$would increase at low galactic latitude, which may prevent telescopes from observing such events. We plot Hammer projection of FRBs in galactic coordinate system in Fig. 8. Hammer projection is an equal-area projection and there are a considerable amount of FRBs located in low galactic latitude area. Furthermore, few of these low-galactic-latitude FRBs have shown extreme values of ${\rm DM_{ISM}}\$even above 200${\rm pc\ cm^{-3}}$. Thus ISM effect could be ignored during our constraint using unlocalized FRBs.