Analysis of birefringence and dispersion effects from spacetime-symmetry breaking in gravitational waves

In this work, we review the effective field theory framework to search for Lorentz and CPT symmetry breaking during the propagation of gravitational waves. The article is written so as to bridge the gap between the theory of spacetime-symmetry breaking and the analysis of gravitational-waves signals detected by ground-based interferometers. The primary physical effects beyond General Relativity that we explore here are dispersion and birefringence of gravitational waves. We discuss their implementation in the open-source LIGO-Virgo algorithm library suite, as well as the statistical method used to perform a Bayesian inference of the posterior probability of the coefficients for symmetry-breaking. We present preliminary results of this work in the form of simulations of modified gravitational waveforms, together with sensitivity studies of the measurements of the coefficients for Lorentz and CPT violation. The findings show the high potential of gravitational wave sources across the sky to probe sensitively for these signals of new physics.


Introduction
Gravitational waves (GWs) are now a ripe testing ground for many aspects of gravitational physics [1][2][3][4]. One of the principle foundations of General Relativity is the Einstein Equivalence Principle, which includes the universality of freefall and the spacetimesymmetry principle of the local Lorentz invariance of physics [5]. The latter principle has seen a boom in tests in the last 20+ years [6], owing primarily to an interesting piece of motivation: that in a fundamental unified theory of physics, local Lorentz invariance may be broken [7][8][9]. The development of an effective field theory framework that describes spacetime-symmetry violations makes comparisons between vastly different kinds of tests possible, generalizing older kinematical test frameworks with a modern viewpoint [10][11][12].
In this article, we discuss the derivation of the polarization-dependent dispersion of GW due to Lorentz and Charge-Parity-Time reversal (CPT) symmetry breaking. We describe the implementation of the modified GW strain in the LIGO-Virgo [2,3] algorithm library, LALSuite [52], as well as the statistical method used to infer the posterior probability of the coefficients for symmetry-breaking. In order to link the theoretical derivation to the analysis of astrophysical signals, we provide detailed explanations of the steps necessary to measure the coefficients for CPT and Lorentz violation, alongside simulations of the modified signals and studies of the sensitivity of current GW interferometers for parameter inference.
The layout of the article is as follows. In Section 2, we describe the theoretical methodology for effective field theory descriptions of local Lorentz violation, including a scalar field example, and an effective quadratic action for the spacetime metric fluctuations. This Section also includes the discussion of the modified plane wave solutions, and the conversion of various expressions to Système International (SI) units. Following this, in Section 3 we describe the implementation of the modified GW signals and the statistical method used for the inference of the coefficients controlling the Lorentz and CPT-breaking effects on propagation. Section 4 includes simulations for a particular subset of the possible forms of Lorentz and CPT violations. A summary and outlook is included in Section 5.
For the bulk of the paper, we work in natural units whereh = c = 1 and Newton's gravitational constant is G N = 1, except when we explicitly write some expressions in SI units. Our convention is to work with Greek letters for spacetime indices and latin letters for spatial indices. The flat spacetime metric signature aligns with the common General Relativity related convention − + ++.

Background and General Relativity
As in the typical gravitational wave scenario, we expand the spacetime metric g µν around a flat Minkowski background η µν as Far from the source at the detectors, GWs are treated as perturbations h µν around the Minkowski metric where h µν << η µν (e.g., components on the order of 10 −21 ). However, one does not assume h µν is small compared to unity in all regions. In particular, in solving for the complete solution in the far radiation zone, one needs to solve in the near zone as well, for example in a post-Newtonian series [53,54]. In standard General Relativity, one solves the Einstein field equations for the metric; The form in (1) is a rewritten form, not yet a specific solution. The full Einstein field equations can be written in the "relaxed" form, as where (T M ) µν is the matter stress-energy tensor, and τ µν is the energy-momentum pseudo tensor [55], and κ = 8πG N . Note that in this equation, (G L ) µν is the Einstein tensor linearized in h µν .
In the wave zone, where the gravitational fields are weak, the equation (2) becomes simply the "vacuum" equations (G L ) µν = 0, which admit wave solutions with two transverse degrees of freedom after choosing a gauge. The transverse-traceless gauge (TT-gauge) is used to describe the propagation of GWs; in this gauge, GR predicts two linearly independent polarizations labelled "+" and "×", with a phase angle difference of π/4, The observable signal comes from the LIGO and Virgo detector responses to the incoming GW, where F A, * are the antenna response patterns of the detectors. The angles θ and φ are the source sky locations, τ is the time delay between detectors receiving the signal, and ψ is the GW frame rotation with respect to the detectors' frame. Note the individual, polarization terms above are not gauge independent, as they depend on ψ, yet the entire observed signal is gauge independent. We return to this point below when discussing Lorentz violation effects for GWs.

Scalar field example
To help understand this framework, and how the action is developed, we consider first a scalar field action in flat spacetime. A free massless scalar field Φ is described by the action When varying this action Φ → Φ + δΦ one obtains, to first order in δΦ, and applying the Leibniz rule, Because of the total derivative, the first term is total 4 divergence and hence is normally considered a surface term, to be evaluated on the 3 dimensional hypersurface Σ bounding the volume of spacetime considered. Since the variational principle in field theory normally assumes the variation δΦ vanishes on the boundary, this term vanishes. What is left is proportional to the arbitrary variation δΦ and therefore if δI = 0 is imposed we obtain the field equations: where = ∂ α ∂ α .
In the effective field theory framework description of Lorentz violation, terms are added to the action (5) that are formed from contractions of general background coefficients with arbitrary numbers of indices k µνλ... and terms involving the scalar field like ∂ µ Φ∂ ν Φ. This is based on the premise that any form of Lorentz violation can be described by the coupling of known matter fields to a fixed background field k µνλ... [10,11]. Under particle Lorentz transformations, the matter fields transform as tensors, while the background field remain fixed. On the other hand, under observer transformations, both background and matter fields transform. The latter condition reflects the idea that physics should be independent of coordinates. These concepts are detailed in the literature. Most notably, see references [58,59] for illustrations in classical mechanics contexts.
There are several treatments of the origin of Lorentz violation that can play a role in the phenomenology of the effective field theory test framework (SME). The Lorentz violation can be explicit, in which the coefficients are prescribed, a priori unknown background fields, unaccompanied by additional dynamical modes. On the other hand, a more elegant mechanism of spontaneous Lorentz-symmetry breaking can be considered. In this latter case, the underlying action for the model is Lorentz invariant but through a dynamical process, nonzero vacuum expectation values for tensor fields can arise [7]. Other scenarios with alternative geometries like Riemann-Finsler geometry have been explored [60][61][62][63][64][65]. Much theoretical discussion of these topics exists in the literature [12,[66][67][68][69][70][71], but we do not delve into details here.
As a simple start, one might consider trying to add a vector coupled to a first derivative of the scalar to (5), as in for arbitrary background vector k ν (we assume the explicit symmetry breaking case for the moment). However, this can be shown to be equivalent to a surface term: where d 3 Σ ν is the hypersurface "area" element. Since the variation δΦ is assumed to vanish on the hypersurface, this contribution will vanish from the field equations. Alternatively, variation of (8) yields a null result more directly: where the last line identically is zero.
To obtain Lorentz-violating terms that yield physical results we modify the action in ( 5) as where k µν are the coefficients for Lorentz violation [11,72], containing 10 independent coefficients describing the degree of Lorentz violation. Note that we assume here that the coefficients are constants in the chosen coordinate system (i.e., the partials vanish, ∂ α k µν = 0). Upon variation, as in (5), we obtain the modified field equations To complete the discussion here we also consider the plane wave solutions to (12). This is achieved by assuming Φ takes the form Φ = Ae ip µ x µ , where x µ is spacetime position and p µ = (ω, p) is the four-momentum for the plane wave. This yields the momentum-space equation Using the definition of the four-momentum we can write this out in a space and time decomposed form: We can solve for the dispersion relation ω( p) and then expand the result to leading order in the coefficients k µν . We obtain This dispersion would modify the propagation of the scalar mode, in particular its speed v = ω/| p| can be written as Note the directional dependence of the speed due to the anisotropic coefficients k 0j and k ij . Even in the case of the isotropic limit, where only k 00 appears, due to the observer Lorentz covariance, this limit is is a special feature of a particular observer frame. For example when viewed by an observer boosted by small β j , anisotropic terms will arise ( e.g., (k ) 0j ∼ −β j k 00 ).
In the typical effective field theory treatment of searches for Lorentz violation, additional, "higher order", terms are also included [73]. Thus the Lagrange density takes the form where now the coefficients are labeled d for the mass dimension of the term in the action, with the scalar field itself having mass dimension 1 and each derivative introducing a mass dimension M 1 . Thus, the result in (11) is the d = 4 limit and the coefficients k µν are dimensionless. In general the coefficients (k (d) ) µνλ... have mass dimension M 4−d .

Gravity sector case
The action from the gravity sector that includes both linearized Lorentz invariant and Lorentz-violating terms, can be described similarly to the scalar case. However, with a multicomponent field, the details of the tensor algebra are more complicated. First we note that linearized General Relativity can be derived from the action where the Einstein tensor is expressed in linearized form with terms of order h 2 and higher discarded. Note that an action quadratic in h µν yields field equations linear in h µν . We now explain in some detail, the construction outlined in Ref. [14]. The starting point for an action that generalizes (18) is whereK (d)µνρσ is an operator given bŷ The operator contains partial derivatives that act on the gravitational field fluctuations h µν ; the K (d)µνρσ 1 ... d−2 are a set of constants in the chosen coordinates. The mass dimension label d refers to the natural units of mass that each term has. At this stage, the nature of these constants is unknown and in what follows we explain the conditions applied to constrain them. One derives the field equations via variation of the action with respect to the fields, similar to the scalar example above. Varying the action (19) with respect to the metric fluctuations h µν , yields In order to completely factor out the variation of the metric field δh µν , integration by parts is performed on the second term. (Note that in doing the integration by parts, we discard surface terms with derivatives of the fluctuations, which is a nontrivial step reflecting the fact that the action contains an arbitrary number of derivatives, going beyond the usual first order derivative form of conventional dynamics.) When d is even, the integration by parts is done an even number of times, creating an overall positive value for the term; if d is odd, the over term is negative in value. We can represent this with (−1) d , and then obtain Since h µν is symmetric, we can indicate the symmetry with parenthesis inK (d)(µν)(αβ) . There are two considerations in (22) to investigate, the first being that only terms contributing to the field equations should survive. Thus we must have The second consideration is the imposition of the linearized gauge symmetry, i.e., h µν → h µν − ∂ µ ξ ν − ∂ ν ξ µ , where ξ µ is an arbitrary vector. 1 If we apply this transformation on the metric within the action (19), i.e., δ ξ h µν = −∂ µ ξ ν − ∂ ν ξ µ , we obtain, from (22), Since ξ µ is arbitrary and derivatives of h ρσ are not necessarily zero, the second condition becomes Under these two conditions (23) and (25), there are three categories of coefficients. These categories are based in part on discrete spacetime symmetry properties of the terms in the action: their behavior under CPT transformations, for which they can be even or odd. Also, the possible tensor index symmetries categorize these coefficients [14]. The three types ofK (d)µνρσ "hat" operators are written aŝ Theŝ operators have even CPT and mass dimension d ≥ 4;q operators have odd CPT and mass dimension d ≥ 5;k operators have even CPT and mass dimension d ≥ 6. The process also reproduces the GR terms. The Lagrange density is then where the first term is an equivalent way of writing the standard GR using the totally antisymmetric Levi-Civita tensor density µρακ (equivalent to (18)). It should be remarked at this point that the Lagrange density in (27) is the most general one constructed purely from the metric fluctuations h µν and taken to quadratic order only. While it includes only constant coefficients in (26), it maintains linearized gauge symmetry. Terms in this Lagrange density can arise in spontaneous-symmetry breaking models, when the additional fluctuations (including possible Nambu-Goldstone and massive modes around the vacuum values have been "integrated out" or "de-coupled" [13,[74][75][76]. 2 On the other hand, examples exist where the quadratic order Lagrange density in (27) can arise from models with explicit symmetry breaking. In either scenario, one is then left with an "effective" Lagrange density, quadratic in the metric fluctuations around a flat background, in which the fluctuations do not appear. Proceeding, the resulting vacuum field equations from (27) are In the absence of Lorentz violation, the field equations (28) reduce to G µν = 0. In the Lorentz gauge, this reduces to h µν = 0, whereh µν = h µν − (1/2)η µν h α α and ∂ µh µν = 0. For plane wave solutionsh µν = A µν e −ip α x α , this yields p 2 = p α p α = 0. This provides the dispersion relation for GR, the equation of motion in energy-momentum space which describes the propagation for GWs. Using the residual gauge freedom in this limit, the number of independent components of the plane wave solutions can be reduced to 2 and will take the form of (3) in the Transverse-Traceless gauge. To find the dispersion relation for the modified equations (28), one again assumes a plane wave form above. There are then at least two approaches to solving the resulting equations, where the components of h µν appear highly coupled with one another due to the extra symmetry-breaking terms in (28). The equations (28) retain the usual gauge freedom, and so one can proceed by choosing a gauge condition and then decomposing the resulting equations into time and space components. For example, using a temporal-type gauge h 0µ = 0 and a helicity basis for the spatial components, one can show that to first order in the coefficients for Lorentz violation, still only 2 degrees of freedom remain [26].
Alternatively, a gauge-independent method for deriving the dispersion relation that uses differential forms exists [73].
Despite the fact that only two physical propagating degrees of freedom remain in the leading order Lorentz violation case, the two modes generally travel at different speeds in the vacuum, resulting in birefringence, and the frequencies of the modes are highly dispersive. (Note that in contrast, for the scalar field example in (15), there is no birefringence effect because there is only one scalar mode whose propagation is modified.) With a helicity basis choice of spatial coordinates, the two propagating modes can be shown to lie in the +2 and −2 helicity projections of the spatial components of the metric fluctuations h ij . The modified dispersion relation can be written as where and All of the derivative factors ∂ µ from (26) are replaced with momenta ∂ µ → ip µ . The plus and minus signs indicate the different dispersion relations for each propagating mode, in vacuum (birefringence). Note that the dispersion and birefringence effects depend on the arrival direction of the plane wavep, revealing this to be a fundamentally anisotropic effect that differs from kinematical isotropic descriptions of symmetry breaking [79].

Gravitational wave signals
Since the terms involving the coefficients in (30) are already at leading order, they can be evaluated with the zeroth-order solution (e.g., p µ = ω(1,p) = | p|(1,p)). This reveals that any effects associated with arriving plane waves should depend on angular functions of the unit vectorp. Further, since LIGO-Virgo analysis use angular sky map coordinate systems, it is advantageous to use the machinery of spherical harmonics and spherical tensors. We can decompose the above coefficients into spherical harmonic form, In these expressions, Y jm (n) are the usual spherical harmonics withn = −p, while ±4 Y jm (n) are spin-weighted spherical harmonics. The coefficients, formerly in cartesian tensor form in (30) are expressed in spherical form k (V)jm , where j = 0, 1, ..., d − 2 and −j ≤ m ≤ j. The meaning of the subscripts I, E, B, V and the relation between the two forms of the coefficients is determined by whether the terms are CPT odd or even and which mass dimensions they encompass, detailed in Refs. [14,73,80].
In GR, there is no difference in the speed between gravitational wave polarizations; both travel at the speed of light (i.e., v = ω/| p| = 1). In the case of a Lorentz violation in the form in (30), the speed of the waves is given by Given enough propagation distance from source to detector, a difference in arrival times may be detectable even for small Lorentz violation, a feature that has been used for photon tests of Lorentz invariance [81][82][83][84][85][86][87]. Using LVC data, we can test for these effects by looking for a phase deviation from GR via polarization comparisons. If Lorentz violation effects are not resolvable given current precision, we can then provide constraints for the LV coefficients.
Modifications in the analysis code use the expressions for the gravitatonal wave strain polarizations. The plane wave solutions will have a phase shift δψ ± due to terms in (36) or (30). Consider first the strain h ∼ e −i(ωt−kl) where l is the distance travelled and k is the wave number. The difference in phase grows in magnitude the farther the gravitational wave travels from the source to detectors. On cosmological scales, it is important to include effects on propagation time from the expanding universe using luminosity distance. Noting k ∼ | p| = ω/v, inputting (36), and including distance and frequency alterations form cosmology, one finds the phase shift expression, where H(z) is the Hubble parameter with redshift z and the observed frequency is related to that emitted via ω obs (1 + z) = ω emit . For each mode, the modified phase shift can be written as where The τ is the effective propagation time due to cosmological redshift z.
It is useful to rewrite the coefficients in terms of effective angles ϑ and ϕ defined by Note that these angles are not the sky location angles θ, and φ. Also using the plus and cross polarizations (3), the modified gravitational wave solutions in terms of the Lorentz-invariant solutions can be written The h LI (+,×) are the Lorentz-invariant gravitational wave for standard GR; one can retrieve GR as a limiting case as β → 0 and δ → 0.
The measured signal at a given detector can be obtained from an equation of the form (4). It is standard in the literature to adopt a Sun-centered Celestial-Equatorial coordinate system (or SCF frame) for reporting measurements of the components of the coefficients for Lorentz violation either in the form s TXY... , ... or in spherical tensor form k (d) (I)10 , ... [6]. Under observer coordinate transformations, the coefficients transform as tensors. In many cases, these transformations can be implemented as global Lorentz transformations on the coefficients. In the present case, we want to ensure the coefficients in the expression for the measured strain are all expressed in terms of the SCF coefficients, thereby leaving any angular, sky location, dependence in the relevant angular variables. Thus, when analyzing data, the signal generically will have extra angular, isotropy breaking, dependence on the sky angles. This will differ significantly from the GR case.

Unit changes and dimension
For applications below, it becomes essential to convert from natural units to SI units when implementing modifications into analysis code. We note here several useful unit substitutions that can used for this and various key equations discussed previously.
Recall natural units are based onh = c = 1. In these units, quantities can have dimensions of energy, typically expressed in terms of electron volts, as GeV= 10 9 eV, for example. For instance, mass dimension d coefficients for Lorentz violation have units of M 4−d . To convert various quantities to SI units, we assume that the starting action has units of joules meters Jm. 3 For instance, the full Einstein Hilbert action in SI units can be written as or for the quadratic action limit of equation (18), simply multiply by c 4 . Units of kg m s −2 come from the factor c 4 G , m 4 from d 4 x, and m −2 from the derivatives contained within the Einstein tensor. Implicit here is the assumption that the metric tensor g µν is dimensionless (the Minkowski metric retains its form η µν = diag (−1, 1, 1, 1)). Likewise, the Lorentzviolating action (19) contains operators with SI units m −2 , and thus from (20), when introducing higher derivatives, the units of the coefficients compensate, thus the coefficients have units m d−4 .
When converting the field equations (28) from position to momentum space, every partial derivative contributes a factor with Planck's constant, i.e., ∂ α → ī h p α . Schematically, the the position space equation has the form where, e.g., for d = 4 a term involving theŝ operators contains coefficients for s (4) coupled to two derivatives that act on h. In momentum space, where the operatorsŝ,q, andk now contain ( ī h ) (d−2) p α 1 ...p α d−2 in place of partials. The units for the coefficients are unchanged, i.e., m d−4 .
One must also keep track of the corrected time-component factors in the four-momenta, p α = (−¯h c ω, p). For instance, the wave speed, via the dispersion relation, becomes where each ζ quantity in (32) inherits a factor of (¯h c ) 2 . To ensure the coefficients, k with similar SI factors for ζ 1 , ζ 2 , ζ 3 .

Analysis method
The coefficients for Lorentz and CPT violation can be measured from the comparison of the speed of gravitational and electromagnetic waves, an analysis that has been performed with gravitational-wave event GW170817 and the associated counterpart gamma-ray burst (GRB) GRB170817A to constrain coefficients of mass dimension 4 with improved accuracy [47]. Using GW signals only, limits on mass dimensions 5 and 6 coefficients have been obtained from the non-observation of a delay between the arrival time of the h + and h × polarizations in the LIGO and Virgo interferometers [14,49,50].
The constraints on the birefringence parameters are obtained from the posterior samples inferred under the assumption of no symmetry breaking, and are limited by the detector resolution to determine the waveform peak frequency, focusing on information from signals at higher frequencies. We aim to compliment prior work by analysing directly the LVC interferometers strain in order to bypass the reliance on posterior parameters inferred under a GR model. Our analysis therefore fully takes into account the correlation between the SME coefficients and the source parameters, including dispersion or birefringence effects, during the inference process.
We have implemented the modification of the GW strain obtained in (42) to estimate the coefficients for symmetry-breaking from the morphology of the signals. As the dispersive and birefringent effects are degenerate with the source properties (e.g. the luminosity distance, due to the additional energy loss during the propagation) we perform a joint estimation of the source parameters and the coefficients for Lorentz and CPT violation taking into account the modifications at all frequencies of the waveform. We implement the Bayesian analysis into a version of the LIGO Algorithm Library suite LALSuite modified for our purposes as described below [52].

Implementation of the modified waveform
The joint measurement of source and beyond-GR constraints have been performed for a variety of new physics parameterizations, including modification of the GW generation and propagation [4]. Following a similar methodology, we implement the modifications of the GW signals derived from the SME framework in the GW simulation package of LALSuite. Such deformations can be anisotropic, as can be infered from the appearance ofn in Eq. (42) via β and δ. Here we focus on the simplest coefficients that produce dispersion and birefringence via Lorentz and CPT violating effects,i.e., those of mass dimension 5. These coefficients are contained within β in (42) and obey the complex There are a priori independent coefficients in this set of terms [26]. We display the first terms within β in SI units: where the sky location of the source (θ, φ) appears. The coefficients are taken as expressed in the SCF in this expression. The general form of the signal observed in the interferometer is: where h (+,×) are the expressions (42), and F (+,×) are the (standard) detector response functions. The expressions for F (+,×) include the rotation angles relating different frames, e.g., the merger frame and the detectors frame as defined in the LALSuite software. The effective propagation time τ parameter of equation (48) is defined in equation (  40) as an integral function of the redshift. Since it needs to be evaluated for every value of the SME coefficients being tested, for computing time feasibility we instead probe the effective coefficient (k (5) (V)jm ) e f f = τ k (5) (V)jm . The value of the SME coefficient is recovered after convergence of the inference process further described in the following Section.
Finally we note that transformations of the coefficients under observer boosts are also computable. This would be important, should it become necessary to include the motion of the Earth, the interferometers, or the motion of a source system's center of mass, relative to the SCF. Currently, it appears the strain measurements are not sensitive to this level of nonrelativistic boosts (e.g. v/c = 10 −4 ).

Bayesian analysis
After implementing the modification of the strain, we include the SME coefficients in LALInference, the parameter estimation package of LALSuite [88]. LALInference performs Bayesian inference of the posterior probability of the GW source parameters with the inclusion of the systematic uncertainties due to the detectors resolutions. The vector set of GR prior parameters, θ GR , includes intrinsic parameters describing the binary system (e.g. the black holes masses and spins) as well as extrinsic parameters placing it in the astrophysical environment (e.g. the sky location, distance, inclination). We add to the preexisting parameters the SME coefficients (k (5) (V)jm ) e f f described in Section 3.1 for the mass dimension 5 case, contained within θ SME .
In order to include the correlation between the GR parameters and the SME coefficients, we perform a simultaneous inference of all the parameters, obtaining the joint posterior probability: P( θ GR , θ SME |d, I) = P(d| θ GR , θ SME ,I) P( θ GR , θ SME |I) where P( θ GR , θ SME |d, I) is the posterior probability, P(d| θ GR , θ SME , I) the likelihood, P( θ GR , θ SME |I) the prior probability and P(d|I) the evidence, and any pertinent background information is included in I. We set a flat prior probability for (k (5) (V)jm ) e f f bounded between |(k (5) (V)jm ) e f f | ∈ [0; 10 −10 ], with maximal value well above the existing constraints on the order of 10 −15 [49]. The likelihood is computed in the frequency domain: whereh i is the frequency-domain template signal,d i are the data observed by the interferometers, T is the duration of the signal, and S n the power spectral density of the detector noise.
Due to the large number of parameters describing the GW emitted by the coalescence of binary systems, the posterior probability is infered with Markov Chain (MC) methods. The chains perform semi-random walks in the parameter space where the recorded steps of the walks are proportional to the quantity in Eq. (50). Different algorithms have been shown to be able to perform parameter inference, of which Markov Chain Monte-Carlo (MCMC) with parallel tempering and nested sampling are implemented in the LVC algorithm library. The method returns joint posterior probabilities of the GR parameters and the SME coefficients. From this we extract the marginalised posterior probability on a subset of parameters by integrating over the distribution of the other variables. The credible intervals are finally obtained by summing the volume of the posterior probability corresponding to the desired fraction of confidence. We present the results of Bayesian inference on simulated signals in Section 4, and will provide the results of the ongoing analysis of LVC detections in a separate publication.

Sensitivity study
As an illustration, we assume for the following, one non-zero coefficient k (5) (V)00 corresponding to isotropic polarization-dependent dispersion. Figure 1 plots the waveforms for both GR and the modified wave form for different values of k (d) (V)00 . We assume a non-spinning binary system that has a luminosity distance of 4 Gpc, and equal masses of m 1 = m 2 = 50M . Note that significant differences in the waveform shape occur for coefficient values as small 10 −13 m, impacting both the amplitude and frequency of the signal. This result can be compared with simulations using analytical template models presented in Ref. [26]. In the latter publication in Figures 1 and 2, simulated waveforms with non-zero coefficients for Lorentz and CPT violation appear to modify the waveform mostly around peak amplitude times, whereas the simulations here in Figure 1 show modification at earlier times. Using the methodology outlined in Section 3.1, we perform a Bayesian inference of the source parameters and the coefficients for Lorentz-violation with simulated dispersed signals in order to study the potential to measure the coefficients with the LVC detections. We simulate a GW emitted by a non-spinning binary system of black holes with symmetric masses of 50M located at 5 Gpc where the dispersion is controlled by one coefficient set to a value of k  These results, obtained with a single event, present encouraging prospects towards the measurement of the coefficients for symmetry-breaking with the current generation of GW interferometers. The second catalog of GW detections, encompassing the two first observing run as well the first half of the third observation run of the LVC, contains 50 events from the coalescence of binary systems of astrophysical compact objects, of which 46 are consistent with black hole systems [89]. Comparing our results with measurements of the mass of the graviton, that also induce a modified dispersion of the GW signal, we note that the constraint has been improved of one order of magnitude from a single event to the analysis of a larger population of GW detections. The constraint from GW159014 was m g ≤ 1.2 · 10 −22 eV/c 2 while it is now m g ≤ 1.76 · 10 −23 eV/c 2 when analysing 33 events from the second GW catalog [2,4]. Based on such results, we can conjecture that the constraints on the SME coefficients from the full catalog of GW detections will provide more stringent measurements than the preliminary sensitivity study shown here. The robustness of such measurements due to the waveform modelling approximant has been explored in [3], showing that those systematic uncertainties do not lead to a large bias nor re-estimation of the constraints at the current detector sensitivity. Other studies show that transient noise may impact the measurement [90] by mimicking a GR deviation, an effect that we palliate by using the LVC-released power spectral densities and frequency ranges that exclude the presence of glitches in the strain data.

Conclusions and Future Work
We describe the implementation of an effective field theory framework for testing Lorentz and CPT symmetry into a version of the LIGO-Virgo Algorithm Library suite LALSuite. The Lorentz-and CPT-violating modifications include the coefficients controlling birefringence and dispersion effects on the gravitational wave polarizations. This work does not rely on posterior results inferred by the LVC that assume no deviations from standard GR; we implement the modifications due to dispersion directly at the level of the templates used for the Bayesian inference of the GW source and propagation parameters in order to incorporate the full information provided by the signal morphology.
Initially, one starts with the action in the effective field theory framework that is quadratic in the metric fluctuations h µν , (19), and after theoretical constraints including gauge invariance, we arrive at the general result in (27). From the field equations (28), a dispersion relation is derived, both in terms of the coefficients from the Lagrange density (30), and in terms of spherical coefficients in a special observer frame (33)- (35). The result shows birefringence and dispersion for the two propagating modes; moreover these effects will vary with sky location of the source. Then considering the expression for propagating and applying a modified phase shift, including cosmological considerations, one can rewrite the expressions for the plus and cross polarizations (42), which are directly implemented within the modified package. Through Bayesian inference, we can perform a parameter estimation to constrain the coefficients for Lorentz violation Samples of visible effects are shown in the sensitivity plots in Section 4.
The theoretical derivations and sensitivity studies presented in this article precede the measurement of SME coefficients with the events detected by the LVC. This computationally intensive analysis is currently ongoing and the results will be reported in a future publication, where we aim to fulfil our analysis for coefficients for Lorentz and CPT violation of mass dimension five and six, with a global analysis. In such a global analysis, the availability of what is now a plethora of GW sources across the sky has the potential to disentangle measurements for a large set of coefficients and thereby obtain an exhaustive search for signals of new physics.