Combined search for anisotropic birefringence in the gravitational-wave transient catalog GWTC-1

The discovery of gravitational waves (GWs) provides an unprecedented arena to test general relativity, including the gravitational Lorentz invariance violation (gLIV). In the propagation of GWs, a generic gLIV leads to anisotropy, dispersion, and birefringence. GW events constrain the anisotropic birefringence particularly well. Kosteleck\'y and Mewes (2016) performed a preliminary analysis for GW150914. We improve their method and extend the analysis systematically to the whole GW transient catalog, GWTC-1. This is the first global analysis of the spacetime anisotropic Lorentzian structure with a catalog of GWs, where multiple events are crucial in breaking the degeneracy among gLIV parameters. With the absence of abnormal propagation, we obtain new limits on 34 coefficients for gLIV in the nonminimal gravity that surpass previous limits by $\sim 10^2$-$10^5$.


I. INTRODUCTION
Lorentz invariance is a cherished symmetry laying to the heart of modern physics. However, motivated by contemporary open questions, there are good reasons to believe that, the Lorentz symmetry breaks at some yet unknown energy scale [1,2]. For example, in the string theory it is perceived to have a Higgs-like spontaneous symmetry breaking of the Lorentz invariance [3,4] which leads to various observable phenomena [5,6]. While the Lorentz invariance violation (LIV) in principle happens in different species sectors [7][8][9], it is the most interesting to study the gravitational LIV (gLIV) due to the following fact [10][11][12][13]. Up to now, the canonical theory of gravitation, namely the general relativity, while being extremely faithful in describing various gravity experiments [14,15], still refuses to embed into the framework of quantum field theory which successfully describes the other three fundamental interactions in a unified way. Therefore, searches for gLIV are closely linked to searches for quantumgravity candidate theories [15][16][17][18].
Nowadays the most popular framework to investigate the possibility of gLIV is the standard-model extension (SME) [10] and the parameterized post-Newtonian formalism [14,19]. We will focus on the former. The SME is an agnostic effective field theory which incorporates all gaugepreserving, Lorentz-violating, energy-momentum-conserving operators that are constructed from the field operators in general relativity and particle physics. New Lorentz-covariant operators are built from contraction of known fields with LIV coefficients. Unless being protected, operators with lower mass dimension are expected to dominate the observables as per the spirit of effective field theories. Here we focus on the subset of the SME that deals with the pure gravity sector; interested readers are referred to Ref. [20] for discussions on the mattergravity couplings. In an effective-field-theoretic framework, to be compatible with the Riemann-Cartan geometry, the symmetry breaking shall be spontaneous instead of explicit [21]. Therefore the SME, while being particle Lorentz-violating, is observer Lorentz-invariant. * lshao@pku.edu.cn The convenience in using the SME to test the Lorentz symmetry has resulted in a flourish of studies during the past decades from different kinds of experiments concerning various particle species [6]. For the gravity sector, constraints come from lunar laser ranging [22,23], atom interferometers [24], cosmic rays [25], precision pulsar timing [26][27][28][29][30][31], planetary orbital dynamics [32], laboratory short-range experiments [33][34][35][36], super-conducting gravimeters [37], and recently gravitational waves (GWs) [38][39][40]; see Hees et al. [11] for a comprehensive review. These constraints are complementary to each other, and in many cases they are individually competent at probing different parts of the parameter space. For example, timing of binary pulsars was shown to be good at systematically probing the lowest-order CPTviolating [29,41], as well as gravitational weak-equivalenceprinciple-violating [30,42], operators via the post-Newtonian dynamics, while laboratory short-range experiments excel in constraining the nonminimal operators in a static Newtonian setup [35].
A newly emerging probe to gLIV phenomena in gravity is the recently discovered GWs by the Advanced LIGO/Virgo detectors [38,43,44]. The propagation of GWs has become one of the central topics in looking for gLIV clues [45][46][47][48][49][50][51][52][53][54]. The existence of gLIV would generally cause anisotropy, dispersion, and birefringence [38,39]. Following the theoretical work by Will [45] and Mirshekari et al. [46], the LIGO/Virgo Collaboration have conducted extensive tests of the GW dispersion relation caused by the gLIV or a massive graviton [48,55]. However, most of the past study focused on the boost breaking, and ignored the breaking of the rotational symmetry. While the commutator of two boost generators gives a rotation generator, it is inevitable to host the rotation breaking when the boost symmetry is broken, unless the object under investigation is exactly at rest with respect to a preferred frame. Such a preferred frame may not even exist in effective field theories [56,57]. This motivates us to look at the anisotropic gLIV and test it thereof in a more generic way.
In this work we follow the preliminary analysis done by Kostelecký and Mewes [38], and extend it systematically with all the GW events detected in the first and second observing runs [44]. It is the first study of this kind, and probes the gLIV with a globally coherent approach. Using the fact that the SME is observer Lorentz-invariant, we combine different GW events from the GW transient catalog, GWTC-1 [44], including ten binary black holes (BBHs) and one binary neutron star (BNS), GW 170817 [58]; see Appendix A for basic parameters of GW events in the GWTC-1. A global analysis with the basis of spin-weighted spherical harmonics provides us new limits on the gLIV coefficients in the nonminimal gravity at mass dimensions 5 and 6. These limits are orders of magnitude tighter than the existing ones. GW observations truly represent precious treasures in studying various gravity theories. While more detections are continuously being made, the investigation in this paper will be improved dramatically and new clues to quantum gravity might be drawn.
The paper is organized as follows. In the next section, we briefly discuss the modified dispersion relation for the linearized gravity in SME, and its effect on the cosmological GW propagation. Anisotropic birefringent effects can be constrained particularly well with GWs. Therefore we focus on the anisotropic birefringence in Sec. III, and lay down the strategies for the practical application to the GWTC-1. In Sec. IV we present detailed Monte Carlo analysis using the posteriors on the GW parameters for the 11 events in the catalog, provided by the LIGO/Virgo Collaboration. Numerical constraints on 34 gLIV coefficients are obtained with the maximal-reach [59] and global approaches. The global analysis, simultaneously with tens of gLIV parameters analyzed, is done for the first time with GWs in a coherent way. It greatly extends previous work done by the LIGO/Virgo Collaboration [48,55,60] and Kostelecký and Mewes [38]. Finally, Sec. V summarizes the paper and discusses future directions for improvement. For readers' convenience, extra information on the GWTC-1 catalog, the spin-weighted spherical harmonics, and a fitting formula to the GW peak frequency are provided respectively in Appendices A, B, and C.
Throughout the paper, we follow conventions used by Misner et al. [61] and Kostelecký [10]. Unless otherwise stated, we use natural units where = G = c = 1.

II. DISPERSION AND PROPAGATION OF GRAVITATIONAL WAVES
The gravity sector in the SME was given fully in the Riemann-Cartan spacetime [10]. However, we restrict ourselves to a 4-dimensional Riemannian spacetime, and only consider the part of spacetime where, after fixing the gauge, the linearized gravity is a sufficiently good approximation. The metric, g µν = η µν +h µν , is expanded around the Minkowski metric, η µν , with a perturbation, h µν . Kostelecký and Mewes [38] constructed the general quadratic Lagrangian density for GWs in the presence of gLIV operators of arbitrary mass dimension d, with the "hat" denoting its operator nature, and K (d)µνρσi 1 i 2 ···i d−2 are normal tensorial components whose mass dimension is 4 − d.
In addition, any contraction of these operators with a derivative is zero. Similar to the previous study in the photon sector of the electrodynamic SME [63,64], the covariant dispersion relation for the two tensor modes of a GW with 4-momentum p µ = (ω, p) is [38], where Here the decomposition is done by the handedness of GWs, instead of the usual "+" and "×" modes. In above equations, ζ 0 and ζ 3 are rotation scalars, ζ 1 and ζ 2 are helicity-4 tensors, and the derivative "∂ µ " in operators is understood to be replaced by "ip µ " [38].
To report experimental constraints on the coefficients for Lorentz/CPT violation in the SME, it is conventional to use the canonical Sun-centered celestial-equatorial frame [63]. Concerning the rotational behavior, it is useful to decompose ζ α with spin-weighted spherical harmonics [38,64], wheren ≡ −p is the direction to the source, and |s| ≤ j ≤ d − 2; see Appendix B for more details on the spin-weighted spherical harmonics. It was shown that [38], (i) anisotropic effects are governed by the coefficients with j 0; (ii) frequency-dependent dispersions are governed by all coefficients except k (4) (I) jm ; (iii) vacuum birefringent effects are governed by all coefficients except k (d) (I) jm .
Besides, d ≥ 4 is even for k (d) (I) jm ; d ≥ 5 is odd for k (d) (V) jm ; and d ≥ 6 is even for k (d) (E) jm and k (d) (B) jm . Assuming that the corrections to Einstein's general relativity are small, we obtain the speed of GWs from Eq. (3), In accordance to the spirit of effective field theories, we further assume that the gLIV happens dominantly at a specific mass dimension d. If it had happened at multiple dimensions, we equivalently take the dimension d, usually the lowest relevant mass dimension where gLIV happens, which introduces the maximum effect. By this assumption we ease the summation over "d" in Eqs. (7)(8)(9), and denote the original ζ α as ζ α (d) (α = 0, 1, 2, 3). Further, we introduce a notation,ζ α where an important fact is that the terms inside the curly bracket is, while being direction-dependent, energyindependent. Now consider two gravitons emitted at t e and t e with energies ω e and ω e respectively. After traveling over a same  (12), for mass dimensions d = 5 and d = 6. The uncertainties are given at the 90% confidence level. Notice that, because the merger part for the BNS was not observed [58,67], in our test we conservatively use f GW = 800 Hz for GW170817. comoving distance [45,46,68], they arrive at a GW detector on the Earth at t a and t a . Following the method developed by Will [45] and Mirshekari et al. [46], we derive where ∆t e ≡ t e − t e , ∆t a ≡ t a − t a , ∆ω d−4 − ω e d−4 , and z is the cosmological redshift. A new distance-like quantity in the above equation is defined as [45], where H 0 = (67.4 ± 0.5) km s −1 Mpc −1 is the Hubble constant, Ω m = 0.315 ± 0.007 is the fraction of matter energy density, and Ω Λ = 0.685 ± 0.007 is the fraction of vacuum energy density in our current Universe [69]. We have used the standard ΛCDM model, which should be rather precise for the GW events with relatively low redshifts. For the 11 GW events in the GWTC-1, their D (d) α for d = 5 and d = 6 are given in Fig. 1 and Table I.

III. CONSTRAINING ANISOTROPIC BIREFRINGENCE
The modified dispersion relation (3) introduces anisotropy, dispersion, and birefringence to the propagation of GWs [38]. As the first application, Kostelecký and Mewes [38] used the observation that there is no indication of mode splitting at the amplitude peak of GW150914 [43]. They took a rough value for the upper limit of the time difference for the arrival of two circular modes, ∆t ≤ 3 ms, and used a central frequency f ∼ 100 Hz, to derive the limit on the difference in the propagation speed between the two circularly polarized CPT-conjugate eigenmodes. They obtained the first constraint on the gLIV in the pure-gravity sector with d = 5, and a competitive limit to existing laboratory bounds on bire- for θ 160 • and φ 120 • [38], where (θ, φ) is the (very) rough sky position of GW150914 in the Sun-centered celestial-equatorial frame. The results are, though heuristic, very encouraging.
Equations (13) and (14) are actually bounds on a set of linear combination of gLIV coefficients. Now we improve the analysis method to a global approach, and extend the study to the whole GWTC-1 catalog. Due to the presence of multiple events, we are privileged to carry out a global analysis that breaks the degeneracy of various gLIV parameters. Such a global analysis was not possible for the time being of Ref. [38] with only GW150914 detected then.
In order to construct the gLIV tests with Eq. (11) in practice, we make the following considerations.
(i) For the same reason as that in the photon sector [63,64,71], in these propagation tests vacuum birefringent phenomena can constrain gLIV parameters more tightly than the dispersive ones. Because ζ 0 (d) , so asζ 0 (d) , introduces polarization-independent time delays, they are comparably loosely bound. Nevertheless, if GW companion particles (photons or neutrinos) are detected, the dispersion can be fairly constrained; for example, see the bound on the SMEs (4) 00 parameter via the simultaneous observation of GW170817 and GRB 170817A [72]. Probablity density Probability density for f GW , generated using the posterior samples provided by the LIGO/Virgo Collaboration [44,65,66] and the fit in Eq. (C2) [70].
To bound the scope of this paper with fair workload, these polarization-independent, dispersive-only delays, encoded in ζ 0 (d) andζ 0 (d) , are omitted in the following analysis. They can be incorporated in future work.
(ii) For the GW frequency at the amplitude peak, we use the dedicated fit for the (2, 2) mode in the Appendix A.3 of Ref. [70]. The fit was obtained by combining catalogs of numerical relativity [73] and test-particle Teukolskycode waveforms. It includes the contribution from the mass ratio of the binary and the orbit-aligned spins of the binary components. The fit is adopted in the socalled SEOBNRv4 waveform family [70], and is being extensively applied in the LIGO/Virgo daily data analysis. An explicit expression for the global fit is given in Appendix C and several examples for different spin combinations are given in Fig. 6 therein. Using the posterior samples provided by the LIGO/Virgo Collaboration [44,65,66], we plot the GW peak frequencies for the 11 GW events in Fig. 2. Notice that, while the fit was obtained from BBH simulations without matter effects, its prediction to the BNS event GW170817, whose nuclear matter effects enter the waveform at the fifth post-Newtonian order [67,74], should be indicative for the analysis in this work. Because the merger part for BNS was not actually observed in GW detectors [58,67], in our tests we conservatively use a rather arbitrary value f GW = 800 Hz for GW170817. It roughly corresponds to the cutoff sensitivity at high-frequency end for LIGO/Virgo detectors at the time of GW170817.
Larger values of f GW would lead to tighter constraints.
(iii) For the time delay between two circularly polarized GW modes, we use a simple order-of-magnitude estimation, |ω GW ∆t| ≤ 2π/ρ, where ∆t is the time difference between two circular modes, and ρ is the network signalto-noise ratio (SNR) of the observed GW events (see Table VI in Appendix A). We expect this estimation to work fairly well at the current stage in constraining the gLIV phenomena, instead of discovering the gLIV phenomena. One direction to improve the investigation would be using the matched-filtering technique with modified/deformed gravitational waveforms [39]. For such an improved study, our results from individual events can be rescaled accordingly. Study along this line can be used to verify our assumption. For now we leave a dedicated analysis to future work.
In a short summary, on one hand we have 11 two-sided constraints from 11 GW events in the GWTC-1, and on the other hand, we need to constrain 16 independent components when d = 5, and 9 + 9 = 18 independent components when d = 6. Therefore, it is an over-constraining system. The reason for the GW propagation constraints being "two sided" is that, we do not expect either circular mode travels faster than the other one; otherwise, the deformation in the waveforms would be quite obvious [39].

IV. NUMERICAL RESUTLS
With all the practical considerations in Secs. II and III taken into account, we present our final numerical constraints on the gLIV coefficients in this section. As mentioned above, it is an over-constraining system. Therefore, we can consequently bound globally all independent coefficients at a specific mass dimension d using the whole GW catalog. However, we will first consider the maximal-reach scenario [6,59] where only one independent component is assumed to be nonzero. It eases comparison with results in Ref. [38], and provides us some insight on the figure of merit for different GW events in testing the anisotropic birefringence.
Because the GW parameters are generally highly correlated, we use the posterior samples provided by the LIGO/Virgo Collaboration [65,66] which include all correlations among parameters and are publicly available. For the ten BBHs, we use the posterior samples tagged with Overall posterior that are combined results derived from the effective-one-body and phenomenological waveform families [44], while for the BNS GW170817, we use the samples tagged with IMRPhenomPv2NRT lowSpin posterior that are derived from an effectively-precessing phenomenological waveform family with the tidal deformability effects incorporated [67]. Other choices do not change our limits in a significant way. The GW parameters we use include, (i) the sky position represented by the right ascension, α, and the declination, δ (see Fig. 5 in Appendix A); (ii) the intrinsic GW parameters including the component masses, m 1 and m 2 , and the (dimensionless) orbital-aligned component spins, χ 1 and χ 2 ; and (iii) the luminosity distance d L from where the redshift, and then the D (d) α in Eq. (12), are derived with the standard ΛCDM model [69].
In the maximal-reach scenario, we assume that only one gLIV parameter is nonzero at a time [59]. For each GW event, we randomly draw samples from the posteriors. We calculate the individual limit for each gLIV parameter for each sample point. The results are stored for statistical inference afterwards. After accumulating enough samples, the 1-σ limits are obtained from their corresponding distributions. In the calculation, all statistical uncertainties are taken into account properly.
In Tables II and III, respectively for mass dimensions d = 5 and d = 6, we list the maximal-reach limits for each gLIV component. We have put the tightest limit from an individual GW event, along with the event name in the fourth and fifth columns. We can see that, for mass dimension-5 operators, the individually tightest limits come from GW151226, GW170608, and GW170817, while for mass dimension-6 operators, the individually tightest limits come from GW170608 and GW170817. These events all have low component masses. Usually, low-mass events are detected closer to the Earth, because their GW strain amplitudes are smaller than highmass ones if they were put at a same distance [75,76]. The finite sensitivity of GW detectors can only pick up those lowmass events relatively nearby [76]. For the three individually best events here (GW151226, GW170608, and GW170817), they all have the cosmological redshift z 0.1 (see Table VI). Although the GW propagation tests benefit from large distances [see Eq. (12)], which usually correspond to high-mass GW events (see Table VI), our anisotropic birefringent tests also depend on the GW frequency with a powerlaw index d−4 [see Eq. (11)], which is positive for d = 5 and 6. High-mass GW events have a lower GW frequency [see Eq. (C2)] [75], which decreases their ability to constrain the gLIV parameters. Our results clearly show that, already at mass dimension 5, low-mass events are preferred to the tests of the gLIV parameters in the SME. At mass dimension 6, the benefit from smaller masses is even more prominent.
For the maximal-reach scenario, we also list the combined limits from multiple GWs by assuming each event mutually independently. Therefore, one gets a combined limit where σ i is the limit from an individual GW event i. As we can see, the combined limit improves the limit from the individually best event by a factor smaller than two. Therefore, in the maximal-reach approach, the limit is dominated by the individually best GW event. Worth to mention that, the limits in Tables II and III   TABLE IV. Global constraints on k (5) (V) jm at the 68% confidence level. Notice that the unit differs from that in Table II. The maximal-reach limits should be considered as quite optimistic, and in reality they should be bounds on some linear combinations of the underlying set of gLIV coefficients [38]. Therefore, it becomes very intriguing to break the degeneracy among these parameters and check the real power of GW events in probing the gLIV when a number of events are available [6]. In the following, we use the simple facts, that (i) different GW events come from different directions [44] and (ii) symmetry breaking in the SME is observer Lorentzinvariant [21], to globally constrain these gLIV parameters. This is possible, because the coefficients of the linear combinations are direction-dependent, depending on the direction n via the spin-weighted spherical harmonics functions; see Eqs. (7)(8)(9) and Appendix B.
We have two sets of global analysis, for mass dimension 5 and mass dimension 6. In each global analysis, we assume that all anisotropic-birefringent gLIV coefficients can be nonzero at that specific d. Therefore, we have in total 16 independent coefficients for d = 5, and 18 independent coefficients for d = 6. Fortunately, the 11 events in GWTC-1 provide us in total 22 useful constraints.
The global analysis for d = 5 is relatively easier, because (i) the square root in Eq. (11) is plainly calculated wheñ ζ 1 (d=5) =ζ 2 (d=5) = 0, and (ii) only the s = 0 spherical harmonics are involved. Similarly to the maximal-reach case, we randomly draw posterior samples, but now simultaneously for 11 events. Their time delays are drawn from zero-mean Gaus- sian distributions with their standard variances determined in Sec. III. Then, for each random draw, we construct the global likelihood as a function of the 16 gLIV coefficients. The likelihood is maximized by the routines in the Scipy.Optimize package [77]. Thus, values of 16 coefficients are calculated in each random draw. The resulted values for the 16 gLIV coefficients are recorded for later inference. After we accumulate enough draws, we extract the constraints on the 16 gLIV coefficients from these recorded distributions.
The 1-dimensional marginalized constraints are listed in Table IV, and the contour plots for these parameters are given in Fig. 3. It is interesting to note that, (i) the global constraints are only about a factor of 10 weaker than the maximal-reach limits, and (ii) these 16 gLIV parameters are hardly correlated. The small correlation between parameters is resulted from the use of multiple events. It shows the advantage of constructing an over-constraining system with multiple events, analogous to the cases of using millisecond pulsars under a similar con- text [26,29,30].
The global analysis for d = 6 is somehow complicated. Now we have the spin-weighted spherical harmonics with s = ±4. To reduce the computational cost, we use s Y * jm = (−1) s+m −s Y j(−m) to calculate −4 Y jm from +4 Y jm . More importantly, in the most general case for d = 6, we haveζ 1 (d=6) 0 andζ 2 (d=6) 0. As a consequence, the square root in Eq. (11) is nontrivial. Therefore, the calculation with these highly nonlinear features takes a much longer computer time, namely, the Scipy.Optimize routines now need much more significant computational time to iterate, in order to locate the maximum of the 18-dimensional likelihood functions. Nevertheless, the principles are similar to the d = 5 case. The results for d = 6 are give in Table V and Fig. 4. Similar conclusions that were made for d = 5 operators in the last paragraph can be made for d = 6 operators as well.

V. DISCUSSION
Gravitational Lorentz invariance violation (gLIV) is a central topic in the program of using the gravitational waves (GWs) to probe the most fundamental principles in modern physics [38, 47-49, 55, 78, 79]. If the Lorentz symmetry is broken in the gravity sector, there might be a preferred frame [14,19] where the breaking is isotropic, but a more generic breaking does not need a preferred frame [10,57]. Isotropic gLIV, being a specific class of gLIV, was studied in details [46,48,49,54,55]. Most of previous work focused on non-birefringent phenomena in searches of the gLIV. Extensions to these studies are needed. On one hand, because the commutator of two boost generators is a generator for rotation, anisotropy is inevitable in a complete search for the gLIV. On the other hand, in effective field theories, the most general gLIV has birefringent behaviors for the two circularly polarized GW eigen-modes in the nonminimal gravity. In this work, we take a leap to systematically study anisotropic birefringent phenomena related to the GW propagation with the GW transient catalog GWTC-1 [44,65,66].
One of the best theoretical frameworks, in carrying out these anisotropic birefringent gLIV tests, is the standardmodel extension (SME) [10,38,39]. We follow the spirit of effective field theories, and assume that the gLIV happens at a specific mass dimension d. The lowest mass dimensions for vacuum birefringence are d = 5 and d = 6. When d = 5, there are in total 16 independent gLIV coefficients in k (5) (V) jm , while d = 6, there are 18 independent gLIV coefficients in k (6) (E) jm and k (6) (B) jm . In our global tests, we simultaneously include all gLIV operators that lead to birefringence at that particular mass dimension. We use the posterior samples for the events in GWTC-1, provided by the LIGO/Virgo Collaboration, to coherently solve for gLIV parameters, and obtain their constraints thereof. Our maximal-reach limits are listed in Tables II and III, and our global limits are presented in Tables IV  and V, as well as in Figs. 3 and 4. No violation of Einstein's general relativity was found, and the constraints on 34 gLIV coefficients are improved by factors ranging from ∼ 10 2 to ∼ 10 5 , with respect to previous limits [38].
The gLIV tests can be improved in multiple directions in the future. (I) A more sophisticated matched-filtering analysis with gLIV-deformed waveforms [39] can be used to validate the assumptions made in this work, though such an analysis could have cost mighty computational time in practice for a catalog of GWs. (II) Another possibility in testing the GW propagation can involve modified cosmological behaviors. For example, Nishizawa [49] considered a generic GW propagation equation, which -besides the cosmological expansion encoded in H ≡ a /a with a being the cosmological scale factor -includes running of the Planck mass M * via ν = H −1 d ln M 2 * /dt, the velocity of GWs c T , the mass of graviton µ, and the anisotropic source term Γγ i j . The philosophy of the approach (15) is different from ours which is based on a modified action of the gravity sector [see Eq. (2)]. Nevertheless, a grander theoretical framework that includes both a modified-gravity action and a modified cosmology might be feasible. It lays beyond the scope of this work however. (III) The final obvious direction to improve the tests in this work is to involve more GW events. Actually, more and yet more accurately-measured GW events are undoubtfully expected [76]. With the ongoing third observing run by the LIGO/Virgo Collaboration, more GW event candidates are already revealed for possible electromagnetic followups [80]. At the time of writing, a second BNS, GW190425, was published [81]. This event is weaker than GW170817, but it will help in constraining the gLIV parameters in the global analysis. With the KAGRA [82] and IndiGO GW detectors coming online in the near future, even better limits will be placed on the gLIV. Ultimately, we hope some positive clues to the long-sought quantum gravity theory might be uncovered via GWs.  In the Advanced LIGO/Virgo's first and second observing runs, respectively taking place from September 12, 2015 to January 19, 2016, and from November 30, 2016 to August 25, 2017, ten confident detections of BBHs and one confident detection of BNSs were reported [44]. In Table VI, we list the basic parameters and their uncertainties for these events [44]. As sky position is important in testing anisotropic birefringence, in Fig. 5 we show the 68% credible regions for sky location of GW events in the GWTC-1 [44], in a Mollweide projection in the equatorial coordinate. Data are taken from the associated data release for the sky maps [83] to the LIGO/Virgo's GWTC-1 [65], hosted by the Gravitational Wave Open Science Center (GWOSC) [84]. The plot has made use of the ligo.skymap package [85], maintained by Leo Singer.
In all of the numerical calculations of this paper, we have used the posterior samples provided by the LIGO/Virgo Collaboration, hosted at the GWOSC [66].
Notice that the uncertainties are rather heterogeneous for different GW events, due to the operation of GW detectors in practice at the detection time. In particular, the sky localization plays an essential role in determining the anisotropic behavior of gLIV, which should be taken into full account, as is done in this work.