Combined limit on the photon mass with nine localized fast radio bursts

The postulate of the "constancy of light speed" is a major pillar of Maxwell's electromagnetism and Einstein's theory of special relativity. This postulate implies that the photon, the fundamental quanta of light, should be massless. Nevertheless, several theories have challenged the concept of the massless photon, starting from the famous de Broglie-Proca theory (de Broglie 1922; Proca 1936) and nowadays proposing models with massive photons for dark energy (e.g., Kouwn et al. 2016). Within that context, determining the rest mass of the photon has historically been an effective way to test the validity of this postulate. However, according to the uncertainty principle, it is impossible to perform any experiment that would fully prove that the photon rest mass is exactly zero. Adopting the age of the universe of about 10¹⁰ years, the ultimate upper bound on the photon rest mass would be m_γ ≤ ℏ/Δ tc² ≈ 10⁻⁶⁹ kg (e.g., Goldhaber & Nieto 1971; Tu et al. 2005). The best experimental strategy is therefore to place ever stricter upper limits on m_γ and push the experimental limits more closely towards the ultimate upper limit.

From a theoretical perspective, a nonzero photon mass can be incorporated into electromagnetism straightforwardly via the Proca equations, which are the relativistic generalization of Maxwell's equations. Using the Proca equations, it is possible to consider some observable effects associated with massive photons. One can then constrain the photon rest mass by searching for such effects (see Goldhaber & Nieto 1971; Lowenthal 1973; Tu et al. 2005; Okun 2006; Goldhaber & Nieto 2010; Spavieri et al. 2011 for reviews). Based on such effects, various experimental methods have been developed to set upper limits on the photon rest mass, such as Coulomb's inverse square law (Williams et al. 1971), Ampère's law (Chernikov et al. 1992), torsion balance (Lakes 1998; Luo et al. 2003b,a), gravitational deflection of electromagnetic waves (Lowenthal 1973; Accioly & Paszko 2004), Jupiter's magnetic field (Davis et al. 1975), magnetohydrodynamic phenomena of the solar wind (Ryutov 1997, 2007; Retinò et al. 2016), mechanical stability of the magnetized gas in galaxies (Chibisov 1976), stability of plasma in Coma cluster (Goldhaber & Nieto 2003), black hole bombs (Pani et al. 2012), pulsar spindown (Yang & Zhang 2017), the frequency dependence in the velocity of light (Lovell et al. 1964; Warner & Nather 1969; Schaefer 1999; Wu et al. 2016; Bonetti et al. 2016, 2017; Zhang et al. 2016; Shao & Zhang 2017; Wei et al. 2017; Wei & Wu 2018; Xing et al. 2019), and so on. Among these methods, the most immediate and robust is to measure the frequency-dependent dispersion of light.

If the photon rest mass is nonzero, then the dispersion relation is given by

$$\begin{eqnarray}&&E=\sqrt{{p}^{2}{c}^{2}+{m}_{\gamma }^{2}{c}^{4}}\,.\end{eqnarray} \tag{ 1 }$$

The massive photon group velocity then takes the form:

$$\begin{eqnarray}&&\upsilon =\displaystyle \frac{\partial E}{\partial p}=c\sqrt{1-\displaystyle \frac{{m}_{\gamma }^{2}{c}^{4}}{{E}^{2}}}\approx c\left(1-\displaystyle \frac{1}{2}\displaystyle \frac{{m}_{\gamma }^{2}{c}^{4}}{{h}^{2}{\nu }^{2}}\right)\,,\end{eqnarray} \tag{ 2 }$$

where ν is the frequency and the last approximation holds when ${m}_{\gamma }\ll h\nu /{c}^{2}\simeq 7\times {10}^{-42}\left(\displaystyle \frac{\nu }{{\rm{GHz}}}\right)\,\,{\rm{kg}}$. Equation (2) means that low frequency photons propagate in a vacuum slower than high frequency ones. Two photons with different frequencies, if emitted simultaneously from the same source, would reach us at different times. Moreover, it is obvious that observations of shorter time structures in lower energy bands from astronomical sources at cosmological distances are particularly sensitive to the photon mass. Thanks to their short time durations, long propagation distances, and low frequency emissions, extragalactic fast radio bursts (FRBs) provide the most promising celestial laboratory so far for constraining the photon mass through the dispersion method (Wu et al. 2016; Bonetti et al. 2016, 2017; Shao & Zhang 2017; Xing et al. 2019). For instance, taking the controversial redshift measurement of FRB 150418 (z = 0.492) ¹ Wu et al. (2016) set a tight upper limit on the photon mass as m_γ ≤ 5.2 × 10⁻⁵⁰ kg, which is three orders of magnitude better than the previous best constraint from other astronomical sources with the same method (see also Bonetti et al. 2016). Later, Bonetti et al. (2017) used the reliable redshift measurement of FRB 121102 to constrain the photon mass to be m_γ ≤ 3.9 × 10⁻⁵⁰ kg. Using the time-frequency structure of subpulses in FRB 121102, Xing et al. (2019) obtained a stricter limit on the photon mass of m_γ ≤ 5.1 × 10⁻⁵¹ kg.

From observations, the pulse arrival times of FRBs strictly follow the 1/ν² law, which is in good agreement with the propagation of radio signals through a cold plasma. However, a similar frequency-dependent dispersion $\propto {m}_{\gamma }^{2}/{\nu }^{2}$ (see Eq. (2)) could also be resulted from a nonzero photon mass. The similarity between the frequency-dependent dispersions induced by the plasma and nonzero photon mass effects complicates the derivation of a constraint on m_γ . Fortunately, since the two dispersions have different dependences on the redshift z, it has been suggested that they could in principle be distinguished when a handful redshift measurements of FRBs become available and making the sensitivity to m_γ to be improved (Bonetti et al. 2016, 2017; Bentum et al. 2017). To date, nine FRBs, including FRB 121102 (Spitler et al. 2016; Chatterjee et al. 2017; Tendulkar et al. 2017), FRB 180916.J0158+65 (Marcote et al. 2020), FRB 180924 (Bannister et al. 2019), FRB 181112 (Prochaska et al. 2019), FRB 190523 (Ravi et al. 2019), FRB 190102, FRB 190608, FRB 190611, and FRB 190711 (Macquart et al. 2020), have already be localized. With the measurements of nine FRB redshifts in hand, it is interesting to investigate what level of photon mass limits can be improved by taking advantage of the different redshift dependences of the dispersions from the plasma and photon mass. In this work, we develop a statistical approach to obtaining a combined limit on m_γ using the nine FRBs with redshift measurements. In addition, such an approach can also simultaneously lead to an estimate of the mean value of the dispersion measure (DM) contributed from the local host galaxy.

We should note that a recent work by Shao & Zhang (2017) extended previous studies to those FRBs without redshift measurement. They used an uninformative prior for the redshift to derive a combined constraint on m_γ from a catalog of 21 FRBs. A comparison with our results will be interesting, because the exact value of the redshift (or distance) plays an important role in constraining the photon mass. It is therefore useful to crosscheck, especially as using the real redshift measurements of FRBs, as we do in this work, avoids potential bias from the prior of redshift.

This paper is arranged as follows. In Section 2, we introduce the analysis method used to disentangle the dispersions from the plasma and photon mass and to constrain the photon mass. A combined limit on the photon mass from nine localized FRBs will be presented in Section 3, followed by our conclusions in Section 4. Throughout this paper, we adopt the flat ΛCDM model with the cosmological parameters recently derived from the latest Planck observations (Planck Collaboration et al. 2018): the Hubble constant H₀ = 67.36 km s⁻¹ Mpc⁻¹, the matter energy density Ω_m = 0.315, the vacuum energy density Ω_Λ = 0.685, and the baryonic matter energy density Ω_b = 0.0493.

2.1. Dispersion from a Nonzero Photon Mass

It is evident from Equation (2) that two massive photons with different frequencies (ν_l < ν_h ), which are emitted simultaneously from a same source at a redshift z, would reach us at different times. The arrival time delay induced by the nonzero photon mass effect can be expressed as

$$\begin{eqnarray}&&\Delta {t}_{{m}_{\gamma }}=\displaystyle \frac{1}{2{H}_{0}}{\left(\displaystyle \frac{{m}_{\gamma }{c}^{2}}{h}\right)}^{2}({\nu }_{l}^{-2}-{\nu }_{h}^{-2}){H}_{\gamma }(z)\,,\end{eqnarray} \tag{ 3 }$$

where H_γ (z) is a dimensionless redshift function,

$$\begin{eqnarray}&&{H}_{\gamma }(z)=\displaystyle {\int }_{0}^{z}\displaystyle \frac{{(1+{z}^{{\prime} })}^{-2}d{z}^{{\prime} }}{\sqrt{{\Omega }_{{\rm{m}}}{(1+{z}^{{\prime} })}^{3}+{\Omega }_{\Lambda }}}\,.\end{eqnarray} \tag{ 4 }$$

2.2. Dispersion from the Plasma Effect

Because of the dispersive nature of plasma, lower frequency radio photons would travel across the ionized median slower than higher frequency ones (Bentum et al. 2017). The arrival time difference between two photons with different frequencies (ν_l < ν_h ), which caused by the plasma effect, is described as

$$\begin{eqnarray}\begin{array}{lll}\Delta {t}_{{\rm{DM}}} & = & \displaystyle \int \displaystyle \frac{{\rm{d}}l}{c}\displaystyle \frac{{\nu }_{p}^{2}}{2}({\nu }_{l}^{-2}-{\nu }_{h}^{-2})\\ & = & \displaystyle \frac{{e}^{2}}{8{\pi }^{2}{m}_{e}{\epsilon }_{0}c}({\nu }_{l}^{-2}-{\nu }_{h}^{-2}){{\rm{DM}}}_{{\rm{astro}}},\end{array}\end{eqnarray} \tag{ 5 }$$

where ν_p = (n_ee²/4π² m_eε₀)^1/2 is the plasma frequency, n_e is the electron number density, m_e and e are the mass and charge of an electron, respectively, and ε₀ is the permittivity. Here DM_astro is defined as the integrated electron number density along the line of sight, i.e., DM_astro ≡ ∫ n_e dl. In a cosmological setting, DM_astro ≡ ∫ n_e,z (1 + z)⁻¹ dl, where n_e,z is the rest-frame number density of electrons and z is the redshift (Deng & Zhang 2014).

For a cosmic source, the DM_astro should include the following components:

$$\begin{eqnarray}\begin{array}{lll}{{\rm{DM}}}_{{\rm{astro}}} & = & {{\rm{DM}}}_{{\rm{MW}}}{{\rm{+DM}}}_{{\rm{MWhalo}}}\\ & & {{\rm{+DM}}}_{{\rm{IGM}}}{\rm{+}}\displaystyle \frac{{{\rm{DM}}}_{{\rm{host}}}}{1+z},\end{array}\end{eqnarray} \tag{ 6 }$$

where DM_MW, DM_MWhalo, DM_IGM, and DM_host represent the DM contributed by the Milky Way ionized interstellar medium, the Milky Way hot gaseous halo, the intergalactic medium (IGM), and the host galaxy (including the host galaxy interstellar medium and source environment), respectively. The cosmological redshift factor, 1 + z, converts the rest-frame DM_host to the measured value in the observer frame (Deng & Zhang 2014).

2.3. Methodology

As described earlier, the radio signals of FRBs are known to arrive with a frequency-dependent delay in time of the 1/ν/² behavior, which are expected from both the line-of-sight free electron content and mass effects on photon propagation (see Eqs. (3) and (5)). In our analysis, we suppose that the observed time delay is attributed to both the plasma and nonzero photon mass effects, i.e.,

$$\begin{eqnarray}&&\Delta {t}_{{\rm{obs}}}=\Delta {t}_{{\rm{DM}}}+\Delta {t}_{{m}_{\gamma }}\,.\end{eqnarray} \tag{ 7 }$$

In practice, the observed DM of an FRB, DM_obs, is directly obtained from the fitting of the ν⁻² form of the observed time delay. In our analysis, it thus equals to

$$\begin{eqnarray}&&{{\rm{DM}}}_{{\rm{obs}}}={{\rm{DM}}}_{{\rm{astro}}}+{{\rm{DM}}}_{\gamma }\,,\end{eqnarray} \tag{ 8 }$$

where DM_γ represents the "effective DM" arises from a nonzero photon mass (Shao & Zhang 2017),

$$\begin{eqnarray}&&{{\rm{DM}}}_{\gamma }\equiv \displaystyle \frac{4{\pi }^{2}{m}_{e}{\epsilon }_{0}{c}^{5}}{{h}^{2}{e}^{2}}\displaystyle \frac{{H}_{\gamma }(z)}{{H}_{0}}{m}_{\gamma }^{2}\,.\end{eqnarray} \tag{ 9 }$$

The measurements of DM_obs and their corresponding uncertainties σ_obs for all nine localized FRBs are presented in Table 1. These DM_obs values were determined by maximizing the peak signal-to-noise ratios ² . In order to identify DM_γ as radical as a massive photon effect from DM_obs, one needs to know the different DM contributions in Equation (6). For a well-localized FRB, the DM due to the Milky Way, DM_MW, can be well estimated from the Galactic electron density models. For FRB 181112, its DM_MW value is derived based on the Galactic electron density model of Prochaska et al. (2019). For the other eight FRBs, we adopt the NE2001 model for DM_MW (Cordes & Lazio 2002). Table 1 also contains the estimated DM_MW of nine localized FRBs, which are available in the FRB catalog (Petroff et al. 2016). The DM_MWhalo term is not well determined, but is estimated to be in the range of approximately 50 – 80 pc cm⁻³ (Prochaska & Zheng 2019). Hereafter we assume DM_MWhalo = 50 pc cm⁻³. However, DM_host is poorly known, which depends on the type of FRB host galaxy, the relative orientations of the host and source, and the near-source plasma (Xu & Han 2015). For a given galaxy type, the scale of DM_host along a certain line of sight is related to the size of the galaxy and electron density, and the average electron density is proportional to the root of the galaxy Hα luminosity (Luo et al. 2018). As the Hα luminosity scales with the star formation rate (SFR; Kennicutt et al. 1994; Madau et al. 1998), we model DM_host as a function of star formation history (Luo et al. 2018), ${{\rm{DM}}}_{{\rm{host}}}(z){{\rm{=DM}}}_{{\rm{host}},0}\sqrt{{\rm{SFR}}(z){\rm{/SFR}}(0)}$, where DM_host,0 is the present value of DM_host(z = 0) and ${\rm{SFR}}(z)=\displaystyle \frac{0.0157+0.118z}{1+{(z/3.23)}^{4.66}}\,{M}_{\odot }{{\rm{yr}}}^{{\rm{-1}}}{{\rm{Mpc}}}^{{\rm{-3}}}$ is the adopted star formation history (Hopkins & Beacom 2006; Li 2008). In our following likelihood estimation, DM_host,0 is treated as a free parameter. The IGM portion of DM is related to the redshift of the source and the fraction of ionized electrons in helium and hydrogen on the propagation path. If both helium and hydrogen are fully ionized (valid below z ∼ 3), it can be written as (Deng & Zhang 2014)

$$\begin{eqnarray}&&{{\rm{DM}}}_{{\rm{IGM}}}=\displaystyle \frac{21c{H}_{0}{\Omega }_{b}{f}_{{\rm{IGM}}}}{64\pi G{m}_{p}}{H}_{e}(z)\,,\end{eqnarray} \tag{ 10 }$$

where m_p is the mass of a proton, f_IGM ≃ 0.83 is the fraction of baryons in the IGM (Fukugita et al. 1998), and H_e (z) is the redshift-dependent dimensionless function,

$$\begin{eqnarray}&&{H}_{e}(z)=\displaystyle {\int }_{0}^{z}\displaystyle \frac{(1+{z}^{{\prime} })d{z}^{{\prime} }}{\sqrt{{\Omega }_{{\rm{m}}}{(1+{z}^{{\prime} })}^{3}+{\Omega }_{\Lambda }}}\,.\end{eqnarray} \tag{ 11 }$$

In Figure 1, we illustrate the dependence of functions H_e (z) and H_γ (z) on the redshift z. It is obvious that the DM contributions from the IGM and a possible photon mass have different redshift dependences. Bonetti et al. (2016, 2017) pointed out that the different behavior of the two redshift functions can be used to break parameter degeneracy and to improve the sensitivity to the photon mass when a few FRB redshifts are available (see also Bentum et al. 2017; Shao & Zhang 2017). For now, nine FRB redshifts have been measured, which might make this point come true.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1. Redshifts and Dispersion Measures of Nine Localized FRBs

FRB Name	z	DMobs (pc cm−3)	DMMW (pc cm−3)	References
FRB 121102	0.19273	558.1 ± 3.3	188	[1], [2], [3]
FRB 180916.J0158+65	0.0337	348.76 ± 0.1	199	[4]
FRB 180924	0.3214	361.42 ± 0.06	40.5	[5]
FRB 181112	0.4755	589.27 ± 0.03	102	[6]
FRB 190102	0.291	363.6 ± 0.3	57.3	[7]
FRB 190523	0.66	760.8 ± 0.6	37	[8]
FRB 190608	0.1178	338.7 ± 0.5	37.2	[7]
FRB 190611	0.378	321.4 ± 0.2	57.83	[7]
FRB 190711	0.522	593.1 ± 0.4	56.4	[7]

Notes: The References are ^[1]Spitler et al. (2016); ^[2]Chatterjee et al. (2017); ^[3]Tendulkar et al. (2017); ^[4]Marcote et al. (2020); ^[5]Bannister et al. (2019); ^[6]Prochaska et al. (2019); ^[7]Macquart et al. (2020); ^[8]Ravi et al. (2019).

For a set of nine localized FRBs, a combined limit on the photon mass m_γ can be obtained by maximizing the likelihood function:

$$\begin{eqnarray}\begin{array}{lll} {\mathcal L} & = & \displaystyle \prod _{i}\displaystyle \frac{1}{\sqrt{2\pi }{\sigma }_{{\rm{tot}},i}}\\ & & \times \exp \left\{-\displaystyle \frac{{[{{\rm{DM}}}_{{\rm{obs}},i}{{\rm{-DM}}}_{{\rm{astro}},i}{{\rm{-DM}}}_{\gamma }({m}_{\gamma },{z}_{i})]}^{2}}{2{\sigma }_{{\rm{tot}},i}^{2}}\right\},\end{array}\end{eqnarray} \tag{ 12 }$$

where ${\sigma }_{{\rm{tot}}}^{2}$ is the total variance of each FRB,

$$\begin{eqnarray}&&{\sigma }_{{\rm{tot}}}^{2}={\sigma }_{{\rm{obs}}}^{2}+{\sigma }_{{\rm{MW}}}^{2}+{\sigma }_{{\rm{MWhalo}}}^{2}+{\sigma }_{{\rm{int}}}^{2}\,,\end{eqnarray} \tag{ 13 }$$

where σ_obs, σ_MW, and σ_MWhalo are the uncertainties of DM_obs, DM_MW, and DM_MWhalo, respectively, and σ_int is the global intrinsic scatter. The ATNF pulsar catalog shows that the average uncertainty of DM_MW for sources at high Galactic latitude (|b| > 10°) is about 10 pc cm⁻³ (Manchester et al. 2005), and we take this value as σ_MW. Recent observations imply that the Galactic halo contributes between 50 – 80 pc cm⁻³ (Prochaska & Zheng 2019). Here we adopt the possible range of DM_MWhalo as its uncertainty, i.e., σ_MWhalo = 30 pc cm⁻³. To be conservative, we introduce an additional free parameter σ_int to account for the intrinsic scatter that might originate from the large IGM fluctuation and the diversity of host galaxy contribution. In this case, the free parameters are DM_host,0, σ_int, and the photon mass m_γ .

We use the Python Markov Chain Monte Carlo module, EMCEE (Foreman-Mackey et al. 2013), to obtain the posterior distributions of the three parameters (m_γ , DM_host,0, and σ_int) from nine FRBs with redshift measurements. The marginalized accumulative probability distribution of the photon mass m_γ is displayed in Figure 2. The 68 % and 95 % confidence-level upper limits on m_γ are

$$\begin{eqnarray}&&{m}_{\gamma }\le 7.1\times {10}^{-51}\,{\rm{kg}}\simeq 4.0\times {10}^{-15}\,{\rm{eV}}/{c}^{2}\end{eqnarray} \tag{ 14 }$$

and

$$\begin{eqnarray}&&{m}_{\gamma }\le 1.3\times {10}^{-50}\,{\rm{kg}}\simeq 7.3\times {10}^{-15}\,{\rm{eV}}/{c}^{2}\,,\end{eqnarray} \tag{ 15 }$$

respectively. These results are comparable with or represent sensitivities improved by a factor of 7-fold over existing photon mass limits from the single FRBs (Wu et al. 2016; Bonetti et al. 2016, 2017; Xing et al. 2019).

In Figure 3, we show the 1D and 2D marginalized posterior probability distributions with 1 – 3σ confidence regions for the parameters. We find that the mean value of DM contributed from the local host galaxy is estimated to be ${{\rm{DM}}}_{{\rm{host}},0}{\rm{=27}}.{0}_{{\rm{-17}}.6}^{{\rm{+23}}.3}\,{{\rm{pccm}}}^{{\rm{-3}}}$ at 68 % confidence level, which is well consistent with those inferred from observations of localized FRBs (Bannister et al. 2019; Ravi et al. 2019). Meanwhile, an estimation for the global intrinsic scatter arises from both the large IGM fluctuation and the diversity of host galaxy contribution, ${\sigma }_{{\rm{int}}}={92.5}_{-23.7}^{+36.4}\,{{\rm{pccm}}}^{{\rm{-3}}}$, is obtained.

Zoom In Zoom Out Reset image size

The rest mass of the photon, m_γ , can be effectively constrained by measuring the frequency-dependent time delays of radio waves from distant astrophysical sources. In particular, extragalactic FRBs are well suited for these studies because they are short-duration radio transients involving long cosmological propagation distances. However, there is a similar frequency-dependent dispersion expected from both the plasma and possible photon mass effects on photon propagation. A key challenge in this method of constraining m_γ , therefore, is to distinguish the dispersions from the plasma and photon mass. Most previous limits are based on the DM of a single FRB, in which the degeneracy of these two dispersions cannot be well handled (Wu et al. 2016; Bonetti et al. 2016, 2017). In view of their different redshift dependences (see Eqs. (4) and (11)), it has been suggested that the DM contributions from the IGM and photon mass could in principle be disentangled by more FRBs with different redshift measurements, enabling the sensitivity to m_γ to be improved (Bonetti et al. 2016, 2017; Bentum et al. 2017).

Recently, nine FRBs with redshift measurements have been reported. These provide a good opportunity not only to disentangle the degeneracy problem but also to place more stringent limits on the photon mass. In this work, we suppose that the observed DM of an FRB, DM_obs, is attributed to both the plasma and nonzero photon mass effects. Using the measurements of DM_obs and z of the nine FRBs, we place a combined limit on the photon mass at 68% (95%) confidence level, i.e., m_γ ≤ 7.1 × 10⁻⁵¹ kg, or equivalently m_γ ≤ 4.0 × 10⁻¹⁵ eV/c² (m_γ ≤ 1.3 × 10⁻⁵⁰ kg, or equivalently m_γ ≤ 7.3 × 10⁻¹⁵ eV/c²), which is as good as or represents a factor of 7 improvement over the results previously obtained from the single FRBs. Moreover, a reasonable estimation for the mean value of DM contributing from host galaxies is simultaneously achieved; i.e., ${{\rm{DM}}}_{{\rm{host}},0}{\rm{=27}}.{0}_{{\rm{-17}}.6}^{{\rm{+23}}.3}\,{{\rm{pccm}}}^{{\rm{-3}}}$. Previously, Shao & Zhang (2017) derived a combined limit of m_γ ≤ 8.7 × 10⁻⁵¹ kg at 68 % confidence level from a catalog of 21 FRBs (20 of them without redshift measurement), where uninformative prior is made to the unknown redshift. With another prior of z that traces the SFR, Shao & Zhang (2017) obtained a combined limit of 7.5 × 10⁻⁵¹ kg at 68 % confidence level, showing that changes in the prior of redshift can lead to slight difference. The precision of our constraint from nine localized FRBs is comparable to these, while avoiding potential bias from the prior of redshift.

More FRBs with redshift measurements are essential for using this method presented here to constrain the photon rest mass. It is encouraging that the Deep Synoptic Array 10-dish prototype (DSA-10) will gradually be replaced with the ongoing 110-dish DSA in 2021 and the ultimate 2000-dish DSA in 2026, which are expected to detect and localize over 100 and 10⁴ FRBs per year (Hallinan et al. 2019; Kocz et al. 2019). With the rapid progress in localizing FRBs, the photon mass limits will be further improved.

We are grateful to the anonymous referee for helpful comments. This work is partially supported by the National Natural Science Foundation of China (Grant Nos. 11673068, 11725314 and U1831122), the Youth Innovation Promotion Association (2017366), the Key Research Program of Frontier Sciences (Grant Nos. QYZDB-SSW-SYS005 and ZDBS-LY-7014), and the Strategic Priority Research Program "Multi-waveband gravitational wave universe" (Grant No. XDB23000000) of Chinese Academy of Sciences.