New Limits on the Lorentz/CPT Symmetry Through 50 Gravitational-wave Events

While the discovery of gravitational waves (GWs; Abbott et al. 2016) provides strong evidence in support of general relativity (GR), it also provides new approaches for testing GR. The classical GR is the most successful theory of gravity to date, but it is not compatible with quantum field theory, which describes the other three fundamental interactions very well. In some candidate theories of quantum gravity, like string theory, the Lorentz invariance (LI), one of the foundations of modern physics, could spontaneously break (Kostelecký & Samuel 1989; Kostelecký & Potting 1991). LI violation (LIV) could happen in many different sectors, such as the electron sector, the photon sector, and the gravity sector. Experimental bounds on these sectors can be found in Kostelecký & Russell (2011).

The standard-model extension (SME; Colladay & Kostelecký 1997; Kostelecký 2004) is a powerful and popular framework to explore LIV. SME is based on a theoretical effective-field approach and contains all possible operators for Lorentz and CPT violations in GR as well as the standard model of particle physics. In the SME, an LIV term in the Lagrangian contains an LIV operator constructed from the contraction of a conventional tensor operator with a violating tensor coefficient. In the framework of SME, the symmetry breaking is spontaneous (Bluhm & Kostelecký 2005) and the violations are observer LI.

In the past years, a number of experiments have been performed to investigate the gravitational LIV (gLIV) in SME (for a review, see, e.g., Hees et al. 2016). Their constraints come from lunar laser ranging (Battat et al. 2007; Bourgoin et al. 2017), atom interferometry (Mueller et al. 2008), planetary ephemerides (Iorio 2012; Hees et al. 2015), Čerenkov radiation of cosmic rays (Kostelecký & Tasson 2015), pulsar timing (Shao 2014a, 2014b, 2019; Jennings et al. 2015; Shao & Wex 2016; Shao & Bailey 2018; Shao 2019), short-range gravity experiments (Bailey et al. 2015; Shao et al. 2016), Very Long Baseline Interferometry (Le Poncin-Lafitte et al. 2016), GWs (Kostelecký & Mewes 2016; Schreck 2017; Shao 2020), and so on.

As shown by Kostelecký & Mewes (2016) and Abbott et al. (2017), GWs provide an opportunity for tests of LIV in the pure gravity sector. In the SME, each LIV term can be cataloged with a specific mass dimension d (Kostelecký & Mewes 2009; Kostelecký & Li 2019). Terms with a higher d are considered to appear at higher energy scales, which means a higher order correction. Similar to the situation of the photon sector (Kostelecký & Mewes 2002, 2009), the anisotropic birefringence phenomenon can be used to restrict the Lorentz/CPT violations tightly in nonminimal gravity with the mass dimension of LIV operators d ≥ 5. At the time of Kostelecký & Mewes (2016), there was only one GW event, GW150914, so at a specific mass dimension d they only had one inequality to constrain a linear combination of coefficients. At d = 5 and 6, they obtained, $\left|{\sum }_{{jm}}{Y}_{{jm}}(\theta ,\phi )\,{k}_{(V){jm}}^{(5)}\right|\leqslant 2\,\times \,{10}^{-14}\,{\rm{m}}$ and $\left|{{\sum }_{{jm}}}_{+4}{Y}_{{jm}}(\theta ,\phi )\left({k}_{(E){jm}}^{(6)}+{{ik}}_{(B){jm}}^{(6)}\right)\right|\leqslant 8\,\times \,{10}^{-9}\,{{\rm{m}}}^{2}$, where θ ≈ 160° and ϕ ≈ 120° are the rough sky position of GW150914 in the celestial-equatorial frame; Y_jm (also denoted as ₀ Y_jm ) and ₊₄ Y_jm are the spin-weighted spherical harmonics that are direction-dependent; ${k}_{(V){jm}}^{(5)}$, ${k}_{(V){jm}}^{(6)}$, and ${k}_{(B){jm}}^{(6)}$ are SME coefficients controlling the gLIV. The constraint at d = 5 is the first constraint on gLIV operators, while the d = 6 constraint is comparable to the limits from short-range laboratory tests (see, e.g., Shao et al. 2016).

Afterwards, 11 events were published in the GW transient catalog GWTC-1 (Abbott et al. 2019a), which enabled a global analysis to constrain all the coefficients simultaneously. Shao (2020) extended the analysis method in Kostelecký & Mewes (2016) to a global approach, using the whole GWTC-1. The analysis broke the degeneracy in SME coefficients and improved the limits by 2–5 orders of magnitude.

Notice that even in the mass dimension 5 there are still 16 independent components of the violation coefficients, but 11 events could only provide 11 inequalities on the linear combination of them. Therefore, these inequalities cannot fully break the parameter degeneracy at d = 5, not to mention the higher mass dimensions. However, with more and more observations of GW events, now with the addition of GWTC-2 (Abbott et al. 2021), we can completely decouple the degenerate parameters.

In this work we will focus on constraining the gLIV coefficients using the whole GW transient catalogs GWTC-1 and GWTC-2 (Abbott et al. 2019a, 2021), which contain 50 GW events in total. Benefiting from the multiple GW events coming from different sky areas, we finally break the degeneracy among all coefficients at d = 5 and d = 6, and obtain limits that are 2 to 7 times tighter than the existing best results.

This article is organized as follows. In Section 2 we introduce the dispersion relation of GWs in the linearized gravity with gLIV, and the effects of anisotropic birefringence during the propagation of GWs. In Section 3 we select the source parameters from the posterior samples of 50 GW events. The constraints on the time delay between the two modes of GWs are given. In Section 4 we carry out a global analysis of the SME coefficients using Monte Carlo methods. The global limits are simultaneously obtained, breaking the degeneracy among various coefficients. Section 5 presents a summary and discusses some possible future improvements. Throughout the article, we use natural units, where ℏ = c = 1.

To study the linearized gravity, we first decompose the metric into the Minkowski metric and a metric perturbation, that is, g_{μ ν} = η_{μ ν} + h_{μ ν} . Then the Lagrangian at dimension d in the pure gravity sector of SME, containing the LI and LIV terms, is expressed as (Kostelecký & Mewes 2016),

$$\begin{eqnarray}&&{{ \mathcal L }}_{{{ \mathcal K }}^{(d)}}=\displaystyle \frac{1}{4}{h}_{\mu \nu }{\hat{{ \mathcal K }}}^{(d)\mu \nu \rho \sigma }{h}_{\rho \sigma },\end{eqnarray} \tag{ 1 }$$

where ${\hat{{ \mathcal K }}}^{(d)\mu \nu \rho \sigma }={{ \mathcal K }}^{(d)\mu \nu \rho \sigma {\alpha }_{1}{\alpha }_{2}\cdots {\alpha }_{d-2}}{\partial }_{{\alpha }_{1}}{\partial }_{{\alpha }_{2}}\cdots {\partial }_{{\alpha }_{d-2}}$ is an operator at mass dimension d ≥ 2 with tensor coefficients ${{ \mathcal K }}^{(d)\mu \nu \rho \sigma {\alpha }_{1}{\alpha }_{2}\cdots {\alpha }_{d-2}}$. Obviously, because the indices "μ ν" and "ρ σ" will be contracted with the symmetric metric perturbation h_{μ ν} , the two pairs of indices are symmetric. Dedicated investigation (Kostelecký & Mewes 2016) showed that there are 14 independent classes of operators that control the behavior of GWs. However, the operators still need to obey the infinitesimal gauge transformation h_{μ ν} → h_{μ ν} + ∂_(μ ξ_ν), where ξ^ν is an arbitrary infinitesimal vector field. It results in only 3 of the 14 classes of operators surviving. Summing the Lagrangian at all mass dimensions, Kostelecký & Mewes (2016) obtained the Lagrangian,

$$\begin{eqnarray}&&{ \mathcal L }={{ \mathcal L }}_{0}+\displaystyle \frac{1}{4}{h}_{\mu \nu }\left({\hat{s}}^{\mu \rho \nu \sigma }+{\hat{q}}^{\mu \rho \nu \sigma }+{\hat{k}}^{\mu \nu \rho \sigma }\right){h}_{\rho \sigma }\,,\end{eqnarray} \tag{ 2 }$$

where ${{ \mathcal L }}_{0}=\tfrac{1}{4}{\epsilon }^{\mu \rho \alpha \kappa }{\epsilon }^{\nu \sigma \beta \lambda }{\eta }_{\kappa \lambda }{h}_{\mu \nu }{\partial }_{\alpha }{\partial }_{\beta }{h}_{\rho \sigma }$ is the linearized Einstein–Hilbert part. The remaining three LIV operators have similar forms like ${\hat{{ \mathcal K }}}^{(d)\mu \nu \rho \sigma }$, but the positions of the contracted indices and their symmetries are different. Their specific forms and properties are discussed in Kostelecký & Mewes (2016).

With the Lagrangian in Equation (2), it is straightforward to obtain equations of motion in vacuum (Mewes 2019). For plane waves, applying Fourier transformation turns the derivative ∂_μ into the 4-momentum, changing these differential equations into algebraic equations of p^μ = (ω, p ). Assuming that the violation coefficients are small, Mewes (2019) found the solution in the temporal gauge and derived the dispersion relation,

$$\begin{eqnarray}&&\omega =\left(1-{\varsigma }^{0}\pm \sqrt{{\left|{\varsigma }_{(+4)}\right|}^{2}+{\left|{\varsigma }_{(0)}\right|}^{2}}\right)p,\end{eqnarray} \tag{ 3 }$$

where the leading part gives the ordinary dispersion relation ω = p, which means that GWs travel at the speed of light. The other terms are combinations of ${\hat{s}}^{\mu \rho \nu \sigma }$, ${\hat{q}}^{\mu \rho \nu \sigma }$, and ${\hat{k}}^{\mu \nu \rho \sigma }$, with the derivatives replaced by 4-momentums, representing the contributions of gLIV.

In the helicity basis, ς⁰ and ς₍₀₎ are real and have zero helicity, while ς₍₊₄₎ is complex and has helicity +4 (Mewes 2019). Taking the ordinary dispersion relation at leading order, we can decompose the frequency ω and the propagation direction $\hat{{\boldsymbol{p}}}$ in the violating terms. Using the spin-weighted spherical harmonics _s Y_jm to expand the direction-dependent part, Mewes (2019) got the explicit series expansion of the violating terms,

$$\begin{eqnarray}\begin{array}{rcl}{\varsigma }^{0} & = & \displaystyle \sum _{{djm}}{\left(-1\right)}^{j}\,{\omega }^{d-4}{\,}_{0}{Y}_{{jm}}(\hat{{\boldsymbol{p}}}){k}_{(I){jm}}^{(d)},\\ {\varsigma }_{(+4)} & = & \displaystyle \sum _{{djm}}{\left(-1\right)}^{j}\,{\omega }^{d-4}{\,}_{-4}{Y}_{{jm}}(\hat{{\boldsymbol{p}}})\\ & & \times \left({k}_{(E){jm}}^{(d)}+{{ik}}_{(B){jm}}^{(d)}\right),\\ {\varsigma }_{(0)} & = & \displaystyle \sum _{{djm}}{\left(-1\right)}^{j}\,{\omega }^{d-4}{\,}_{0}{Y}_{{jm}}(\hat{{\boldsymbol{p}}}){k}_{(V){jm}}^{(d)}.\end{array}\end{eqnarray} \tag{ 4 }$$

Coefficients ${k}_{{jm}}^{(d)}$ are the so-called SME coefficients that control the gLIV. All coefficients are complex numbers and obey $\ {k}_{j(-m)}^{(d)}={\left(-1\right)}^{m}{k}_{{jm}}^{(d)* }$. For the indices, m takes − j, ⋯ , j, while j takes ∣s∣, ⋯ , d − 2, with s being the spin weight in _s Y_jm . The valid values of the mass dimension d, are more complicated: (i) ${k}_{(I){jm}}^{(d)}$ only exists at even dimensions with d ≥ 4; (ii) ${k}_{(V){jm}}^{(d)}$ only exists at odd dimensions with d ≥ 5; and (iii) ${k}_{(E){jm}}^{(d)}$ and ${k}_{(B){jm}}^{(d)}$ only exist at even dimensions with d ≥ 6.

From Equation (3), one obtains the phase speed of GWs,

$$\begin{eqnarray}&&{v}_{\pm }=1-{\varsigma }^{0}\pm | \vec{\varsigma }\,| ,\end{eqnarray} \tag{ 5 }$$

where $| \vec{\varsigma }\,| \equiv \sqrt{{\left|{\varsigma }_{(+4)}\right|}^{2}+{\left|{\varsigma }_{(0)}\right|}^{2}}.$ Like Shao (2020), we assume that gLIV mainly occurs at a specific dimension, so it is convenient to define a new parameter ${\varsigma }^{(d)}(\hat{{\boldsymbol{p}}})\equiv \varsigma (\hat{{\boldsymbol{p}}}){\omega }^{4-d}$, which only depends on the propagation direction of GWs. Then we can separate the frequency and direction parts, and Equation (5) simplifies to

$$\begin{eqnarray}&&{v}_{\pm }=1-{\omega }^{d-4}\left({\varsigma }^{(d)0}\mp {\left|\vec{\varsigma }\,\right|}^{(d)}\right).\end{eqnarray} \tag{ 6 }$$

Note that the terms in brackets are only dependent on direction now.

The focus of our study lies in the birefringence phenomenon, occurring when ${\left|\vec{\varsigma }\,\right|}^{(d)}\ne 0$. This can be divided into two simplified situations: (i) ${\varsigma }_{(0)}^{(d)}\ne 0$, ${\varsigma }_{(+4)}^{(d)}=0$, and (ii) ${\varsigma }_{(+4)}^{(d)}\ne 0$, ${\varsigma }_{(0)}^{(d)}=0$. They correspond to odd dimensions with d ≥ 5 and even dimensions with d ≥ 6 respectively. When birefringence occurs, the two modes propagate at different speeds v₊ and v₋, which will cause an extra time difference between them after traveling from source to detector. In an expanding universe, Kostelecký & Mewes (2016) derived the time difference between two modes,

$$\begin{eqnarray}&&{\rm{\Delta }}t=2\,{\left|\vec{\varsigma }\,\right|}^{(d)}{\left(\displaystyle \frac{\omega }{1+z}\right)}^{d-4}{\int }_{0}^{z}\displaystyle \frac{{\left(1+z^{\prime} \right)}^{d-4}}{H(z^{\prime} )}{\rm{d}}z^{\prime} ,\end{eqnarray} \tag{ 7 }$$

where ${\left|\vec{\varsigma }\,\right|}^{(d)}$ is a function of the SME coefficients, ω is the frequency when GWs were emitted, z is the cosmological redshift, and H(z) is the Hubble parameter. Using the ΛCDM model, we have $H(z)={H}_{0}\sqrt{{{\rm{\Omega }}}_{{\rm{m}}}{\left(1+z\right)}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}}$, where H₀ ≈ 67.4 km s⁻¹ Mpc⁻¹ is the Hubble constant, Ω_m ≈ 0.315 is the fraction of matter density, and Ω_Λ ≈ 0.685 is the fraction of vacuum energy density. Once the SME coefficients and the information of GWs are given, with Equation (7), we can calculate the time difference between two modes for every GW event.

As discussed in the previous section, a general gLIV will result in anisotropy, dispersion, and birefringence for GWs, wherein the anisotropic birefringence phenomenon can be used to constrain the SME coefficients very well. Kostelecký & Mewes (2016) first attempted to impose a restriction on SME coefficients at d = 5 and 6 using GW150914. At the amplitude peak, they found no indication of birefringent splitting. Therefore, the time difference must be less than the width of the peak, which is approximately 3 ms. They conservatively estimated that the frequency f was 100 Hz near the peak. Substituting these values into Equation (7), they obtained the first set of constraints. For GWs in GWTC-1, Shao (2020) carried out the first global analysis of the SME coefficients for GWs and significantly improved the bounds. Following them, here we use 50 GW events and constrain fully the coefficient space at mass dimensions d = 5 and 6.

For the source parameters of GWs in GWTC-1, we follow the considerations in Shao (2020) and use the posterior samples provided by the LIGO/Virgo Collaboration (Abbott et al. 2019a). We use the Overall_posterior samples for the 10 binary black holes (BBHs), and the IMRPhenomPv2NRT_lowSpin_posterior samples for the binary neutron star (BNS), GW170817. For events in GWTC-2, we use the PublicationSamples data set in the posterior files, which were actually used to create results in Abbott et al. (2021). We have checked that, in practice, different waveforms do not have significant impacts on our final results (Shao 2020).

For the frequency at amplitude peak of the chirp signal, we use the fitting formula in Appendix A.3 in Bohé et al. (2017). However, the formula was obtained from BBH simulations, which means that when the GW source is a BNS or a neutron star–black hole (NS-BH) binary, the matter effects of NSs are not taken into account. Therefore, we choose a lower frequency, 800 Hz, for GW170817 as it was roughly the sensitivity upper edge (Shao 2020). For other sources whose lighter component mass is less than 3 M_⊙, GW190425 (Abbott et al. 2020a), GW190426_152155, and GW190814 (Abbott et al. 2020b), which may be of BNS or NS-BH origins, we give rough estimates through the last few cycles near merger of the template waveforms or the real strain data. We use 800 Hz for GW190425, 700 Hz for GW190426_152155, and 280 Hz for GW190814. These values are rather conservative, because higher frequency will give better limits according to Equation (7).

For the time difference between two helicity modes, we adopt the same criterion as in Shao (2020). For a definite GW event, we assume that the time delay between two modes satisfies ∣ωΔt∣ ≤ 2π/ρ, or equivalently,

$$\begin{eqnarray}&&\left|{\rm{\Delta }}t\right|\leqslant \displaystyle \frac{1}{\rho f},\end{eqnarray} \tag{ 8 }$$

where f is the GW frequency at the amplitude peak in the source frame and ρ is the network signal-to-noise ratio that can be obtained in Abbott et al. (2019a, 2021). We defer a refined criterion, possibly with matched filtering, to future work.

For the birefringence phenomenon, we will discuss the two lowest dimensions 5 and 6 separately for nonminimal gravity in the SME. According to Section 2, the indices of these SME coefficients are obtained as follows. For d = 5, one has j = 0, 1, 2, 3, thus there are 10 independent ${k}_{(V){jm}}^{(5)}$. However, ${k}_{(V){jm}}^{(5)}$ is real only when m = 0, so there are 16 independent components for ${k}_{(V){jm}}^{(5)}$ in total. For d = 6, note that because ∣s∣ ≤ j ≤ d − 2 and − j ≤ m ≤ j, j only takes 4. For the same reason as d = 5, there are nine independent components for ${k}_{(E){jm}}^{(6)}$ and ${k}_{(B){jm}}^{(6)}$, respectively, thus one has 18 independent components in total.

First, we take d = 5 to show how 50 GW events can constrain the coefficient space. For every GW event, based on Equations (7) and (8) we have an inequality of the 16 components, and each inequality gives an area between two hypersurfaces symmetrical about the origin. ⁴ It is easy to see that Δt is actually a linear combination of ${k}_{(V){jm}}^{(5)}$, and the combination coefficients are the spin-weighted spherical harmonics, ₀ Y_jm , which are direction-dependent. Based on the fact that GW events scatter in different sky areas, these linear combinations are expected to be linearly independent. Consequently, the 50 pairs of hypersurfaces can enclose a region containing the origin, where all SME coefficients must be. In this way we have completed a simple but global limitation of the coefficient space. In this closed area every coefficient is bounded. Then we can obtain the value range of each parameter in principle.

To incorporate all uncertainties and statistically obtain the distribution of these SME coefficients, one requires the distribution of the "measured" time delay. As a conservative consideration, we can take the upper bound of Δt in Equation (7) as its standard deviation. Because the violations in SME are expected to be tiny, we take the central value of Δt to be zero. We assume that the "measured" time delay obeys a Gaussian normal distribution with μ = 0 and σ = 1/(ρ f). Based on this, we carry out a global analysis following the steps below.

• 1.  
  For each GW event, we randomly draw the source parameters in the posterior samples, and calculate the frequency at the amplitude peak.
• 2.  
  We randomly draw the time delays of the 50 GW events from their respective distributions and take these time delays as the result of one measurement in the Monte Carlo processes.
• 3.  
  Given the SME coefficients, we calculate the theoretical time delay ${\rm{\Delta }}{t}_{n}^{(\mathrm{th})}$ for the n-th GW event where n takes 1, 2, ⋯ , 50.
• 4.  
  Considering that the measured time delay obeys a Gaussian distribution with $\mu ={\rm{\Delta }}{t}_{n}^{(\mathrm{th})}$ and σ = 1/(ρ_n f_n ), we construct the multi-Gaussian likelihood as a function of the SME coefficients for all events. After maximizing the likelihood, the corresponding SME coefficients at the extreme point are what we extract from each draw.
• 5.  
  We repeat steps (1)–(5) enough times until a stable distribution is obtained.

For the 16 components at d = 5, the violin plots of the one-dimensional marginalized distributions are given in Figure 1, while their confidence intervals are given in Table 1. Because the dependence of Δt on the SME coefficients is linear, the marginalized distributions are very close to Gaussian distributions. Similar to Shao (2020), no obvious correlation is found among the components due to the different properties of the GW events, including sky locations, peak GW frequencies, and so on.

Zoom In Zoom Out Reset image size

Table 1. The 68% Confidence Intervals of the 16 Independent Components of ${k}_{(V){jm}}^{(5)}$

j	m	Component	Confidence Interval [10−16 m]
0	0	${k}_{(V)00}^{(5)}$	(−4.00, 4.60)
1	0	${k}_{(V)10}^{(5)}$	(−3.81, 4.14)
1	1	$\mathrm{Re}\,{k}_{(V)11}^{(5)}$	(−3.37, 3.16)
		$\mathrm{Im}\,{k}_{(V)11}^{(5)}$	(−3.38, 3.26)
2	0	${k}_{(V)20}^{(5)}$	(−4.34, 4.19)
2	1	$\mathrm{Re}\,{k}_{(V)21}^{(5)}$	(−2.54, 2.70)
		$\mathrm{Im}\,{k}_{(V)21}^{(5)}$	(−2.78, 2.93)
2	2	$\mathrm{Re}\,{k}_{(V)22}^{(5)}$	(−3.64, 3.70)
		$\mathrm{Im}\,{k}_{(V)22}^{(5)}$	(−2.86, 3.05)
3	0	${k}_{(V)30}^{(5)}$	(−4.21, 4.33)
3	1	$\mathrm{Re}\,{k}_{(V)31}^{(5)}$	(−3.06, 3.02)
		$\mathrm{Im}\,{k}_{(V)31}^{(5)}$	(−3.19, 3.09)
3	2	$\mathrm{Re}\,{k}_{(V)32}^{(5)}$	(−3.13, 3.02)
		$\mathrm{Im}\,{k}_{(V)32}^{(5)}$	(−2.71, 2.64)
3	3	$\mathrm{Re}\,{k}_{(V)33}^{(5)}$	(−3.22, 3.16)
		$\mathrm{Im}\,{k}_{(V)33}^{(5)}$	(−3.49, 3.49)

Download table as:  ASCIITypeset image

For the 18 components at mass dimension d = 6, we have to analyze with more care. Since ς₍₊₄₎ is a complex number, the calculation of $\sqrt{{\left|{\varsigma }_{(+4)}\right|}^{2}}$ is much more difficult than d = 5, whereas ς₍₀₎ is real and the square root is trivial. Due to this reason, the dependence of Δt on the SME coefficients is highly nonlinear, leading the marginalized distributions to be more peaked than the Gaussian distribution in the end. Despite the complications above, correlations between the 18 independent components are still small. Their one-dimensional marginalized distributions are given in Figure 2 and the 68% confidence intervals are illustrated in Table 2.

Zoom In Zoom Out Reset image size

Table 2. Same as Table 1, for ${k}_{(E){jm}}^{(6)}$ and ${k}_{(B){jm}}^{(6)}$

j	m	Component	Confidence Interval [10−12 m2]
4	0	${k}_{(E)40}^{(6)}$	(−1.70, 1.73)
4	1	$\mathrm{Re}\,{k}_{(E)41}^{(6)}$	(−1.76, 1.90)
		$\mathrm{Im}\,{k}_{(E)41}^{(6)}$	(−2.05, 2.16)
4	2	$\mathrm{Re}\,{k}_{(E)42}^{(6)}$	(−2.08, 2.15)
		$\mathrm{Im}\,{k}_{(E)42}^{(6)}$	(−1.75, 1.78)
4	3	$\mathrm{Re}\,{k}_{(E)43}^{(6)}$	(−2.24, 2.09)
		$\mathrm{Im}\,{k}_{(E)43}^{(6)}$	(−2.19, 2.11)
4	4	$\mathrm{Re}\,{k}_{(E)44}^{(6)}$	(−2.12, 2.09)
		$\mathrm{Im}\,{k}_{(E)44}^{(6)}$	(−2.10, 2.01)
4	0	${k}_{(B)40}^{(6)}$	(−1.97, 2.00)
4	1	$\mathrm{Re}\,{k}_{(B)41}^{(6)}$	(−1.70, 1.67)
		$\mathrm{Im}\,{k}_{(B)41}^{(6)}$	(−1.95, 1.91)
4	2	$\mathrm{Re}\,{k}_{(B)42}^{(6)}$	(−2.01, 2.37)
		$\mathrm{Im}\,{k}_{(B)42}^{(6)}$	(−1.78, 1.81)
4	3	$\mathrm{Re}\,{k}_{(B)43}^{(6)}$	(−2.16, 2.03)
		$\mathrm{Im}\,{k}_{(B)43}^{(6)}$	(−2.02, 2.18)
4	4	$\mathrm{Re}\,{k}_{(B)44}^{(6)}$	(−2.09, 2.12)
		$\mathrm{Im}\,{k}_{(B)44}^{(6)}$	(−2.16, 2.09)

Download table as:  ASCIITypeset image

GWs provide excellent opportunities to make a profound study of the gravitational Lorentz invariance violation (gLIV). Many studies (Nishizawa 2018; Abbott et al. 2019b; Mirshekari et al. 2012; Wang & Zhao 2020; Wang et al. 2021) have been performed to test the Lorentz symmetry in general relativity (GR) using GW events. While a general gLIV includes anisotropy, dispersion, and birefringence, previous work mostly emphasized the isotropic birefringent phenomenon.

In our work, we focus on the anisotropic birefringence of GWs in the popular theoretical framework called the standard-model extension (SME; Colladay & Kostelecký 1997; Kostelecký 2004). Assuming that gLIV mainly occurs at a particular mass dimension, for the two lowest dimensions that can produce birefringence, d = 5 and 6, we separately carry out global analysis on them and constrain the whole coefficient space fully. Thanks to the multiple GW events in the GW transient catalogs GWTC-1 and GWTC-2, we completely break the degeneracy among the 16 independent components of the SME coefficients ${k}_{(V){jm}}^{(5)}$ at d = 5, and 18 independent components at d = 6. This is reflected in the results of Tables 1 and 2, which are obtained by simultaneously constraining all the coefficients, instead of setting only one coefficient to be nonzero each time (i.e., the "maximum reach" approach; Tasson 2019). As can be seen, there is no obvious violation of GR in our results, and our limits on the SME coefficients are 2 to 7 times tighter than that in the earlier work (Shao 2020).

While in this paper, we have only made use of one type of phenomena in the birefringent gLIV, namely the differences of the arrival times between two modes, there are many more approaches to improve our tests. First, more information can be extracted from the birefringence phenomena, such as changes of the polarizations. Similar analysis has been conducted in the photon sector (Kostelecký & Mewes 2008, 2013). Second, Equation (7) is a relatively simple estimation of the time delay (Shao 2020), and in fact Mewes (2019) has worked out the whole deformed waveform analytics and we can use these modified templates to implement the matched-filtering analysis. For mass dimensions 5 and 6, the computational cost is expected to be significant (O'Neal-Ault et al. 2021). Third, in a very unlikely scenario, the pattern function of GW detectors might hide one of the two polarizations. We think it is impossible for this to have happened in all 50 GW events that we used here. However, to entirely rule out this possibility, one needs to incorporate the polarization of the events, probably in a matched-filtering analysis (Wang et al. 2021; O'Neal-Ault et al. 2021). Finally, advanced detectors like the KAGRA (Aso et al. 2013) have recently joined the global efforts to detect GWs. With the ground-based detector network, more and more GWs in a wider frequency range will obviously provide direct improvements to our results.

In conclusion, the GWs are undoubtedly fantastic phenomena for testing the fundamental principles of modern physics, and we expect a brighter future in the coming years for fundamental physics thanks to the use of these natural treasures.

We are grateful to Yi-Fan Wang for discussions and the LIGO/Virgo Collaboration for providing the posterior samples of their parameter-estimation studies. This work was supported by the National Natural Science Foundation of China (11975027, 11991053, 11721303), the National SKA Program of China (2020SKA0120300), the Young Elite Scientists Sponsorship Program by the China Association for Science and Technology (2018QNRC001), the Max Planck Partner Group Program funded by the Max Planck Society, and the High-performance Computing Platform of Peking University. Z.W. is supported by the Principal's Fund for the Undergraduate Student Research Study at Peking University.