A Data-driven Technique Using Millisecond Transients to Measure the Milky Way Halo

In the early universe the majority of baryons resided in a cool, diffuse plasma, which is predicted to have collapsed into sheetlike and filamentary structures that make up the intergalactic medium (IGM). Around the time of structure formation, dark matter collapses into halos, pulling baryons with it. As the gas falls inward, it is shock-heated to form a hot, diffuse plasma, known as halo gas or the circumgalactic medium (CGM). Approximately 10% of the gas cools and falls into the center of the halo to form stars and the interstellar medium (ISM; e.g., White & Rees 1978).

Comparing the baryonic mass fraction detected for galaxies (${M}_{b}/{M}_{\mathrm{halo}}$) to the cosmic mean (${{\rm{\Omega }}}_{b}/{{\rm{\Omega }}}_{m}$), however, reveals a baryonic deficit (e.g., Dai et al. 2010). The missing baryons may have been ejected back into the IGM before forming stars or perhaps have yet to be detected (e.g., Prochaska et al. 2011; Booth et al. 2012). In the latter scenario, the CGM presents itself as a possible refuge.

This issue holds for the CGM of our Galaxy. While it is evident that its stars and ISM correspond to ≲25% of the baryonic mass available to a halo with mass ${M}_{\mathrm{halo}}={10}^{12.2}{M}_{\odot }$ (the current estimate; Boylan-Kolchin et al. 2013), the mass and distribution of gas within our Galactic halo are not well determined even despite our close proximity. The key observables that constrain the Galactic CGM include soft X-ray emission from the plasma (Henley et al. 2010), X-ray, and UV absorption lines of oxygen ions (Faerman et al. 2017; Kovács et al. 2019), density constraints from ram pressure stripping of the Large Magellanic Cloud (LMC; Salem et al. 2015), and dispersion measure (DM) observations from pulsars toward the LMC (Manchester et al. 2006). These have provided valuable constraints for models of the Galactic halo, but still allow for large variations in the mass and spatial extent of the gas (Fang et al. 2013; Faerman et al. 2013; Bregman et al. 2018; Prochaska & Zheng 2019).

A primary challenge to assessing the Galactic CGM is that the gas is too diffuse (especially at large radii) to be imaged directly. Furthermore, the absorption-line measurements (e.g., O vi and O vii) require substantial ionization and/or metallicity corrections to infer the total gas. In this respect, the DM measurements toward the LMC provide the most direct probe of the ionized gas, yet it lies at only ≈1/4 the virial radius ${r}_{200}$ of the Galaxy. Ideally, one would prefer to record DM measurements to ${r}_{200}$ and also across the sky to search for asymmetries in the halo gas distribution. Just such an opportunity is now afforded (albeit with caveats, as we will discuss) by the transients known as fast radio bursts (FRBs).

FRBs are the population of ∼millisecond chirps of bright radio emission at approximately GHz frequencies discovered serendipitously (Lorimer et al. 2007) and now pursued in earnest with dedicated projects and facilities (Caleb et al. 2016; CHIME/FRB Collaboration et al. 2018; Law et al. 2018; Kocz et al. 2019). Recorded in each FRB event is its DM value ${\mathrm{DM}}_{\mathrm{FRB}}$. The majority greatly exceed estimates for our Galactic ISM and CGM, lending strong statistical support that FRBs have an extragalactic origin (Cordes & Chatterjee 2019; Petroff et al. 2019). This inference has been confirmed by a small but growing set of FRBs localized to $\approx 1^{\prime\prime}$ and then shown to reside in a distant galaxy (Tendulkar et al. 2017; Bannister et al. 2019; Prochaska et al. 2019; Ravi et al. 2019; Marcote et al. 2020). As a result, the community now recognizes FRBs as a viable tool to probe ionized gas across the universe, e.g., to conclusively detect the so-called "missing" baryons of the present-day universe (Fukugita et al. 1998; Macquart 2018).

Owing to its integral nature, ${\mathrm{DM}}_{\mathrm{FRB}}$ includes contributions from all of the electrons along the sight line: the IGM, gas in distant Galactic halos, the ionized gas of the system hosting the FRB, and our Milky Way (MW). Indeed, the host and Galaxy contributions (${\mathrm{DM}}_{\mathrm{host}}$, DM${}_{\mathrm{MW}}$) are frequently considered a "nuisance" to proposed analyses of the cosmic web. In this Letter, however, we view them as a highly desired signal, i.e., a new opportunity to constrain the Galactic CGM.

There are two primary challenges that this Letter addresses: how to use pulsars and FRBs to probe the DM of Galactic halos, and how to do so with a limited data set. The first problem is addressed by constraining the DM contribution of the MW halo to the total observed DM of pulsars and FRBs.

For the second challenge, only ∼100 FRBs have been observed to date; this necessitates techniques that are well suited to dealing with small data sets. We propose the use of standard kernel density estimation (KDE; Silverman 1986) and asymmetric, variable-bandwidth KDE (Chen 2000; Hoffmann & Jones 2015) to find probability density functions (PDFs) of the DM distribution of pulsars and of FRBs, respectively. Other density estimation techniques are explored—namely, density estimation using field theory (DEFT; Kinney 2014, 2015; Chen et al. 2018) and a generalized extreme value (GEV), but prove to be insufficient (see Appendices B and C for details). From the PDFs one can estimate the maximum MW halo DM given by pulsars, and the minimum MW halo and host halo DM given by FRBs. This infers constraints on the DM of the MW CGM and part of the host CGM.

We measure a MW halo DM of ${63}_{-21}^{+27}\ (\mathrm{stat})\pm 9\ (\mathrm{sys})\mathrm{pc}\,{\mathrm{cm}}^{-3}$, corresponding to a 1σ confidence detection. The precision of this measurement is limited by the FRB sample size and we predict a robust detection of the MW halo with the incorporation of FRB detections anticipated in the coming year. The techniques presented here will make the best precision and least ambiguous measurement of the MW halo in several years with samples of 10⁴ FRBs.

The Letter is structured as follows. Section 2 outlines the core concepts of this work. Section 3 details the density estimation techniques used in the analysis. The methodology and results are presented in Section 4, where Section 4.2 provides constraints based on observed data and Section 4.3 provides an analysis based on simulations. The results and implications are discussed in Section 5, and conclusions are summarized in Section 6.

Pulsars and FRBs are both millisecond radio transients. The former lie in the disk of the MW galaxy and the latter are extragalactic. Since the group velocity of the electromagnetic wave depends on the free electron density (n_e) along the path of propagation, the arrival time of the transient signal is extended. This spread is described by the DM:

$$\begin{eqnarray}&&\mathrm{DM}=\int \displaystyle \frac{{n}_{e}{ds}}{1+z}.\end{eqnarray} \tag{ 1 }$$

DMs can therefore be used to study the distribution of baryons along the line of sight between a transient source and an observer.

Figure 1 shows a schematic of how electrons are distributed relative to pulsars and FRBs. Galactic halos are assumed to be devoid of radio transients, but contain a significant column density of electrons. Pulsars have been detected predominantly in the Galactic disk or nearby globular clusters⁵ (Manchester et al. 2005). Those with a known distance have been used to create detailed models of the electron density distribution of the Milky Way disk (Cordes & Lazio 2002, 2003; Gaensler et al. 2008; Yao et al. 2017). In the following we adopt both the NE2001⁶ and YMW16⁷ algorithms.

Zoom In Zoom Out Reset image size

If we assume FRBs are distributed throughout their host galaxies and throughout space, then the lowest ${\mathrm{DM}}_{\mathrm{FRB}}$ values set a bound on the electron column density associated with the halos of the Milky Way and the typical host galaxy. This measurement is the focus of this Letter. Table 1 provides a summary of the notation used in this Letter.

Table 1.  Notation

Quantity	Description
${\mathrm{DM}}_{\mathrm{pulsar}}$	The total DM measurement of a pulsar
${\mathrm{DM}}_{\mathrm{FRB}}$	The total DM measurement of an FRB
${\mathrm{DM}}_{\mathrm{ISM}}^{\delta }$	DM from a fraction of the Galactic ISM
${\mathrm{DM}}_{\mathrm{ISM}}$	Total sight line DM for the Galactic ISM
${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}$	DM of all gas in our Galactic halo
${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}^{\delta }$	DM from a fraction of gas in our Galactic halo
${\mathrm{DM}}_{\mathrm{IGM}}$	DM from the IGM (gas between halos)
${\mathrm{DM}}_{\mathrm{cosmic}}$	DM from all cosmic gas (IGM + halos)
$\langle {\mathrm{DM}}_{\mathrm{cosmic}}\rangle$	Average DM from all cosmic gas
${\mathrm{DM}}_{\mathrm{host}}$	DM from FRB host galaxy halo

Download table as:  ASCIITypeset image

2.1. Constraints from Pulsars

We consider

$$\begin{eqnarray}&&{\mathrm{DM}}_{\mathrm{pulsar}}={\mathrm{DM}}_{\mathrm{ISM}}^{\delta }+{\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}^{\delta },\end{eqnarray} \tag{ 2 }$$

with ${\mathrm{DM}}_{\mathrm{ISM}}^{\delta }$ as the ISM contribution and ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}^{\delta }$ as the halo contribution. We then define an ISM-corrected quantity ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{pulsar}}$, which subtracts the total ISM contribution along the pulsar sight line,

$$\begin{eqnarray}&&{\rm{\Delta }}{\mathrm{DM}}_{\mathrm{pulsar}}={\mathrm{DM}}_{\mathrm{pulsar}}-{\mathrm{DM}}_{\mathrm{ISM}}.\end{eqnarray} \tag{ 3 }$$

Most pulsars have unknown distances yet are expected to lie predominantly in the Galactic disk, with a scale height of 100 pc (Faucher-Giguère & Kaspi 2006). Therefore, ${\mathrm{DM}}_{\mathrm{ISM}}$ is generally larger than ${\mathrm{DM}}_{\mathrm{pulsar}}$ and the majority of ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{pulsar}}$ values will be negative. Any positive values could be attributed to the halo, and therefore the maximum ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{pulsar}}$ yields a lower limit:

$$\begin{eqnarray}&&{\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}\gt \max \left[{\rm{\Delta }}{\mathrm{DM}}_{\mathrm{pulsar}}\right].\end{eqnarray} \tag{ 4 }$$

Such an analysis must allow for uncertainties in the modeling of ${\mathrm{DM}}_{\mathrm{ISM}}$, but for high Galactic latitudes these uncertainties are expected to be less than $10\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$.

2.2. Constraints from FRBs

${\mathrm{DM}}_{\mathrm{FRB}}$ has contributions from the ISM, the MW halo, cosmic gas, and the FRB host galaxy,

$$\begin{eqnarray}&&{\mathrm{DM}}_{\mathrm{FRB}}={\mathrm{DM}}_{\mathrm{ISM}}+{\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}+{\mathrm{DM}}_{\mathrm{cosmic}}+{\mathrm{DM}}_{\mathrm{host}}.\end{eqnarray} \tag{ 5 }$$

Similar to the pulsars, we define an ISM-corrected measure:

$$\begin{eqnarray}&&{\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}={\mathrm{DM}}_{\mathrm{FRB}}-{\mathrm{DM}}_{\mathrm{ISM}}.\end{eqnarray} \tag{ 6 }$$

From the full distribution of ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}$, we will examine the lowest values on the expectation that these have lower ${\mathrm{DM}}_{\mathrm{cosmic}}$ contributions. For reference, an FRB at z = 0.03 (e.g., Marcote et al. 2020) has an average $\langle {\mathrm{DM}}_{\mathrm{cosmic}}\rangle \approx 25\mathrm{pc}\,{\mathrm{cm}}^{-3}$.

The lowest values of ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}$ should also reflect the lowest combinations of ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}$ and ${\mathrm{DM}}_{\mathrm{host}}$. We expect significant variations in the latter both due to the distribution of host galaxy masses and also from variations in the FRB location within the galaxy. We express ${\mathrm{DM}}_{\mathrm{host}}^{\min }$ as the minimum of this distribution, which may be 10 to several tens of $\mathrm{pc}\,{\mathrm{cm}}^{-3}$.

Regarding variations in ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}$, galaxy formation models tend to predict a nearly spherical distribution of gas, especially beyond the inner halo (but see Yamasaki & Totani 2020, which includes a nonspherical component). Spherically symmetric models of our Galaxy yield less than $10\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$ variations in ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}$ even though the Sun is located off-center (Prochaska & Zheng 2019). In the following, we will assume a single ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}$ unless otherwise discussed. One recovers

$$\begin{eqnarray}&&{\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}+{\mathrm{DM}}_{\mathrm{host}}^{\min }=\min \left[{\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}\right],\end{eqnarray} \tag{ 7 }$$

and therefore

$$\begin{eqnarray}&&{\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}\lt \min \left[{\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}\right].\end{eqnarray} \tag{ 8 }$$

KDE is a nonparametric technique that estimates an unknown density by constructing a kernel at each data point and summing their contributions. Owing to their shapes, the distributions of ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{pulsar}}$ and ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}$ are each suited to a different class of KDE. ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{pulsar}}$ has smooth edges and can be adequately modeled with a Gaussian kernel and a fixed bandwidth. The sharp edge of ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}$, however, necessitates a varying bandwidth and a kernel with a steep cutoff.

In Section 3.1 we outline standard KDE and in Section 3.2 we describe the modifications for asymmetric, bandwidth-varying KDE.

3.1. Standard KDEs

Consider an independent and identically distributed sample $\{{X}_{i}:i=1,...,n\}$ drawn from some unknown distribution f(x). We wish to obtain an estimate $\hat{f}(x)$ of this distribution using KDE:

$$\begin{eqnarray}&&\hat{f}(x)=\displaystyle \frac{1}{n}\sum _{i=1}^{n}{K}_{h}\left({X}_{i}-x\right)=\displaystyle \frac{1}{{nh}}\sum _{i=1}^{n}K\left(\displaystyle \frac{{X}_{i}-x}{h}\right),\end{eqnarray} \tag{ 9 }$$

where K is the kernel and $h\gt 0$ is the bandwidth. The kernel is the underlying distribution function and the bandwidth is a smoothing parameter. In standard KDE symmetric kernels are used, such as Gaussian, triangular, cosine, biweight, triweight, or Epanechnikov. While an Epanechnikov kernel is most optimal in terms of the mean squared error, a Gaussian kernel is the most widely used: the loss of efficiency is marginal ($\sim 5 \%$) and the distribution offers convenient mathematical properties. As such, a Gaussian kernel is used in our analysis of ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{pulsar}}$. Bandwidth selection is a trade-off between the bias of the KDE and its variance. Often the bandwidth is chosen to minimize the mean integrated squared error (MISE),

$$\begin{eqnarray}&&\mathrm{MISE}(h)={\rm{E}}\left[\int {\left(\hat{f}(x)-f(x\right)}^{2}{dx}\right],\end{eqnarray} \tag{ 10 }$$

which is equivalent to the expected L₂ risk function. f(x) is unknown; however, it can be approximated through various techniques (see Jones et al. 1996). One can also use rule-of-thumb bandwidth estimators, such as Silverman's (1986) and Scott's (1979); however, these assume the underlying distribution is Gaussian. In our analysis we use scikit-learn to select the optimal bandwidth via cross-validation.

The KernelDensity() function invokes a nearest neighbors based approach: instead of using the full data set to estimate the density at each point, a number of neighboring points are selected based on the bandwidth. This improves the algorithm efficiency by ignoring distant points that have a negligible effect. KDEs are generated for a range of bandwidths, and GridSearchCV() is used to find the optimal bandwidth. Here n-fold cross-validation is performed. The pulsar data are divided into n subsets, a KDE is generated using the data from $n-1$ subsets (training data), and the performance of the KDE is evaluated on the remaining subset (test data) by calculating the log-likelihood, $\sum \mathrm{log}\hat{p}({x}_{i})$. This process is repeated n times, using a different subset as the test set each time, to give a final (averaged) log-likelihood score. In this manner, scores are calculated for a range of bandwidths. The bandwidth with the maximum log-likelihood is selected for the analysis ($h\approx 10$).

3.2. Asymmetric KDEs

Standard KDE performs well when the underlying distribution is unbounded and the density of data is relatively uniform. We will show, however, that the ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}$ distribution has data concentrated toward the front of the distribution and is bounded on $[0,\infty )$. This presents two problems that standard KDE cannot resolve. First, a fixed bandwidth h entails a trade-off between large- and small-scale structure: overdense regions will be oversmoothed by a large h, and underdense regions will be overfitted if h is too small. Second, symmetric kernels have significant bias at or near a boundary, known as edge or boundary effects. A fixed and symmetric kernel will allocate weight outside of the density region when smoothing the distribution.

Various techniques have been developed that attempt to resolve this issue, e.g., data reflection (Schuster 1985), boundary kernels (Müller 1991, 1993, 1994), the hybrid method (Hall & Wehrly 1991), generating pseudo-data (Cowling & Hall 1996), data binning and local polynomial fitting (Cheng et al. 1997), and others. One can also invoke asymmetric kernels (such as gamma, lognormal, and inverse Gaussian) and variable bandwidths. In this work we use gamma estimators developed by Chen (2000) and expanded upon by Jeon & Kim (2013) and Hoffmann & Jones (2015).

The gamma PDF with standard gamma function ${\rm{\Gamma }}(\cdot )$ is given by

$$\begin{eqnarray}&&{K}_{k,\theta }(x)=\displaystyle \frac{{x}^{k-1}\exp \left(-\tfrac{x}{\theta }\right)}{{\theta }^{k}{\rm{\Gamma }}(k)},\end{eqnarray} \tag{ 11 }$$

with scale parameter k and shape parameter θ. Chen (2000) take $k={\rho }_{h}(x)$ and $\theta =h$ with random gamma variables X_i to obtain

$$\begin{eqnarray}&&{K}_{{\rho }_{h}(x),h}({X}_{i})=\displaystyle \frac{{X}_{i}^{{\rho }_{h}(x)-1}\exp \left(-\tfrac{{X}_{i}}{h}\right)}{{h}^{{\rho }_{h}(x)}{\rm{\Gamma }}({\rho }_{h}(x))},\end{eqnarray} \tag{ 12 }$$

with

$$\begin{eqnarray*}{\rho }_{h}(x)=\left\{\begin{array}{ll}\tfrac{x}{h}, & \mathrm{if}\ x\geqslant 2h,\\ {\left(\tfrac{x}{2h}\right)}^{2}+1, & \mathrm{if}\ x\in [0,2h).\end{array}\right.\end{eqnarray*}$$

The resulting gamma estimator is given by

$$\begin{eqnarray}&&\hat{f}(x)=\displaystyle \frac{1}{n}\sum _{i=1}^{n}{K}_{{\rho }_{h}(x),h}({X}_{i}).\end{eqnarray} \tag{ 13 }$$

The shape of gamma kernels vary naturally, allowing for different smoothness at different points of the distribution. Further, because gamma kernels are nonnegative, the gamma estimator itself is unlikely to deviate below zero. The bandwidths h depend either on the point of estimation ($h(x);$ a balloon estimator), or on the sample associated with a kernel ($h({X}_{i});$ sample-smoothing estimator). In this analysis we consider the former.

Another challenge for standard KDEs is that regions with few samples have overestimated densities and regions with many are underestimated. Shifted KDEs minimize this bias by moving samples from higher- to lower-density regions. Combining this with balloon estimators (Hoffmann & Jones 2015), one has

$$\begin{eqnarray}&&\hat{f}(x)=\displaystyle \frac{1}{n}\sum _{i=1}^{n}{K}_{{\rho }_{h}(x),h(x)}\left({X}_{i}-{h}^{p}(x)\delta (x\right),\end{eqnarray} \tag{ 14 }$$

where p is the order of the kernel and $\delta (x)$ is the shift. The kernel is shifted by ${h}^{p}(x)\delta (x)$, which vanishes for small bandwidths. For our analyses, we use Python code by Hoffmann & Jones (2015),⁸ where the optimal bandwidth for each kernel is chosen by minimizing the MISE.

4.1. Bounding the DM Distributions

4.2. Observed Sample

To define our sample of pulsars and FRBs, we use the largest aggregation sites for each type of object. For pulsars, we downloaded the Australia Telescope National Facility (ATNF) pulsar catalog (version 1.61; Manchester et al. 2005). For FRBs, we downloaded the FRBCat (downloaded 2020 February 25, verified events only; Petroff et al. 2016).

The Milky Way electron distribution is more complex at low Galactic latitudes owing to contributions from spiral arms, HII regions, and supernova remnants. Electron density models are most complex on size scales smaller than 200 pc and within 1 kpc of the Sun (Cordes & Lazio 2003). To minimize systematic error introduce by the model, we only consider sources more than 200/1000 ≈ 20° from the galactic plane; we also compare the results with a second, more conservative cut to estimate systematic error. We also remove all pulsars within 5° of the Magellanic clouds. For a latitude limit of $| b| \gt 20$°, the samples include 371 pulsars and 83 FRBs. For a latitude limit of $| b| \gt 30$°, the samples include 215 pulsars and 64 FRBs. Owing to the significant decrease in FRB data for $| b| \gt 30$°, the final results presented in this Letter use a Galactic cut of $| b| \gt 20$°.

This analysis requires correcting by the total ${\mathrm{DM}}_{\mathrm{ISM}}$ contribution estimated from the Milky Way. Even at high Galactic latitudes, the electron density models have systematic uncertainties on the order of tens of percent due to modeling errors (Schnitzeler 2012). We estimate ${\mathrm{DM}}_{\mathrm{ISM}}$ with both the NE2001 (Cordes & Lazio 2002, 2003) and YMW16 (Yao et al. 2017) models as a way of estimating potential systematic errors.

We then generated distributions of ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{pulsar}}$ and ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}$, as given by Equations (3) and (6). These are shown in Figures 2(a) and (b). As expected, the majority of ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{pulsar}}$ values are negative with a small tail to positive values. In contrast, the ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}$ distribution is exclusively positive and rises sharply at ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}\approx 64\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$.

Zoom In Zoom Out Reset image size

We applied KDE (with Gaussian and gamma kernels, respectively) to the observed ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{pulsar}}$ and ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}$ distributions to derive PDFs for each. The dark, thick curves in Figures 2(a) and (b) show the results. Also overlaid on the figures are a series of distributions derived from 1000 resampled data sets (100 shown). Table 2 reports the final results for both models on both Galactic latitude samples. In general, we find that the uncertainty on ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}$ values are dominated by the size of the FRB sample. However, the uncertainty on the two distributions is largely insensitive to the Galactic latitude cut. The YMW16 model tends to have slightly smaller ${\mathrm{DM}}_{\mathrm{ISM}}$ values for this sample, which yields larger $\max \left[{\rm{\Delta }}{\mathrm{DM}}_{\mathrm{pulsar}}\right]$ and $\min \left[{\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}\right]$ estimates. However, the separation of these distributions is not sensitive to the Galactic electron density model.

Table 2.  Constraints Derived from (a) Pulsar and (b) FRB Observations

	Latitude	$\max \left[{\rm{\Delta }}{\mathrm{DM}}_{\mathrm{pulsar}}\right]$	${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}$
NE2001	$| b| \gt 20$°	$-2\pm 2\,(\mathrm{stat})\pm 9\,(\mathrm{sys})\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$	$\gt -11\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$
	$| b| \gt 30$°	$-4\pm 3\,(\mathrm{stat})\pm 8\,(\mathrm{sys})\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$	$\gt -13\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$
YMW16	$| b| \gt 20$°	$7\pm 2\,(\mathrm{stat})\pm 9\,(\mathrm{sys})\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$	$\gt -2\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$
	$| b| \gt 30$°	$4\pm 2\,(\mathrm{stat})\pm 8\,(\mathrm{sys})\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$	$\gt -5\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$

	Latitude	$\min \left[{\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}\right]$	${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}$
NE2001	$| b| \gt 20$°	${54}_{-19}^{+40}\,(\mathrm{stat})\pm 9\,(\mathrm{sys})\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$	$\lt 127\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$
	$| b| \gt 30$°	${45}_{-9}^{+39}\,(\mathrm{stat})\pm 7\,(\mathrm{sys})\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$	$\lt 110\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$
YMW16	$| b| \gt 20$°	${63}_{-21}^{+27}\,(\mathrm{stat})\pm 9\,(\mathrm{sys})\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$	$\lt 123\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$
	$| b| \gt 30$°	${52}_{-11}^{+37}\,(\mathrm{stat})\pm 7\,(\mathrm{sys})\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$	$\lt 113\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$

Note. NE2001 and YMW16 are used to model ${\mathrm{DM}}_{\mathrm{ISM}}$ with $| b| \gt 20$° and $| b| \gt 30$°. $\max \left[{\rm{\Delta }}{\mathrm{DM}}_{\mathrm{pulsar}}\right]$ and $\min \left[{\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}\right]$ are calculated at 1σ, and upper and lower limits for ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}$ at 95% c.l. Systematic errors are taken to be the difference between NE2001 and YMW16 estimates. KDE with Gaussian kernels and fixed bandwidths are used to model ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{pulsar}}$, and KDE with gamma kernels and varying bandwidths are used to model ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}$.

Download table as:  ASCIITypeset image

4.3. Simulated Sample

We now simulate ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}$ to explore how the estimation of min[${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}]$ is likely to improve as more FRB data become available and to assess our choice of metric for $\min \left[{\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}\right]$. From Equation (5),

$$\begin{eqnarray}&&{\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}={\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}+{\mathrm{DM}}_{\mathrm{cosmic}}+{\mathrm{DM}}_{\mathrm{host}}.\end{eqnarray} \tag{ 15 }$$

${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}$ has a positive minimum, whereas ${\mathrm{DM}}_{\mathrm{cosmic}}$ and ${\mathrm{DM}}_{\mathrm{host}}$—in principle—have minimums of zero. As such, ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}$ provides a zero-point offset for ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}$, i.e., $\min \left[{\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}\right]\gt 0$. For the following simulation, ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}$ is chosen to be a delta function at 30 $\mathrm{pc}\,{\mathrm{cm}}^{-3}$ and ${\mathrm{DM}}_{\mathrm{host}}$ is approximated by a lognormal distribution with a mean of $\mu =40$ $\mathrm{pc}\,{\mathrm{cm}}^{-3}$ and a standard deviation of $\sigma =0.5$. Other models for these quantities are explored in Appendix A.

To generate a cosmic DM contribution to the simulation, we must adopt a distribution of redshifts for the FRBs. We choose to estimate it from the observed ${\mathrm{DM}}_{\mathrm{FRB}}$ values. Specifically, we adopt a DM–z relation⁹ to convert the observed sample of ${\mathrm{DM}}_{\mathrm{FRB}}$ values to a set of redshifts. Here the observed sample set has $| b| \gt 20$° and ${\mathrm{DM}}_{\mathrm{ISM}}$ is subtracted off with NE2001. We then applied a standard KDE with a Gaussian kernel to build a PDF of the z values from which random draws may be taken. The draws are fed back into the DM–z relationship to obtain the average cosmic contribution to the DM,

$$\begin{eqnarray}&&\langle {\mathrm{DM}}_{\mathrm{cosmic}}(z)\rangle =\int \displaystyle \frac{{\bar{n}}_{e}{ds}}{1+z},\end{eqnarray} \tag{ 16 }$$

where ${\bar{n}}_{e}={f}_{d}(z){\rho }_{b}(z){\mu }_{e}/{\mu }_{m}{m}_{p}$ is the average electron density, f_d is the fraction of cosmic baryons in diffuse ionized gas, ${\rho }_{b}\equiv {{\rm{\Omega }}}_{b}{\rho }_{c}$ is the cosmic baryonic mass density, and ${\mu }_{m}$ and ${\mu }_{e}$ describe properties of helium.

We allow for deviations of ${\mathrm{DM}}_{\mathrm{cosmic}}$ from the average value following the formalism presented in Macquart & Ekers (2018). Our treatment is simpler than theirs; specifically, we assume that the fractional standard deviation of $\langle {\mathrm{DM}}_{\mathrm{cosmic}}\rangle$ is ${\sigma }_{\mathrm{DM}}={{Fz}}^{-1/2}$ with F = 0.2. We may then generate a simulated ${\mathrm{DM}}_{\mathrm{cosmic}}$ distribution based on the z distribution and random draws from a Gaussian characterized by ${\sigma }_{\mathrm{DM}}=1$ and truncated at $\pm 1\sigma$. Throughout, we enforce ${\mathrm{DM}}_{\mathrm{cosmic}}\gt 0$. The resultant ${\mathrm{DM}}_{\mathrm{cosmic}}$ values are added to ${\mathrm{DM}}_{\mathrm{halo}}$ and ${\mathrm{DM}}_{\mathrm{host}}$ to give the simulated PDF of ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}$.

Figure 3(a) shows a realization of this simulated PDF for n = 10,000 draws. This realization has an absolute minimum of ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}=30\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$ and rises sharply due to the host and ${\mathrm{DM}}_{\mathrm{cosmic}}$ contributions. The dark-red curve is the KDE (gamma) using the original data set and the other red curves are distributions generated with resampled data.

Zoom In Zoom Out Reset image size

We explore the sensitivity of the analysis and results to samples size n as follows. For n = 100, 1000, and 10,000, we draw a random set of ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}$ values and model the distributions with KDE (gamma). We then estimate a minimum value from the gradients of the PDFs, i.e., $\min \left[{\mathrm{DM}}_{\mathrm{FRB}}\right]$ is the value that maximizes the slope of the KDE. Since each n PDF is complemented by 1000 PDFs resampled from the original data set, 1000 minima are available for error estimation. The distribution of min[${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}$] values are shown in Table 3. As n increases, the dispersion in min[${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}$] decreases and the central values approach $\approx 34\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$ (Figure 3(b)). Adding more than 10,000 samples has no notable effect on the results.

Table 3.  Simulation Estimates for Different Sample Sizes with $\min \left[{\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB},\mathrm{sim}}\right]=30\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$

No. FRBs	$\min \left[{\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}\right]$	${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}$
100	$37\pm 24\,(\mathrm{stat})\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$	$\lt 114\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$
1000	$35\pm 7\,(\mathrm{stat})\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$	$\lt 55\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$
10000	$34\pm 2\,(\mathrm{stat})\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$	$\lt 44\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$

Note. The second column gives the recovered measurements for $\min \left[{\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}\right]$ at 1σ and the last column gives an upper limit for ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}$ (95% c.l.).

Download table as:  ASCIITypeset image

The simulation estimates are skewed to the left for small n and approach a Gaussian distribution with increased confidence as n increases (Figure 3(b)). While the mean values of the distributions are similar (Table 3), a sample size of n = 100 is inadequate to place a constraint with reasonable confidence. The confidence level does however improve significantly as n approaches 10,000.

Other choices for ${\mathrm{DM}}_{\mathrm{host}}$ are explored to ensure the metric $\min [{\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}]=\max \left[f^{\prime} ({\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}})\right]$ is reasonably robust to changes in the FRB simulation. Results are consistent, as detailed in Appendix A. The smoother the leading edge of ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}$, i.e., the smoother ${\mathrm{DM}}_{\mathrm{host}}$, the more conservative the limits become, and a very sharp edge for ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}$, i.e., a delta function for ${\mathrm{DM}}_{\mathrm{host}}$, is described well by the metric. These cases represent extreme examples of possible host galaxy DM distributions.

The principle empirical result of our work is a conservative upper limit on the DM contribution of the Milky Way halo. At 1σ, ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}={63}_{-21}^{+27}\,(\mathrm{stat})\pm 9\,(\mathrm{sys})\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$ ($| b| \gt 20$°, YMW16). This can be converted to a conservative upper limit of ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}\lt 123\mathrm{pc}\,{\mathrm{cm}}^{-3}$ (95% c.l.). This includes the ISM and halo, and potentially a nonzero contribution from the FRB host galaxy, which is plausibly several tens of $\mathrm{pc}\,{\mathrm{cm}}^{-3}$ (see below). This limit also includes a nonzero contribution from the cosmic web (${\mathrm{DM}}_{\mathrm{cosmic}}$). That contribution is difficult to estimate at present but we note that the lowest-redshift FRB ($z=0.03;$ Marcote et al. 2020) would yield an average ${\mathrm{DM}}_{\mathrm{cosmic}}$ of $\approx 25\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$. A more realistic, yet speculative, upper limit to ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}$ may therefore be $\approx 50\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$.

The results presented include two measurements of uncertainty: systematic uncertainties related to ISM models and statistical uncertainties related to the estimation techniques. Another point to consider is the effect that Galactic latitude has on results. Owing to the complexity of the electron distribution at lower Galactic latitudes, we consider cuts of $| b| \gt 20$° and $| b| \gt 30$°. Results are largely insensitive to this cut; however, the loss of data at $| b| \gt 30$° (371 to 215 pulsars, and 83 to 64 FRBs), motivates a cut of $| b| \gt 20$° for our final analysis.

Pulsar constraints are dominated by uncertainties in modeling ${\mathrm{DM}}_{\mathrm{ISM}}$. We find that, on average, ${\mathrm{DM}}_{\mathrm{ISM}}$ values recovered from NE2001 are $\approx 10\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$ lower than those from YMW16. Given the expectation that ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}\gt 0$, we use YMW16 in our analysis (see Table 2(b)). This gives a final result of ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}\gt -2\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$ (95% c.l.). We note that characterizing the line of sight to MW pulsars may help find HII regions that bias the ${\mathrm{DM}}_{\mathrm{ISM}}$ estimate, allowing for improvement in the pulsar sample.

FRB constraints are predominantly limited by sample size n, i.e., our simulations show a significant improvement as n increases. For an absolute value of ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}=30\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$, limits for $n=100,1000$, and 10,000 are ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}\lt 114\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$, ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}\lt 55\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$, and ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}\lt 44\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$, respectively (95% c.l.). This suggests that once thousands of FRBs have been observed, the constraints will greatly improve.

Even the conservative limit of ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}\lt 123\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$ offers a valuable bound to models of the Galactic halo and the Local Group that our Galaxy resides within. Scenarios that adopt a Galactic halo mass ${M}_{\mathrm{halo}}\approx {10}^{12.2}\,{M}_{\odot }$ that has retained all of its cosmic average of baryons estimate ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}\gt 50\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$ (Prochaska & Zheng 2019; but see Keating& Pen 2020). Furthermore, models that would predict the gas traces the dark matter profile would yield ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}\gt 200\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$ (Figure 4); these are ruled out by our FRB analysis, and also their overestimated X-ray emission (e.g., Fang et al. 2015). Our results also place an upper bound on the average contribution from the Local Group medium, consistent with current estimates (Prochaska & Zheng 2019). Clearly, as the observed FRB sample increases—one expects a dramatic leap from the CHIME survey (CHIME/FRB Collaboration et al. 2018)—the resultant limits may well distinguish between models where the Galaxy has retained the majority of its baryons from those where they have been expelled.

Zoom In Zoom Out Reset image size

To illustrate the potential constraints, Figure 4 shows a model-based estimate for ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}$ for a dark matter halo with mass ${M}_{\mathrm{halo}}={10}^{12.2}{M}_{\odot }$, baryonic mass ${M}_{b}={{\rm{\Omega }}}_{b}/{{\rm{\Omega }}}_{m}{M}_{\mathrm{halo}}\approx 2.4\,\times {10}^{11}{M}_{\odot }$, and that 75% of those baryons are in a diffuse, ionized halo. The density profile is assumed to follow a modified Navarro–Frenk–White (NFW) profile parameterized by y₀ and α (see Mathews & Prochaska 2017; Prochaska & Zheng 2019). The upper limit to ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}$ estimated from our analysis prefers larger $\alpha ,{y}_{0}$ with a strict NFW profile ($\alpha =0,{y}_{0}=1)$ ruled out at high confidence unless ${M}_{b}\ll {{\rm{\Omega }}}_{b}/{{\rm{\Omega }}}_{m}{M}_{\mathrm{halo}}$. Larger $\alpha ,{y}_{0}$ are inferred for our Galaxy and external ones from absorption-line analyses (e.g., Faerman et al. 2017; Mathews & Prochaska 2017).

We emphasize that ongoing FRB projects will offer complementary constraints on the magnitude and distribution of contributions from the host and the cosmic web to the upper limit on ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}$. In particular, well-localized FRBs reveal the host galaxy population and the redshift distribution of FRB events. From follow-up observations of the hosts, one may estimate the DM contribution from the host galaxy ISM through measurements of the Balmer line emission (e.g., Tendulkar et al. 2017; Chittidi et al. 2020). The two systems analyzed thus far yield ${\mathrm{DM}}_{\mathrm{host},\mathrm{ISM}}\approx 50\mbox{--}200\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$. There are other FRBs (e.g., FRB 180924; Bannister et al. 2019) where the Balmer emission is low or even negligible at the FRB location and we infer ${\mathrm{DM}}_{\mathrm{host},\mathrm{ISM}}\lt 50\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$. Within the next year, we expect to have a sample of ∼20 hosts to derive the distribution.

One may additionally translate the estimated stellar mass of the host galaxy into a model-based estimate for the DM contribution from the halo gas of the host (Bannister et al. 2019; Prochaska & Zheng 2019). Current estimates range from $\approx 50\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$ for the most massive hosts (Bannister et al. 2019) to $\lt 20\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$ for FRB 181112 (Prochaska & Zheng 2019). From the redshift distribution of the localized FRBs, one may estimate the minimum typical contribution of ${\mathrm{DM}}_{\mathrm{cosmic}}$ to the ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}$ limit. This bears an important caveat that the selection biases of the localized sample will not match those of the larger ensemble (e.g., due to differences in the radio frequencies and/or flux limit). One will need to account for these differences. Alternatively, one may focus on the analysis of the a localized sample alone once it grows to a sufficient sample size.

Last, we emphasize that other, future observations will also offer constraints on ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}$ independent of FRB analyses. We anticipate high-precision X-ray absorption-line spectroscopy of the Galactic halo from the upcoming Japanese XRISM mission. With a spectral resolution that will greatly exceed current X-ray satellites, the data will yield much more reliable estimates of ${{\rm{O}}}^{+5}$ and ${{\rm{O}}}^{+6}$ column densities across the sky. At the very least, these yield conservative lower limits to ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}$. Another promising yet still unrealized opportunity is to discover pulsars in Andromeda or any other Local Group galaxy. These would offer a strict upper bound on ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}$ or even a well-informed value along that sight line.

We have demonstrated how density estimation techniques can be used to probe the DM—i.e., the line-of-sight electron column density—of the MW Galactic halo. For the corrected ${\mathrm{DM}}_{\mathrm{pulsar}}$ and ${\mathrm{DM}}_{\mathrm{FRB}}$ distributions, we recover $\max \left[{\mathrm{DM}}_{\mathrm{pulsar}}\right]\approx 7\,\pm 2\,(\mathrm{stat})\pm 9\,(\mathrm{sys})\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$ and $\min \left[{\mathrm{DM}}_{\mathrm{FRB}}\right]\approx {63}_{-21}^{+27}\,(\mathrm{stat})\,\pm 9\,(\mathrm{sys})\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$ (1σ uncertainty). Conservative upper and lower limits on the Galactic halo DM are also derived: ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}\gt -2\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$ and ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}\lt 123\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$ (95% c.l.). Here the lower bound given by pulsars reflects only a fraction of the MW halo DM, and the upper bound given by FRBs includes a nominal contribution from the FRB host galaxy and IGM. In the latter case, the localization of FRBs at very low distances and/or on the outskirts of galaxies would establish that the minimum DM would be more representative of the MW halo. Scenarios consistent with this include the collapse of compact objects (e.g., Falcke & Rezzolla 2013) that have been expelled from a host galaxy, as well as more exotic theories such as tiny electromagnetic explosions (which may occur in dark matter halos; Thompson 2017a, 2017b) and cosmic strings (e.g., Vachaspati 2008; Yu et al. 2014; Zadorozhna 2015; Brandenberger et al. 2017).

We do not consider how ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}$ may vary as a function of Galactic latitude. It may be possible with a sample of a couple thousand FRBs per region of sky, but it is left to future work.

Our current estimates cannot yet discern whether the Milky Way has retained its cosmic average of baryons (${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}\,\gt 50\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$); however, in the near future, as more FRBs are reported, results may offer a valuable complement to other analyses. In the least, the methodology provides a reasonable—albeit conservative—estimate of ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}$ and a minimum contribution from ${\mathrm{DM}}_{\mathrm{host}}$. This may discern the viability of Galactic halo models and aid in the search for missing baryons.

We would like to thank the anonymous referee for insightful, thorough, and valuable input. E.P. and J.X.P., as members of the Fast and Fortunate for FRB Follow-up team F⁴ (http://www.ucolick.org/f-40), acknowledge support from NSF grant AST-1911140. C.J.L. acknowledges support under NSF grant 2022546. This work was initiated as a project for the Kavli Summer Program in Astrophysics held at the University of California, Santa Cruz in 2019. The program was funded by the Kavli Foundation, The National Science Foundation, UC Santa Cruz, and the Simons Foundation. We thank them for their generous support. E.P. is supported by a L'Oréal-UNESCO For Women in Science Young Talents Fellowship, by a PhD fellowship from the South African National Institute for Theoretical Physics (NITheP), and by a top-up bursary from the South African Research Chairs Initiative of the Department of Science and Technology (SARChI) and the National Research Foundation (NRF) of South Africa. Any opinion, finding and conclusion or recommendation expressed in this material is that of the authors and the NRF does not accept any liability in this regard.

We postulate that the minimum of the ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}$ distribution can be approximated by $\min [{\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}]=\max \left[f^{\prime} ({\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}})\right]$. This metric is based on the prior that the underlying distribution has a sharp leading edge and is motivated by simulations. To bear weight, the metric must hold for a wide range of reasonable ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}$ distributions.

The MW can be given by a delta function (as its DM is thought to vary by 10 $\mathrm{pc}\,{\mathrm{cm}}^{-3}$) and the cosmic DM distribution can be modeled theoretically. The distribution of host galaxy DMs, however, is unknown. In the main analysis we consider a lognormal distribution with μ = 40 $\mathrm{pc}\,{\mathrm{cm}}^{-3}$ and a standard deviation of σ = 0.5. Here we consider two extreme variations: a delta function at $30\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$ and a broad Gaussian distribution with $\mu =60\,\mathrm{pc}\,{\mathrm{cm}}^{-3}$ and σ = 0.5. The former distribution makes the edge of ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}$ sharper and the latter makes it smoother. The metric is a reasonable approximation for the combined ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}+{\mathrm{DM}}_{\mathrm{host}}$ contribution when each distribution is sharp (Figure A1(a)). When ${\mathrm{DM}}_{\mathrm{host}}$ has a smooth edge, the estimates are more conservative (Figure A1(b)). Thus, provided the leading edge of ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}$ is sufficiently sharp, the metric for determining the distribution minimum can be considered reasonably robust.

Looking at Figure A1, a sample size of n = 1000 appears sufficient to provide an estimate consistent with that of n = 10,000. For n = 100, distributions are wide and skewed to the left, providing results that are clearly premature.

Zoom In Zoom Out Reset image size

Density estimation using field theory (DEFT; Kinney 2014, 2015; Chen et al. 2018) is a newly developed technique specifically developed for the small data regime. It takes a Bayesian field theory approach to density estimation in small data sets using a Laplace approximation of the Bayesian posterior (also see Riihimaki & Vehtari 2014). An advantage of DEFT over standard density estimation methods is that the method does not require the manual identification of critical parameters nor does it require the specification of boundary conditions. The DEFT simulations in this Letter use the Python package Statistics Using Field Theory (SUFTware) by Chen et al. (2018).

Consider n data points (${x}_{1},{x}_{2},\,\ldots ,\,{x}_{n}$) drawn from a known probability distribution ${Q}_{\mathrm{true}}(x)$ with x intervals of length L. We wish to find the best estimate ${Q}^{* }(x)$ of this distribution and the accompanying ensemble of other plausible estimates. Each distribution Q(x) is parameterized by a real field $\phi (x)$, ensuring that Q(x) is positive and normalized:

$$\begin{eqnarray}&&Q(x)=\displaystyle \frac{{e}^{-\phi (x)}}{\int {dx}^{\prime} {e}^{-\phi (x^{\prime} )}}.\end{eqnarray} \tag{ B1 }$$

Using scalar field theory, a prior $p(\phi | {\ell })$ is formulated that favors smooth probability densities. Specifically, Kinney (2015) consider priors of the form

$$\begin{eqnarray}&&p(\phi | {\ell })=\displaystyle \frac{{e}^{-{S}_{{\ell }}^{0}\left[\phi \right]}}{{Z}_{{\ell }}^{0}},\end{eqnarray} \tag{ B2 }$$

with action

$$\begin{eqnarray}&&{S}_{{\ell }}^{0}\left[\phi \right]=\int \displaystyle \frac{{dx}}{L}\displaystyle \frac{{{\ell }}^{2\alpha }}{2}{\left({\partial }^{\alpha }\phi \right)}^{2},\end{eqnarray} \tag{ B3 }$$

and partition function

$$\begin{eqnarray}&&{Z}_{{\ell }}^{0}=\int { \mathcal D }\phi {e}^{-{S}_{{\ell }}^{0}[\phi ]}.\end{eqnarray} \tag{ B4 }$$

Here, ℓ gives the length scale below which ϕ fluctuations are strongly damped and $\alpha \gt 0$ is an integer in the range [1, ..., 4] that determines the smoothness. The resultant posterior is given by

$$\begin{eqnarray}&&p(\phi | \mathrm{data},{\ell })=\displaystyle \frac{{e}^{-{S}_{{\ell }}\left[\phi \right]}}{{Z}_{{\ell }}},\end{eqnarray} \tag{ B5 }$$

with nonlinear action

$$\begin{eqnarray}&&{S}_{{\ell }}\left[\phi \right]=\int \displaystyle \frac{{dx}}{L}\left\{\displaystyle \frac{{{\ell }}^{2}\alpha }{2}{\left({\partial }^{\alpha }\phi \right)}^{2}+{nLR}\phi +{{ne}}^{-\phi }\right\},\end{eqnarray} \tag{ B6 }$$

and partition function

$$\begin{eqnarray}&&{Z}_{{\ell }}=\int { \mathcal D }\phi {e}^{-{S}_{{\ell }}[\phi ]}.\end{eqnarray} \tag{ B7 }$$

$R(x)=\tfrac{1}{n}{\sum }_{i=1}^{n}\partial (x-{x}_{i})$ is a histogram that summarizes the data.

Maximum a posteriori (MAP) density estimation approximates the posterior $p(\phi | \mathrm{data},{\ell })$ as a δ function given by the mode of the posterior, at which the action ${S}_{{\ell }}\left[\phi \right]$ is then minimized. It has been shown that even without imposing boundary conditions on ϕ, ${S}_{{\ell }}\left[\phi \right]$ has a unique minimum (Kinney 2015). The optimal length scale ${{\ell }}^{* }$ is identified by maximizing the Bayesian evidence $p(\mathrm{data}| {\ell })$.

The uncertainty in the DEFT estimate Q* is determined by sampling the Bayesian posterior,

$$\begin{eqnarray}&&p(Q| \mathrm{data})=\int {dlp}({\ell }| \mathrm{data})p(Q| \mathrm{data},{\ell }),\end{eqnarray} \tag{ B8 }$$

by first drawing ℓfrom $p({\ell }| \mathrm{data})$ and then drawing Q from $p(Q| \mathrm{data},{\ell })$. Laplace approximation is used to estimate $p(Q| \mathrm{data},{\ell })$ by constructing a Gaussian centered at its MAP value. This gives the Laplace posterior,

$$\begin{eqnarray}&&{p}_{\mathrm{Lap}}(Q| \mathrm{data})=\int {dlp}({\ell }| \mathrm{data}){p}_{\mathrm{Lap}}(Q| \mathrm{data},{\ell }),\end{eqnarray} \tag{ B9 }$$

from which an ensemble of distributions Q can be sampled. Some of the Qs, however, are clearly not representative of the underlying distribution. Importance resampling is thus used to remove unfavorable distributions, where each ϕ is given a weight,

$$\begin{eqnarray}&&{w}_{{\ell }}[\phi ]=\exp \left({S}_{{\ell }}^{\mathrm{Lap}}[\phi ]-{S}_{{\ell }}[\phi ]\right),\end{eqnarray} \tag{ B10 }$$

proportional to its probability of being drawn (Chen et al. 2018). DEFT uses importance resampling with replacement; however, for this work we invoke importance resampling without replacement.

When a posterior turns out to be a poor approximation of the target distribution, a few of the sampled distributions are given very large weights and the majority are given small weights (Skare et al. 2003; Gelman et al. 2014). When resampling with replacement, the heavily weighted distributions become significantly over represented. In our case, ∼60%–70% of the sampled distributions were duplicates, which lead to notable bias when calculating the upper and lower bounds of ${\mathrm{DM}}_{\mathrm{MW},\mathrm{halo}}$. As such, we use a set of the most probable distributions, with limited replications. Specifically, we select 500 out of 1000 distributions via importance sampling without replacement. This lowered the duplication rate to ∼10%.

We approximate the FRB distribution described in Section 4.3 using DEFT for n = 100, n = 1000, and n = 10,000. Even for large n DEFT is unable to adequately describe the sharp edge of the simulated distribution. In Figure B1(a), the PDF tail extends below zero, violating the physical condition that ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}\gt 0$. Further, the PDF cuts straight through the front of the simulated distribution and so bypasses the structure of the distribution's edge.

Zoom In Zoom Out Reset image size

A standard statistical technique for estimating the maximum values of an ensemble is to fit it with a GEV PDF (e.g., Coles 2001). This technique, however, is most applicable for assessing the upper limit of a distribution with a long tail. For ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}$, this holds for the largest values but the lowest values rise sharply as one may expect from the MW and host contributions.

Nevertheless, we attempted to estimate the minimum of ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}$ following the standard practice of assessing the maximum of the negative of the distribution (Coles 2001). The results reported a minimum value at effectively infinite confidence at the lowest ${\rm{\Delta }}{\mathrm{DM}}_{\mathrm{FRB}}$ in the distribution and we found the results were unstable to random sampling.