Probing cosmic anisotropy with GW/FRB as upgraded standard sirens

TL;DR

This paper demonstrates that combined GW/FRB observations as upgraded standard sirens significantly improve the sensitivity to cosmic anisotropy testing compared to using GW data alone, both theoretically and through simulations.

Contribution

It introduces a novel approach of using combined GW/FRB data as upgraded standard sirens to enhance tests of cosmic anisotropy.

Findings

Combined GW/FRB measurements double the sensitivity to cosmic anisotropy.

Theoretical analysis confirms improved anisotropy constraints.

Simulations support the effectiveness of the combined method.

Abstract

Recently it was shown that cosmic anisotropy can be well tested using either standard siren measurement of luminosity distance $d_\mathrm{L}(z)$ from gravitational-wave (GW) observation or dispersion measure ($\mathrm{DM}(z)$) from fast radio burst (FRB). It was also observed that the combined measurement of $d_\mathrm{L}(z)\cdot\mathrm{DM}(z)$ from the GW/FRB association system as suggested in some of FRB models is more effective to constrain cosmological parameters than $d_\mathrm{L}(z)$ or $\mathrm{DM}(z)$ separately due to its independence from Hubble constant. In this paper, we will show both theoretically and with simulation that, this upgraded sirens from combined GW/FRB observations could test cosmic anisotropy with a double relative sensitivity compared to the usual standard siren from GW observation alone.

Equations

$$d_{\mathrm{L}}(z)=\frac{1+z}{H_0}{\int }_0^z\frac{d{z}^{\prime }}{E(z^{\prime })},$$

$$\mathcal{A}=$$

$$\mathcal{H}(f)=\mathcal{A}{f}^{-7/6}{e}^{i\Psi (f)},$$

$$h(t)=F_{+}(\theta ,\varphi ,\psi )h_{+}(t)+F_{\times }(\theta ,\varphi ,\psi )h_{\times }(t).$$

$$\mathrm{DM}=\int \frac{n_{e,z}\mathrm{d}l}{1+z},$$

$$\Delta t=\frac{e^2}{2\pi m_{e}c}\left(\frac{1}{{\nu }_1^2}-\frac{1}{{\nu }_2^2}\right)\mathrm{DM}.$$

$${\mathrm{DM}}_{\mathrm{obs}}={\mathrm{DM}}_{\mathrm{MW}}+{\mathrm{DM}}_{\mathrm{IGM}}+\frac{{\mathrm{DM}}_{\mathrm{HG}}}{1+z},$$

$$⟨{\mathrm{DM}}_{\mathrm{IGM}}(z)⟩=\frac{3H_0{\Omega }_b{f}_{\mathrm{IGM}}}{8\pi m_p}{\int }_0^z\frac{\chi (z^{\prime })(1+z^{\prime })}{E(z^{\prime })}d{z}^{\prime },$$

$$\displaystyle\chi_{\rm e,He}(z)=\left\{\begin{array}[]{lcl}1,&&{z<3};\\ 0.025z^{3}-0.244z^{2}+0.513z+1.006,&&{z>3}.\end{array}\right.$$

$$d_{\mathrm{L}}^{\mathrm{fid}}(\hat{z})=d_{\mathrm{L}}^{\mathrm{iso}}(z)[1+g(\hat{n}\cdot \hat{z})],$$

$${\mathrm{DM}}_{\mathrm{IGM}}^{\mathrm{fid}}(\hat{z})={\mathrm{DM}}_{\mathrm{IGM}}^{\mathrm{iso}}(z)[1+g(\hat{n}\cdot \hat{z})].$$

$$\hat{n}=(\cos \varphi \sin \theta ,\sin \varphi \sin \theta ,\cos \theta ),$$

$$d_{\mathrm{L}}^{\mathrm{fid}}(\hat{z})\cdot {\mathrm{DM}}_{\mathrm{IGM}}^{\mathrm{fid}}(\hat{z})=d_{\mathrm{L}}^{\mathrm{iso}}(z)\cdot {\mathrm{DM}}_{\mathrm{IGM}}^{\mathrm{iso}}(z)[1+g(\hat{n}\cdot \hat{z}){]}^2.$$

$$\frac{1}{d_{\mathrm{L}}^{\mathrm{fid}}(\hat{z})\cdot {\mathrm{DM}}_{\mathrm{IGM}}^{\mathrm{fid}}(\hat{z})}\frac{\mathrm{d}[d_{\mathrm{L}}^{\mathrm{fid}}(\hat{z})\cdot {\mathrm{DM}}_{\mathrm{IGM}}^{\mathrm{fid}}(\hat{z})]}{\mathrm{d}g}$$

$$=$$

$$F_{+}^{(1)}(\theta ,\varphi ,\psi )=$$

$$F_{\times }^{(1)}(\theta ,\varphi ,\psi )=$$

$$\displaystyle\sigma_{d_{\rm L}}^{\rm inst}\simeq\sqrt{\bigg{\langle}\frac{\partial\mathcal{H}}{\partial d_{\rm L}},\frac{\partial\mathcal{H}}{\partial d_{\rm L}}\bigg{\rangle}^{-1}}.$$

$${\sigma }_{d_{\mathrm{L}}}^{\mathrm{inst}}\simeq \frac{2d_{\mathrm{L}}}{\rho }$$

$${\sigma }_{d_{\mathrm{L}}}$$

$$\displaystyle=\sqrt{\bigg{(}\frac{2d_{\rm L}}{\rho}\bigg{)}^{2}+(0.05zd_{\rm L})^{2}}.$$

$${\sigma }_{d_{\mathrm{L}}\cdot {\mathrm{DM}}_{\mathrm{IGM}}}=\sqrt{({\mathrm{DM}}_{\mathrm{IGM}}\cdot {\sigma }_{d_{\mathrm{L}}}{)}^2+(d_{\mathrm{L}}\cdot {\sigma }_{{\mathrm{DM}}_{\mathrm{IGM}}}{)}^2}.$$

$$\rho =\sqrt{\sum _{i=1}^N⟨{\mathcal{H}}^{(i)},{\mathcal{H}}^{(i)}⟩},$$

$$⟨a,b⟩\equiv 4{\int }_{f_{\min }}^{f_{\max }}\frac{\tilde{a}(f){\tilde{b}}^{\ast }(f)+{\tilde{a}}^{\ast }(f)\tilde{b}(f)}{2}\frac{df}{S_h(f)}$$

$$S_h(f)=$$

$$h=0.678,{\Omega }_{\mathrm{m}}=0.308,{\Omega }_{\mathrm{b}}=0.049$$

$$P(z)\propto \frac{4\pi d_{\mathrm{c}}^2(z)R(z)}{(1+z)E(z)},$$

$$\displaystyle R(z)=\left\{\begin{array}[]{lcl}1+2z,&&{z\leq 1};\\ \frac{3}{4}(5-z),&&{1<z<5};\\ 0,&&{z\geq 5}.\end{array}\right.$$

$$\displaystyle\chi^{2}=\sum_{i=1}^{N}\bigg{[}\frac{d_{\rm L}^{i}\cdot\mathrm{DM}_{\rm IGM}^{i}-d_{\rm L}^{\rm fid}(\hat{z})\cdot\mathrm{DM}^{\rm fid}_{\rm IGM}(\hat{z})}{\sigma_{d_{\rm L}\cdot\mathrm{DM}_{\rm IGM}}^{i}}\bigg{]}^{2},$$

Full text

1 Introduction

The power of the standard siren [1, 2] starts to be appreciated since the first discovery of binary neutron star (BNS) merger [3, 4]. The luminosity distance can be directly determined from the gravitational waveform of coalescing binaries and the redshift can be inferred from the associated electromagnetic (EM) counterpart. In a previous work [5], we have demonstrated the promiseful perspective of constraint ability on cosmic anisotropy using gravitational wave (GW) as standard siren. We found that with a few hundreds standard siren events, the cosmic isotropy can be ruled out at $3\sigma$ confidence level (C.L.) and the dipole direction can be constrained within $10\sim 20\%$ at 3$\sigma$ C.L. if the dipole amplitude is of order $\mathcal{O}(10^{-2})$. See [6, 7, 8, 9, 10] for more recent studies on the standard siren and its application in cosmology and [11] for a review and references therein.

On the other hand, it was indicated in some models [12, 13] of fast radio burst (FRB) [14] that, the observed millisecond-duration radio transients could originate from the mergers of binary neutron stars (BNS) [15, 16, 17], black hole-neutron star (BHNS) [18], double charged black holes (BHs) [19] or charged and rotating BHs [20] (See also [21, 22, 23] for the blitzar model of FRB from a single collapsing magnetised neutron star). A typical FRB could exhibit high dispersion measure (DM) of order $\mathcal{O}(10^{2}\sim 10^{3})\,\mathrm{pc}\cdot\mathrm{cm}^{-3}$ with negligible uncertainties of order $\mathcal{O}(10^{-2}\sim 10^{-1})\,\mathrm{pc}\cdot\mathrm{cm}^{-3}$ [12]. Recently, the cosmic anisotropy was tested in [24] using dispersion measure for some simulated FRB events with known redshifts. With a few hundreds FRB events, the cosmic anisotropy can be found if the dipole amplitude is of order $\mathcal{O}(10^{-2})$. Other applications of FRBs on cosmology can be found in [25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36].

In particular, an ungraded standard siren was proposed in [31] when the binary system exhibits both GW and FRB signals. Combining the luminosity distance $d_{\rm L}(z)$ from GW and dispersion measures $\mathrm{DM}_{\mathrm{IGM}}(z)$ from FRB through intergalactic medium (IGM), the authors found that it could constrain the cosmological parameters more strongly than that using separately either $d_{\rm L}(z)$ or $\mathrm{DM}_{\mathrm{IGM}}(z)$. It was also discussed that the FRB event rate is around $1.4\times 10^{3}\rm Gpc^{-3}yr^{-1}$ [19], meanwhile the BNS merger rate is estimated to be $1540_{-1220}^{+3200}\rm Gpc^{-3}yr^{-1}$[3], which are at the same order of magnitude. It is worth noting that, the identified system with GW/FRB association is nevertheless idealistic, but still possible for the next generation of GW detectors, and there could be at least a few hundreds such events detectable per year [31].

In the present work, we will use this upgraded siren to test the constraint ability on the cosmic anisotropy, assuming that GWs and FRBs have the same progenitor system like BNS or BHNS. The paper is organized as follows. In Sec. 2, the upgraded siren with GW/FRB association is reviewed and compared to the usual standard siren when the dipole modulation is concerned. In Sec. 3, we present the error estimation of the jointed measurements of luminosity distance and dispersion measure. In Sec. 4, the GW/FRB standard siren events are simulated, and the presumed dipole modulation is recovered by adopting the Markov Chain Monte Carlo (MCMC) approach. In Sec. 5, the results are summarized in conclusions. Throughout this paper, a flat universe is assumed for simplicity, and the geometric unit $c=G=1$ is adopted.

2 GW/FRB standard siren

In this section, we will briefly summarize some key points of the GW/FRB standard siren, more details can be found in [5] and [31]. A sensitivity comparison to GW standard siren is given at the end of this section.

For a GW standard siren event at redshift $z$, the isotropic luminosity distance in a spatially flat FLRW universe

[TABLE]

can be directly inferred from the amplitude

[TABLE]

of Fourier transform of GW signal

[TABLE]

whose detector response in the transverse-traceless gauge is

[TABLE]

Here $E(z)\equiv H(z)/H_{0}=\sqrt{(1-\Omega_{m})+\Omega_{m}(1+z)^{3}}$ is dimensionless Hubble parameter, $H_{0}=100h\,\rm km\,s^{-1}\,Mpc^{-1}$ is the Hubble constant, $\Omega_{m}$ is the matter density parameter today, $\iota$ is the inclination angle between the orbit of the GW source and the line-of-sight of the observer, $\mathcal{M}_{c}$ is the observed chirp mass $\mathcal{M}_{c}=(1+z)M\eta^{3/5}$ with total mass $M=m_{1}+m_{2}$ and symmetric mass ratio $\eta=m_{1}m_{2}/M^{2}$. The phase in Eq. (2.3) is computed in the post-Newtonian formalism up to 3.5 PN and the specific expression can be found in [37]. $F_{+,\times}$ are the antenna pattern functions for the two polarizations $h_{+,\times}$, ($\theta,\phi$) are the directional angles of the source in the detector frame, and $\psi$ is the polarization angle. Note that, one can take $\iota=0^{\circ}$ as argued in [5] and references therein.

On the other hand, the dispersion measure is defined as the observed column density of the free electrons along the line-of-sight [38, 39, 25, 26],

[TABLE]

which can be measured from the observed time delay between two frequencies of emitted electromagnetic signal [38, 39, 25, 26],

[TABLE]

When applied to FRB, the observed DM consists of three parts:

[TABLE]

which come from the Milky Way (MW), intergalactic medium (IGM), and the FRB host galaxy (HG), respectively. Once $\mathrm{DM}_{\rm obs}$, $\mathrm{DM}_{\rm MW}$ [40, 41, 42] and $\mathrm{DM}_{\rm HG}$ are determined from observations, one can obtain the observed $\mathrm{DM}_{\rm IGM}$, which contributes most among all with most relevance to cosmological study. The mean value of $\mathrm{DM}_{\rm IGM}$ can be written as [38, 39, 25, 26]

[TABLE]

where the fraction of baryon mass in the IGM $f_{\rm IGM}\simeq 0.83$, $\Omega_{b}$ is the current baryon mass fraction of the universe , and $m_{p}$ is the mass of proton. $\chi(z)=\frac{3}{4}y_{1}\chi_{\rm e,H}(z)+\frac{1}{8}y_{2}\chi_{\rm e,He}(z)$, with $y_{1}\sim 1,y_{2}\sim 1$. $\chi_{\rm e,H}(z)$ and $\chi_{\rm e,He}(z)$ are the ionization fractions for H and He, respectively. Since H is essentially fully ionized at $z<6$, we could take $\chi_{\rm e,H}(z)=1$. For He, however, it is fully ionized at $z<3$, so we have the approximate expression for $\chi_{\rm e,He}(z)$, which is given by [27]

[TABLE]

In the rest part of the paper, we will use the notation $\mathrm{DM}_{\rm IGM}$ instead of $\langle\mathrm{DM}_{\rm IGM}\rangle$. It was observed in [31] that the constraint ability of $d_{\mathrm{L}}$ or $\mathrm{DM}_{\rm IGM}$ method alone on the cosmological parameters is limited due to the fact that they both depend on the Hubble constant $H_{0}$, while $d_{\mathrm{L}}\cdot\mathrm{DM}_{\rm IGM}$ is independent of $H_{0}$, as can be seen from Eqs. (2.1) and (2.8). Note the fact that the local measurement of $H_{0}$ is different from that given by the Planck with 4.4 $\sigma$ confidence, which is the so-called “Hubble Tension” [43, 44, 45]. Thus the inclusion of $\mathrm{DM}_{\rm IGM}$ to $d_{\mathrm{L}}$ can avoid this problem and exhibit better performance.

One of ways to characterize the cosmic anisotropy is to introduce a dipole modulation with amplitude $g$ and direction $\hat{n}$ on the isotropic luminosity distance,

[TABLE]

where the isotropic luminosity distance $d_{\mathrm{L}}^{\rm iso}(z)$ is calculated according to Eq. (2.1). Since the definition of dispersion measure also contains the distance element $\mathrm{d}l$ along the line-of-sight, it can be also modified in a similar way as [24]

[TABLE]

where the isotropic dispersion measure $\mathrm{DM}_{\mathrm{IGM}}^{\rm iso}(z)$ is computed according to (2.8). The dipole direction is given by

[TABLE]

where $\theta\in[0,\pi)$ and $\phi\in[0,2\pi)$. Therefore, the fiducial combination takes the form of

[TABLE]

To see that the upgraded siren from GW/FRB performs better than the standard siren from GW alone to test the cosmic anisotropy, one can calculate their relative sensitivity, namely the relative change of measured quantity with respect to the change of the dipole amplitude,

[TABLE]

The relative sensitivity of the upgraded siren from GW/FRB is twice as large as that of the standard siren from GW alone. Thus we can conclude that $d_{\rm L}\cdot\mathrm{DM}_{\rm IGM}$ is more sensitive than $d_{\rm L}$ with respect to the presumed anisotropy.

3 Error estimation

ET is a proposed third-generation ground-based GW detector, and the frequency it will cover ranges from 1 to $10^{4}$ Hz. The exact forms of pattern functions for ET are given by [46]

[TABLE]

and the rest of the pattern functions are $F_{+,\times}^{(2)}(\theta,\phi,\psi)=F_{+,\times}^{(1)}(\theta,\phi+2\pi/3,\psi)$ and $F_{+,\times}^{(3)}(\theta,\phi,\psi)=F_{+,\times}^{(1)}(\theta,\phi+4\pi/3,\psi)$, respectively, since the three interferometers align with an angle $60^{\circ}$ with each other. The error estimations are presented below.

Assuming that the error on the luminosity distance is uncorrelated with errors on the remaining GW parameters, the instrumental error can be estimated with Fisher matrix by

[TABLE]

Since $\mathcal{H}\propto d_{\rm L}^{-1}$, one could approximate $\sigma_{d_{\rm L}}^{\rm inst}\simeq d_{\rm L}/\rho$. Moreover, a factor of 2 is added to the instrumental error

[TABLE]

for a conservative estimation [47] when the correlation between $d_{\rm L}$ and $\iota$ is taken into account. Another error to be considered is $\sigma_{d_{\rm L}}^{\rm lens}$ due to the effect of weak lensing, and we assume $\sigma_{d_{\rm L}}^{\rm lens}/d_{\rm L}=0.05z$ as that in [46]. We therefore take the total uncertainty on the luminosity distance as

[TABLE]

Thus, the corresponding error for the product $d_{\mathrm{L}}\cdot\mathrm{DM}_{\rm IGM}$ is given by [31]

[TABLE]

In the following section, we will use $d_{\mathrm{L}}\cdot\mathrm{DM}_{\rm IGM}$ and Eq. (3.5) to simulate the corresponding measurement through Gaussian distribution.

4 Simulation test

Next we simulate a set of GW/FRB events. A GW event is usually claimed when the signal-to-noise ratio (SNR) of the detector network reaches above 8. The combined SNR for the network of $N$ ($N=3$ for ET) independent interferometers is given by

[TABLE]

where the scalar product is defined as

[TABLE]

for given $\tilde{a}(f)$ and $\tilde{b}(f)$ as the Fourier transforms of some functions $a(t)$ and $b(t)$, and $S_{h}(f)$ denotes the one-side noise power spectral density (PSD) characterizing the performance of GW detector. The noise PSD for ET is [5, 46].

[TABLE]

where $x=f/f_{p}$ with $f_{p}\equiv 100\rm\,Hz$. The lower and upper cutoff frequencies $f_{\rm min}=1\,\rm Hz$ and $f_{\rm max}=2f_{\rm LSO}$, where the orbit frequency at the last stable orbit $f_{\rm LSO}=1/6^{3/2}2\pi\mathcal{M}_{\rm obs}$ with the observed total mass $M_{\rm obs}=(1+z)M$.

We take the standard $\Lambda$CDM model as our isotropic cosmological model with fiducial values

[TABLE]

from the Planck 2015 data [48]. The redshift distribution of the observable sources takes the form of [46]

[TABLE]

where $d_{\rm c}$ is the comoving distance defined as $d_{\rm c}(z)\equiv\int_{0}^{z}1/E(z^{\prime})dz^{\prime}$, and $R(z)$ describes the NS-NS merger rate, which is given by [49]

[TABLE]

In order to recover the presumed anisotropy with MCMC method, we calculate $\chi^{2}$ as

[TABLE]

where $N$ denotes the number of data sets, $d_{\rm L}^{i}\cdot\mathrm{DM}_{\rm IGM}^{i},\ \sigma_{d_{\rm L}\cdot\mathrm{DM}_{\rm IGM}}^{i}$ are the $i$-th combination of luminosity distance and dispersion measure, and the corresponding error of the simulated data, respectively.

The specific steps to simulate the mock data of standard siren events are similar with those in [5, 6], which are listed as follows.

Simulate $N$ GW/FRB datasets, according to the redshift distribution in Eq. (4.5), and the angles $\theta$ and $\phi$ are randomly sampled within the intervals $[0,\pi)$ and $[0,2\pi)$. 2. 2.

Calculate the isotropic luminosity distance $d_{\mathrm{L}}^{\rm iso}(z)$ according to Eq. (2.1), isotropic dispersion measure $\mathrm{DM}_{\mathrm{IGM}}^{\rm iso}(z)$ from (2.8) and anisotropic fiducial combined $d_{\mathrm{L}}^{\mathrm{fid}}(\hat{z})\cdot\mathrm{DM}_{\rm IGM}^{\rm fid}(\hat{z})$ according to Eq. (2.15). 3. 3.

Randomly sample the mass of neutron star and black hole within $[1,2]$ $M_{\odot}$ and $[3,10]$ $M_{\odot}$, respectively. Although the current black holes detected by LIGO/Virgo are commonly heavier than $10M_{\odot}$ (e.g., $30M_{\odot}$ in GW 150914 [50]), we still assume a conservative mass distribution for solar-mass black holes given by [51], as done by many other authors [6, 46, 47]222In fact we also did a similar analysis by assuming the mass distribution of black holes in $[10,100]$ $M_{\odot}$, we found that the results are qualitatively similar and do not change much. . Another reason to assume the black holes in region $[3,10]$ $M_{\odot}$ is that the FRBs originate from the mergers of binary neutron stars (BNS) or black hole-neutron stars (BHNS). Assuming the ratio between the number of BHNS and that of BNS is roughly 0.03 [52]. Evaluate the SNR and the error $\sigma_{d_{L}\cdot\mathrm{DM}_{\rm IGM}}$. 4. 4.

Simulate $d_{\rm L}^{\rm mea}\cdot\mathrm{DM}_{\rm IGM}^{\rm mea}$ with Gaussian distribution $\mathcal{N}(d_{\rm L}^{\rm fid}\cdot\mathrm{DM}_{\rm IGM}^{\rm fid},\sigma_{d_{\rm L}\cdot\mathrm{DM}_{\rm IGM}})$ and calculate the $\chi^{2}$ in Eq (4.10). Note that here we do not sample the measurement of $d_{\rm L}^{\rm mea}$ (or $\mathrm{DM}_{\rm IGM}^{\rm mea})$ from $d_{\rm L}^{\rm mea}=\mathcal{N}(d_{\rm L}^{\rm fid},\sigma_{d_{\rm L}})$ (or $\mathrm{DM}_{\rm IGM}^{\rm mea}=\mathcal{N}(\mathrm{DM}_{\rm IGM}^{\rm fid},\sigma_{\mathrm{DM}_{\rm IGM}})$) separately, as done in [31]. 5. 5.

Apply the MCMC method to calculate the likelihood function of ($g,\theta,\phi$) and find out the constrained dipole modulation $(g^{c},\theta^{c},\phi^{c})$, which will be compared with the presumed fiducial dipole modulation $(g^{f},\theta^{f},\phi^{f})$.

5 Results and discussions

A representative result is illustrated in Fig. 1 to recover the presumed fiducial value of dipole amplitude $g^{f}=0.1$ (left column) and the dipole direction angles (medium and right columns). The top-row panels are obtained using GW siren, while the medium-row panels are obtained using GW/FRB siren, of which the error bars are compared in the bottom-row panels with GW siren (dashed) and GW/FRB siren (solid), respectively. As expected theoretically from (2), the error bars are all reduced when DM information is included, indicating a better performance of the GW/FRB siren. From Fig. 1, we can clearly see that GW sirens are not able to constrain dipole amplitude very well, while the inclusion of FRB could improve the constraint ability significantly. As for dipole directions, GW sirens have already constrained them quite well, so that GW/FRB sirens do not show considerable improvements, especially when increasing $N$ to a larger value.

To summary, we use the combination of GWs and FRBs as the upgraded sirens to investigate the constraint ability on the anisotropy of the universe expansion. This combination takes the advantage of being independent of $H_{0}$, which reduces the corresponding errors. With the presumed dipole anisotropy, we construct the simulated data of GW events from BNS and BHNS for ET, assuming GWs and FRBs share the same progenitor system. Compared with the result presented in our previous work [5], we found a significant improvement on recovering the dipole amplitude. Although there do exist some challenges for this upgraded siren, as discussed in [31], we do believe that it can become a quite promising approach for studying the constraint ability on the cosmic anisotropy when more and more BNS or BHNS events are detected in the future.

Acknowledgments

We want to thank Xue-Wen Liu for helpful discussions and Jun-Jie Wei, He Gao for useful correspondence. This work is supported by the National Natural Science Foundation of China Grants No. 11690022, No. 11435006 and No. 11647601, and by the Strategic Priority Research Program of CAS Grant No. XDB23030100 and by the Key Research Program of Frontier Sciences of CAS. SJW is supported by the postdoctoral scholarship of Tufts University.