Bounding the photon mass with cosmological propagation of fast radio bursts

Photon is the fundamental quantum of electromagnetic fields, whose mass, $m_{\gamma}$, should be strictly zero in Maxwell's theory. But not all theories adopt this hypothesis. If the rest mass of the photon is not zero, there will be an additional time delay between photons of different frequencies after they travel through a fixed distance. By analyzing the time delay, we can measure or constrain the photon mass. Fast radio bursts (FRBs) -- transient radio bursts characterized by millisecond duration and cosmological propagation -- are excellent astrophysical laboratories to constrain $m_{\gamma}$. In this work we use a catalog of 129 FRBs in a Bayesian framework to constrain $m_{\gamma}$. As a result, we obtain a new bound on the photon mass, $m_{\gamma} \leq 3.1\times 10^{-51}\rm\,kg\simeq 1.7 \times 10^{-15}\,eV/c^2$ ($m_{\gamma} \leq 3.9\times 10^{-51}\rm\,kg \simeq 2.2 \times 10^{-15}\,eV/c^2$) at the $68\%$ $(95\%$) confidence level. The result represents the best limit purely from kinematic analysis of light propagation. The bound on the photon mass will be tighter in the near future with increment in the number of FRBs, more accurate measurement of the redshift for FRBs, and refinement in the knowledge about the origin of dispersion measures (DMs).


Introduction
"Constancy of light speed" is an important postulate of Einstein's theory of special relativity [1], and it is adopted in the general relativity and quantum field theories. In quantum mechanics, the particle-wave duality translates the constancy of light speed into the "masslessness of photons", indicating that the rest mass of the photon should be strictly zero. Nevertheless, there are still a number of theories which challenge this postulate, including the famous Maxwell-de Broglie-Proca theory where photons are massive [2,3]. In the Standard Model of particle physics, the gauge symmetry of quantum electrodynamics prevents photons from acquiring a non-zero mass. Therefore, confirmation of the masslessness of photons can also be viewed as a need for a gauge theory of quantum electrodynamics.
In the viewpoint of testing fundamental principles and pushing the boundary of contemporary physical theories to a larger extent, we need to tighten the empirical limits on the photon mass from various experiments. A non-zero mass of the photon can cause a series of experimental or observational consequences. There are already many efforts to constrain the rest mass of the photon [4,5,6,7,8], such as early experiment in testing the Coulomb law [9], toroid Cavendish balance [10,11], gravitational deflection of electromagnetic waves [4], magneto-hydro-dynamic phenomena of the Solar wind [12,13,14], spindown of pulsar [15], Jupiter's magnetic field [16], the mechanical stability of magnetized gas in galaxies [17], "black hole bombs" [18,19] and so on. These are based on the dynamics to test the photon mass, and they depend on specific massive photon theories. Compared with the dynamic tests, kinematic tests [20,21,22,23,24,25,26,27,28,29] are cleaner since they do not involve any dynamic hypothesis thus are largely independent of the underlying massive photon theory. In this paper, we constrain the photon mass from the propagation of electromagnetic waves of fast radio bursts (FRBs) [30], which is a kinematic method. Due to the uncertainty principle of quantum mechanics and the finite age of our Universe (∼ 10 10 years), there is an ultimate lower limit on the photon rest mass that we in principle could probe. It equals to m γ ≈ /∆t 0 c 2 0 ≈ 10 −69 kg, where is the reduced Planck constant, ∆t 0 is the age of the Universe, c 0 is the limiting velocity for photons of infinite frequency [5,31].
Assuming that photons have a non-zero rest mass, m γ , according to Einstein's special relativity, the dispersion equation can be written as, where p is the momentum. We can get the group velocity for a photon with energy E = hν as, where h is the Planck constant and ν is the photon frequency. The last term only holds when m γ hν/c 2 0 7 × 10 −42 (ν/GHz) kg. It is suggested by Eq. (2) that when the energy is smaller, the relative modification is larger. We constrain the photon mass by the accumulative time delay from propagation of light. Since a shorter emitting time-scale, a longer propagation distance, and a lower energy of the photons can provide a larger time delay, FRBs have some features that deserve our attention in a test of such kind [21,22,23,25,28,29]: (i) FRBs belong to a kind of instantaneous events with a time scale in milliseconds; (ii) FRBs are widely considered to be cosmological objects and their redshifts can be significant in general; (iii) FRBs are observed at radio frequency which means a lower energy for photons, so they can provide a tight bound on m γ via Eq. (2).
The above features make FRBs a great celestial laboratory to limit m γ . Previous bounds on the photon mass from FRBs are collected in Table 1 [21,23,25,27,28]. With considerably increased radio data from telescopes such as the Canadian Hydrogen Intensity Mapping Experiment (CHIME), the Australian Square Kilometre Array Pathfinder (ASKAP), the Molonglo Observatory Synthesis Telescope (UTMOST), the Five-hundred-meter Aperture Spherical radio Telescope (FAST), and the Apertif Radio Transient System installed in the Westerbork Synthesis Radio Telescope (WSRT), the observation of FRBs has developed a lot in recent years [32]. Therefore, it is timely to use the most recent data to put an updated bound on the photon mass.
In this paper, we collect a catalog of FRBs to bound the photon mass. Figure 1 shows the sky distribution of FRBs which are used in this work [33]. 1 Based on previous work, we use the observed dispersion measure (DM) to bound the photon mass. In order to do so, we estimate the redshifts of FRBs through package FRUITBAT [34] 2 and estimate the DM of host galaxies in a Bayesian framework. Also with the help of a Bayesian method, we obtain a new photon mass bound, m γ ≤ 3.1 × 10 −51 kg 1.7 × 10 −15 eV/c 2 (m γ ≤ 3.9 × 10 −51 kg 2.2 × 10 −15 eV/c 2 ) at the 68% (95%) confidence level (C.L.), which improves the previous bound of the photon mass from similar methods by a factor of 2-3.
The paper is organized as follows. In the next section, we introduce a theoretical framework for data analysis which includes a hypothesis for the ν −2 -behaved time delay, and a Bayesian framework to estimate the DM of host galaxies and the photon mass. In Sec. 3, we present our results on the estimation of the DM of host galaxies and the bound of the photon mass from a combination of 129 FRBs. In Sec. 4 we present some discussions on the results.

Theoretical framework
In Sec. 2.1 we introduce a hypothesis for the ν −2behaved time delay of FRBs in observation, and then utilize a Bayesian framework in Sec. 2.2 to analyze the observed FRB data. The analysis includes the estimation of the DM of host galaxies and the bound of the photon mass.

A hypothesis on the time delay
If photons have a non-zero rest mass, there will be a time delay in the arrival times of FRB photons observed at different frequencies. Observations of FRBs show an indisputable ν −2 -dependent relation in the time delay, ∆t ∝ ν −2 [33]. In our calculation, we attribute the time delay to two causes: (i) the time delay of electromagnetic waves caused by the propagation of photons through the ionized medium, and (ii) the time delay caused by a non-zero rest mass of photons.
On one hand, the astrophysical term, ∆t DM , is caused by the interaction between the ionizing medium and the  electromagnetic waves. For a photon with energy, E = hν, relative to a photon with an infinite energy, the time delay ∆t DM is [35], where the plasma frequency ν p ≡ n e e 2 /4π 2 m e 0 with e the charge of an electron, n e the number density of electrons, m e the mass of an electron, and 0 the permittivity of free space. The definition of DM astro is the integral of the electron number density along the propagation path, DM astro ≡ n e dl [35]. Considering a non-negligible cosmological redshift, one has DM astro ≡ n (0) e (1 + z) −1 dl [36], where z is the redshift, and n (0) e is the electron number density in the rest frame.
There are multiple sources of DM astro , mainly including contributions from the Milky Way (MW) galaxy, DM MW , from the intergalactic medium (IGM), DM IGM , and from the host galaxy, DM host . Therefore, the total DM contributed from the electromagnetic wave propa-gating through the ionized medium can be written as, where we have included contributions from the nearsource plasma (e.g. supernova remnant, pulsars' wind nebula, HII region, etc. [37]) collectively in DM host . Generally, people consider that the halo of the MW is mainly composed of dark matter and there is no other evident interaction but gravity between dark matter and electromagnetic waves. So we ignore the contribution from the MW halo to the total DM in our work [38,39]. On the other hand, comparing fiducial photons of infinite frequency with Eq. (2), the time delay caused by the mass of the photon, denoted as ∆t m γ , is, where the Hubble constant H 0 = 67.66 ± 0.42 km s −1 Mpc −1 [40], and H γ (z) is a dimension-less function of the source redshifit z, where the vacuum energy density Ω Λ = 0.6889 ± 0.0056, and the matter energy density Ω m = 0.3111 ± 0.0056 [40].
In our hypothesis, the total time delay can be written as, ∆t = ∆t DM + ∆t m γ , and these two time delays are both ν −2 -behaved which means that the total observed DM, DM obs , can be written in a similar form, DM obs = DM astro +DM γ , where the equivalent DM corresponding to the time delay caused by the mass of photons is [25],

Bayesian framework
We work in a Bayesian framework developed in Ref. [25]. In Bayesian analysis, the posterior distribution of parameter θ is given as, where D is data, H is the hypothesis given in Sec. 2.1, I is all other relevant background knowledge. P (D|θ, H, I) is the likelihood for the data, P (θ|H, I) is the prior on the parameter θ, and P(D|H, I) is model evidence. With the prior and the likelihood, we can get the posterior distribution of parameter θ from Eq. (8). In our work, we apply Bayesian analysis on the parameter estimation of DM host and m γ as follows. The model evidence merely plays a normalization factor in our problem.
We have separated DM obs into DM astro and DM γ , and will analyze each term in DM astro to deduce the likelihood of DM host .
Firstly, there are some electron distribution models with due uncertainties for the MW from combining different observational results. The most widely used models are the NE2001 model [41] and the YMW16 model [42]. We mainly use the NE2001 model in our work. As for the uncertainty estimation, we take the difference between the above two models as a standard deviation. In principle, there are potential contributions from the Galactic halo. As we discussed in Sec. 2.1, considering that the contribution from the Galactic halo is small and negligible in most cases, we do not count it into DM obs . Even if there is a certain contribution, it can be covered by our estimation of the error in DM MW . Excluding DM MW , we can get the excess DM of the astrophysical term, DM FRB ≡ DM astro − DM MW = DM host + DM IGM .
Secondly, for the DM from IGM, DM IGM , we use a one-parameter model to account for the scattering in the electrons from foreground structures [43,44]. Its probability distribution is, where the normalization constant A guarantees the integral of the probability function to 1. The fractional fluctuation of the dispersion, s = Fz −0.5 , is caused by the cosmological large-scale structures which may exist in the foreground, and F is used to quantify the strength of the baryonic feedback. In Ref. [43], the authors selected α = 3 and β = 3. The parameter C 0 can be derived from the demand that the mean of the distribution of ∆ satisfies ∆ = 1. Taking the average with respect to the sky location, DM IGM at redshift z is, where m p is the proton mass, Ω b = 0.0487 ± 0.0007 is the baryonic matter energy density [40], f IGM 0.83 is the baryonic fraction of IGM [45]. In the electron mass fraction, g(z) ≡ 3 4 y 1 χ e,H (z) + 1 8 y 2 χ e,He (z), the hydrogen and helium mass fractions are normalized to their typical values 3/4 and 1/4 respectively. We have used y 1 ∼ 1, y 2 ∼ 1, and χ e is the ionized fraction in IGM, which for hydrogen χ e,H (z) 1 at z < 6, and for helium, χ e,He (z) 1 at z < 2. Therefore, we will only consider FRBs with z < 2, and we have [46], where we have rewritten the Hubble constant H 0 in a dimensionless form, using h 70 = H 0 /(70 km s −1 Mpc −1 ). For the last term in DM astro , we adopt a model for the probability distribution of DM host . What is different from the above two terms is that, unlikely the NE2001 model or Eq. (12), there is no persuasive astrophysical model yet for DM host . So we take advantage of Bayesian analysis, and use a lognormal distribution model in our work [43]. There are two main considerations: (i) the lognormal distribution is a positive definite distribution, and (ii) it will have asymmetric tails at high values.
These high values have many potential origins like the gas local to the FRB (such as the HII region and the circumstellar medium) or the interstellar medium. These two characteristics agree with the physical meaning of DM host . Specifically, we use the probability distribution of DM host as [43], where µ, σ host are the estimated parameters which are used to calculate the expectation and variance of DM host respectively. Note that in this equation we are assuming the standard "pc cm −3 " unit for DM. Strictly speaking, the DM host obtained in Ref. [43] uses the assumption m γ = 0. It would be ideal to allow for a non-zero photon mass. However, compared with DM IGM , DM host is not that large to make great difference on the final result. So in this part, we temporally ignore the contribution of the photon mass term, and use DM obs as DM astro . This treatment is tentatively feasible if we are putting a bound on the photon mass, as we are doing it in this work, but will need to be refined if we are to measure a non-zero m γ .
We proceed to estimate the model likelihood by computing the likelihoods of all FRBs and combining the P IGM in Eq. (9) and P host in Eq. (14), where N FRBs is the number of FRBs with the index i running over FRBs. In the above equation, the likelihood for each FRB is Finally, we get a likelihood in a four-parameter space, (Ω b h 70 , F, µ, σ host ), among which µ, σ host are the (dimensionless) parameters that we need to calculate the DM host . We choose the priors of each parameter to be uniform in the following ranges: Ω b h 70 ∈ [0.015, 0.095], F ∈ [0.01, 0.5], e µ ∈ [20,200], and σ host ∈ [0.2, 2]. The priors that we choose are based on reasonable consideration and they were used in other literature as well. If one has chosen other priors, as long as that are not extreme, we do not expect significant changes to our results.
From discussions so far, we see that once we get DM host , with above DM IGM and DM MW , we can get DM astro and eventually deduce the value of DM γ to bound m γ . Now we proceed to obtain the likelihood L 2 for m γ . As we discussed, we can get DM MW in the NE2001 model and the YMW16 model. We take the result of the NE2001 model as the central value and the difference between the two models as the standard deviation. We take the average DM IGM as the expectation value using Eq. (12), and the contribution of foreground structures is covered in the uncertainty of DM IGM . We associate a 20% uncertainty to Ω b f IGM in the hope to cover the inhomogeneity of IGM along different lines of sight. We use our result in the parameter estimation of DM host described earlier. To take the cosmological evolution into account, we have multiplied it by a factor of (1 + z) −1 .
After estimating DM astro , we deduce DM γ which is related directly to m γ . We apply a Bayesian framework to bound the photon mass. Similarly, we adopt the posterior distribution, . (17) Assuming independence of observations of different FRBs, we construct the logarithm of likelihood as [25], where DM γ is given in Eq. (7), DM astro is composed of the above listed terms in our Markov-chain Monte Carlo (MCMC) runs, σ includes all uncertainties added in quadratic (including uncertainties in DM obs , DM MW , DM IGM , DM host , and the redshift z).
We adopt a uniform prior on m γ in the range [10 −69 , 10 −42 ] kg. The uncertainty principle of quantum mechanics gives the lower bound since the age of our Universe is finite, while the upper bound is chosen to satisfy the demands in the approximation in Eq. (2). This prior covers a very wide range which we view as a conservative choice. In practice, for a uniform prior, because the measure from 0 to 10 −69 kg is totally negligible in any sense, if we have chosen the lower end at m γ = 0 which encloses the zero mass case in the Maxwell's theory, our results in the next section stay practically the same. Table 2: Repeating FRBs that are used to estimate DM host . Sky position is given in Galactic longitude, g l , and Galactic latitude, g b . DM is in unit of pc cm −3 , where DM obs is from the fitting of the ν −2 -behavior in the frequency-dependent time delay, and DM MW is based on the NE2001 Galactic electron density model [41]. Redshift is given by the localization method.

FRB
Telescope g l (deg)

Results
As mentioned, we use MCMC techniques to calculate the posterior of DM host and m γ . We use PYTHON implementation of an affine-invariant MCMC ensemble sampler [47], emcee, 3 to perform simulations.
We use a PYTHON package, FRUITBAT [34], to estimate the redshifts of FRBs in the catalog. FRUIT-BAT is aimed to assist estimation of redshifts, energies and the Galactic DM contributions of FRBs. Since FRUITBAT only considers the line-of-sight DM MW and DM IGM , without the host galaxy and foreground structures that contribute to the DM of many lines of sight, its result could be close to, but larger than, what is given in the catalog already. Comparing the redshift determined by direct measurements (e.g. from localization of the host galaxy for the repeating FRBs) with that obtained by FRUITBAT [33,48,49,50], we choose 30% uncertainty for these redshifts. For the repeating FRBs in the catalog, we only use the first burst of each repeating FRB, because our calculation require that every single burst is independent. Meanwhile, we use FRUIT-BAT to estimate DM MW in the NE2001 model and the YMW16 model.
There are several repeating FRBs observed in recent years. They represent precious sources in studying the origin of FRBs. Before calculating DM host , we first select sources with redshift measurement. We excluded FRB 190614 [51] and FRB 20200120E [52] since their redshift measurements are not that accurate and need more follow-up observations. We also ignore FRB 121102 [48] and FRB 180916 [50] owing to their low Galactic latitude. The rest 7 repeating FRBs that we use are listed in Table 2.
As stated in the last section, we use the NE2001 model and Eq. (12) to estimate DM MW and DM IGM re-spectively. In our MCMC runs in estimating DM host , there are 30 chains with 15000 samples per chain after discarding the "burn-in" samples [47]. In Fig. 2, we show the samples returned by the MCMC sampler after discarding the "burn-in" samples. We choose our model for DM host to follow a log-normal distribution, characterized by an expectation e µ+σ 2 /2 and a median e µ . The standard deviation of the distribution is (e σ 2 − 1) 1/2 e µ+σ 2 /2 . Since the standard deviation given by MCMC is larger, we use the error of the expectation of DM host as the standard deviation in the subsequent calculation but not the one calculated by the log-normal model. Our calculation gives rather conservative results with DM host = 124 pc cm −3 , and σ = 90 pc cm −3 .
Up to March 15, 2021, there are 140 available sources in the FRB catalog [33]. We ignore FRBs whose errors of the observed DMs are not given, and ones whose redshifts are larger than 2 in order to satisfy the approximation in Eq. (11). There are 129 FRBs left in our calculation. Ideally, we would calculate the loglikelihood in the form of Eq. (18) which means simultaneously analyzing all FRB data in one run, but the computational cost is quite high due to the large dimensionality of the parameter space. Considering that FRBs are uncorrelated with each other, we combine the individual bounds on m γ from the MCMC run of every single FRB. Such an approach significantly reduces the computational burden. It is worth noting that we discard weak bounds whose average values are larger than 10 −50.25 kg because they give little contribution on the final combined bound of m γ from all FRBs.
In our MCMC runs in estimating m γ , there are 26 chains with 15000 samples per chain after discarding the burn-in samples. In Fig. 3 at 68% C.L. and 95% C.L. respectively. It should be noted that we have ignored the contribution of DM γ when we estimate the DM host since the degeneracy of DM host and DM γ is very high. They are completely degenerate for a single FRB. However we want to stress that we obtained DM host through a Bayesian estimation of more than one hundred sources, and the degeneracy of DM host and DM γ has been reduced to some degree, since it is very unlikely for the Nature to hide a common value of the photon mass for all FRBs through tuning the DM associated with each FRB. That is the advantage of using a catalog of FRBs [25]. As we discussed in Sec. 2.2, this method is statistically feasible because we are bounding the photon mass instead of measuring a non-zero photon mass. In the future, if one is at the stage to measure a non-zero photon mass, one needs to design more sophisticated methods to overcome or reduce the degenecy between DM host and DM γ .
However, for a conservative consideration, we calculate the photon mass bounds with selected FRBs by Eq. (7) under the situation that DM astro = 0, in which the calculation is completely model-independent and most conservative. The histogram of the obtained photon mass bounds is shown in Fig. 4. Also, we read from the histogram at 95% C.L. a very conservative bound for the photon mass, m γ ≤ 1.9 × 10 −50 kg.

Discussion
Though Maxwell's theory is a well established theory in physics, it is still intriguing to study its fundamentals to unprecedented precision due to its tremendous importance in our everyday life and physical theories' integrity. During past decades, there are various bounds on the photon mass from dynamic tests. We give a brief overview for comparison with ours in Eqs. (19)(20). In terrestrial laboratories, physicists have tested the Coulomb law to great precision, and obtained a photon mass, m γ < 1.6 × 10 −44 kg, in the Proca theory [9]. The Cavendish balance, which generates a magnetic dipole vector potential moment with a suspended toroidal coil and let the photon mass exhibit as the torque, gives a bound, m γ < 1.2 × 10 −48 kg [5,10,11]. As astrophysical observations reveal their potential in testing fundamental physics [5], a few intriguing tests were conducted. Magneto-hydro-dynamic phenomena of the Solar wind, using observations of the flow pattern, give a bound m γ < 1.5 × 10 −48 kg [14]. Another interesting bound, m γ < 7.1 × 10 −56 kg, comes from the so-called "black hole bombs", where massive fields around rotating black holes trigger a superradiant instability [18]. We can see that the constraints given by different tests vary greatly, spanning more than ten orders of magnitude. As we stressed, these bounds are in general dependent on the underlying massive photon theory.
In contrast, our limit in this paper is gained from a purely kinematic test which is cleaner since it is generally quite independent of the underlying massive photon theory. In the history of bounding the photon mass by kinematic tests, pulsars were the earlier celestial objects that were used [53]. But even today, pulsars are mostly observed in the MW which means that the distance of the propagation is limited. With the discovery of γ-ray burst (GRB), it was found that GRBs are at a cosmological distance, which makes them a better celestial ob-ject than pulsars. In Refs. [24,54], GRB were used to bound m γ , and the tightest bound is m γ < 1.1×10 −47 kg. Generically speaking, γ-ray photons are of high energy, and GRBs are more suitable to test ultra-violet aspects of theories, e.g. in testing the Lorentz symmetry of spacetime [55,56,57]. In contrast, FRBs have low energy photons, and as emphasized earlier, they are more suitable to test the photon mass.
In recent years, the ability in detecting FRBs has been greatly improved, so that the potential in FRBs to constrain the photon mass has been continuously explored (see Table 1). From the study of the single burst such as FRB 121102 to the study of 9 repeating bursts [29] or the study of subbursts of FRB 121102 [28], the research on FRBs is versatile in probing the lower limit of the photon mass. In our work, we give the tightest bound on the photon mass in kinematic tests so far with cosmological propagation of FRBs in Eqs. (19)(20).
Our improvement is due to that, (a) more FRBs are used to bound the photon mass, and (b) the redshift estimation done by the FRUITBAT package [34].
In the near future, the bound given by FRBs will be even tighter for several reasons: (i) expected increment in the number of FRBs; (ii) more accurate measurement of the redshifts of FRBs; (iii) refinement in the knowledge about the origin of DMs, such as a better model for the electron distribution of the MW, a better understanding on the various astrophysical contributions to DM IGM like the evolution of f IGM which may correlate to the redshift [58]; (iv) more information on the host galaxy and FRBs themselves with improved analysis on the signal of FRBs and their Faraday rotation, and so on; (v) and detection at lower frequency, which is advantageous in constraining m γ [see Eq. (5)].
With the increment in the number of FRBs, the population research comes to reality. Through PYTHON packages such as frbpoppy [59], we hope to make predictions on the detection capabilities of current and upcoming telescopes, and give projected prediction on the bound of the photon mass. Such a study lays beyond the scope of current paper.