Analysis of Birefringence and Dispersion Effects from Spacetime-Symmetry Breaking in Gravitational Waves

Abstract

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-wave 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, and we discuss 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 sensitively probe for these signals of new physics.

1. 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 spacetime-symmetry 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].

The specific consequences of local Lorentz-symmetry breaking for GWs have been studied in several works, within a general effective field theory framework [13,14,15,16,17,18], and in specific models [19,20,21,22,23,24,25]. In particular, the effects on propagation have been determined for generic Lorentz-violating terms in the linearized gravity limit [26], which is the focus in this work.

Examples of searches for Lorentz violation in gravity include table-top tests like gravimetry [27,28,29,30,31,32,33], short-range gravity tests [34,35,36], near-Earth tests [37,38,39,40], solar system planetary tests [41,42,43], and astrophysical tests with pulsars [44,45,46]. Measurements of simultaneous gravitational and electromagnetic radiation have yielded limits on certain types of Lorentz violation in gravity versus light [47]. Recent work has begun to look at the available GW catalog to place constraints on coefficients describing Lorentz violation for gravity [48,49,50]. Additionally, closely related searches for parameterizations of deviations from General Relativity have been completed [3,4,24,51].

In this article, we discuss the derivation of the polarization-dependent dispersion of GWs 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], and we describe 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 a 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 where $\hslash =c=1$, and Newton’s gravitational constant is $G_N\ne 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 $-+++$.

2. Theoretical Framework

2.1. Background and General Relativity

As in the typical gravitational wave scenario, we expand the spacetime metric $g_{\mu \nu }$ around a flat Minkowski background ${\eta }_{\mu \nu }$ as

$$g_{\mu \nu }={\eta }_{\mu \nu }+h_{\mu \nu }.$$

(1)

Far from the source at the detectors, GWs are treated as perturbations $h_{\mu \nu }$ around the Minkowski metric where $h_{\mu \nu }<<{\eta }_{\mu \nu }$ (e.g., components on the order of ${10}^{-21}$). However, one does not assume $h_{\mu \nu }$ 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

$${\left(G_L\right)}^{\mu \nu }=\kappa [{\left(T_M\right)}^{\mu \nu }+{\tau }^{\mu \nu }],$$

(2)

where ${\left(T_M\right)}^{\mu \nu }$ is the matter stress-energy tensor, ${\tau }^{\mu \nu }$ is the energy-momentum pseudo tensor [55], and $\kappa =8\pi G_N$. Note that, in this equation, ${\left(G_L\right)}^{\mu \nu }$ is the Einstein tensor linearized in $h_{\mu \nu }$.

In the wave zone, where the gravitational fields are weak, Equation (2) becomes simply the “vacuum” equations ${\left(G_L\right)}^{\mu \nu }=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 $\pi /4$:

$$h_{\mu \nu }=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& h_{+}& h_{\times }& 0\\ 0& h_{\times }& -h_{+}& 0\\ 0& 0& 0& 0\end{array}\right).$$

(3)

The observable signal comes from the LIGO and Virgo detector responses to the incoming GW,

$$S_A(t,\theta ,\varphi ,\psi )=F_{A,+}(\theta ,\varphi ,\psi )h_{+}(t,\theta ,\varphi ,\psi ,\tau )+F_{A,\times }(\theta ,\varphi ,\psi )h_{\times }(t,\theta ,\varphi ,\psi ,\tau ),$$

(4)

where $F_{A,*}$ are the antenna response patterns of the detectors. The angles $\theta$ and $\varphi$ are the source sky locations, $\tau$ is the time delay between detectors receiving the signal, and $\psi$ is the GW frame rotation with respect to the detectors’ frame. Note that the individual, polarization terms above are not gauge-independent, as they depend on $\psi$, yet the entire observed signal is gauge-independent. We return to this point below when discussing Lorentz violation effects for GWs.

2.2. Spacetime-Symmetry Breaking Scenario

We consider the observable effects on gravitational wave propagation subject to Lorentz- and CPT-breaking terms in an effective field theory framework known as the Standard-Model Extension (SME) [10,11,12,13,14,56,57].

2.2.1. 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 $\Phi$ is described by the action

$$I_{\mathrm{sc}}=-\frac{1}{2}\int d^{4}x{\eta }^{\mu \nu }({\partial }_{\mu }\Phi )({\partial }_{\nu }\Phi ).$$

(5)

When varying this action $\Phi \to \Phi +\delta \Phi$, one obtains, to the first order in $\delta \Phi$, and applying the Leibniz rule,

$$\begin{array}{ccc}\hfill \delta I_{\mathrm{sc}}& =& -\int d^{4}x{\eta }^{\mu \nu }({\partial }_{\mu }\delta \Phi )({\partial }_{\nu }\Phi )\hfill \\ & =& -\int d^{4}x[{\partial }_{\mu }\left(\delta \Phi ({\partial }^{\mu }\Phi )\right)-\delta \Phi {\partial }_{\mu }{\partial }^{\mu }\Phi ].\hfill \end{array}$$

(6)

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 three-dimensional hypersurface $\Sigma$ bounding the volume of spacetime considered. Since the variational principle in field theory normally assumes that the variation $\delta \Phi$ vanishes on the boundary, this term vanishes. What is left is proportional to the arbitrary variation $\delta \Phi$; therefore, if $\delta I=0$ is imposed, we obtain the field equations:

$$\square \Phi =0,$$

(7)

where $\square ={\partial }_{\alpha }{\partial }^{\alpha }$.

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^{\mu \nu \lambda \dots }$ and terms involving the scalar field such as ${\partial }_{\mu }\Phi {\partial }_{\nu }\Phi$. 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^{\mu \nu \lambda \dots }$ [10,11]. Under particle Lorentz transformations, the matter fields transform as tensors, while the background field remains 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 Refs. [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, a priori unknown background fields, are prescribed and 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 such as 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 Equation (5), as in

$$\Delta I_{\mathrm{sc}}=\int d^{4}x{k}_{\nu }(\Phi {\partial }^{\nu }\Phi )$$

(8)

for an arbitrary background vector $k_{\nu }$ (we assume the explicit symmetry breaking case for the moment). However, this can be shown to be equivalent to a surface term:

$$\begin{array}{ccc}\hfill \Delta I_{\mathrm{sc}}& =& \frac{1}{2}\int d^{4}x{k}_{\nu }{\partial }^{\nu }\left({\Phi }^2\right)\hfill \\ & =& \frac{1}{2}\int d^3{\Sigma }_{\nu }{k}^{\nu }{\Phi }^2,\hfill \end{array}$$

(9)

where $d^3{\Sigma }_{\nu }$ is the hypersurface “area” element. Since the variation $\delta \Phi$ 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:

$$\begin{array}{ccc}\hfill \delta \Delta I_{\mathrm{sc}}& =& \int d^{4}x{k}_{\nu }\left[\delta \Phi {\partial }^{\nu }\Phi +\Phi \delta ({\partial }^{\nu }\Phi )\right]\hfill \\ & =& \int d^{4}x{k}_{\nu }({\partial }^{\nu }\Phi -{\partial }^{\nu }\Phi )\delta \Phi ,\hfill \end{array}$$

(10)

where the last line identically is zero.

To obtain Lorentz-violating terms that yield physical results, we modify the action in (5) as

$$I_{\mathrm{sc}}=-\frac{1}{2}\int d^{4}x\left({\eta }^{\mu \nu }({\partial }_{\mu }\Phi )({\partial }_{\nu }\Phi )+({\partial }_{\mu }\Phi )k^{\mu \nu }{\partial }_{\nu }\Phi )\right),$$

(11)

where $k^{\mu \nu }$ 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, ${\partial }_{\alpha }{k}^{\mu \nu }=0$). Upon variation, as in (5), we obtain the modified field equations

$$\square \Phi +k^{\mu \nu }{\partial }_{\mu }{\partial }_{\nu }\Phi =0.$$

(12)

To complete the discussion here, we also consider the plane wave solutions to (12). This is achieved by assuming that $\Phi$ takes the form $\Phi =A{e}^{i{p}_{\mu }{x}^{\mu }}$, where $x^{\mu }$ is the spacetime position and $p^{\mu }=(\omega ,\overrightarrow{p})$ is the four-momentum for the plane wave. This yields the momentum-space equation

$$p_{\mu }{p}^{\mu }+k_{\mu \nu }{p}^{\mu }{p}^{\nu }=0.$$

(13)

Using the definition of the four-momentum, we can write this out in a space and time decomposed form:

$${\omega }^2(1-k_{00})-2k_{0j}{p}^j\omega -k_{ij}{p}^i{p}^j-{\overrightarrow{p}}^2=0.$$

(14)

We can solve for the dispersion relation $\omega \left(\overrightarrow{p}\right)$ and then expand the result to the leading order in the coefficients $k_{\mu \nu }$. We obtain

$$\omega \approx |\overrightarrow{p}|\left(1+\frac{1}{2}(k_{00}+2k_{0j}{\widehat{p}}^j+k_{ij}{\widehat{p}}^i{\widehat{p}}^j)\right).$$

(15)

This dispersion would modify the propagation of the scalar mode; in particular, its speed $v=\omega /|\overrightarrow{p}|$ can be written as

$$v\approx 1+\frac{1}{2}(k_{00}+2k_{0j}{\widehat{p}}^j+k_{ij}{\widehat{p}}^i{\widehat{p}}^j).$$

(16)

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 ${\beta }^j$, anisotropic terms will arise (e.g., ${\left(k^{\prime }\right)}_{0j}\sim -{\beta }^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

$$I_{\mathrm{sc}}=-\frac{1}{2}\int d^{4}x\left({\eta }^{\mu \nu }({\partial }_{\mu }\Phi )({\partial }_{\nu }\Phi )+({\partial }_{\mu }\Phi )\sum _d{\left(k^{\left(d\right)}\right)}^{\mu \nu \lambda \dots }({\partial }_{\nu }{\partial }_{\lambda }\dots \Phi )\right),$$

(17)

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 mass dimension $M^1$. Thus, the result in (11) is the $d=4$ limit, and the coefficients $k^{\mu \nu }$ are dimensionless. In general, the coefficients ${\left(k^{\left(d\right)}\right)}^{\mu \nu \lambda \dots }$ have mass dimension $M^{4-d}$.

2.2.2. 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 the linearized General Relativity can be derived from the action

$$I_{GR}=-\frac{1}{4\kappa }\int d^{4}x{h}_{\mu \nu }{G}^{\mu \nu },$$

(18)

where the Einstein tensor is expressed in a linearized form with terms of order $h^2$ and higher discarded. Note that an action quadratic in $h_{\mu \nu }$ yields field equations that are linear in $h_{\mu \nu }$.

We now explain in some detail the construction outlined in [14]. The starting point for an action that generalizes (18) is

$$I=\frac{1}{8\kappa }\int d^{4}x{h}_{\mu \nu }{\widehat{K}}^{\left(d\right)\mu \nu \rho \sigma }{h}_{\rho \sigma },$$

(19)

where ${\hat{K}}^{\left(d\right)\mu \nu \rho \sigma }$ is an operator given by

$${\widehat{K}}^{\left(d\right)\mu \nu \rho \sigma }=K^{\left(d\right)\mu \nu \rho \sigma {ϵ}_1\dots {ϵ}_{d-2}}{\partial }_{{ϵ}_1}\dots {\partial }_{{ϵ}_{d-2}}.$$

(20)

The operator contains partial derivatives that act on the gravitational field fluctuations $h_{\mu \nu }$; $K^{\left(d\right)\mu \nu \rho \sigma {ϵ}_1\dots {ϵ}_{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_{\mu \nu }$ yields

$$\begin{array}{ccc}\hfill \delta I& =& \frac{1}{8\kappa }\int d^{4}x[\delta h_{\mu \nu }{K}^{\left(d\right)\mu \nu \rho \sigma {ϵ}_1\dots {ϵ}_{d-2}}{\partial }_{{ϵ}_1}\dots {\partial }_{{ϵ}_{d-2}}{h}_{\rho \sigma }\hfill \\ & & \phantom{\frac{1}{8}\int d^{4}x}+h_{\mu \nu }{K}^{\left(d\right)\mu \nu \rho \sigma {ϵ}_1\dots {ϵ}_{d-2}}{\partial }_{{ϵ}_1}\dots {\partial }_{{ϵ}_{d-2}}\delta h_{\rho \sigma }].\hfill \end{array}$$

(21)

In order to completely factor out the variation of the metric field $\delta h_{\mu \nu }$, 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

$$\delta I=\frac{1}{8\kappa }\int d^{4}x\delta h_{\alpha \beta }\left[K^{\left(d\right)\left(\alpha \beta \right)\left(\mu \nu \right){ϵ}_1\dots {ϵ}_{d-2}}+{(-1)}^d{K}^{\left(d\right)\left(\mu \nu \right)\left(\alpha \beta \right){ϵ}_1\dots {ϵ}_{d-2}}\right]{\partial }_{{ϵ}_1}\dots {\partial }_{{ϵ}_{d-2}}{h}_{\mu \nu }.$$

(22)

Since $h_{\mu \nu }$ is symmetric, we can indicate the symmetry with parenthesis in ${\widehat{K}}^{\left(d\right)\left(\mu \nu \right)\left(\alpha \beta \right)}$.

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

$$K^{\left(d\right)\left(\alpha \beta \right)\left(\mu \nu \right){ϵ}_1\dots {ϵ}_{d-2}}+{(-1)}^d{K}^{\left(d\right)\left(\mu \nu \right)\left(\alpha \beta \right){ϵ}_1\dots {ϵ}_{d-2}}\ne 0.$$

(23)

The second consideration is the imposition of the linearized gauge symmetry, i.e., $h_{\mu \nu }\to h_{\mu \nu }-{\partial }_{\mu }{\xi }_{\nu }-{\partial }_{\nu }{\xi }_{\mu }$, where ${\xi }^{\mu }$ is an arbitrary vector. General gauge-breaking terms are considered in Ref. [15]. If we apply this transformation on the metric within the action (19), i.e., ${\delta }_{\xi }{h}_{\mu \nu }=-{\partial }_{\mu }{\xi }_{\nu }-{\partial }_{\nu }{\xi }_{\mu }$, we obtain, from (22),

$$\begin{array}{ccc}\hfill {\delta }_{\xi }I& =& \frac{1}{8\kappa }\int d^{4}x{\partial }_{\alpha }{\xi }_{\beta }\left[{(-1)}^d{\widehat{K}}^{\left(d\right)\left(\mu \nu \right)\left(\alpha \beta \right)}+{\widehat{K}}^{\left(d\right)\left(\alpha \beta \right)\left(\mu \nu \right)}\right]h_{\mu \nu },\hfill \\ & =& -\frac{1}{8\kappa }\int d^{4}x{\xi }_{\nu }\left[{(-1)}^d{\widehat{K}}^{\left(d\right)\left(\rho \sigma \right)\left(\mu \nu \right)}+{\widehat{K}}^{\left(d\right)\left(\mu \nu \right)\left(\rho \sigma \right)}\right]{\partial }_{\mu }{h}_{\rho \sigma }.\hfill \end{array}$$

(24)

Since ${\xi }_{\mu }$ is arbitrary, and derivatives of $h_{\rho \sigma }$ are not necessarily zero, the second condition becomes

$$\left[{(-1)}^d{\widehat{K}}^{\left(d\right)\left(\rho \sigma \right)\left(\mu \nu \right)}+{\widehat{K}}^{\left(d\right)\left(\rho \sigma \right)\left(\mu \nu \right)}\right]{\partial }_{\mu }=0.$$

(25)

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. Additionally, the possible tensor index symmetries categorize these coefficients [14]. The three types of ${\widehat{K}}^{\left(d\right)\mu \nu \rho \sigma }$ “hat" operators are written as

$$\begin{array}{ccc}\hfill {\widehat{s}}^{\mu \rho \nu \sigma }& =& s^{\left(d\right)\mu \rho {ϵ}_1\nu \sigma {ϵ}_2\dots {ϵ}_{d-2}}{\partial }_{{ϵ}_1}\dots {\partial }_{{ϵ}_{d-2}},\hfill \\ \hfill {\widehat{q}}^{\mu \rho \nu \sigma }& =& q^{\left(d\right)\mu \rho {ϵ}_1\nu {ϵ}_2\sigma {ϵ}_3\dots {ϵ}_{d-2}}{\partial }_{{ϵ}_1}\dots {\partial }_{{ϵ}_{d-2}}\hfill \\ \hfill {\widehat{k}}^{\mu \nu \rho \sigma }& =& k^{\left(d\right)\mu {ϵ}_1\nu {ϵ}_2\rho {ϵ}_3\sigma {ϵ}_4\dots {ϵ}_{d-2}}{\partial }_{{ϵ}_1}\dots {\partial }_{{ϵ}_{d-2}}.\hfill \end{array}$$

(26)

The $\widehat{s}$ operators have an even CPT and mass dimension $d\ge 4$; $\widehat{q}$ operators have an odd CPT and mass dimension $d\ge 5$; $\widehat{k}$ operators have an even CPT and mass dimension $d\ge 6$. The process also reproduces the GR terms.

The Lagrange density is then

$$\begin{array}{ccc}\hfill \mathcal{L}& =& \frac{1}{8\kappa }{ϵ}^{\mu \rho \alpha \kappa }{ϵ}^{\nu \sigma \beta \lambda }{\eta }_{\kappa \lambda }{h}_{\mu \nu }{\partial }_{\alpha }{\partial }_{\beta }{h}_{\rho \sigma }\hfill \\ & & +\frac{1}{8\kappa }{h}_{\mu \nu }({\widehat{s}}^{\mu \rho \nu \sigma }+{\widehat{q}}^{\mu \rho \nu \sigma }+{\widehat{k}}^{\mu \rho \nu \sigma })h_{\rho \sigma },\hfill \end{array}$$

(27)

where the first term is an equivalent way of writing the standard GR using the totally antisymmetric Levi-Civita tensor density ${ϵ}^{\mu \rho \alpha \kappa }$ (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_{\mu \nu }$ and taken to the quadratic order only. While it includes only constant coefficients in (26), it maintains a 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].1 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

$$0=G^{\mu \nu }-[\frac{1}{4}({\widehat{s}}^{\mu \rho \nu \sigma }+{\widehat{s}}^{\mu \sigma \nu \rho })+\frac{1}{2}{\widehat{k}}^{\mu \nu \rho \sigma }+\frac{1}{8}({\widehat{q}}^{\mu \rho \nu \sigma }+{\widehat{q}}^{\nu \rho \mu \sigma }+{\widehat{q}}^{\mu \sigma \nu \rho }+{\widehat{q}}^{\nu \sigma \mu \rho })]h_{\rho \sigma }.$$

(28)

In the absence of Lorentz violation, the field Equation (28) reduce to $G^{\mu \nu }=0$. In the Lorentz gauge, this reduces to $\square {\overline{h}}^{\mu \nu }=0$, where ${\overline{h}}^{\mu \nu }=h^{\mu \nu }-(1/2){\eta }^{\mu \nu }{h}_{\phantom{\alpha }\alpha }^{\alpha }$ and ${\partial }_{\mu }{\overline{h}}^{\mu \nu }=0$. For plane wave solutions ${\overline{h}}_{\mu \nu }=A_{\mu \nu }{e}^{-i{p}_{\alpha }{x}^{\alpha }}$, this yields $p^2=p^{\alpha }{p}_{\alpha }=0$. This provides the dispersion relation for GR,

$$\omega =|\overrightarrow{p}|,$$

(29)

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 Equation (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_{\mu \nu }$ appear highly coupled with one another due to the extra symmetry-breaking terms in (28). The Equation (28) retain the usual gauge freedom, 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\mu }=0$ and a helicity basis for the spatial components, one can show that, to the 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

$$\omega =|\overrightarrow{p}|\left(1-{\zeta }^0\pm \left|\overrightarrow{\zeta }\right|\right),$$

(30)

where

$$|\overrightarrow{\zeta }|=\sqrt{{\left({\zeta }^1\right)}^2+{\left({\zeta }^2\right)}^2+{\left({\zeta }^3\right)}^2}$$

(31)

and

$$\begin{array}{ccc}\hfill {\zeta }^0& =& \frac{1}{4|\overrightarrow{p}{|}^2}\left(-{\widehat{s}}^{\mu \nu }{}_{\mu \nu }+\frac{1}{2}{\widehat{k}}^{\mu \nu }{}_{\mu \nu }\right),\hfill \\ \hfill {\left({\zeta }^1\right)}^2+{\left({\zeta }^2\right)}^2& =& \frac{1}{8|\overrightarrow{p}{|}^4}\left({\widehat{k}}^{\mu \nu \rho \sigma }{\widehat{k}}_{\mu \nu \rho \sigma }-{\widehat{k}}^{\mu \rho }{}_{\nu \rho }{\widehat{k}}_{\mu \sigma }{}^{\nu \sigma }+\frac{1}{8}{\widehat{k}}^{\mu \nu }{}_{\mu \nu }{\widehat{k}}^{\rho \sigma }{}_{\rho \sigma }\right),\hfill \\ \hfill {\left({\zeta }^3\right)}^2& =& \frac{1}{16|\overrightarrow{p}{|}^4}\left(-\frac{1}{2}{\widehat{q}}^{\mu \rho \nu \sigma }{\widehat{q}}_{\mu \rho \nu \sigma }-{\widehat{q}}^{\mu \nu \rho \sigma }{\widehat{q}}_{\mu \nu \rho \sigma }+({\widehat{q}}^{\mu \rho \nu }{}_{\rho }+{\widehat{q}}^{\nu \rho \mu }{}_{\rho }){\widehat{q}}_{\mu \sigma \nu }{}^{\sigma }\right).\hfill \end{array}$$

(32)

All of the derivative factors ${\partial }_{\mu }$ from (26) are replaced with momenta ${\partial }_{\mu }\to i{p}_{\mu }$. 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 wave $\widehat{p}$, revealing this to be a fundamentally anisotropic effect that differs from kinematical isotropic descriptions of symmetry breaking [79].

2.2.3. Gravitational Wave Signals

Since the terms involving the coefficients in (30) are already at the leading order, they can be evaluated with the zeroth-order solution (e.g., $p^{\mu }=\omega (1,\widehat{p})=\left|\overrightarrow{p}\right|(1,\widehat{p})$). This reveals that any effects associated with arriving plane waves should depend on angular functions of the unit vector $\widehat{p}$. Further, since LIGO-Virgo analysis uses 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 a spherical harmonic form,

$$\begin{array}{ccc}\hfill {\zeta }^0& =& \sum _{djm}{\omega }^{d-4}{Y}_{jm}\left(\widehat{\mathbf{n}}\right)k_{\left(I\right)jm}^{\left(d\right)},\hfill \end{array}$$

(33)

$$\begin{array}{ccc}\hfill {\zeta }^1\mp i{\zeta }^2& =& \sum _{djm}{\omega }^{d-4}{}_{\pm 4}{Y}_{jm}\left(\widehat{\mathbf{n}}\right)\left(k_{\left(E\right)jm}^{\left(d\right)}\pm i{k}_{\left(B\right)jm}^{\left(d\right)}\right),\hfill \end{array}$$

(34)

$$\begin{array}{ccc}\hfill {\zeta }^3& =& \sum _{djm}{\omega }^{d-4}{Y}_{jm}\left(\widehat{\mathbf{n}}\right)k_{\left(V\right)jm}^{\left(d\right)}.\hfill \end{array}$$

(35)

In these expressions, $Y_{jm}\left(\widehat{\mathbf{n}}\right)$ are the usual spherical harmonics with $\widehat{n}=-\widehat{p}$, while ${}_{\pm 4}{Y}_{jm}\left(\widehat{\mathbf{n}}\right)$ are spin-weighted spherical harmonics. The coefficients, formerly in Cartesian tensor form in (30), are expressed in spherical form $k_{\left(I\right)jm}^{\left(d\right)}$, $k_{\left(E\right)jm}^{\left(d\right)}$, $k_{\left(B\right)jm}^{\left(d\right)}$, and $k_{\left(V\right)jm}^{\left(d\right)}$, where $j=0,1,\dots ,d-2$ and $-j\le m\le j$. The meaning of the subscripts $I,E,B,V$ and the relation between the two forms of the coefficients are determined by whether the terms are CPT odd or even and by 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=\omega /|\overrightarrow{p}|=1$). In the case of a Lorentz violation in the form in (30), the speed of the waves is given by

$$v=1-{\zeta }^0\pm \left|\overrightarrow{\zeta }\right|.$$

(36)

Given enough propagation distance from source to detector, a difference in arrival times may be detectable even for a 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 the 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 $\delta {\psi }_{\pm }$ due to terms in (36) or (30). Consider first the strain

$$h\sim e^{-i(\omega t-kl)}$$

(37)

where l is the distance travelled, and k is the wave number. The difference in phase grows in magnitude as the gravitational wave travels from the source to the detectors. On cosmological scales, it is important to include effects on propagation time from the expanding universe using luminosity distance. Noting $k\sim |\overrightarrow{p}|=\omega /v$, inputting (36), and including distance and frequency alterations form cosmology, one finds the phase shift expression

$${\delta }_{\pm }={\omega }_{obs}{\int }_0^{z}d{z}^{\prime }\frac{(-{\zeta }^0\pm |\overrightarrow{\zeta }\left|\right)}{H\left(z^{\prime }\right)},$$

(38)

where $H\left(z\right)$ is the Hubble parameter with redshift z, and the observed frequency is related to that emitted via ${\omega }_{obs}(1+z)={\omega }_{emit}$.

For each mode, the modified phase shift can be written as

$$\delta {\psi }_{\pm }=-\delta \pm \beta ,$$

(39)

where

$$\begin{array}{cc}\hfill \delta & ={\omega }^{d-3}\tau {\zeta }^{\left(d\right)0},\hfill \\ \hfill \beta & ={\omega }^{d-3}\tau \left|{\overrightarrow{\zeta }}^{\left(d\right)}\right|,\hfill \\ \hfill \tau & ={\int }_0^{z}dz\frac{{(1+z)}^{d-4}}{H\left(z\right)}\hfill \end{array}$$

(40)

and $|\overrightarrow{\zeta }|={\omega }^{d-4}\left|{\overrightarrow{\zeta }}^{\left(d\right)}\right|$ and ${\zeta }^0={\omega }^{d-4}{\zeta }^{\left(d\right)0}$. The $\tau$ is the effective propagation time due to cosmological redshift z.

It is useful to rewrite the coefficients in terms of effective angles $\vartheta$ and $\phi$ defined by

$$sin\vartheta =\frac{|{\zeta }^1\mp i{\zeta }^2|}{|\overrightarrow{\zeta }|},\phantom{30}cos\vartheta =\frac{{\zeta }^3}{|\overrightarrow{\zeta }|},\phantom{30}{e}^{\mp i\phi }=\frac{{\zeta }^1\mp i{\zeta }^2}{\sqrt{{\left({\zeta }^1\right)}^2+{\left({\zeta }^2\right)}^2}}.$$

(41)

Note that these angles are not the sky location angles $\theta$ and $\varphi$. Using the plus and cross polarizations (3) in addition, the modified gravitational wave solutions in terms of the Lorentz-invariant solutions can be written

$$\begin{array}{ccc}\hfill h_{(+)}& =& e^{i\delta }(cos\beta -isin\vartheta cos\phi sin\beta )h_{(+)}^{LI}\hfill \\ & & -e^{i\delta }sin\beta (cos\vartheta +isin\vartheta sin\phi )h_{(\times )}^{LI}\hfill \\ \hfill h_{(\times )}& =& e^{i\delta }(cos\beta +isin\vartheta cos\phi sin\beta )h_{(\times )}^{LI}\hfill \\ & & +e^{i\delta }sin\beta (cos\vartheta -isin\vartheta sin\phi )h_{(+)}^{LI}.\hfill \end{array}$$

(42)

The $h_{(+,\times )}^{LI}$ are the Lorentz-invariant gravitational wave for standard GR; one can retrieve GR as a limiting case as $\beta \to 0$ and $\delta \to 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\dots },\dots$ or in spherical tensor form $k_{\left(I\right)10}^{\left(d\right)},\dots$ [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.

2.3. Unit Changes and Dimension

For applications below, it becomes essential to convert from natural units to SI units when implementing modifications into the 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 on $\hslash =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 $\mathrm{Jm}$.2 For instance, the full Einstein Hilbert action in SI units can be written as

$$I_{EH}=\frac{c^4}{16\pi G}\int d^{4}x\sqrt{-g}R,$$

(43)

or, for the quadratic action limit of Equation (18), simply multiply by $c^4$. Units of $\mathrm{kg}\mathrm{m}{\mathrm{s}}^{-2}$ come from the factor $\frac{c^4}{G}$, ${\mathrm{m}}^4$ comes from $d^{4}x$, and ${\mathrm{m}}^{-2}$ comes from the derivatives contained within the Einstein tensor. Implicit here is the assumption that the metric tensor $g_{\mu \nu }$ is dimensionless (the Minkowski metric retains its form ${\eta }_{\mu \nu }=\mathrm{diag}(-1,1,1,1)$). Likewise, the Lorentz-violating action (19) contains operators with SI units ${\mathrm{m}}^{-2}$; thus, from (20), when introducing higher derivatives, the units of the coefficients compensate, so the coefficients have units ${\mathrm{m}}^{\mathrm{d}-4}$.

When converting the field Equation (28) from position to momentum space, every partial derivative contributes a factor with Planck’s constant, i.e., ${\partial }_{\alpha }\to \frac{i}{\hslash }{p}_{\alpha }$. Schematically, the position space equation has the form

$$\partial \partial h+s^{\left(4\right)}\partial \partial h+q^{\left(5\right)}\partial \partial \partial h+\dots =0,$$

(44)

where, e.g., for $d=4$ a term involving the $\widehat{s}$ operators contains coefficients for $s^{\left(4\right)}$ coupled to two derivatives that act on h. In momentum space,

$${\left(\frac{i}{\hslash }\right)}^{2}pph+{\left(\frac{i}{\hslash }\right)}^2{s}^{\left(4\right)}pph+{\left(\frac{i}{\hslash }\right)}^3{q}^{\left(5\right)}ppph+\dots =0,$$

(45)

where the operators $\widehat{s}$, $\widehat{q}$, and $\widehat{k}$ now contain ${\left(\frac{i}{\hslash }\right)}^{(d-2)}{p}_{{\alpha }_1}\dots p_{{\alpha }_{d-2}}$ in the place of partials. The units for the coefficients are unchanged, i.e., ${\mathrm{m}}^{\mathrm{d}-4}$.

One must also keep track of the corrected time-component factors in the four momenta, $p_{\alpha }=(-\frac{\hslash }{c}\omega ,\overrightarrow{p})$. For instance, the wave speed, via the dispersion relation, becomes

$$v_{\pm }=\hslash \omega /\left|\overrightarrow{p}\right|=c\left(1+c^2(-{\zeta }^0\pm |\overrightarrow{\zeta }\left|\right)\right)$$

(46)

where each $\zeta$ quantity in (32) inherits a factor of ${\left(\frac{\hslash }{c}\right)}^2$. To ensure the coefficients, $k_{\left(I\right)jm}^{\left(d\right)}$, $k_{\left(E\right)jm}^{\left(d\right)}$, $k_{\left(B\right)jm}^{\left(d\right)}$, and $k_{\left(V\right)jm}^{\left(d\right)}$ have SI units of ${\mathrm{m}}^{\mathrm{d}-4}$, we redefine Equations (33)–(35) by implementing a factor of $c^{(2-d)}$, e.g.,

$${\zeta }^0=c^{(2-d)}\sum _{djm}{\omega }^{d-4}{Y}_{jm}\left(\widehat{\mathbf{n}}\right)k_{\left(I\right)jm}^{\left(d\right)},$$

(47)

with similar SI factors for ${\zeta }^1,{\zeta }^2$, and ${\zeta }^3$.

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 the mass dimension 5 and 6 coefficients have been obtained from the non-observation of a delay between the arrival time of the $h_{+}$ and $h_{\times }$ 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 complement prior work by directly analyzing 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 performed 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 implemented the Bayesian analysis into a version of the LIGO Algorithm Library suite LALSuite modified for our purposes as described below [52].

3.1. 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 modifications of the GW generation and propagation [4]. Following a similar methodology, we implemented 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 inferred from the appearance of $\widehat{n}$ in Equation (42) via $\beta$ and $\delta$. 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 $\beta$ in (42) and obey the complex conjugate relation $k_{jm}^{\left(d\right)*}={(-1)}^m{k}_{j(-m)}^{\left(d\right)}$, for $j=0,1,2,3$, $-j\le m\le j$. There are a priori independent coefficients in this set of terms [26]. We display the first terms within $\beta$ in SI units:

$${\beta }^{\left(5\right)}=\frac{{\omega }^2{\tau }^{\left(5\right)}}{2\sqrt{\pi }c}\left|k_{\left(V\right)00}^{\left(5\right)}-\sqrt{\frac{3}{2}}sin\theta \left(e^{i\varphi }{k}_{\left(V\right)11}^{\left(5\right)}+e^{-i\varphi }{k}_{\left(V\right)11}^{\left(5\right)*}\right)+\sqrt{3}cos\theta k_{\left(V\right)10}^{\left(5\right)}+\dots \right|.$$

(48)

where the sky location of the source ($\theta$, $\varphi$) appears. The coefficients are taken as expressed in the SCF in this expression.

The general form of the signal observed in the interferometer is:

$$S_A=F_{(+)}{h}_{(+)}+F_{(\times )}{h}_{(\times )},$$

(49)

where $h_{(+,\times )}$ are the expressions (42), and $F_{(+,\times )}$ are the (standard) detector response functions. The expressions for $F_{(+,\times )}$ 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 $\tau$ 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 ${\left(k_{\left(V\right)jm}^{\left(5\right)}\right)}_{eff}=\tau k_{\left(V\right)jm}^{\left(5\right)}$. 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}$).

3.2. Bayesian Analysis

After implementing the modification of the strain, we included 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, ${\overrightarrow{\theta }}_{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, and inclination). We add to the preexisting parameters the SME coefficients ${\left(k_{\left(V\right)jm}^{\left(5\right)}\right)}_{eff}$ described in Section 3.1 for the mass dimension 5 case, contained within ${\overrightarrow{\theta }}_{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:

$$\begin{array}{c}\hfill P({\overrightarrow{\theta }}_{GR},{\overrightarrow{\theta }}_{SME}|d,I)=\frac{P\left(d\right|{\overrightarrow{\theta }}_{GR},{\overrightarrow{\theta }}_{SME},I)P({\overrightarrow{\theta }}_{GR},{\overrightarrow{\theta }}_{SME}|I)}{P\left(d\right|I)},\end{array}$$

(50)

where $P({\overrightarrow{\theta }}_{GR},{\overrightarrow{\theta }}_{SME}|d,I)$ is the posterior probability, $P\left(d\right|{\overrightarrow{\theta }}_{GR},{\overrightarrow{\theta }}_{SME},I)$ the likelihood, $P({\overrightarrow{\theta }}_{GR},{\overrightarrow{\theta }}_{SME}|I)$ the prior probability, and $P\left(d\right|I)$ the evidence, and any pertinent background information is included in I. We set a flat prior probability for ${\left(k_{\left(V\right)jm}^{\left(5\right)}\right)}_{eff}$ bounded between $|{\left(k_{\left(V\right)jm}^{\left(5\right)}\right)}_{eff}|\in [0;{10}^{-10}]$, with a maximal value well above the existing constraints on the order of ${10}^{-15}$ [49]. The likelihood is computed in the frequency domain:

$$\begin{array}{c}\hfill P\left(d\right|{\overrightarrow{\theta }}_{GR},{\overrightarrow{\theta }}_{SME},I)=exp\sum _i\left[-\frac{2|{\tilde{h}}_i({\overrightarrow{\theta }}_{GR},{\overrightarrow{\theta }}_{SME})-{\tilde{d}}_i{|}^2}{T{S}_n\left(f_i\right)}-\frac{1}{2}log\left(\frac{\pi T{S}_n\left(f_i\right)}{2}\right)\right],\end{array}$$

(51)

where ${\tilde{h}}_i$ is the frequency-domain template signal, ${\tilde{d}}_i$ represents the data observed by the interferometers, T is the duration of the signal, and $S_n$ is 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 inferred 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 Equation (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 marginalized posterior probability on a subset of parameters by integrating it 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.

4. Sensitivity Study

As an illustration, we assume, for the following, one non-zero coefficient $k_{\left(V\right)00}^{\left(5\right)}$ corresponding to isotropic polarization-dependent dispersion. Figure 1 plots the waveforms for both GR and the modified wave form for different values of $k_{\left(V\right)00}^{\left(d\right)}$. We assume a non-spinning binary system that has a luminosity distance of 4 Gpc, and equal masses of $m_1=m_2=50M_{\odot }$. Note that significant differences in the waveform shape occur for coefficient values as small as ${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 Figure 1 and Figure 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 performed 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 simulated a GW emitted by a non-spinning binary system of black holes with symmetric masses of $50M_{\odot }$ located at 5 Gpc where the dispersion is controlled by one coefficient set to a value of $k_{\left(V\right)00}^{\left(d\right)}={10}^{-14}$. Figure 2 shows the posterior probability on the luminosity distance and the coefficient, where both are recovered around the simulated values. The $1\sigma$ credible interval shows a constraint on $k_{\left(V\right)00}^{\left(d\right)}$ where the zero value is excluded, showing that the coefficient can be measured with a single event provided that it is relatively large. The $k_{\left(V\right)00}^{\left(d\right)}$ posterior probability density marginalized over the source and systematic uncertainties is shown in the violin plot.

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 first two observing runs as well as 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, which also induce a modified dispersion of the GW signal, we note that the constraint has been improved by one order of magnitude from a single event in the analysis of a larger population of GW detections. The constraint from GW159014 was $m_g\le 1.2\times {10}^{-22}{\mathrm{eV}/\mathrm{c}}^2$, while it is now $m_g\le 1.76\times {10}^{-23}{\mathrm{eV}/\mathrm{c}}^2$ when analyzing 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 a 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.

5. Conclusions and Future Work

We described 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 implemented 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_{\mu \nu }$, (19), and after theoretical constraints including gauge invariance, we arrive at the general result in (27). From the field Equation (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 the sky location of the source. Thus, 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 complete our analysis for coefficients for Lorentz and CPT violation of mass dimension 5 and 6, with a global analysis. In 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.