# Gravitational-wave lensing in Einstein-aether theory

Julius Streibert  Fachbereich Physik, Freie Universität Berlin, D-14195
Berlin, Germany  Max Planck Institute for Gravitational Physics (Albert
Einstein Institute), D-14476 Potsdam, Germany  Hector O. Silva  Max Planck
Institute for Gravitational Physics (Albert Einstein Institute), D-14476
Potsdam, Germany  Miguel Zumalacárregui  Max Planck Institute for
Gravitational Physics (Albert Einstein Institute), D-14476 Potsdam, Germany

###### Abstract

Einstein-aether theory provides a model to test the validity of local Lorentz
invariance in gravitational interactions. The speed of gravitational waves as
measured from the binary neutron star event GW170817 sets stringent limits on
Einstein-aether theory, but only on a combination of the theory’s free
parameters. For this reason, a significant part of the theory’s parameter
space remains unconstrained by observations. Motivated by this, we explore the
propagation of gravitational waves in Einstein-aether theory over an
inhomogeneous background (i. e., gravitational wave lensing) as a potential
mechanism to break the degeneracies between the theory’s free parameters, and
hence enable new constraints on the theory to be obtained. By bringing the
field equations into the form of the so-called kinetic matrix and applying a
formalism known as the propagation eigenstate framework, we find that the
speed of gravitational waves is modified by inhomogeneities in the aether
field. However, the modification is common to both gravitational polarizations
and vanishes in the limit in which gravitational waves propagate with luminal
speed. This lens-dependent gravitational wave speed contrasts with the lens-
induced birefringence observed in other theories beyond general relativity,
like Horndeski’s theory. While the potential to improve tests based on
gravitational-wave speed is limited, our formalism sets the basis to fully
describe signal propagation over inhomogeneous spacetimes in Einstein-aether
theory and other extensions of general relativity.

##  I Introduction

Theories of gravity other than Einstein’s general relativity (GR), termed
“beyond GR”, are motivated, among other reasons, by unexplained phenomena like
dark matter and dark energy as well as by attempts to construct a theory of
quantum gravity [1]. Some theories beyond GR can be categorized by which types
of fundamental fields they admit in addition to the metric tensor. Of these,
theories with additional scalar fields, such as Horndeski’s theory, are
amongst the most well studied [2, 1]. Other theories of modified gravity have
received comparatively less attention [3]. This includes vector-tensor
theories.

The introduction of a vector field to a theory can break local Lorentz
invariance (LLI) by introducing a preferred direction at each point in
spacetime [4]. Such violations of LLI are of theoretical interest. In
particular, the motivation for investigating theories violating LLI stems from
theories of quantum gravity, some of which utilize or even require a violation
of LLI [5]. LLI violations are expected if spacetime is discrete at the Planck
scale, as predicted by certain realizations of quantum gravity [6]. Note that
LLI can only be tested up to the maximum energy and velocity scale that is
experimentally accessible: the possibility of violations at even higher
energies remains [4].

Einstein-aether theory is a vector-tensor theory that violates LLI by means of
an additional unit timelike vector field, the aether. The theory was
originally conceived by Jacobson and Mattingly as “gravity with a dynamical
preferred frame” [4] and is well-posed in all viable regions of its parameter
space [7]. The theory admits spherically symmetric solutions [8] including
analogues to stars [9] and black holes [10]. By making some additional
assumptions Einstein-aether theory becomes the IR limit of an extended form of
Hořava–Lifshitz theory [11], linking it to the quantum gravity motivations for
violations of LLI.

A wide variety of both theoretical and observational constraints can be
imposed on Einstein-aether theory. While some of its parameters (for example
the Eddington-Robertson-Schiff parameters) match those of GR [8], others
deviate. For example, Einstein-aether theory admits two forms of gravitational
constant, local and cosmological, both of which can in principle deviate from
the GR case [12]. Other parameters that can differ from GR are the post-
Newtonian preferred frame parameters both in the weak field [13, 14] and in
the strong field [15].

Each of the five fluctuations propagates with a different speed, the value of
which depends on the theory’s coupling constants [16]. The presence of
additional polarizations can be probed by gravitational interferometers [17,
18, 19, 20], although the coincident arrival of signals with different
helicity is not expected, due to differences in speed accumulating over the
signal’s travel time.

The (necessarily positive) mode energy densities of Einstein-aether theory
were first calculated by Eling [21]. The same year it was shown that Einstein-
aether mode speeds can not be significantly subluminal by argument of
gravitational Cherenkov radiation [22]. Radiation damping due to Einstein-
aether radiation was first discussed in Ref. [23], this opened the door to
constraints derived from observations of binary pulsars [24, 15].

Stringent limits on Einstein-aether can be obtained by comparing the arrival
times of gravitational and electromagnetic radiation emitted by the same
source [25, 26, 27, 28, 29, 30]. This was possible after the observation of a
neutron-star merger, GW170817 [31], the prompt emission of a short gamma-ray
burst [32] and other electromagnetic counterparts [33]. GW170817 was shown to
be only consistent with standard GW polarizations [34, 35]: therefore, the
short time delay between gravitational and electromagnetic signals at $\sim
40$ Mpc translated into an exquisite limit on the speed of “standard” GWs. The
theory parameters that affect this quantity were thus exquisitely constrained
[36, 6].

GW170817 shifted the landscape of constraints on Einstein-aether theory. Since
tensor wave speeds were constrained to a precision of $10^{-15}$, previously
used series expansion became invalid [6]. In particular, this made old bounds
obtained from binary pulsar observations obsolete and allowed new ones to be
derived [37]. There have also been attempts to constrain Einstein-aether
theory with waveform modeling and parameter estimation based on both GW170817
and GW190425 [38].

All tests of GW propagation in Einstein-aether theory have assumed a perfectly
homogeneous background. While a reasonable first step to model GW propagation
on the average universe [39], this approach neglects inhomogeneities in the
metric and Einstein-aether field. Inhomogeneities in the metric are
responsible for gravitational lensing, the deflection, magnification, and
delay of signals propagating through the universe. In some cases,
gravitational lensing can split a source into multiple images and profoundly
distort their shape. Among many applications, gravitational lensing of
electromagnetic sources has provided tests of gravity on extragalactic and
cosmological scales [40, 41]. Observations of lensed GWs can be used in a
similar way [42] and will also improve our ability to test the presence of
additional polarizations [43, 44], modified speeds [45, 46] and modified
propagation [47, 48, 49]. In addition, it can be used to constrain the mass of
gravitons [50] and the charged hair of black holes [51].

Besides lensing, novel propagation effects can be caused by inhomogeneities in
the additional fields, e. g. the Einstein-aether. A non-homogeneous
configuration allows interactions between excitations with different helicity
(scalar, vector, and tensor fluctuations). In addition to modifying the
propagation speed, nonsymmetric configurations can introduce birefringence, a
speed difference between polarizations with the same helicity. Lens-induced GW
birefringence happens in scalar-tensor theories, leading to different
propagation speed for the $+$ and $\times$ polarization [52]. Beyond-Einstein
birefringence can be tested without the need of an electromagnetic
counterpart, producing limits comparable to those of the neutron-star merger
[53]. Other classes of theories are predict birefringence between the left-
and right-polarized GWs [54, 55, 56, 57].

The purpose of this paper is to examine the role that inhomogeneities in the
aether field can have on the propagation of GWs, a phenomenon we will
generically refer to as “GW lensing” (although inhomogeneities do not need to
produce multiple images or even be located between the source and observer).
Our goal is to establish how the propagation of standard and additional
polarizations depends on the theory parameters and whether they can lead to
novel tests, analogous to GW birefringence.

To study GW lensing in Einstein-aether theory we will first review the basics
of the theory itself: we discuss Einstein-aether theory and its GWs without
lensing and with a flat Minkowski background in Sec. II. We introduce lensing
by switching to a non-Minkowski background in Sec. III. We utilize certain
approximations and initially obtain the new propagating modes to the lowest
order in these approximations. In Sec. IV, we then turn to consecutively
higher orders. Finally we take stock of our results and of what future steps
can be taken in Sec. V.

We use the mostly-plus metric signature and use geometrical units with
$c=G=1$, unless stated otherwise.

##  II Einstein-Aether theory

In this section, we review the basics of Einstein-aether theory. In Sec. II.1,
we present the theory’s action and field equations. We then linearize the
theory around a Minkowski background and identify the theory’s degrees of
freedom in Sec. II.2. Finally, we review the current observational constraints
on the theory in Sec. II.3.

###  II.1 Action and field equations

Einstein-aether theory is the most general vector-tensor theory of gravity
that satisfies the following conditions [58, 4]:

  1. 1.

the vector field $u^{\alpha}$, called the aether, is constrained to be unit
timelike, $u^{\alpha}u_{\alpha}=-1$.

  2. 2.

the equations of motion are of at most second order in derivatives.

  3. 3.

the equations of motion are of at most second order in the aether.

The aether is introduced in order to break local Lorentz invariance (LLI) by
introducing a preferred direction at each point in spacetime [4]. Since a
static aether would violate general covariance, the aether must instead by
dynamic [4], itself subject to field equations. A dynamic aether with an
unconstrained modulus typically evolves towards zero, returning the GR limit
and restoring LLI. For this reason, we demand the aether to be unit timelike.

We now construct the theory’s action. In principle, it can contain Lagrangians
for gravity, aether, and matter fields, as well as gravity-aether and matter-
aether couplings. Here, we will assume that aether and matter do not couple
directly. This leaves us with the Einstein-Hilbert term for gravity, $R$, the
standard matter Lagrangian $\mathcal{L}_{\rm m}$ as well as a new Einstein-
aether Lagrangian density $\mathcal{L}_{\text{ae}}$ containing the gravity-
aether coupling:

$S=\frac{1}{2\kappa_{\text{ae}}}\int{\rm
d}^{4}x\,\sqrt{-g}\,\left(R+\mathcal{L}_{\text{ae}}+2\kappa_{\text{ae}}\mathcal{L}_{m}\right).$
(1)

Here, $\kappa_{\text{ae}}$ is the Einstein-aether equivalent of the Einstein
gravitational constant, which may in principle differ from the GR value of
$8\pi$.111To be precise, $\kappa_{\text{ae}}$ can be shown to be related to
two of the four coupling constants of Einstein-aether theory, $c_{1}$ and
$c_{4}$: $\kappa_{\text{ae}}=4\pi(2-c_{1}-c_{4})$ [59]. The Einstein-aether
Lagrangian density contains two terms:

$\displaystyle\mathcal{L}_{\text{ae}}$
$\displaystyle=-M^{\alpha\beta}_{\phantom{\alpha\beta}\mu\nu}u^{\mu}_{;\alpha}u^{\nu}_{;\beta}+\lambda\left(g_{\alpha\beta}u^{\alpha}u^{\beta}+1\right).$
(2)

The first contains the possible gravity-aether couplings while the second
enforces the unit timelike constraint for the aether via a Lagrange
multiplier.

The tensor $M^{\alpha\beta}_{\phantom{\alpha\beta}\mu\nu}$ contains the four
dimensionless coupling constants, $c_{1}$ through $c_{4}$, of Einstein-aether
theory:

$\displaystyle M^{\alpha\beta}_{\phantom{\alpha\beta}\mu\nu}$
$\displaystyle=c_{1}g^{\alpha\beta}g_{\mu\nu}+c_{2}\delta^{\alpha}_{\mu}\delta^{\beta}_{\nu}+c_{3}\delta^{\alpha}_{\nu}\delta^{\beta}_{\mu}-c_{4}u^{\alpha}u^{\beta}g_{\mu\nu}.$
(3)

We also use the following notation for linear combinations of the coupling
constants:

$\displaystyle c_{ijk\dots}$ $\displaystyle=c_{i}+c_{j}+c_{k}+\dots,$ (4)
$\displaystyle c_{+}$ $\displaystyle=c_{1}+c_{3},$ (5) $\displaystyle c_{-}$
$\displaystyle=c_{1}-c_{3}.$ (6)

General relativity is recovered when all coupling constants, $c_{1}$ through
$c_{4}$, vanish.

The theory’s field equations follow from varying Einstein-aether action with
respect to the Lagrange multiplier, the metric tensor and the aether. The
Lagrange multiplier variation enforces the constraint:

$\displaystyle g_{\alpha\beta}u^{\alpha}u^{\beta}$ $\displaystyle=-1.$ (7)

The metric variation gives us the Einstein-aether equivalent of the Einstein
field equations:

$\displaystyle G^{\alpha\beta}-T^{\alpha\beta}_{\text{ae}}$
$\displaystyle=\kappa_{\text{ae}}T_{m}^{\alpha\beta},$ (8)

where $T^{\alpha\beta}_{\text{ae}}$ is called the aether stress-energy tensor.
It can be defined in terms of the aether current
$J^{\alpha}_{\phantom{\alpha}\mu}$ and the aether acceleration $a^{\mu}$:

$\displaystyle T^{\alpha\beta}_{\text{ae}}$
$\displaystyle=-D_{\mu}\left[u^{(\beta}J^{\alpha)\mu}-J^{\mu(\alpha}u^{\beta)}-J^{(\alpha\beta)}u^{\mu}\right]$
$\displaystyle\mathrel{\phantom{=}}-c_{1}\left[\left(D_{\mu}u^{\alpha}\right)\left(D^{\mu}u^{\beta}\right)-\left(D^{\alpha}u_{\mu}\right)\left(D^{\beta}u^{\mu}\right)\right]$
$\displaystyle\mathrel{\phantom{=}}+c_{4}a^{\alpha}a^{\beta}+\lambda
u^{\alpha}u^{\beta}-\frac{1}{2}g^{\alpha\beta}J^{\delta}_{\phantom{\delta}\sigma}D_{\delta}u^{\sigma},$
(9) $\displaystyle J^{\alpha}_{\phantom{\alpha}\mu}$
$\displaystyle=M^{\alpha\beta}_{\phantom{\alpha\beta}\mu\nu}D_{\beta}u^{\nu},$
(10) $\displaystyle a^{\mu}$ $\displaystyle=u^{\alpha}D_{\alpha}u^{\mu}.$ (11)

The variation with respect to the aether yields another set of field
equations:

$\displaystyle\text{\AE}_{\mu}$
$\displaystyle=D_{\alpha}J^{\alpha}_{\phantom{\alpha}\mu}+c_{4}a_{\alpha}D_{\mu}u^{\alpha}+\lambda
u_{\mu}=0$ (12)

At this stage, $\lambda$ is usually immediately eliminated [8, 23, 6]. For
this, it will prove convenient to split the four components of the aether
field equation across two new equations by multiplying with the aether and
taking symmetric and antisymmetric part:

$\displaystyle 0$ $\displaystyle=\text{\AE}^{(\alpha}u^{\beta)}$
$\displaystyle=u^{(\beta}D_{\mu}J^{\mu|\alpha)}+c_{4}a_{\mu}u^{(\beta}D^{\alpha)}u^{\mu}+\lambda
u^{\alpha}u^{\beta},$ (13) $\displaystyle 0$
$\displaystyle=\text{\AE}^{[\alpha}u^{\beta]}$
$\displaystyle=u^{[\beta}D_{\mu}J^{\mu|\alpha]}+c_{4}a_{\mu}u^{[\beta}D^{\alpha]}u^{\mu}.$
(14)

Solving Eq. (13) for $D_{\mu}\left(J^{\mu(\alpha}u^{\beta)}\right)$ and
substituting the result in Eq. (9), we find:

$\displaystyle T^{\alpha\beta}_{\text{ae}}$
$\displaystyle=-D_{\mu}\left[u^{(\beta}J^{\alpha)\mu}-J^{(\alpha\beta)}u^{\mu}\right]$
$\displaystyle\mathrel{\phantom{=}}-\,c_{1}\left[\left(D_{\mu}u^{\alpha}\right)\left(D^{\mu}u^{\beta}\right)-\left(D^{\alpha}u_{\mu}\right)\left(D^{\beta}u^{\mu}\right)\right]$
$\displaystyle\mathrel{\phantom{=}}+\,c_{4}a^{\alpha}a^{\beta}+J^{\mu(\alpha}D_{\mu}u^{\beta)}-c_{4}a_{\mu}u^{(\beta}D^{\alpha)}u^{\mu}$
$\displaystyle\mathrel{\phantom{=}}-\frac{1}{2}g^{\alpha\beta}J^{\delta}_{\phantom{\delta}\sigma}D_{\delta}u^{\sigma}.$
(15)

$\text{\AE}_{\mu}$ has four components and one variable was eliminated by
using the symmetric part of its field equation. It then follows that the
antisymmetric part has three remaining components. We choose the three entries
$\text{\AE}^{[0}u^{1]}$, $\text{\AE}^{[0}u^{2]}$ and $\text{\AE}^{[0}u^{3]}$
as these independent components. All remaining entries of
$\text{\AE}^{[\alpha}u^{\beta]}$ can be expressed as a linear combination of
the $\text{\AE}^{[0}u^{i]}$. For example, $\text{\AE}^{[1}u^{2]}$ is equal to
$u^{1}\text{\AE}^{[0}u^{2]}/u^{0}-u^{2}\text{\AE}^{[0}u^{1]}/u^{0}$. Our
remaining set of field equations then is:

$\displaystyle g_{\alpha\beta}u^{\alpha}u^{\beta}$ $\displaystyle=-1,$ (16)
$\displaystyle G^{\alpha\beta}-T^{\alpha\beta}_{\text{ae}}$
$\displaystyle=\kappa_{\text{ae}}T_{m}^{\alpha\beta},$ (17)
$\displaystyle\text{\AE}^{[0}u^{i]}$ $\displaystyle=0.$ (18)

###  II.2 Linearized theory and propagating degrees of freedom

We now linearize our field equations by expanding the metric tensor around the
flat Minkowski background, up to first order in a perturbation $h_{\mu\nu}$.

$\displaystyle\overline{g}_{\mu\nu}$ $\displaystyle=\eta_{\mu\nu}+h_{\mu\nu},$
(19)

and we denote linearization with an overline. By solving the background field
equations it can be shown that the aether corresponding to a flat Minkowski
background is $\delta^{\mu}_{0}$ [16]. Calling the aether perturbation
$w^{\mu}$, the aether is then linearized as follows:

$\displaystyle\overline{u}^{\mu}$ $\displaystyle=\delta^{\mu}_{0}+w^{\mu}.$
(20)

Note that using covariant linearization is also possible [6]. This means we
could also have set $\overline{u}_{\mu}^{(0)}=\delta_{\mu}^{0}$, where the
superscript $(0)$ denotes zeroth-order in the perturbation.

We now go to combined first order in $h_{\mu\nu}$ and $w^{\mu}$ in the field
equations. The linearization of the constraint equation (16) immediately
eliminates one component:

$\displaystyle
0=\overline{g_{\alpha\beta}u^{\alpha}u^{\beta}}+1=h_{00}-2w^{0}\Leftrightarrow
w^{0}=\frac{h_{00}}{2}.$ (21)

To linearize the metric field equation (17) and the asymmetric aether field
equation (18) we need the linearization of $D_{\beta}u^{\nu}$. We note that
Christoffel symbols of the Minkowski background in Cartesian coordinates have
no zeroth order in the perturbation:

$\displaystyle\Gamma^{\nu(0)}_{\beta\mu}$ $\displaystyle=0.$ (22)

Since the partial derivative of a constant is zero
$\left(D_{\beta}u^{\nu}\right)^{(0)}$ vanishes. Quantities derived from
covariant derivatives of the aether, like $a^{\mu}$ and
$J^{\alpha}_{\phantom{\alpha}\mu}$, also vanish at zeroth order. Therefore, we
see that covariant derivatives of the aether, aether acceleration and aether
current must be of first order at minimum, which means that products of any of
these are of higher order , and any term multiplied with them must be of
zeroth order.

In addition to linearizing gravity, we are right now interested in the
propagation of GWs, not their emission. Therefore, we consider free GWs, set
$T^{\alpha\beta}_{\rm m}=0$ and work with the unsourced, linearized field
equations:

$\displaystyle\overline{G}^{\alpha\beta}-\overline{T}^{\alpha\beta}_{\text{ae}}$
$\displaystyle=0,$ (23) $\displaystyle\overline{\text{\AE}^{[0}u^{i]}}$
$\displaystyle=0.$ (24)

To linear order in perturbations, the quantities that appear in the field
equations are:

$\displaystyle\overline{T}^{\alpha\beta}_{\text{ae}}$
$\displaystyle=\partial_{0}\overline{J}^{(\alpha\beta)}-\delta^{(\beta}_{0}\partial_{\mu}\overline{J}^{\alpha)\mu},$
(25) $\displaystyle\overline{\text{\AE}^{[0}u^{i]}}$
$\displaystyle=-\frac{1}{2}\partial_{\mu}\overline{J}^{\mu i},$ (26)
$\displaystyle\overline{J}^{\alpha}_{\phantom{\alpha}\mu}$
$\displaystyle=M^{\alpha\beta(0)}_{\phantom{\alpha\beta}\mu\nu}\overline{D_{\beta}u^{\nu}},$
(27) $\displaystyle M^{\alpha\beta(0)}_{\phantom{\alpha\beta}\mu\nu}$
$\displaystyle=c_{1}\eta^{\alpha\beta}\eta_{\mu\nu}+c_{2}\delta^{\alpha}_{\mu}\delta^{\beta}_{\nu}+c_{3}\delta^{\alpha}_{\nu}\delta^{\beta}_{\mu}-c_{4}\delta^{\alpha}_{0}\delta^{\beta}_{0}\eta_{\mu\nu}.$
(28)

Instead of solving the ten components of the metric field equations and the
remaining three components of the aether field equations, we will combine the
two into an equation for what we call the effective tensor $X^{\alpha\beta}$:

$\displaystyle X^{\alpha\beta}$ $\displaystyle=0,$ (29) $\displaystyle X^{00}$
$\displaystyle=\overline{G}^{00}-\overline{T}^{00}_{\text{ae}},$ (30)
$\displaystyle X^{0i}$
$\displaystyle=\overline{G}^{0i}-\overline{T}^{0i}_{\text{ae}}+\overline{\text{\AE}^{[0}u^{i]}},$
(31) $\displaystyle X^{i0}$
$\displaystyle=\overline{G}^{i0}-\overline{T}^{i0}_{\text{ae}}-\overline{\text{\AE}^{[i}u^{0]}},$
(32) $\displaystyle X^{ij}$
$\displaystyle=\overline{G}^{ij}-\overline{T}^{ij}_{\text{ae}}.$ (33)

This approach mostly matches that of Ref. [23], the key difference is
discussed in Appendix A. Since the metric field equation contains purely
symmetric terms, we do not lose any information by combining it with the
equation for the purely antisymmetric $\overline{\text{\AE}^{[0}u^{i]}}$.

In order to solve the linearized field equations, we perform a scalar-vector-
tensor (SVT) decomposition of metric and aether perturbations. Specifically,
we use a decomposition introduced by Foster [23]:

$\displaystyle h_{0i}$ $\displaystyle=\gamma_{i}+\gamma_{,i},$ (34)
$\displaystyle h_{ij}$
$\displaystyle=\phi_{ij}+\frac{1}{2}P_{ij}f+2\phi_{(i,j)}+\phi_{,ij},$ (35)
$\displaystyle w_{i}$ $\displaystyle=v_{i}+v_{,i}.$ (36)

where $P_{ij}$ is the following differential operator,

$\displaystyle P_{ij}x$ $\displaystyle=\delta_{ij}\Delta x-x_{,ij}.$ (37)

and $\Delta=\partial^{k}\partial_{k}$ is the Laplacian operator.

In Eqs. (34) through (36), $\gamma$, $f$, $\phi$ and $v$ are scalars;
$\gamma_{i}$, $\phi_{i}$ and $v_{i}$ are transverse vectors, meaning that
their divergence vanishes, whereas $\phi_{ij}$ is a symmetric transverse
traceless matrix, meaning that
$\phi^{ij}_{\phantom{ij},j}=\phi^{ij}_{\phantom{ij},i}=\phi^{i}_{\phantom{i}i}=0$.
In essence, we decompose the vector parts into transverse and longitudinal
parts and the tensor parts into transverse and longitudinal as well as
traceful and traceless parts.

Note that we can freely raise and lower spatial indices since after
linearization, and the raising and lowering operations are done with
$\eta_{\alpha\beta}$, the spatial part of which is $\delta_{ij}$. This means
that from now on we also apply the Einstein sum convention to terms like
$x_{kk}$. We can easily calculate the trace $h=\eta^{\mu\nu}h_{\mu\nu}$ in
terms of our decomposition:

$\displaystyle h$ $\displaystyle=F+\Delta\phi-h_{00}.$ (38)

Here we adopt the notation $F=\Delta f$ used by Foster [23] since $F$ will at
some points turn out to be a more intuitive variable than $f$.

The field equations can be simplified by eliminating some components. Since
$h_{\mu\nu}$ is a symmetric 4×4 tensor it has ten independent components. In
terms of our decomposition, there are four scalar components $h_{00}$,
$\gamma$, $f$ and $\phi$. As transverse vectors, $\gamma_{i}$ and $\phi_{i}$
have two independent components each, as does the symmetric transverse
traceless tensor $\phi_{ij}$ (three components are constrained by
transverseness and one by tracelessness). We can remove some of these degrees
of freedom (DOF) by using gauge invariance. Einstein-aether theory is
invariant under general coordinate transformations $\varepsilon^{\mu}$:

$\displaystyle x^{\prime\mu}$ $\displaystyle=x^{\mu}+\varepsilon^{\mu}.$ (39)

Because $\varepsilon_{\mu}$ has four entries, so we can eliminate four
components with this gauge invariance. We will choose our gauge so that

$\gamma=\phi_{i}=v=0\,.$ (40)

We show that this choice is possible in Appendix B.

We can now calculate the linearized covariant derivatives
$\overline{D_{\beta}u^{\nu}}$, linearized aether current
$\overline{J}^{\alpha}_{\phantom{\alpha}\mu}$ and, finally, the effective
tensor $X^{\alpha\beta}$ in terms of the SVT decomposition of the
perturbations. We find:

$\displaystyle X_{00}$ $\displaystyle=-\frac{1}{2}\Delta
F+\frac{c_{14}}{2}\Delta h_{00},$ (41) $\displaystyle X_{0i}$
$\displaystyle=-\frac{1}{2}\Delta\gamma_{i}-\frac{1}{2}\dot{F}_{,i}-c_{14}\left(\ddot{v}_{i}+\ddot{\gamma}_{i}\right)$
$\displaystyle\mathrel{\phantom{=}}+\frac{c_{-}}{2}\Delta\left(v_{i}+\gamma_{i}\right)+\frac{c_{14}}{2}\dot{h}_{00,i},$
(42) $\displaystyle X_{i0}$
$\displaystyle=-\frac{1}{2}\Delta\gamma_{i}-\frac{1}{2}\dot{F}_{,i}$
$\displaystyle\mathrel{\phantom{=}}-\frac{1}{2}\Delta\left\\{c_{+}v_{i}+\partial_{i}\left[\left(c_{+}+c_{2}\right)\dot{\phi}+c_{2}\dot{f}\right]\right\\},$
(43) $\displaystyle X_{ij}$
$\displaystyle=-\frac{1}{2}\left(\Delta\phi_{ij}-\ddot{\phi}_{ij}\right)-\dot{\gamma}_{(i,j)}$
$\displaystyle\mathrel{\phantom{=}}+\frac{1}{4}P_{ij}\left(F-\ddot{f}-2h_{00}-2\ddot{\phi}\right)-\frac{1}{2}\ddot{f}_{,ij}$
$\displaystyle\mathrel{\phantom{=}}-\frac{c_{+}}{2}\ddot{\phi}_{ij}-c_{+}\dot{v}_{(i,j)}$
$\displaystyle\mathrel{\phantom{=}}-\frac{1}{2}P_{ij}\left[c_{2}\left(\ddot{\phi}+\ddot{f}\right)+\frac{c_{+}}{2}\ddot{f}\right]$
$\displaystyle\mathrel{\phantom{=}}-\frac{1}{2}\partial_{i}\partial_{j}\left[\left(c_{2}+c_{+}\right)\ddot{\phi}+c_{2}\ddot{f}\right],$
(44)

where an overdot denotes a derivative with respect to time, and $c_{\pm}$ was
defined in Eqs. (5) and (6).

In addition to the SVT decomposition of the perturbations we now also perform
a SVT decomposition of the effective tensor $X^{\alpha\beta}$ itself. By
setting the result to zero we arrive at the new form of the field equations.
We also drop prefactors and global derivatives since we want no homogeneous or
static terms in the solutions. We obtain four redundant components and four
constraints:

$\displaystyle h_{00}$ $\displaystyle=\frac{F}{c_{14}},$ $\displaystyle v_{i}$
$\displaystyle=-\frac{\gamma_{i}}{c_{+}},$ $\displaystyle\phi$
$\displaystyle=-\frac{1+c_{2}}{c_{123}}f.$ (45)

Substituting these constraints into the remaining components yields the
gravitational-wave equations of Einstein-aether theory. We obtain five wave
equations for the single scalar, two vector, and two tensorial DOF, namely,

$\displaystyle 0$ $\displaystyle=\Delta f-c_{S}^{-2}\,\ddot{f},$ (46)
$\displaystyle 0$
$\displaystyle=\Delta\gamma_{i}-{c_{V}^{-2}}\,\ddot{\gamma}_{i},$ (47)
$\displaystyle 0$
$\displaystyle=\Delta\phi_{ij}-{c_{T}^{-2}}\,\ddot{\phi}_{ij},$ (48)

that propagate with the speeds,

$\displaystyle\frac{1}{c_{S}^{2}}$
$\displaystyle=\frac{(1-c_{+})(2+3c_{2}+c_{+})c_{14}}{c_{123}(2-c_{14})},$
(49) $\displaystyle\frac{1}{c_{V}^{2}}$
$\displaystyle=\frac{2c_{14}(1-c_{+})}{2c_{1}-c_{+}c_{-}},$ (50)
$\displaystyle\frac{1}{c_{T}^{2}}$ $\displaystyle=1-c_{+},$ (51)

respectively.

During this derivation we have divided by certain combinations of the theory
constants: $c_{14}$, $c_{+}$ and $c_{123}$ for the constraint equations as
well as $2-c_{14}$ and $2c_{1}-c_{+}c_{-}$ for the wave equations. We have
therefore implicitly excluded cases where these quantities are equal to zero.
If we were to fine-tune these combinations so that they vanish then the
derivation of certain modes falls apart. Since we demand that our solutions
have no static or homogeneous terms, these modes are most often constrained to
zero instead of being described by a wave equation.222 There is only one
exception to this: if $c_{+}=0$ the constraint equation for $v_{i}$ falls
apart and $\gamma_{i}$ is constrained to zero, but $v_{i}$ still obeys a wave
equation with mode speed $c_{V}$. However, since we do not consider direct
couplings between aether and matter in this paper, this kind of wave has no
measurable effect.  Similar issues arise in the cases $1-c_{+}=0$ as well as
$2+3c_{2}+c_{+}=0$ and $2c_{1}-c_{+}c_{-}=0$. All three cases result in
vanishing inverse mode speeds, changing the character of the equation from
hyperbolic to elliptic. Therefore we now explicitly exclude these cases as
well.

A possible solution to each wave equation is given by a monochromatic plane
wave that propagates along the $z$ direction:

$\psi=\hat{\psi}\exp\left[i\omega\left(\frac{z}{c_{p}}-t\right)-\frac{i\pi}{2}\right]$
(52)

where the subscript $p$ indicates that the propagation speed corresponds to
either a scalar, vector or tensorial polarization. Here, $\hat{\psi}$ is the
wave amplitude, $\omega$ is the angular frequency and $-i\pi/2$ is a phase
chosen such that the real part of Eq. (52) is a sine of $\omega(z-t)$. This
means that $\psi=0$ for $z=t$ and, specifically, $z=t=0$. We choose this
convention so that the coordinates of test particles in the $xy$ plane are
unaffected at $z=t=0$.

Using Eq. (52), we can construct the full metric and aether perturbations,
which expressed in terms of the SVT decomposition reads:

$\displaystyle w^{\mu}=\left(\begin{array}[]{c}\dfrac{h_{00}}{2}\\\ \\\
v_{i}\end{array}\right),\;$ $\displaystyle
h_{\mu\nu}=\left(\begin{array}[]{ccc}h_{00}&&\gamma_{i}\\\ \\\
\gamma_{i}&&\phi_{ij}+\dfrac{1}{2}P_{ij}f+\phi_{,ij}\end{array}\right).$ (59)

If we specify the $z$ direction as propagation direction, then there are six
gauge invariant polarizations for GWs in theories beyond GR:

  1. 1.

the scalar breathing mode $h_{b}$,

  2. 2.

the scalar longitudinal mode $h_{l}$,

  3. 3.

the vector $x$ mode $h_{x}$,

  4. 4.

the vector $y$ mode $h_{y}$,

  5. 5.

the tensor plus mode $h_{+}$,

  6. 6.

the tensor cross mode $h_{\times}$.

The effect of all six polarizations on a ring of test particles is shown in
Fig. 1.

Figure 1: Qualitative effect of GWs beyond GR on circles and spheres composed
of test particles. We assume the GW wavelengths to be much larger than the
diameter of any such circle or sphere. Tensor and vector waves deform circles
into ellipses in different planes and along different axes. These modes
preserve the circle areas. The scalar modes instead expand and contract
circles and spheres, changing their area and volume, respectively : in the
$xy$ plane for the breathing mode and in the $z$ direction for the
longitudinal mode. Figure modified from Ref. [60], compare also Ref. [61].

In Einstein-aether theory the vector modes are proportional to $\gamma_{1}$
and $\gamma_{2}$ and the tensor modes are proportional to $\phi_{11}$ and
$\phi_{12}$. The breathing mode is proportional to $F$ and the longitudinal
mode is proportional to both $F$ and a linear combination of theory constants
that we call $\alpha$:

$\displaystyle\alpha=\frac{(1-c_{+})(2+3c_{2}+c_{+})}{c_{123}(2-c_{14})}-\frac{1+c_{2}}{c_{123}}.$
(61)

This means that $h_{b}$ and $h_{l}$ are combined into a joint scalar
polarization. Whether or not longitudinal effects are in phase with breathing
effects, are out of phase or vanish completely depends on the value of
$\alpha$. In Fig. 2 we show the effect of the joint breathing and longitudinal
scalar polarization on a sphere of test particles for $\alpha=0$ and $\pm 1$.

Figure 2: Effect of scalar waves with amplitude 0.3 and different values of
$\alpha$ on an a sphere of test particles with diameter $1$ and a scalar
wavelength much larger than $1$. Time is measured in fractions of the period
$T$, and different curves correspond to $\alpha=0$ and $\pm 1$. The sphere
expands and contracts both in the $xy$ plane and in the $z$ direction. The
effect differs based on the value of $\alpha$. For example, at $t=T/4$ the
sphere expands in the $z$ direction for $\alpha=1$, while it contracts for
$\alpha=-1$.

A more detailed derivation of the effects of Einstein-aether waves can be
found in Appendix C.

###  II.3 Theory constraints

Theory and observations impose a variety of constraints on Einstein-aether
theory. The theory allows for scalar, vector and tensor propagating modes with
speeds that differ from $1$. To avoid imaginary frequency instabilities, the
squared speeds must be nonnegative, $c_{S,V,T}^{2}\geq 0$ [16, 59]. To avoid
the existence of gravitational Cherenkov radiation, they can also not be
smaller than $1$ up to a precision of about $10^{-15}$,
$c_{S,V,T}^{2}>1-10^{-15}$ [22].333This is equivalent to $c_{S,V,T}>1-5\times
10^{-16}$. This constraint can equivalently be obtained by considerations on
the Weak Cosmic Censorship conjecture under the test particle approximation:
speeds $c_{T}^{2}<1$ would allow for the overcharging of extremal charged
Einstein-aether black holes [62].444In Ref. [62] the constraint is stated to
be $0\leq c_{+}<1$. This corresponds to $c_{T}^{2}\geq 1$, without the window
of width $10^{-15}$ below $1$ allowed for by the Cherenkov constraint.
However, since said reference utilizes the test particle approximation, we
will assume the result to be not perfectly precise so that the small window is
still open. In addition to the constraints on the mode speeds, the Einstein-
aether mode energy densities must be positive [21] and the coefficients of the
time kinetic term must be positive so as to avoid ghosts [6].

Multimessenger observations from the GW signal GW170817 and the gamma-ray
burst GRB 170817A give us the best current observational constraint on the
theory. They constrain the speed of tensor waves to be above $1-3\times
10^{-15}$ and below $1+7\times 10^{-16}$ [36, 6]. However, $c_{T}$ depends
only on $c_{+}$, which means that the “parameter plane” spanned by
$c_{1}+c_{3}=0$ remains unconstrained. Part of the motivation for studying GW
lensing in Einstein-aether theory is a hope of breaking this degeneracy and
restricting parameters directly.

Further observational constraints can be imposed on the Newtonian
gravitational constant $G_{N}$ and cosmological gravitational constant
$G_{\text{cosmo}}$, which in Einstein-aether theory no longer have the same
value [12, 59]. The simple observation that gravity is attractive demands
$G_{N}>0$, although it can be shown that this constraint is automatically
included in the constraint from gravitational Cherenkov radiation [59].
Observations of primordial 4He abundance give rise to a so-called
nucleosynthesis constraint that demands $\left\lvert
G_{\text{cosmo}}/G_{N}\right\rvert-1\lesssim 1/8$ [12, 59].

Another group of observational constraints stems from parametrized post-
Newtonian analysis [63, 64, 65, 66]. The parametrized post-Newtonian formalism
contains a set of parameters called preferred frame parameters, $\alpha_{1}$
and $\alpha_{2}$. They have been constrained to
$\left\lvert\alpha_{1}\right\rvert\leq 10^{-4}$ and
$\left\lvert\alpha_{2}\right\rvert\leq 10^{-7}$ by using lunar laser ranging
and the solar alignment with the ecliptic [58, 67, 13]. These parameters are
zero in GR, but, in general, are nonvanishing in Einstein-aether theory. In
the strong field regime there exist equivalent parameters $\hat{\alpha}_{1}$
and $\hat{\alpha}_{2}$ which can be constrained to
$\left\lvert\hat{\alpha}_{1}\right\rvert\leq 3.5\times 10^{-5}$ and
$\left\lvert\hat{\alpha}_{2}\right\rvert\leq 1.6\times 10^{-9}$ at 95%
confidence using observations of isolated millisecond pulsars [68, 69, 6].

The final constraint is provided by observations of the orbital decay of
binary and trinary pulsars. These observations can provide constraints for
$c_{+}$ and $c_{-}$ [24, 15], but only under the assumption of small
$\alpha_{1,2}$ compared to $c_{+}$ and $c_{-}$. Now that the best constraint
provided by multimessenger observations is much stricter than the conditions
on the preferred frame parameters these constraints have become obsolete [6,
7]. However, $\alpha_{1}$ can still be constrained with pulsar observations.
The constraint obtained from trinary pulsars is about an order of magnitude
stricter than that from lunar laser ranging:
$\left\lvert\alpha_{1}\right\rvert\leq 10^{-5}$ [37].

##  III Gravitational wave lensing

We now study Einstein-aether theory with a non-Minkowski background. First, in
Sec. III.1), we apply some approximations and obtain the new field equations.
Solving the field equations is not trivial hence, in Sec. III.2, we apply to
the propagation eigenstate framework [52]. This framework allows us to study
the mixing of modes resulting from GW propagation on general backgrounds. The
central object of the framework is the so-called kinetic matrix. Eigenvalues
and eigenvectors of the kinetic matrix will be crucial, so we perform an
eigendecomposition in Sec. III.3. Lastly, in Sec. III.4, we apply the
framework to obtain the speeds of the lensed and mixed propagating states.

###  III.1 Approximations and field equations

After reviewing the standard case, i. e., GWs in Einstein-aether without
lensing and with a flat background, we want to focus on the extension where GW
lensing leads to a mixing between modes. To introduce lensing, instead of
using an altogether different background, we leave the Minkowski background
intact, but no longer assume that the background aether is $\delta^{\mu}_{0}$,
instead allowing it to have spatial components. This is an approximation, the
exact meaning of which is best understood through an expansion of the metric
tensor in Riemann normal coordinates around a point $0$ [70]:

$\displaystyle g_{\alpha\beta}\left(x^{\mu}\right)$
$\displaystyle=g_{\alpha\beta}(0)+\frac{1}{3}R_{\alpha\mu\nu\beta}(0)x^{\mu}x^{\nu}+\mathcal{O}\left(x^{3}\right).$
(62)

We find that the leading order correction to spacetime being flat is given by
the curvature at $0$ and the distance $x^{\mu}$ we have traveled from $0$.
Roughly speaking, our approximation holds if curvature is small on the scale
of lensing. Physically, this corresponds to passing far from a point lens
rather than close by it. We focus on a perturbation of the background aether
because, while perturbations of the metric due to lensing have been studied
before, perturbations of the aether field have not. Since the correction to
the aether will then also be small, we go only to first order in the spatial
components of the aether.

We will call the new, nontrivial background aether $\mu^{\alpha}$. We again
linearize around the contravariant background:

$\displaystyle\overline{u}^{\alpha}$ $\displaystyle=\mu^{\alpha}+w^{\alpha}.$
(63)

We will make one more approximation: that derivatives of $\mu^{\alpha}$ are
small compared to derivatives of $w^{\alpha}$. Formally, this means we are
restricting ourselves to the limit of geometric optics or ray optics.
Fluctuations of the background happen on scale of the lens. Fluctuations of
the perturbation happen on scale of the wavelength. Therefore, the geometric
optics limit corresponds to small wavelengths compared to the scale of the
lens. The frequency range of LIGO is between 10 Hz and 10 kHz [71]. Assuming
that GWs propagate roughly with the speed of light,
$\lambda_{\text{GW}}=c/f_{\text{GW}}$ with wavelength $\lambda_{\text{GW}}$
and frequency $f_{\text{GW}}$ gives us a wavelength range of between 30 and
30,000 km, with typical wavelengths in the range of 1,000 km. While black
holes, planets or stars may be too small for our approximation to hold,
galaxies serve as ideal lenses for this approximation. Even when viewed
orthogonal to the disk, the Milky Way is about $10^{13}$ times bigger than the
characteristic wavelength, $\sim$ 1,000 km.555According to Bland-Hawthorn et
al. [72], the exponential scaleheight of the dominant old thin disk of the
Milky Way galaxy at the sun’s location is between 220 and 450 parsecs, which
corresponds to approximately $10^{19}$ m.

The validity of the high-frequency expansion also requires that the time
delays associated to the lens are larger than the GW period, that is,

$f\gg\frac{1}{8\pi GM_{L}(1+z_{L})}\sim
100\text{Hz}\left(\frac{100M_{\odot}}{(1+z_{L})M_{L}}\right)\,,$ (64)

where $(1+z_{L})M_{L}$ is the redshifted effective mass of the lens [73, Sec.
IIA]. For extended lenses with size $\gg GM$ this is a more restrictive
requirement. If it is not fulfilled, a wave-optic treatment is necessary to
describe GW lensing [74, 75, 76].

With this in mind, we can already transfer some results from the unlensed
case. Consider calculating a linearized covariant derivative of the aether,
$\overline{D_{\gamma}u^{\rho}}$. The derivative operator $D_{\gamma}$ has
zeroth and first order contributions: partial derivatives and linearized
Christoffel symbols, respectively. Similarly, the aether has its background
configuration $\mu^{\rho}$ and a perturbation $w^{\rho}$ to it. The
combination of both first order terms would be of second order, leaving us
with three nonvanishing terms:

$\displaystyle\overline{D_{\gamma}u^{\rho}}$
$\displaystyle=\partial_{\gamma}\mu^{\rho}+\partial_{\gamma}w^{\rho}+\overline{\Gamma}^{\rho}_{\gamma\sigma}\mu^{\sigma}.$
(65)

Since we assume $\partial_{\gamma}\mu^{\rho}\ll\partial_{\gamma}w^{\rho}$ the
first term drops out, leaving us only with contributions of first order in the
perturbation. This means that just as in the unlensed case products of aether,
aether acceleration and aether current must drop out and any term multiplied
with them must be of zeroth order:

$\displaystyle\overline{T}^{\alpha\beta}_{\text{ae}}$
$\displaystyle=\mu^{\mu}\partial_{\mu}\overline{J}^{(\alpha\beta)}-\mu^{(\beta}\partial_{\mu}\overline{J}^{\alpha)\mu},$
(66) $\displaystyle\overline{\text{\AE}^{[0}u^{i]}}$
$\displaystyle=\mu^{[i}\partial_{\mu}\overline{J}^{\mu|0]},$ (67)

that, by inserting the aether current, can be rewritten as

$\displaystyle\overline{T}^{\alpha\beta}_{\text{ae}}$
$\displaystyle=M^{(\alpha|\gamma(0)}_{\phantom{(\alpha|\gamma}\nu\rho}\left[\eta^{\beta)\nu}\mu^{\mu}-\eta^{\mu\nu}\mu^{\beta)}\right]\partial_{\mu}\overline{D_{\gamma}u^{\rho}},$
(68) $\displaystyle\overline{\text{\AE}^{[0}u^{i]}}$
$\displaystyle=M^{\mu\gamma(0)}_{\phantom{\mu\gamma}\nu\rho}\mu^{[i}\eta^{0]\nu}\partial_{\mu}\overline{D_{\gamma}u^{\rho}}.$
(69)

To summarize, we make the following three physical approximations:

  1. 1.

weak field approximation: the GW amplitudes are small,

  2. 2.

weak lensing approximation: the GWs pass far from the lens,

  3. 3.

geometric optics approximation: the lens is much larger than the GW
wavelength.

Mathematically, these three approximations mean that:

  1. 1.

we linearize in $h_{\mu\nu}$ and $w_{\mu}$,

  2. 2.

we keep the Minkowski background and linearize in $\mu_{i}$,

  3. 3.

we neglect $\partial_{\gamma}\mu^{\rho}$ compared to
$\partial_{\gamma}w^{\rho}$,

respectively.

Once again our first step, after linearization, is to eliminate $w^{0}$ from
the field equations by using the constraint equation. The linearized
constraint equation takes the following form:

$\displaystyle 0$
$\displaystyle=\overline{g_{\alpha\beta}u^{\alpha}u^{\beta}}+1$
$\displaystyle=(\mu_{k})^{2}-(\mu_{0})^{2}+h_{00}(\mu_{0})^{2}+2\gamma_{k}\mu_{k}\mu^{0}$
$\displaystyle\mathrel{\phantom{=}}+\,\phi_{kl}\mu_{k}\mu_{l}+\frac{1}{2}(\mu_{k})^{2}F-\frac{1}{2}U^{2}f+U^{2}\phi$
$\displaystyle\mathrel{\phantom{=}}-\,2w^{0}\mu^{0}+2v_{k}\mu_{k}+1.$ (70)

Here we defined $U\psi=\mu_{k}\psi_{,k}$, for an arbitrary scalar function
$\psi$. We now go to first order in $\mu_{k}$, which means $(\mu_{0})^{2}=1$
and $(\mu_{k})^{2}=0$. Equation (70) then simplifies to,

$\displaystyle 0$
$\displaystyle=h_{00}+2\gamma_{k}\mu_{k}-2w^{0}+2v_{k}\mu_{k},$ (71)

which we can solve for $w^{0}$:

$\displaystyle w^{0}$
$\displaystyle=\frac{1}{2}h_{00}+\mu_{k}\gamma_{k}+\mu_{k}v_{k}.$ (72)

In the trivial case $\mu_{k}=0$ and we reclaim the trivial constraint equation
$w^{0}=h_{00}/2$; cf. Eq. (LABEL:eq:wmu_and_hmunu).

We now proceed as we did previously: we calculate the linearized aether
stress-energy tensor [Eq. (68)], $\overline{\text{\AE}^{[0}u^{i]}}$ [Eq. (69)]
and, finally, the rather lengthy components of the effective tensor [Eqs. (30)
through (33)]. We find:

$\displaystyle X_{00}$ $\displaystyle=\frac{1}{2}c_{14}\Delta
h_{00}-\frac{1}{2}\Delta F+\frac{1}{2}c_{2}U\dot{F}$
$\displaystyle\mathrel{\phantom{=}}-\frac{1}{2}(c_{-}-c_{2}+2c_{4})U\Delta\dot{\phi}$
$\displaystyle\mathrel{\phantom{=}}+\frac{1}{2}(c_{+}+2c_{4})\mu_{k}\Delta\gamma_{k}+c_{3}\mu_{k}\Delta
v_{k},$ (73) $\displaystyle X_{0i}$
$\displaystyle=\frac{1}{2}c_{14}\dot{h}_{00,i}-\frac{1}{4}\left(c_{+}+2c_{4}\right)\mu_{i}\Delta
h_{00}$
$\displaystyle\mathrel{\phantom{=}}+\frac{1}{4}(c_{+}+2c_{4})Uh_{00,i}-\frac{1}{2}\dot{F}_{,i}+\frac{1}{4}c_{-}\mu_{i}\Delta
F$
$\displaystyle\mathrel{\phantom{=}}-\frac{1}{4}(c_{-}-2c_{2}+2c_{4})\mu_{i}\ddot{F}-\frac{1}{4}c_{-}UF_{,i}$
$\displaystyle\mathrel{\phantom{=}}+\frac{1}{4}(c_{-}+2c_{4})U\ddot{f}_{,i}+\frac{1}{2}c_{2}\mu_{i}\Delta\ddot{\phi}$
$\displaystyle\mathrel{\phantom{=}}-\frac{1}{2}(c_{-}+2c_{4})U\ddot{\phi}_{,i}-\frac{1}{2}(1-c_{-})\Delta\gamma_{i}-c_{14}\ddot{\gamma}_{i}$
$\displaystyle\mathrel{\phantom{=}}+\frac{1}{2}(c_{+}+2c_{4})\mu_{k}\dot{\gamma}_{k,i}-(c_{+}+2c_{4})U\dot{\gamma}_{i}$
$\displaystyle\mathrel{\phantom{=}}+\frac{1}{2}c_{-}\Delta
v_{i}-c_{14}\ddot{v}_{i}+c_{3}\mu_{k}\dot{v}_{k,i}$
$\displaystyle\mathrel{\phantom{=}}-\frac{1}{2}(c_{+}+4c_{4})U\dot{v}_{i}+\frac{1}{2}c_{-}\mu_{k}\Delta\phi_{ik}$
$\displaystyle\mathrel{\phantom{=}}-\frac{1}{2}(c_{-}+2c_{4})\mu_{k}\ddot{\phi}_{ik},$
(74) $\displaystyle X_{i0}$
$\displaystyle=-\frac{1}{4}\left(c_{-}+2c_{4}\right)\mu_{i}\Delta h_{00}$
$\displaystyle\mathrel{\phantom{=}}+\frac{1}{4}(c_{+}+2c_{2})Uh_{00,i}-\frac{1}{2}(1+c_{2})\dot{F}_{,i}$
$\displaystyle\mathrel{\phantom{=}}-\frac{1}{2}c_{2}UF_{,i}-\frac{1}{2}c_{123}\Delta\dot{\phi}_{,i}-\frac{1}{2}c_{123}U\Delta\phi_{,i}$
$\displaystyle\mathrel{\phantom{=}}-\frac{1}{2}\Delta\gamma_{i}-c_{2}\mu_{k}\dot{\gamma}_{k,i}-\frac{1}{2}c_{+}U\dot{\gamma}_{i}-\frac{1}{2}c_{+}\Delta
v_{i}$
$\displaystyle\mathrel{\phantom{=}}-c_{2}\mu_{k}\dot{v}_{k,i}-\frac{1}{2}c_{+}U\dot{v}_{i},$
(75) $\displaystyle X_{ij}$
$\displaystyle=-\frac{1}{2}P_{ij}h_{00}+\frac{1}{4}(c_{3}-c_{4})\mu_{i}\dot{h}_{00,j}$
$\displaystyle\mathrel{\phantom{=}}+\frac{1}{4}(c_{3}-c_{4})\mu_{j}\dot{h}_{00,i}+\frac{1}{2}c_{2}\delta_{ij}U\dot{h}_{00}$
$\displaystyle\mathrel{\phantom{=}}+\frac{1}{4}P_{ij}F-\frac{1}{4}(c_{+}+2c_{2})\delta_{ij}\ddot{F}-\frac{1}{4}P_{ij}\ddot{f}$
$\displaystyle\mathrel{\phantom{=}}+\frac{1}{4}(c_{+}-2)\ddot{f}_{,ij}+\frac{1}{4}c_{2}\mu_{i}\dot{F}_{,j}+\frac{1}{4}c_{2}\mu_{j}\dot{F}_{,i}$
$\displaystyle\mathrel{\phantom{=}}-\frac{1}{2}(c_{+}+2c_{2})\delta_{ij}U\dot{F}+\frac{1}{2}c_{+}U\dot{f}_{,ij}-\frac{1}{2}P_{ij}\ddot{\phi}$
$\displaystyle\mathrel{\phantom{=}}-\frac{1}{2}c_{2}\delta_{ij}\Delta\ddot{\phi}-\frac{1}{2}c_{+}\ddot{\phi}_{,ij}+\frac{1}{4}c_{123}\mu_{i}\Delta\dot{\phi}_{,j}$
$\displaystyle\mathrel{\phantom{=}}+\frac{1}{4}c_{123}\mu_{j}\Delta\dot{\phi}_{,i}-c_{2}U\delta_{ij}\Delta\dot{\phi}-c_{+}U\dot{\phi}_{,ij}$
$\displaystyle\mathrel{\phantom{=}}-\frac{1}{2}\dot{\gamma}_{i,j}-\frac{1}{2}\dot{\gamma}_{j,i}-\frac{1}{4}c_{-}\mu_{i}\Delta\gamma_{j}-\frac{1}{4}c_{-}\mu_{j}\Delta\gamma_{i}$
$\displaystyle\mathrel{\phantom{=}}-\frac{1}{2}(c_{3}-c_{4})\mu_{i}\ddot{\gamma}_{j}-\frac{1}{2}(c_{3}-c_{4})\mu_{j}\ddot{\gamma}_{i}$
$\displaystyle\mathrel{\phantom{=}}-c_{2}\delta_{ij}\mu_{k}\ddot{\gamma}_{k}-\frac{1}{2}c_{+}\dot{v}_{i,j}-\frac{1}{2}c_{+}\dot{v}_{j,i}$
$\displaystyle\mathrel{\phantom{=}}+\frac{1}{2}c_{3}\mu_{i}\Delta
v_{j}+\frac{1}{2}c_{3}\mu_{j}\Delta
v_{i}-\frac{1}{2}(c_{3}-c_{4})\mu_{i}\ddot{v}_{j}$
$\displaystyle\mathrel{\phantom{=}}-\frac{1}{2}(c_{3}-c_{4})\mu_{j}\ddot{v}_{i}-c_{2}\delta_{ij}\mu_{k}\ddot{v}_{k}-\frac{1}{2}c_{+}Uv_{i,j}$
$\displaystyle\mathrel{\phantom{=}}-\frac{1}{2}c_{+}Uv_{j,i}-\frac{1}{2}\Delta\phi_{ij}+\frac{1}{2}(1-c_{+})\ddot{\phi}_{ij}$
$\displaystyle\mathrel{\phantom{=}}-c_{+}U\dot{\phi}_{ij},$ (76)

The Helmholtz decomposition of $X_{0i}$ and $X_{i0}$ requires the
decomposition of $\mu_{i}$ into transverse and longitudinal parts, but is
otherwise simple. The SVT decomposition of $X_{ij}$ is less straightforward:
some terms simply do not fit into our scheme of having every term be
transverse traceless, a symmetrized derivative of a transverse vector or
expressed in terms of $\partial_{i}\partial_{j}$ or $P_{ij}$. However, while
the SVT decomposition is difficult if not impossible in real space it becomes
possible in Fourier space.666Specifically, the Fourier space SVT decomposition
contains terms proportional to $k^{-2}$. These terms are the problem with a
real space SVT decomposition: they are nonlocal and can only be described with
inverse operators. We work in Fourier space to proceed. We define the Fourier
transform $\mathcal{F}$ and the inverse Fourier transform $\mathcal{F}^{-1}$
of a function $\psi$ as,

$\displaystyle\mathcal{F}[\psi]$ $\displaystyle=\int{\rm d}^{4}x\,\psi
e^{-ix_{\mu}k^{\mu}},$ (77) $\displaystyle\mathcal{F}^{-1}[\psi]$
$\displaystyle=\frac{1}{2\pi}\int{\rm d}^{4}k\,\psi e^{ix_{\mu}k^{\mu}},$ (78)

and obtain the same four constraints we utilized previously up to first order
in $\mu_{k}$:

$\displaystyle h_{00}$
$\displaystyle=\frac{F}{c_{14}}+\frac{2c_{14}(1+c_{2})-c_{123}}{c_{14}c_{123}}\frac{\omega\mu_{l}k_{l}}{k^{2}}F$
$\displaystyle\mathrel{\phantom{=}}+\frac{2c_{3}-(c_{+}+2c_{4})c_{+}}{c_{14}c_{+}}\mu_{k}\gamma_{k}$
(79) $\displaystyle\phi$
$\displaystyle=\frac{1+c_{2}}{c_{123}}\frac{F}{k^{2}}+\frac{1}{c_{14}}\frac{\mu_{l}k_{l}}{\omega
k^{2}}F-\frac{2c_{2}(1-c_{+})}{c_{123}c_{+}}\mu_{k}\gamma_{k},$ (80)
$\displaystyle v_{i}$
$\displaystyle=-\frac{c_{-}+2c_{4}}{2c_{+}c_{14}}\mu_{i}^{T}F-\frac{\gamma_{i}}{c_{+}}-\frac{1-c_{+}}{c_{+}}\frac{\omega\mu_{l}k_{l}}{k^{2}}\gamma_{i}.$
(81)

The remaining equations will serve as the gateway to utilizing what is known
as the propagation eigenstate framework [52] and finding the mixing between
the modes. The transverse part of $X_{0i}=0$ carries the first two DOF:

$\displaystyle 0$
$\displaystyle=-\left(\frac{c_{+}+2c_{4}}{2c_{14}}-\frac{c_{-}}{2}+c_{-}\frac{c_{-}+2c_{4}}{2c_{+}c_{14}}\right)\mu_{i}^{T}k^{2}F$
$\displaystyle\mathrel{\phantom{=}}-\left(\frac{c_{-}-2c_{2}+2c_{4}}{2}+c_{2}\frac{1+c_{2}}{c_{123}}-\frac{c_{-}+2c_{4}}{c_{+}}\right)\mu_{i}^{T}\omega^{2}F$
$\displaystyle\mathrel{\phantom{=}}-\frac{2c_{1}-c_{+}c_{-}}{c_{+}}k^{2}\gamma_{i}+\frac{2c_{14}(1-c_{+})}{c_{+}}\left(1+\frac{\omega\mu_{l}k_{l}}{k^{2}}\right)\omega^{2}\gamma_{i}$
$\displaystyle\mathrel{\phantom{=}}+\left[2(c_{+}+2c_{4})-\frac{c_{+}+4c_{4}}{c_{+}}-c_{-}\frac{1-c_{+}}{c_{+}}\right]\omega\mu_{l}k_{l}\gamma_{i}$
$\displaystyle\mathrel{\phantom{=}}+c_{-}\mu_{k}k^{2}\phi_{ik}-(c_{-}+2c_{4})\omega^{2}\mu_{k}\phi_{ik}.$
(82)

Inserting the constraints into the transverse traceless part of $X_{ij}=0$
yields equations for two more DOF:

$\displaystyle 0$
$\displaystyle=-\frac{2c_{3}+c_{+}c_{-}}{2c_{+}}k^{2}\left(\mu_{i}^{T}\gamma_{j}+\mu_{j}^{T}\gamma_{i}\right)$
$\displaystyle\mathrel{\phantom{=}}+\frac{2c_{3}+c_{+}c_{-}}{2c_{+}}\left(\delta_{ij}k^{2}-k_{i}k_{j}\right)\mu_{k}\gamma_{k}$
$\displaystyle\mathrel{\phantom{=}}+\frac{(1-c_{+})(c_{3}-c_{4})}{c_{+}}\omega^{2}\left(\mu_{i}^{T}\gamma_{j}+\mu_{j}^{T}\gamma_{i}\right)$
$\displaystyle\mathrel{\phantom{=}}-\frac{(1-c_{+})(c_{3}-c_{4})}{c_{+}}\frac{\omega^{2}}{k^{2}}\left(\delta_{ij}k^{2}-k_{i}k_{j}\right)\mu_{k}\gamma_{k}$
$\displaystyle\mathrel{\phantom{=}}-k^{2}\phi_{ij}+(1-c_{+})\omega^{2}\phi_{ij}+2c_{+}\omega\mu_{l}k_{l}\phi_{ij}.$
(83)

Finally, we can insert our constraints into the transverse trace part of
$X_{ij}=0$:

$\displaystyle 0$
$\displaystyle=\left(\frac{1}{c_{14}}-\frac{1}{2}\right)F-\frac{(1-c_{+})(2+c_{+}+3c_{2})}{2c_{123}}\frac{\omega^{2}}{k^{2}}F$
$\displaystyle\mathrel{\phantom{=}}-\left[c_{+}+2c_{2}-\frac{2(1+c_{2})^{2}}{c_{123}}+\frac{2}{c_{14}}\right]\frac{\omega\mu_{l}k_{l}}{k^{2}}F$
$\displaystyle\mathrel{\phantom{=}}+\left[\frac{2c_{3}-(c_{+}+2c_{4})c_{+}}{c_{14}c_{+}}+\frac{c_{-}}{2}+\frac{c_{3}}{c_{+}}\right]\mu_{k}\gamma_{k}$
$\displaystyle\mathrel{\phantom{=}}+\bigg{[}\frac{2c_{2}(1-c_{+})(1+c_{2})}{c_{123}c_{+}}$
$\displaystyle\mathrel{\phantom{=+\bigg{[}}}+(2c_{2}+c_{3}-c_{4})\left(1-\frac{1}{c_{+}}\right)\bigg{]}\frac{\omega^{2}}{k^{2}}\mu_{k}\gamma_{k}.$
(84)

###  III.2 Propagation eigenstate framework

To solve our field equations with the propagation eigenstate framework, we
must first

  1. 1.

normalize our equations so the trivial limit is of the form
$0=c_{p}^{2}\Delta\psi-\ddot{\psi}$

  2. 2.

convert our equations into a matrix equation

We now explicitly consider the $z$ direction as propagation direction. In this
case, all transverse vectors have components only in $x$ or $y$ direction.
Applying this process, we arrive at a matrix equation containing a 5×5 matrix:

$\displaystyle\left(\begin{array}[]{ccccccccc}K_{SS}&&K_{SV}\mu_{1}&&K_{SV}\mu_{2}&&0&&0\\\
\\\ K_{VS}\mu_{1}&&K_{VV}&&0&&K_{VT}\mu_{1}&&K_{VT}\mu_{2}\\\ \\\
K_{VS}\mu_{2}&&0&&K_{VV}&&-K_{VT}\mu_{2}&&K_{VT}\mu_{1}\\\ \\\
0&&K_{TV}\mu_{1}&&-K_{TV}\mu_{2}&&K_{TT}&&0\\\ \\\
0&&K_{TV}\mu_{2}&&K_{TV}\mu_{1}&&0&&K_{TT}\end{array}\right)\left(\begin{array}[]{c}F\\\
\\\ \gamma_{1}\\\ \\\ \gamma_{2}\\\ \\\ \phi_{11}\\\ \\\
\phi_{12}\end{array}\right)$ $\displaystyle=0.$ (103)

The entries of the matrix are second-order derivative operators:

$\displaystyle K_{SS}$ $\displaystyle=\omega^{2}-c_{S}^{2}k^{2}+d_{1}\omega
k\mu_{3},$ (104) $\displaystyle K_{SV}$
$\displaystyle=-d_{2}k^{2}-d_{3}\omega^{2},$ (105) $\displaystyle K_{VS}$
$\displaystyle=-a_{1}k^{2}-a_{2}\omega^{2},$ (106) $\displaystyle K_{VV}$
$\displaystyle=\omega^{2}-c_{V}^{2}k^{2}+\frac{\omega^{3}}{k}\mu_{3}+a_{3}\omega
k\mu_{3},$ (107) $\displaystyle K_{VT}$
$\displaystyle=-a_{4}k^{2}-a_{5}\omega^{2},$ (108) $\displaystyle K_{TV}$
$\displaystyle=-b_{1}k^{2}-b_{2}\omega^{2},$ (109) $\displaystyle K_{TT}$
$\displaystyle=\omega^{2}-c_{T}^{2}k^{2}+b_{3}\omega k\mu_{3}.$ (110)

Each coefficient is given by the theory parameters as follows:

$\displaystyle a_{1}$
$\displaystyle=\frac{c_{+}}{2c_{14}(1-c_{+})}\bigg{(}\frac{c_{+}+2c_{4}}{2c_{14}}-\frac{c_{-}}{2}+c_{-}\frac{c_{-}+2c_{4}}{2c_{+}c_{14}}\bigg{)},$
(111) $\displaystyle a_{2}$ $\displaystyle=\frac{c_{+}}{2c_{14}(1-c_{+})}$
$\displaystyle\mathrel{\phantom{=}}\times\left(\frac{c_{-}-2c_{2}+2c_{4}}{2}+c_{2}\frac{1+c_{2}}{c_{123}}-\frac{c_{-}+2c_{4}}{c_{+}}\right)$
(112) $\displaystyle a_{3}$ $\displaystyle=\frac{c_{+}}{2c_{14}(1-c_{+})}$
$\displaystyle\mathrel{\phantom{=}}\times\left[2(c_{+}+2c_{4})-\frac{c_{+}+4c_{4}}{c_{+}}-c_{-}\frac{1-c_{+}}{c_{+}}\right],$
(113) $\displaystyle a_{4}$
$\displaystyle=-c_{-}\frac{c_{+}}{2c_{14}(1-c_{+})},$ (114) $\displaystyle
a_{5}$ $\displaystyle=(c_{-}+2c_{4})\frac{c_{+}}{2c_{14}(1-c_{+})},$ (115)
$\displaystyle b_{1}$
$\displaystyle=\frac{2c_{3}+c_{+}c_{-}}{2c_{+}(1-c_{+})},$ (116)
$\displaystyle b_{2}$ $\displaystyle=-\frac{c_{3}-c_{4}}{c_{+}},$ (117)
$\displaystyle b_{3}$ $\displaystyle=\frac{2c_{+}}{1-c_{+}},$ (118)
$\displaystyle d_{1}$
$\displaystyle=\frac{2c_{123}}{(1-c_{+})(2+c_{+}+3c_{2})}$
$\displaystyle\mathrel{\phantom{=}}\times\left[c_{+}+2c_{2}-\frac{2(1+c_{2})^{2}}{c_{123}}+\frac{2}{c_{14}}\right],$
(119) $\displaystyle d_{2}$
$\displaystyle=\frac{2c_{123}}{(1-c_{+})(2+c_{+}+3c_{2})}$
$\displaystyle\mathrel{\phantom{=}}\times\left[\frac{2c_{3}-c_{+}(c_{+}+2c_{4})}{c_{14}c_{+}}+\frac{c_{-}}{2}+\frac{c_{3}}{c_{+}}\right],$
(120) $\displaystyle d_{3}$
$\displaystyle=\frac{2c_{123}}{(1-c_{+})(2+c_{+}+3c_{2})}$
$\displaystyle\mathrel{\phantom{=}}\times\bigg{[}\frac{2c_{2}(1-c_{+})(1+c_{2})}{c_{123}c_{+}}$
$\displaystyle\mathrel{\phantom{=\times\bigg{[}}}+(2c_{2}+c_{3}-c_{4})\left(1-\frac{1}{c_{+}}\right)\bigg{]}.$
(121)

For this matrix equation and related expressions we switch from index notation
for tensors to using a more classical vector arrow and circumflex for vectors
and matrices, respectively. We refer to the matrix as the kinetic matrix
$\hat{K}$. The vector $(F,\gamma_{1},\gamma_{2},\phi_{11},\phi_{12})$ is
called $\vec{h}$. Its five entries are called the interaction eigenstates.

Note that, if we were to go beyond leading order in the geometric-optics
approximation , other matrices would appear in the total matrix equation: the
amplitude matrix $\hat{A}$ containing first-order derivative operators and the
mass matrix $\hat{M}$ containing no derivatives. In that case, the
considerations that we apply to the kinetic matrix need to be applied to the
sum of all three matrices $\hat{D}$ [77].

With $\hat{K}$ and $\vec{h}$ as a starting point we can now formally introduce
the propagation eigenstate framework. Since $\hat{K}$ is not diagonal the
interaction eigenstates do no longer obey wave equations; they themselves do
no longer propagate. We now want to find the propagation eigenstates $\vec{H}$
that obey decoupled wave equations, the solutions of which describe
propagation. To that end we seek to diagonalize $\hat{K}$. We do this by
applying a matrix $\hat{M}$ to it so that $\hat{A}=\hat{M}\hat{K}\hat{M}^{-1}$
is diagonal. We can now modify the equation $\hat{K}\vec{h}=0$. We begin by
multiplying with $\hat{M}$ from the left:

$\displaystyle\hat{M}\hat{K}\vec{h}=0.$ (122)

Now, $\hat{M}^{-1}\hat{M}$ is the 5×5 identity matrix, which we can freely
insert after $\hat{K}$:

$\displaystyle\hat{M}\hat{K}\hat{M}^{-1}\hat{M}\vec{h}=0.$ (123)

We now define $\vec{H}=\hat{M}\vec{h}$:

$\displaystyle\hat{A}\vec{H}=0.$ (124)

Therefore, we find that the propagation eigenstates are given by
$\hat{M}\vec{h}$. The matrix $\hat{M}$ is called the mixing matrix because it
describes how the interaction eigenstates mix into the propagation
eigenstates. Since $\hat{M}$ diagonalizes $\hat{K}$ it contains the five
eigenvectors of $\hat{K}$ as row vectors. The $n$th propagation eigenstate can
then be calculated by taking the scalar product of the $n$th eigenvector of
$\hat{K}$ with $\hat{h}$.

We also want to calculate with what speed the propagation eigenstates
propagate. To that end we consider the entries of $\hat{A}$. They are the
eigenvalues of $\hat{K}$. Since they describe wave equations it must be
possible to bring them into the form $\omega^{2}-c_{X}^{2}k^{2}$ with some
mode speed $c_{X}$. By setting this to zero and solving for $\omega/k$ we can
recover $c_{X}$. In other words, setting the $n$th eigenvalue to zero and
solving for $\omega/k$ yields the speed of the $n$th propagation eigenstate.

###  III.3 Eigendecomposition of the kinetic matrix

There are two equations that govern the eigendecomposition of a matrix. The
first one is the equation defining a pair of eigenvalue $E$ and eigenvectors
$\vec{V}$:

$\displaystyle E\vec{V}$ $\displaystyle=\hat{K}\vec{V}.$ (125)

The second one is the condition that the eigenvectors are normalized. This
means that the scalar product of an eigenvector with itself must be $1$:

$\displaystyle\vec{V}\cdot\vec{V}$ $\displaystyle=1.$ (126)

We now split eigenvector, eigenvalue and kinetic matrix into zeroth and first
order in $\mu_{k}$. We will denote the zeroth order in $\mu_{k}$ by
superscript $0\mu$, the first order by $1\mu$ and so on. The two equations we
must solve can then also be split by order:

$\displaystyle E^{0\mu}\vec{V}^{0\mu}$
$\displaystyle=\hat{K}^{0\mu}\vec{V}^{0\mu},$ (127) $\displaystyle
E^{0\mu}\vec{V}^{1\mu}+E^{1\mu}\vec{V}^{0\mu}$
$\displaystyle=\hat{K}^{0\mu}\vec{V}^{1\mu}+\hat{K}^{1\mu}\vec{V}^{0\mu},$
(128) $\displaystyle\vec{V}^{0\mu}\cdot\vec{V}^{0\mu}$ $\displaystyle=1,$
(129) $\displaystyle\vec{V}^{0\mu}\cdot\vec{V}^{1\mu}$ $\displaystyle=0.$
(130)

We note that $\hat{K}^{0\mu}$ only has diagonal elements. It follows from the
zeroth-order equations that the zeroth-order eigenvalues $E^{0\mu}$ are its
diagonal elements: These are: $\omega^{2}-c_{S}^{2}k^{2}$, two instances of
$\omega^{2}-c_{V}^{2}k^{2}$, and two instances of $\omega^{2}-c_{T}^{2}k^{2}$.
The corresponding zeroth-order eigenvectors are $\left(1,0,0,0,0\right)$,
$\left(0,1,0,0,0\right)$ and so forth. As expected, without the background
aether for mixing we recover the interaction eigenstates $F$, $\gamma_{1}$,
$\gamma_{2}$, $\phi_{11}$ and $\phi_{12}$ with their respective speeds; cf.
Eqs. (49), (50) and (51). We are left with the first order equations which ,
as an example, we now solve for $\vec{V}^{0\mu}=\left(0,1,0,0,0\right)$. We
first insert the zeroth-order eigenvalue and eigenvector into the first order
equations:

$\displaystyle K_{VV}^{0\mu}\left(\begin{array}[]{c}V_{1}\\\ \\\ V_{2}\\\ \\\
V_{3}\\\ \\\ V_{4}\\\ \\\
V_{5}\end{array}\right)^{1\mu}+E^{1\mu}\left(\begin{array}[]{c}0\\\ \\\ 1\\\
\\\ 0\\\ \\\ 0\\\ \\\ 0\end{array}\right)$
$\displaystyle=\left(\begin{array}[]{ccccccccc}K_{SS}&&0&&0&&0&&0\\\ \\\
0&&K_{VV}&&0&&0&&0\\\ \\\ 0&&0&&K_{VV}&&0&&0\\\ \\\ 0&&0&&0&&K_{TT}&&0\\\ \\\
0&&0&&0&&0&&K_{TT}\end{array}\right)^{0\mu}\left(\begin{array}[]{c}V_{1}\\\
\\\ V_{2}\\\ \\\ V_{3}\\\ \\\ V_{4}\\\ \\\ V_{5}\end{array}\right)^{1\mu}$
(167)
$\displaystyle\mathrel{\phantom{=}}+\left(\begin{array}[]{ccccccccc}K_{SS}^{1\mu}&K_{SV}\mu_{1}&K_{SV}\mu_{2}&0&0\\\
\\\ K_{VS}\mu_{1}&K_{VV}^{1\mu}&0&K_{VT}\mu_{1}&K_{VT}\mu_{2}\\\ \\\
K_{VS}\mu_{2}&0&K_{VV}^{1\mu}&-K_{VT}\mu_{2}&K_{VT}\mu_{1}\\\ \\\
0&K_{TV}\mu_{1}&-K_{TV}\mu_{2}&K_{TT}^{1\mu}&0\\\ \\\
0&K_{TV}\mu_{2}&K_{TV}\mu_{1}&0&K_{TT}^{1\mu}\end{array}\right)\left(\begin{array}[]{c}0\\\
\\\ 1\\\ \\\ 0\\\ \\\ 0\\\ \\\ 0\end{array}\right),$ (186)

$\displaystyle\left(\begin{array}[]{c}0\\\ \\\ 1\\\ \\\ 0\\\ \\\ 0\\\ \\\
0\end{array}\right)\cdot\left(\begin{array}[]{c}V_{1}\\\ \\\ V_{2}\\\ \\\
V_{3}\\\ \\\ V_{4}\\\ \\\ V_{5}\end{array}\right)^{1\mu}$ $\displaystyle=0.$
(205)

The second equation simply tells us that $V_{2}^{1\mu}=0$ , which we then
insert into Eq. (186) and simplify. These steps lead us to:

$\displaystyle\left(\begin{array}[]{c}K_{VV}^{0\mu}V_{1}^{1\mu}\\\ \\\ 0\\\
\\\ K_{VV}^{0\mu}V_{3}^{1\mu}\\\ \\\ K_{VV}^{0\mu}V_{4}^{1\mu}\\\ \\\
K_{VV}^{0\mu}V_{5}^{1\mu}\end{array}\right)+\left(\begin{array}[]{c}0\\\ \\\
E^{1\mu}\\\ \\\ 0\\\ \\\ 0\\\ \\\ 0\end{array}\right)$
$\displaystyle=\left(\begin{array}[]{c}K_{SS}^{0\mu}V_{1}^{1\mu}\\\ \\\ 0\\\
\\\ K_{VV}^{0\mu}V_{3}^{1\mu}\\\ \\\ K_{TT}^{0\mu}V_{4}^{1\mu}\\\ \\\
K_{TT}^{0\mu}V_{5}^{1\mu}\end{array}\right)+\left(\begin{array}[]{c}K_{SV}\mu_{1}\\\
\\\ K_{VV}^{1\mu}\\\ \\\ 0\\\ \\\ K_{TV}\mu_{1}\\\ \\\
K_{TV}\mu_{2}\end{array}\right).$ (242)

The second row of the foregoing equation tells us that
$E^{1\mu}=K_{TT}^{1\mu}$. The third row is trivially true, meaning that we
have no information about $V_{2}^{1\mu}$. The remaining rows can be solved for
$V_{1}^{1\mu}$, $V_{4}^{1\mu}$ and $V_{5}^{1\mu}$:

$\displaystyle V_{1}^{1\mu}$
$\displaystyle=\frac{K_{SV}}{K_{VV}^{0\mu}-K_{SS}^{0\mu}}\mu_{1},$ (243)
$\displaystyle V_{4}^{1\mu}$
$\displaystyle=\frac{K_{TV}}{K_{VV}^{0\mu}-K_{TT}^{0\mu}}\mu_{1},$ (244)
$\displaystyle V_{5}^{1\mu}$
$\displaystyle=\frac{K_{TV}}{K_{VV}^{0\mu}-K_{TT}^{0\mu}}\mu_{2}.$ (245)

These expressions are defined in terms of the elements of the kinetic matrix
$\hat{K}$, which are defined in terms of $a_{i}$, $b_{i}$ and $d_{i}$ which in
turn are defined in terms of the $c_{i}$; recall Eqs. (104) through (121). We
can simplify Eqs. (243) through (245) by inserting all these definitions into
each other. Additionally, all fractions are multiplied with terms of first
order in $\mu_{k}$. This means we can approximate the fractions to zeroth
order in $\mu_{k}$ and replace $\omega^{2}/k^{2}$ with the correct mode speed.
Since we are discussing a vector mode we replace it with $c_{V}^{2}$. Carrying
out these steps we arrive at simple results, e. g., $V_{4}^{1\mu}=-\mu_{1}$,
and we can calculate all five eigenvalues and eigenvectors:

$\displaystyle E\in\bigg{\\{}$
$\displaystyle\omega^{2}-c_{S}^{2}k^{2}+d_{1}\omega
k\mu_{3},\omega^{2}-c_{V}^{2}k^{2}+\frac{\omega^{3}}{k}\mu_{3}+a_{3}\omega
k\mu_{3},\omega^{2}-c_{V}^{2}k^{2}+\frac{\omega^{3}}{k}\mu_{3}+a_{3}\omega
k\mu_{3},$ $\displaystyle\omega^{2}-c_{T}^{2}k^{2}+b_{3}\omega
k\mu_{3},\omega^{2}-c_{T}^{2}k^{2}+b_{3}\omega k\mu_{3}\bigg{\\}},$ (246)
$\displaystyle\vec{V}\in$ $\displaystyle\left\\{\left(\begin{array}[]{c}1\\\
\\\ \dfrac{c_{14}(c_{+}-2)+2c_{+}}{2c_{14}(1-c_{+})}\mu_{1}\\\ \\\
\dfrac{c_{14}(c_{+}-2)+2c_{+}}{2c_{14}(1-c_{+})}\mu_{2}\\\ \\\ 0\\\ \\\
0\end{array}\right),\left(\begin{array}[]{c}-2\mu_{1}\\\ \\\ 1\\\ \\\
V_{3}^{1\mu}\\\ \\\ -\mu_{1}\\\ \\\
-\mu_{2}\end{array}\right),\left(\begin{array}[]{c}-2\mu_{2}\\\ \\\
V_{2}^{1\mu}\\\ \\\ 1\\\ \\\ \mu_{2}\\\ \\\
-\mu_{1}\end{array}\right),\left(\begin{array}[]{c}0\\\ \\\
\dfrac{c_{+}}{1-c_{+}}\mu_{1}\\\ \\\ -\dfrac{c_{+}}{1-c_{+}}\mu_{2}\\\ \\\
1\\\ \\\ V_{5}^{1\mu}\end{array}\right),\left(\begin{array}[]{c}0\\\ \\\
\dfrac{c_{+}}{1-c_{+}}\mu_{2}\\\ \\\ \dfrac{c_{+}}{1-c_{+}}\mu_{1}\\\ \\\
V_{4}^{1\mu}\\\ \\\ 1\end{array}\right)\right\\}.$ (292)

The propagation eigenstates are:

$\displaystyle H_{1}$
$\displaystyle=F+\frac{c_{14}(c_{+}-2)+2c_{+}}{2c_{14}(1-c_{+})}(\mu_{1}\gamma_{1}+\mu_{2}\gamma_{2}),$
(293) $\displaystyle H_{2}$
$\displaystyle=-2\mu_{1}F+\gamma_{1}+V_{3}^{1\mu}\gamma_{2}-\mu_{1}\phi_{11}-\mu_{2}\phi_{12},$
(294) $\displaystyle H_{3}$
$\displaystyle=-2\mu_{2}F+V_{2}^{1\mu}\gamma_{1}+\gamma_{2}+\mu_{2}\phi_{11}-\mu_{1}\phi_{12},$
(295) $\displaystyle H_{4}$
$\displaystyle=\frac{c_{+}}{1-c_{+}}(\mu_{1}\gamma_{1}-\mu_{2}\gamma_{2})+\phi_{11}+V_{5}^{1\mu}\phi_{12},$
(296) $\displaystyle H_{5}$
$\displaystyle=\frac{c_{+}}{1-c_{+}}(\mu_{2}\gamma_{1}+\mu_{1}\gamma_{2})+V_{4}^{1\mu}\phi_{11}+\phi_{12}.$
(297)

We find that the scalar mode and tensor modes both mix with the vector modes
but not with each other. Mixing between the two vector modes and the two
tensor modes is determined by constants of first order in $\mu_{k}$. We cannot
determine the values of the $V_{i}^{1\mu}$ without going to higher orders in
$\mu_{k}$,777The approach we use for determining eigenvalues and eigenvectors
is essentially perturbation theory that has both degenerate eigenvalues and a
nonsymmetric perturbation. At first order in $\mu_{k}$ we cannot determine the
$V_{i}^{1\mu}$ because the degeneracies are not lifted. At any higher order
$n$ we can no longer assume the $n-1$ order contribution to the eigenvectors
to be orthogonal. Due to these two issues standard perturbation theory methods
fail. but we also do not need to: they do not affect the eigenvalues and the
mode speeds that will lead to observable time lags. These coefficients amount
to a rotation of the two polarizations of equal spin within the propagation
eigenstate to which they dominate: that is, $V_{4}^{1\mu},V_{5}^{1\mu}$ change
the relative contribution of $\phi_{11},\phi_{12}$ in the mostly-tensor
states, $H_{4},H_{5}$ respectively (and similarly for vector polarizations
$\gamma_{1,2}$ within the mostly-vector eigenstates $H_{2},H_{3}$). This
suggests that the effect of the $V_{i}^{1\mu}$ is degenerate with the
polarization angle, and would require additional information to be probed.

###  III.4 Propagation eigenstate speeds

To find the propagation eigenstate speeds we set the eigenvalues of $\hat{K}$
to zero. This is equivalent to solving the wave equation or to finding the
dispersion relation in the high frequency limit. Similarly to the previous
section we now expand $\omega$ as $\omega^{0\mu}+\omega^{1\mu}$. The new
zeroth and first orders of the three nondegenerate eigenvalues are:

$\displaystyle E^{0\mu}$
$\displaystyle\in\bigg{\\{}\left(\omega^{0\mu}\right)^{2}-c_{S}^{2}k^{2},$
$\displaystyle\mathrel{\phantom{\in\bigg{\\{}}}\left(\omega^{0\mu}\right)^{2}-c_{V}^{2}k^{2},$
$\displaystyle\mathrel{\phantom{\in\bigg{\\{}}}\left(\omega^{0\mu}\right)^{2}-c_{T}^{2}k^{2}\bigg{\\}},$
(298) $\displaystyle E^{1\mu}$
$\displaystyle\in\bigg{\\{}2\omega^{0\mu}\omega^{1\mu}+d_{1}\omega^{0\mu}k\mu_{3},$
$\displaystyle\mathrel{\phantom{\in\bigg{\\{}}}2\omega^{0\mu}\omega^{1\mu}+\frac{\left(\omega^{0\mu}\right)^{3}}{k}\mu_{3}+a_{3}\omega^{0\mu}k\mu_{3},$
$\displaystyle\mathrel{\phantom{\in\bigg{\\{}}}2\omega^{0\mu}\omega^{1\mu}+b_{3}\omega^{0\mu}k\mu_{3}\bigg{\\}}.$
(299)

Notably this expansion does not change our results. We would arrive at the
same speeds even if we had expanded $\omega$ at an earlier step. Setting the
zeroth order eigenvalues to zero we, unsurprisingly, arrive at
$\omega^{0\mu}/k=c_{S}$, etc. Inserting the zeroth order solutions into the
first order yields

$\displaystyle 0$ $\displaystyle=2\omega^{1\mu}+d_{1}k\mu_{3},$ (300)
$\displaystyle 0$
$\displaystyle=2\omega^{1\mu}+c_{V}^{2}k\mu_{3}+a_{3}k\mu_{3},$ (301)
$\displaystyle 0$ $\displaystyle=2\omega^{1\mu}+b_{3}k\mu_{3}.$ (302)

Rearranging, we find:

$\displaystyle\frac{\omega^{1\mu}}{k}$
$\displaystyle=-\frac{d_{1}}{2}d_{1}\mu_{3},$ (303)
$\displaystyle\frac{\omega^{1\mu}}{k}$
$\displaystyle=-\frac{a_{3}+c_{V}^{2}}{2}\mu_{3},$ (304)
$\displaystyle\frac{\omega^{1\mu}}{k}$
$\displaystyle=-\frac{b_{3}}{2}\mu_{3}.$ (305)

We can now combine zeroth and first orders to obtain the modified speeds
$c_{S}^{\prime}$, $c_{V}^{\prime}$ and $c_{T}^{\prime}$. We also calculate and
insert the three coefficients $d_{1}/2$, $b_{3}/2$ and $(a_{3}+c_{V}^{2})/2$:

$\displaystyle c_{S}^{\prime}$
$\displaystyle=c_{S}-\frac{c_{123}}{(1-c_{+})(2+c_{+}+3c_{2})}\ $
$\displaystyle\mathrel{\phantom{=c_{S}-}}\times\left[c_{+}+2c_{2}-\frac{2(1+c_{2})^{2}}{c_{123}}+\frac{2}{c_{14}}\right]\mu_{3},$
(306) $\displaystyle c_{V}^{\prime}$
$\displaystyle=c_{V}-\frac{c_{+}(c_{+}+2c_{4})-2c_{4}}{2c_{14}(1-c_{+})}\mu_{3},$
(307) $\displaystyle c_{T}^{\prime}$
$\displaystyle=c_{T}-\frac{c_{+}}{1-c_{+}}\mu_{3}.$ (308)

At this point we take stock of how we can apply our results to observations of
lensed GW signals. We are primarily interested in how the additional modes
arising in Einstein-aether theory affect the speed of tensor waves. This means
we assume that scalar and vector modes are not emitted and that we cannot
measure $c_{S}^{\prime}$ and $c_{V}^{\prime}$.

We find that to first order in $\mu_{k}$ there is only one modified tensor
mode speed $c_{T}^{\prime}$. There is no birefringence – that would require
two different speeds $c_{T,1}^{\prime}$ and $c_{T,2}^{\prime}$. Additionally
$c_{T}^{\prime}$ depends only on $c_{+}$ and $\mu_{3}$:

$\displaystyle
c_{T}^{\prime}=\sqrt{\frac{1}{1-c_{+}}}-\frac{c_{+}}{1-c_{+}}\mu_{3}.$ (309)

This means that even if we were to observe one or more lensed signal (and
model their $\mu_{3}$) we can at best constrain $c_{+}$. Since $c_{+}$ can
already be constrained by observing unlensed signals this is not particularly
useful. In light of this we now turn to higher orders in $\mu_{k}$ , with the
hope that either birefringence arises or that the higher order corrections to
$c_{T}$ depend on parameters other than $c_{+}$.

##  IV Higher orders in µk

In this section we investigate higher orders in $\mu_{k}$. To that end we
calculate the full kinetic matrix in Sec. IV.1. In Sec. IV.2, we use this
result to determine the modified mode speeds up to second order in $\mu_{k}$.
After this example calculation, we obtain and discuss the modified mode speeds
up to fifth order in $\mu_{k}$ in Sec. IV.3.

###  IV.1 Full kinetic matrix

At higher orders in $\mu_{k}$ many intermediate results become quite complex.
We therefore perform most of the remaining calculations with Mathematica and
do not include all of these intermediate results. Some Mathematica-related
details of the calculations are discussed in Appendix D.

We first calculate the full kinetic matrix to all orders in $\mu_{k}$. To make
calculations simpler, we will assume propagation in the $z$ direction from the
beginning. Essentially, this will mean that partial derivatives in the $x$ and
$y$ direction vanish: $\partial_{1}=\partial_{2}=0$. As a consequence, the
operator $P_{ij}$ only has the entries
$P_{11}=P_{22}=\Delta=\partial_{3}^{2}$. The nonvanishing entries of $h_{ij}$
then take the following form:

$\displaystyle h_{01}$ $\displaystyle=\gamma_{1},$ $\displaystyle h_{11}$
$\displaystyle=\frac{F}{2}+\phi_{11},$ $\displaystyle h_{12}$
$\displaystyle=\phi_{12},$ (310) $\displaystyle h_{02}$
$\displaystyle=\gamma_{2},$ $\displaystyle h_{22}$
$\displaystyle=\frac{F}{2}-\phi_{11},$ $\displaystyle h_{33}$
$\displaystyle=\Delta\phi.$ (311)

We can now calculate the full kinetic matrix. As usual, the first step is the
constraint equation, now in its exact form. It separates into parts of zeroth
and first order in the perturbation:

$\displaystyle\mu_{0}$
$\displaystyle=\sqrt{1+\mu_{1}^{2}+\mu_{2}^{2}+\mu_{3}^{2}},$ (312)
$\displaystyle w_{0}$
$\displaystyle=-\frac{\mu_{0}}{2}h_{00}-\frac{\mu_{1}^{2}+\mu_{2}^{2}}{4\mu_{0}}F-(\gamma_{1}\mu_{1}+\gamma_{2}\mu_{2})$
$\displaystyle\mathrel{\phantom{=}}-\frac{1}{\mu_{0}}(v_{1}\mu_{1}+v_{2}\mu_{2})-\frac{\mu_{3}^{2}}{2\mu_{0}}\Delta\phi$
$\displaystyle\mathrel{\phantom{=}}-\frac{\mu_{1}^{2}-\mu_{2}^{2}}{2\mu_{0}}\phi_{11}-\frac{\mu_{1}\mu_{2}}{\mu_{0}}\phi_{12}.$
(313)

With this, we now calculate the linearized derivatives, $M$ tensor, aether
current, aether stress-energy tensor, $\overline{\text{\AE}^{[0}u^{i]}}$ and
effective tensor, as done previously. We do not include their exact forms for
brevity. Since we assume propagation in $z$ direction the SVT decomposition is
now simple:

$\displaystyle X_{10}^{T}$ $\displaystyle=X_{10},$ $\displaystyle X_{20}^{T}$
$\displaystyle=X_{20},$ (314) $\displaystyle X_{30}^{L}$
$\displaystyle=X_{30},$ $\displaystyle X_{01}^{T}$ $\displaystyle=X_{01},$
(315) $\displaystyle X_{02}^{T}$ $\displaystyle=X_{02},$ $\displaystyle
X_{03}^{L}$ $\displaystyle=X_{03},$ (316) $\displaystyle X_{11}^{TTL}$
$\displaystyle=\frac{X_{11}-X_{22}}{2},$ $\displaystyle X_{12}^{TTL}$
$\displaystyle=X_{12},$ (317) $\displaystyle X_{11}^{TTF}$
$\displaystyle=\frac{X_{11}+X_{22}}{2},$ $\displaystyle X_{13}^{LTL}$
$\displaystyle=X_{13},$ (318) $\displaystyle X_{23}^{LTL}$
$\displaystyle=X_{23},$ $\displaystyle X_{33}^{LTF}$ $\displaystyle=X_{33},$
(319)

where $TTL$ means the transverse traceless part, $TTF$ means the transverse
traceful part, $LTL$ means the longitudinal traceless part and $LTF$ means the
longitudinal traceful part.

We use $X_{00}$, $X_{i0}^{T}$ and $X_{i0}^{L}$ to constrain $h_{00}$, $v_{1}$,
$v_{2}$ and $\phi$. We then apply this to the five equations we need for the
kinetic matrix, already properly normalized:

$\displaystyle-\frac{4c_{123}c_{14}X_{11}^{TTF}}{(1-c_{+})(2+c_{+}+3c_{2})}$
$\displaystyle=0,$ (320)
$\displaystyle-\frac{c_{+}X_{01}^{T}}{(1-c_{+})c_{14}}$ $\displaystyle=0,$
(321) $\displaystyle-\frac{c_{+}X_{02}^{T}}{(1-c_{+})c_{14}}$
$\displaystyle=0,$ (322)
$\displaystyle-\frac{2X_{11}^{TTL}}{1-c_{+}\mu_{0}^{2}}$ $\displaystyle=0,$
(323) $\displaystyle-\frac{2X_{12}^{TTL}}{1-c_{+}\mu_{0}^{2}}$
$\displaystyle=0.$ (324)

The prefactors of the $F$ terms of the left-hand sides of Eqs. (320) through
(324) give us the first row of the kinetic matrix. The $\gamma_{1}$ terms give
us the second row and so on, until we have obtained the full kinetic matrix.
The result is complex enough that we opted not to calculate the exact
eigenvalues and mode speeds. Still, we can now easily obtain the entries of
the kinetic matrix and the mode speeds to higher orders in $\mu_{k}$.

###  IV.2 Second order in µk

As an example of how we obtain the higher order mode speed corrections we now
explicitly calculate the second order eigenvalues. The second order of the
equation defining the eigenvalues and eigenvectors is:

$\displaystyle
E^{0\mu}\vec{V}^{2\mu}+E^{1\mu}\vec{V}^{1\mu}+E^{2\mu}\vec{V}^{0\mu}$
$\displaystyle=\hat{K}^{0\mu}\vec{V}^{2\mu}+\hat{K}^{1\mu}\vec{V}^{1\mu}+\hat{K}^{2\mu}\vec{V}^{0\mu}.$
(325)

The only entries in $\hat{K}^{2\mu}$ that are relevant to us are the fourth
and fifth diagonal entries, both of which are equal:

$\displaystyle K_{44}^{2\mu}=K_{55}^{2\mu}$
$\displaystyle=-(2-c_{+})c_{+}c_{T}^{4}k^{2}\mu_{3}^{2}$
$\displaystyle\mathrel{\phantom{=}}-\frac{1}{2}\left(2c_{+}-2c_{4}-c_{+}c_{-}\right)c_{T}^{4}k^{2}\left(\mu_{1}^{2}+\mu_{2}^{2}\right).$
(326)

We use this and insert $E^{0\mu}$, $E^{1\mu}$, $\vec{V}^{0\mu}$ and
$\vec{V}^{1\mu}$ for the first modified tensor polarization.

$\displaystyle\left(\omega^{2}-c_{T}^{2}k^{2}\right)\left(\begin{array}[]{c}V_{1}\\\\[5.69054pt]
V_{2}\\\\[5.69054pt] V_{3}\\\\[5.69054pt] V_{4}\\\\[5.69054pt]
V_{5}\end{array}\right)^{2\mu}+\frac{2c_{+}}{1-c_{+}}\omega
k\mu_{3}\left(\begin{array}[]{c}0\\\\[5.69054pt]
\dfrac{c_{+}}{1-c_{+}}\mu_{1}\\\\[5.69054pt]
-\dfrac{c_{+}}{1-c_{+}}\mu_{2}\\\\[5.69054pt] 0\\\\[5.69054pt]
V_{5}^{1\mu}\end{array}\right)+E^{2\mu}\left(\begin{array}[]{c}0\\\ \\\ 0\\\
\\\ 0\\\ \\\ 1\\\ \\\ 0\end{array}\right)$ (346) $\displaystyle=$
$\displaystyle\left(\begin{array}[]{ccccccccc}K_{SS}&&0&&0&&0&&0\\\\[5.69054pt]
0&&K_{VV}&&0&&0&&0\\\\[5.69054pt] 0&&0&&K_{VV}&&0&&0\\\\[5.69054pt]
0&&0&&0&&K_{TT}&&0\\\\[5.69054pt]
0&&0&&0&&0&&K_{TT}\end{array}\right)^{0\mu}\left(\begin{array}[]{c}V_{1}\\\\[5.69054pt]
V_{2}\\\\[5.69054pt] V_{3}\\\\[5.69054pt] V_{4}\\\\[5.69054pt]
V_{5}\end{array}\right)^{2\mu}$ (357)
$\displaystyle+\left(\begin{array}[]{ccccccccc}K_{SS}^{1\mu}&&K_{SV}\mu_{1}&&K_{SV}\mu_{2}&&0&&0\\\\[5.69054pt]
K_{VS}\mu_{1}&&K_{VV}^{1\mu}&&0&&K_{VT}\mu_{1}&&K_{VT}\mu_{2}\\\\[5.69054pt]
K_{VS}\mu_{2}&&0&&K_{VV}^{1\mu}&&-K_{VT}\mu_{2}&&K_{VT}\mu_{1}\\\\[5.69054pt]
0&&K_{TV}\mu_{1}&&-K_{TV}\mu_{2}&&K_{TT}^{1\mu}&&0\\\\[5.69054pt]
0&&K_{TV}\mu_{2}&&K_{TV}\mu_{1}&&0&&K_{TT}^{1\mu}\end{array}\right)\left(\begin{array}[]{c}0\\\\[5.69054pt]
\dfrac{c_{+}}{1-c_{+}}\mu_{1}\\\\[5.69054pt]
-\dfrac{c_{+}}{1-c_{+}}\mu_{2}\\\\[5.69054pt] 0\\\\[5.69054pt]
V_{5}^{1\mu}\end{array}\right)+\left(\begin{array}[]{ccccc}\ddots&&\vdots&&\vdots\\\\[5.69054pt]
\dots&&K_{44}&&\vdots\\\\[5.69054pt]
\dots&&\dots&&K_{55}\end{array}\right)^{2\mu}\left(\begin{array}[]{c}0\\\\[5.69054pt]
0\\\\[5.69054pt] 0\\\\[5.69054pt] 1\\\\[5.69054pt] 0\end{array}\right).$ (376)

Since we are only interested in the eigenvalues and speeds, we focus on the
fourth row:

$\displaystyle\left(\omega^{2}-c_{T}^{2}k^{2}\right)V_{4}^{2\mu}+E^{2\mu}$
$\displaystyle=$ $\displaystyle
K_{TT}^{0\mu}V_{4}^{2\mu}+\frac{K_{TV}c_{+}}{1-c_{+}}\left(\mu_{1}^{2}+\mu_{2}^{2}\right)+K_{44}^{2\mu}.$
(377)

The first terms on both sides are identical:

$\displaystyle E^{2\mu}$
$\displaystyle=\frac{K_{TV}c_{+}}{1-c_{+}}\left(\mu_{1}^{2}+\mu_{2}^{2}\right)+K_{44}^{2\mu}.$
(378)

We insert $K_{TV}$ and $K_{44}^{2\mu}$ that, remember, are given by Eqs. (109)
and (326), respectively. We find:

$\displaystyle E^{2\mu}$
$\displaystyle=\left(-b_{1}k^{2}-b_{2}\omega^{2}\right)\frac{c_{+}}{1-c_{+}}\left(\mu_{1}^{2}+\mu_{2}^{2}\right)$
$\displaystyle\mathrel{\phantom{=}}-(2-c_{+})c_{+}c_{T}^{4}k^{2}\mu_{3}^{2}$
$\displaystyle\mathrel{\phantom{=}}-\frac{1}{2}\left(2c_{+}-2c_{4}-c_{+}c_{-}\right)c_{T}^{4}k^{2}\left(\mu_{1}^{2}+\mu_{2}^{2}\right).$
(379)

Finally, we insert $b_{1}$ and $b_{2}$:

$\displaystyle E^{2\mu}$
$\displaystyle=\left[-\frac{2c_{3}+c_{+}c_{-}}{2c_{+}(1-c_{+})}k^{2}+\frac{c_{3}-c_{4}}{c_{+}}\omega^{2}\right]\frac{c_{+}}{1-c_{+}}\left(\mu_{1}^{2}+\mu_{2}^{2}\right)$
$\displaystyle\mathrel{\phantom{=}}-(2-c_{+})c_{+}c_{T}^{4}k^{2}\mu_{3}^{2}$
$\displaystyle\mathrel{\phantom{=}}-\frac{1}{2}\left(2c_{+}-2c_{4}-c_{+}c_{-}\right)c_{T}^{4}k^{2}\left(\mu_{1}^{2}+\mu_{2}^{2}\right).$
(380)

To obtain the second-order speed correction we now expand $\omega$ up to
$\omega^{2\mu}$. This again leads to a shuffling of the orders in $\mu_{k}$:

$\displaystyle E^{2\mu}$
$\displaystyle=\left(\omega^{1\mu}\right)^{2}+2\omega^{0\mu}\omega^{2\mu}+\frac{2c_{+}}{1-c_{+}}\omega^{1\mu}k\mu_{3}$
$\displaystyle\mathrel{\phantom{=}}+\left[\frac{c_{3}-c_{4}}{1-c_{+}}\left(\omega^{0\mu}\right)^{2}-\frac{2c_{3}+c_{+}c_{-}}{2(1-c_{+})^{2}}k^{2}\right]\left(\mu_{1}^{2}+\mu_{2}^{2}\right)$
$\displaystyle\mathrel{\phantom{=}}-(2-c_{+})c_{+}c_{T}^{4}k^{2}\mu_{3}^{2}$
$\displaystyle\mathrel{\phantom{=}}-\frac{1}{2}\left(2c_{+}-2c_{4}-c_{+}c_{-}\right)c_{T}^{4}k^{2}\left(\mu_{1}^{2}+\mu_{2}^{2}\right).$
(381)

We set this expression to zero and insert the values for $\omega^{0\mu}$ and
$\omega^{1\mu}$ obtained from the lower orders:

$\displaystyle 0$
$\displaystyle=\frac{c_{+}^{2}}{(1-c_{+})^{2}}k^{2}\mu_{3}^{2}+2c_{T}\omega^{2\mu}k-\frac{2c_{+}}{1-c_{+}}\frac{c_{+}}{1-c_{+}}k^{2}\mu_{3}^{2}$
$\displaystyle\mathrel{\phantom{=}}+\left[\frac{c_{3}-c_{4}}{1-c_{+}}c_{T}^{2}-\frac{2c_{3}+c_{+}c_{-}}{2(1-c_{+})^{2}}\right]\left(\mu_{1}^{2}+\mu_{2}^{2}\right)k^{2}$
$\displaystyle\mathrel{\phantom{=}}-(2-c_{+})c_{+}c_{T}^{4}k^{2}\mu_{3}^{2}$
$\displaystyle\mathrel{\phantom{=}}-\frac{1}{2}\left(2c_{+}-2c_{4}-c_{+}c_{-}\right)c_{T}^{4}k^{2}\left(\mu_{1}^{2}+\mu_{2}^{2}\right).$
(382)

We now combine and cancel terms using $1-c_{+}=c_{T}^{-2}$:

$\displaystyle 0$
$\displaystyle=2c_{T}\omega^{2\mu}k-2c_{+}c_{T}^{4}k^{2}\mu_{3}^{2}-c_{+}c_{T}^{4}k^{2}\left(\mu_{1}^{2}+\mu_{2}^{2}\right),$
(383)

and then solve for $\omega^{2\mu}/k$:

$\displaystyle\frac{\omega^{2\mu}}{k}$
$\displaystyle=c_{+}c_{T}^{3}\mu_{3}^{2}+\frac{1}{2}c_{+}c_{T}^{3}\left(\mu_{1}^{2}+\mu_{2}^{2}\right).$
(384)

Repeating the entire calculation for the second modified tensor mode we find
that we once again only have one modified mode speed:

$\displaystyle c_{T}^{\prime}$ $\displaystyle=c_{T}-c_{+}c_{T}^{2}\mu_{3}$
$\displaystyle\mathrel{\phantom{=}}+c_{+}c_{T}^{3}\mu_{3}^{2}+\frac{1}{2}c_{+}c_{T}^{3}\left(\mu_{1}^{2}+\mu_{2}^{2}\right).$
(385)

Once again there is no birefringence and the modified speed depends only on
$c_{+}$ and the $\mu_{i}$. We can now keep repeating the eigenvalue
calculation in the same manner for increasingly higher orders in $\mu_{k}$.

###  IV.3 Fifth Order in µk

For this paper, we calculated the modified tensor speed up to
$\mathcal{O}\left(\mu_{k}^{5}\right)$:

$\displaystyle c_{T}^{\prime}$
$\displaystyle=c_{T}\bigg{[}1-c_{+}c_{T}\mu_{L}+\frac{1}{2}c_{+}c_{T}^{2}\left(2\mu_{L}^{2}+\mu_{T}^{2}\right)$
$\displaystyle\mathrel{\phantom{=c_{T}\bigg{[}}}-\frac{1}{2}c_{+}c_{T}^{3}\left(1-c_{+}\right)\mu_{L}\left(\mu_{L}^{2}+\mu_{T}^{2}\right)$
$\displaystyle\mathrel{\phantom{=c_{T}\bigg{[}}}+\frac{1}{8}c_{+}^{2}c_{T}^{4}\left(8\mu_{L}^{4}+12\mu_{L}^{2}\mu_{T}^{2}+3\mu_{T}^{4}\right)$
$\displaystyle\mathrel{\phantom{=c_{T}\bigg{[}}}+\frac{1}{8}c_{+}c_{T}^{5}\left(1-6c_{+}-3c_{+}^{2}\right)\mu_{L}\left(\mu_{L}^{2}+\mu_{T}^{2}\right)^{2}\bigg{]}.$
(386)

Here we switched notation to longitudinal ($\mu_{L}=\mu_{3}$) and transverse
($\mu_{T}=\sqrt{\mu_{1}^{2}+\mu_{2}^{2}}$) aether components. We again find
the absence of birefringence and a modified speed depending only on $c_{+}$.
We suspect that this trend will continue at higher orders and that,
unfortunately, we can only constrain $c_{+}$ with our methods.

To make the results more intuitive, Fig. 3 shows the effect of the aether on
the tensor mode speed by plotting $c_{T}^{\prime}-c_{T}$ (up to fifth order in
$\mu_{k}$) as a function of $\mu_{L}$, $\mu_{T}$ and $c_{+}$. The magnitudes
$\left\lvert\mu_{T}\right\rvert$ and $\left\lvert\mu_{L}\right\rvert$ should
be smaller than $1$, so we (somewhat arbitrarily) let them run from $-1/2$ to
$1/2$. For $c_{+}$, we recall that observations constrain $\Delta
c_{T}=c_{T}-1$ to between $-5\times 10^{-16}$ [22] and $7\times 10^{-16}$ [6].
We can express $c_{+}$ as a function of $\Delta c_{T}$:

$\displaystyle c_{+}=1-\frac{1}{(1+\Delta c_{T})^{2}}\approx 2\Delta c_{T}.$
(387)

Therefore, the two extreme values of $c_{+}$ permitted by observations are
$-10^{-15}$ and $1.4\times 10^{-15}$. We use these values for Fig. 3.

Figure 3: Deviation of the fifth order modified tensor mode speed
$c_{T}^{\prime}$ from the background value $c_{T}$ as a function of the
longitudinal aether component $\mu_{L}$ and transverse aether component
$\mu_{T}$. Top: Upper limit of $c_{+}$ permitted by observations. Bottom:
Lower limit of $c_{+}$. A faster speed compared to the background value is
indicated by reds, a slower speed is indicated by blues. Higher values of
$\mu_{L}$ slow tensor waves down for the upper limit of $c_{+}$ and speed them
up for the lower limit. Higher magnitudes of $\mu_{T}$ increase speed for the
upper limit and decrease it for the lower limit.

We find that for positive $c_{+}$ the longitudinal component of the aether,
$\mu_{L}$, slows down tensor waves. Larger values of $\mu_{L}$ correspond to a
stronger slowdown, while negative values of $\mu_{L}$ correspond to negative
slowdown, that is, to a speedup. For negative values of $c_{+}$ the
longitudinal component of the aether instead speeds tensor waves up, acting as
a sort of tailwind. The transverse component of the aether, $\mu_{T}$, has the
same effect independent of its sign: it increases $c_{T}^{\prime}$ for
positive values of $c_{+}$ and decreases it for negative values of $c_{+}$.

The contrast between Einstein-aether theory and Horndeski’s theory, as
investigated in Ref. [52], is notable:

  1. 1.

In Einstein-aether theory there is only one modified speed $c_{T}^{\prime}$
and no birefringence. However, different lensed signals have different speeds.
We refer to this as lens-dependent GW speed. Observing lensed signals allows
us to constrain the parameter combinations that appear in the singular
modified speed $c_{T}^{\prime}$.

  2. 2.

In Horndeski’s theory there are two different modified mode speeds
$c_{T,1}^{\prime}$ and $c_{T,2}^{\prime}$. In this case we have lens-induced
birefringence. We can measure the delay between the two polarizations of a
signal and constrain parameter combinations that appear in
$c_{T,1}^{\prime}-c_{T,2}^{\prime}$. Note that both speeds depend on the
properties of the lens and differently lensed signals may exhibit different
time delays.

Testing these situations requires comparing the arrival time of different
signals. This happens because the total travel time is subject to very large
uncertainties, e. g., from gravitational potentials of intervening galaxies
[78]. Testing lens-dependent GW speed therefore requires an electromagnetic or
neutrino counterpart [32, 79]. This reduces the number of testable events and
their redshift, limiting the chances of a close lens-source alignment.
However, in the case of lens-induced birefringence, each GW polarization acts
as an independent messenger. This allows all GW events to be considered,
regardless of non-gravitational counterparts [53]. Figure 4 showcases the
differences between the two lens-induced effects.

Figure 4: Differences between lens-induced effects in different theories. The
propagation of a wave packet from left to right is shown. Time is measured in
arbitrary units. Top and middle: lens-dependent GW speed as observed in
Einstein-aether theory. There is no difference in speed between the two tensor
modes. There is no birefringence, any one signal does not get scrambled.
However, signals lensed by different lenses propagate with different speeds.
Bottom: Lens-induced birefringence as observed in Horndeski’s theory. The
modified cross and plus mode propagate with different speeds. Birefringence
scrambles the signal.

Regarding observations, the main signature is the dependence of
$c_{T}^{\prime}$ on $\mu_{L}$, $\mu_{T}$ and $c_{+}$. To test it one needs to
complete the following steps:

  1. 1.

measure $c_{T}^{\prime}$ for a lensed signal

  2. 2.

model $\mu_{L}$ and $\mu_{T}$ for the signal

  3. 3.

calculate $c_{+}$ by using Eq. (386)

  4. 4.

compare with other constraints imposed on $c_{+}$

With this result we now move to conclude this paper and to give an outlook to
the future.

##  V Summary and outlook

In this paper, we studied GW lensing in Einstein-aether theory. We discussed
Einstein-aether theory with its unit timelike vector field, the aether. We
briefly summarized previous results on the propagation of Einstein-aether
excitations over Minkowski spacetime with a timelike aether field. Solving the
linearized Einstein-aether field equations shows the existence of three new
propagating modes. There is the scalar mode $F$ that contracts and expands an
object both longitudinally and transversally and two vector modes $\gamma_{1}$
and $\gamma_{2}$. The propagation speeds of all fluctuations, including the
standard GW polarizations, are dependent on the theory parameters.

The novelty of our analysis is to combine Einstein-aether theory and GW
lensing, that is, describe the evolution of fluctuations in a non-homogeneous
background (e. g. induced by a gravitational lens). We found that
approximations are necessary to allow analytical calculations. We made three
such approximations: small GW amplitudes, small background aether
inhomogeneities and high GW frequency. This allowed us to bring the Einstein-
aether field equations into the form of a matrix equation containing the
kinetic matrix $\hat{K}$. The propagation eigenstate framework [52] (see also
Ref. [77]) allowed us to use eigenvalues and eigenvectors of $\hat{K}$ to
study the mixing of modes and the resulting modified mode speeds.

Implementing the small background aether inhomogeneities approximation proved
challenging. By expanding to only the first order of the spatial components of
the background aether $\mu_{k}$, calculations were simplified dramatically.

Background aether inhomogeneities introduce a linear correction to the GW
speed (Eq. (308)). Two features of this correction stand out

  1. 1.

It is common to both GW polarizations. Therefore, Einstein-aether theory
predicts no birefringence, and deviations on the GW speed can only be tested
by multimessenger observations, using the arrival time relative to another
signal (e. g. an electromagnetic counterpart).

  2. 2.

It only depends on $c_{+}$, the same parameter that controls the speed of GWs
on a homogeneous background. Therefore, existing constraints limit the
magnitude of lensing-induced corrections and their potential to test
additional sectors of Einstein-aether theory.

The situation is substantially different from scalar-tensor gravity, where the
existence of birefringence [52] allows stringent tests using catalogues of
black hole mergers [53]. We verified that these results hold up to fifth order
in the background aether inhomogeneities
$\mathcal{O}\left(\mu_{k}^{5}\right)$, plausibly discarding the possibility
that they are due to our approximations. This difference may be due to the
lack of coupling between tensor and scalar perturbations mediated by
nonhomogeneous aether configurations, at leading order of the high-frequency
expansion. In Horndeski, this coupling is mediated by non-linear derivative
interactions $\mathcal{L}_{H}\supset
R(\partial\phi)^{2},(\Box\phi)^{2},\dots$, which are absent in Einstein-aether
theory.

There are many aspects in which GW lensing in Einstein-aether theory should be
further studied. An important next step is to consider the next orders in the
frequency expansion: corrections $\sim f^{0}$ will describe the amplitude of
the propagating degrees of freedom. At this order, it is possible that
gravitational lenses introduce observable effects, in the form of amplitude
oscillations [80], but between tensor, vector and scalar polarizations. These
effects can be tested by the presence of multiple polarizations at the
detector [19]. Oscillations between states with different helicities are
forbidden by Lorentz invariance on a symmetric background and can only be
mediated by inhomogeneities in the aether field. An additional order in the
expansion, $\sim f^{-2}$, describes dispersive effects in GW propagation. GW
dispersion affects the waveform, and can be thus tested on all signals,
without an electromagnetic counterpart [18]. Dispersion is a distinct
signature of massive gravity [81, 82]. In GR by strong gravitational fields
and even be birefringent [83, 84, 85]. However, GW dispersion and additional
polarizations would be a clear signature of lens-induced and beyond-GR
effects.

Further studies should overcome the perturbative expansion on the spatial
components of the aether field $\mu_{k}\ll 1$. As a first attempt, this
sufficed to characterize the lensing effect on the GW speed. An analysis
without assuming small inhomogeneities could be done in general or relying on
spherically symmetric solutions of the theory [9]. Moreover, a nonperturbative
analysis will shed light on the many cancellations leading to the properties
of modified GW speed: a lack of birefringence in contrast to scalar-tensor
theories [52] and the fact that lens-induced corrections only depend on
$c_{+}$, the one-parameter combination responsible for the GW speed on
homogeneous backgrounds. Identifying the structural property of the theory
responsible for this feature will likely shed light on how birefringence and
anomalous propagation may emerge in other classes of deviations from
Einstein’s theory.

In this paper we used the propagation eigenstate framework to study Einstein-
aether theory. However, we can in principle apply it to many other theories
beyond GR, with Lorentz-violating scenarios being especially compelling.
Hořava–Lifshitz gravity [86] is a theory of quantum gravity built with
concepts from condensed matter physics. Just like with Einstein-aether theory,
it breaks local Lorentz invariance. However, while Einstein-aether theory is
equally anisotropic at all energy ranges, Hořava–Lifshitz gravity is strongly
anisotropic in the quantum regime: at high energies, high frequencies, short
wavelengths and short length scales. At low energies and long length scales,
in the regime of classical gravity, the anisotropy is much weaker [86].
Blas–Pujolàs–Sibiryakov–Hořava (BPSH) theory [87] is another interesting
extension. While Hořava–Lifshitz gravity has a static preferred direction,
BPSH theory has a dynamic preferred direction described by a vector field,
analogous to Einstein-aether theory. Compared to Einstein-aether theory, there
is an additional demand on the vector field in BPSH theory: it must be
hypersurface orthogonal. This means that it must be possible to split
spacetime into infinitely many spacelike, three-dimensional hypersurfaces so
that the vector field is always orthogonal to these hypersurfaces. In fact,
Einstein-aether theory is the hypersurface orthogonal low energy limit of BPSH
theory [11].

Our work is a first step towards understanding GW propagation in Einstein-
aether theories with inhomogeneous field configurations, that is, GW lensing.
The leading order results reflect a remarkable simplicity, as lens-induced
corrections are common to both $+,\times$ polarizations and depend only on the
parameters that modify the GW speed on a spatially-homogeneous background.
This feature, likely anchored in the inner structure of the theory, calls out
for a simpler explanation. It is nevertheless plausible that additional
effects in the amplitude and phase of GWs will enable novel tests, not only of
Einstein-aether, but of Lorentz invariance in general. Given the vast numbers
of signals available to next-generation ground detectors [88, 89, 90], space-
interferometers pushing towards lower frequencies [91, 92, 93] and growing
evidence for nanohertz GWs from pulsar-timing arrays [94, 95, 96, 97], further
developments in GW lensing will vastly enhance our capacity to test Einstein’s
theory and its underlying assumptions.

## Acknowledgements

We thank Guilherme Brando, Sumanta Chakraborty, Lyla Choi, Alessandro Longo,
Stefano Savastano, Hiroki Takeda, Giovanni Tambalo and Héctor Villarrubia-Rojo
for discussions. HOS acknowledges funding from the Deutsche
Forschungsgemeinschaft (DFG) - Project No.: 386119226.

##  Appendix A Effective Tensor Subtleties

Our definition of the effective tensor deviates from the one used by Foster
[23] in one key point: He includes the term $\overline{\text{\AE}^{[i}u^{j]}}$
in the definition of $X^{ij}$. Physically, this ensures the conservation of
source terms. We will refer to his definition of the effective tensor as
$\tilde{X}_{\alpha\beta}$. Assuming a sufficiently slowly varying background
aether $\mu^{\alpha}$, the two sets of definitions are related as follows:

$\displaystyle\tilde{X}^{00}$ $\displaystyle=X^{00},$ (388)
$\displaystyle\tilde{X}^{0i}$ $\displaystyle=X^{0i},$ (389)
$\displaystyle\tilde{X}^{i0}$ $\displaystyle=X^{i0},$ (390)
$\displaystyle\tilde{X}^{ij}$
$\displaystyle=X^{ij}+\frac{\mu^{i}}{2\mu^{0}}\cdot
X^{0j}-\frac{\mu^{i}}{2\mu^{0}}\cdot X^{j0}$
$\displaystyle\mathrel{\phantom{=}}-\frac{\mu^{j}}{2\mu^{0}}\cdot
X^{0i}+\frac{\mu^{j}}{2\mu^{0}}\cdot X^{i0},$ (391) $\displaystyle X^{ij}$
$\displaystyle=\tilde{X}^{ij}-\frac{\mu^{i}}{2\mu^{0}}\cdot\tilde{X}^{0j}+\frac{\mu^{i}}{2\mu^{0}}\cdot\tilde{X}^{j0}$
$\displaystyle\mathrel{\phantom{=}}+\frac{\mu^{j}}{2\mu^{0}}\cdot\tilde{X}^{0i}-\frac{\mu^{j}}{2\mu^{0}}\cdot\tilde{X}^{i0}.$
(392)

As $X^{\alpha\beta}$ and $\tilde{X}^{\alpha\beta}$ can be written as linear
combinations of each other’s components, so can the field equations
$X^{\alpha\beta}=0$ and $\tilde{X}^{\alpha\beta}=0$. Therefore, both sets of
field equations are mathematically equivalent.

Foster’s definition works well enough in the exact case, but becomes
problematic when combined with series expansions in $\mu^{k}$ as used in Sec.
III. $\tilde{X}^{\alpha\beta}$ is subject to the conservation law
$\tilde{X}^{\alpha\beta}_{\phantom{\alpha\beta},\beta}=0$, which is crucial
for this definition to yield the correct number of wave equations. Similarly
important is the condition

$\displaystyle\tilde{X}^{[ij]}=\frac{\mu^{j}}{\mu^{0}}\tilde{X}^{[i0]}+\frac{\mu^{i}}{\mu^{0}}\tilde{X}^{[0j]},$
(393)

which stems from $\text{\AE}^{[\alpha}u^{\beta]}$ only having three
independent components. Now consider the antisymmetric part of the $i$
component of the conservation law:

$\displaystyle\partial_{0}\tilde{X}^{[i0]}=\partial_{j}\tilde{X}^{[ij]}.$
(394)

If we apply $\partial_{j}$ to Eq. (393), we can combine it with Eq. (394) into
a condition that must be upheld if we want to obtain the correct number of
wave equations:

$\displaystyle\partial_{0}\tilde{X}^{[i0]}=\frac{\mu^{j}}{\mu^{0}}\partial_{j}\tilde{X}^{[i0]}+\frac{\mu^{i}}{\mu^{0}}\partial_{j}\tilde{X}^{[0j]}.$
(395)

If we violate Eq. (395) then we lose one of the two important redundancies and
will, in the end, obtain superfluous false wave equations. This means that if
we consider a series expansion in $\mu^{k}$, Eq. (395) must hold at all orders
and all terms on both sides must always be expanded to the same order in
$\mu^{k}$. However, this is clearly impossible: If
$\partial_{0}\tilde{X}^{[i0]}$ requires $\tilde{X}^{[i0]}$ to be expanded to
the $n$th order in $\mu^{k}$, then
$\mu^{j}\partial_{j}\tilde{X}^{[i0]}/\mu^{0}$, which contains the additional
factor $\mu^{j}$, requires $\tilde{X}^{[i0]}$ to be expanded to $n-1$st order
in $\mu^{k}$. As we cannot expand $\tilde{X}^{[i0]}$ to different orders
simultaneously, we find that Foster’s definition is unsuitable for us and opt
to trade the conservation law for mathematical simplicity.

##  Appendix B Gauge considerations

To prove that $\gamma=\phi_{i}=v=0$ is a possible gauge, we consider how a
gauge transformation affects perturbations on a Minkowski background:

$\displaystyle h^{\prime}_{\mu\nu}$
$\displaystyle=h_{\mu\nu}+\varepsilon_{\mu,\nu}+\varepsilon_{\nu,\mu},$ (396)
$\displaystyle w^{\prime}_{\mu}$
$\displaystyle=w_{\mu}-\mu^{\alpha}\partial_{\alpha}\varepsilon_{\mu}+\varepsilon^{\alpha}\partial_{\alpha}\mu_{\mu}.$
(397)

Note that some authors switch primed and unprimed quantities, or, conversely,
switch the sign for the $\varepsilon$ terms. Performing another SVT
decomposition we split $\varepsilon_{\mu}$ into temporal part $\xi_{0}$,
spatial transverse part $\xi_{i}$ and spatial longitudinal part $\xi_{,i}$. We
then find the following transformation rules for the metric perturbation:

$\displaystyle h_{00}$ $\displaystyle\rightarrow h_{00}+2\dot{\xi}_{0},$ (398)
$\displaystyle\gamma_{i}$ $\displaystyle\rightarrow\gamma_{i}+\dot{\xi}_{i},$
(399) $\displaystyle\gamma$
$\displaystyle\rightarrow\gamma+\xi_{0}+\dot{\xi},$ (400)
$\displaystyle\phi_{ij}$ $\displaystyle\rightarrow\phi_{ij},$ (401)
$\displaystyle f$ $\displaystyle\rightarrow f,$ (402) $\displaystyle\phi_{i}$
$\displaystyle\rightarrow\phi_{i}+\xi_{i},$ (403) $\displaystyle\phi$
$\displaystyle\rightarrow\phi+2\xi.$ (404)

For the aether perturbation we will assume a more general background
$\mu_{\alpha}$ instead of $\delta^{0}_{\alpha}$. We then also need to
decompose the spatial background aether $\mu_{i}$ into transverse part
$\nu_{i}$ and longitudinal part $\nu_{,i}$:

$\displaystyle w_{0}$ $\displaystyle\rightarrow
w_{0}-\mu^{\alpha}\partial_{\alpha}\xi_{0}-\xi_{0}\dot{\mu}_{0}+\xi_{k}\mu_{0,k}+\xi_{,k}\mu_{0,k},$
(405) $\displaystyle v_{i}$ $\displaystyle\rightarrow
v_{i}-\mu^{\alpha}\partial_{\alpha}\xi_{i}-\xi_{0}\dot{\nu}_{i}+\xi_{k}\nu_{i,k}+\xi_{,k}\nu_{i,k},$
(406) $\displaystyle v$ $\displaystyle\rightarrow
v-\mu^{\alpha}\partial_{\alpha}\xi-\xi_{0}\dot{\nu}+\xi_{k}\nu_{,k}+\xi_{,k}\nu_{,k}.$
(407)

Whether or not a gauge is possible to reach depends on whether or not
$\xi_{0}$, $\xi_{i}$ and $\xi$ can be found for it. For the gauge
$\gamma=\phi_{i}=v=0$, $\xi_{0}$ and $\xi_{i}$ can be specified while $\xi$ is
the solution of a partial differential equation:

$\displaystyle\xi_{0}$ $\displaystyle=-\gamma-\dot{\xi},$ (408)
$\displaystyle\xi_{i}$ $\displaystyle=-\phi_{i},$ (409)
$\displaystyle\dot{\xi}$ $\displaystyle=A_{k}\xi_{,k}+B,$ (410) $\displaystyle
A_{i}$ $\displaystyle=\frac{\nu_{i}}{\mu_{0}+\dot{\nu}},$ (411) $\displaystyle
B$
$\displaystyle=-\frac{v-\phi_{k}\nu_{,k}+\gamma\dot{\nu}}{\mu_{0}+\dot{\nu}}.$
(412)

We are unable to solve this equation, but to enable our gauge transformation
we only need to prove the existence of a solution. This existence is
guaranteed by the Cauchy–Kowalevski theorem. We will adopt it from a textbook
by Folland [98] to our specific situation:

  * •

one single partial differential equation

  * •

$A_{i}$, $B$ and $\xi$ depend on spacetime coordinates $x^{\mu}$ only

  * •

$A_{i}$ and $B$ are defined on all of spacetime (due to $\mu_{0}\approx 1$ and
$\nu\approx 0$ the denominators cause no issues)

The theorem then takes the following form:

“Suppose that $A_{k}$ and $B$ are analytic $\mathbb{R}$-valued functions
defined on all of spacetime. Then the Cauchy problem
$\dot{\xi}=A_{k}\xi_{,k}+B$ with initial conditions $\xi(t=0,x^{i})=0$ has a
unique analytic solution on all of spacetime.”

This theorem ensures that our choice of gauge is possible.

##  Appendix C Derivation of the effects of Einstein-aether waves

We recall our starting points: monochromatic plane waves in $z$ direction and
the metric and aether perturbations $h_{\mu\nu}$ and $w^{\mu}$.

$\displaystyle\psi=\hat{\psi}\exp\left[i\omega\left(\frac{z}{c_{p}}-t\right)-\frac{i\pi}{2}\right],$
(413) $\displaystyle w^{\mu}=\left(\begin{array}[]{c}\dfrac{h_{00}}{2}\\\ \\\
v_{i}\end{array}\right),\;h_{\mu\nu}=\left(\begin{array}[]{ccc}h_{00}&&\gamma_{i}\\\
\\\ \gamma_{i}&&\phi_{ij}+\dfrac{1}{2}P_{ij}f+\phi_{,ij}\end{array}\right).$
(420)

When acting on plane waves, $P_{ij}$ simplifies to:

$\displaystyle P_{11}\psi$ $\displaystyle=P_{22}\psi$
$\displaystyle=\Delta\psi,$ (421) $\displaystyle P_{i\neq j}\psi$
$\displaystyle=P_{33}\psi$ $\displaystyle=0.$ (422)

Similarly, $\partial_{i}\partial_{j}$ becomes $\delta_{i3}\delta_{j3}\Delta$
if we assume our solution. We combine this with the constraint equation to
obtain metric and aether perturbations:

$\displaystyle w^{\mu}$
$\displaystyle=\left(\begin{array}[]{c}\dfrac{F}{2c_{14}}\\\ \\\
-\dfrac{\gamma_{1}}{c_{+}}\\\ \\\ -\dfrac{\gamma_{2}}{c_{+}}\\\ \\\
0\end{array}\right),\;$ (430) $\displaystyle h_{\mu\nu}$
$\displaystyle=\left(\begin{array}[]{ccccccc}\dfrac{F}{c_{14}}&&\gamma_{1}&&\gamma_{2}&&0\\\
\\\ \gamma_{1}&&\phi_{11}+\dfrac{F}{2}&&\phi_{12}&&0\\\ \\\
\gamma_{2}&&\phi_{12}&&-\phi_{11}+\dfrac{F}{2}&&0\\\ \\\
0&&0&&0&&-\dfrac{1+c_{2}}{c_{123}}F\end{array}\right).$ (438)

While our result for the metric perturbation is correct, it is difficult to
describe the effects of GWs in this gauge since the time coordinate is
affected. The synchronous gauge, defined by $h_{\alpha 0}=0$, is better suited
for this purpose. Therefore, we make use of a result by Hiroki Takeda [60]. In
theories beyond GR with all modes speeds equal to $1$ the six gauge invariant
polarizations can be expressed in terms of the components of $h_{\mu\nu}$:

$\displaystyle h_{b}$ $\displaystyle=\frac{1}{2}\left(h_{11}+h_{22}\right),$
$\displaystyle h_{l}$ $\displaystyle=h_{00}+2h_{03}+h_{33},$ (439)
$\displaystyle h_{x}$ $\displaystyle=h_{01}+h_{13},$ $\displaystyle h_{y}$
$\displaystyle=h_{02}+h_{23},$ (440) $\displaystyle h_{+}$
$\displaystyle=\frac{1}{2}\left(h_{11}-h_{22}\right),$ $\displaystyle
h_{\times}$ $\displaystyle=h_{12}.$ (441)

We now need to correct for the fact that Einstein-aether theory has mode
speeds different from $1$. To see what new forms these polarizations take we
use the (gauge invariant) Ricci tensor. By calculating its components and
forming linear combinations, we can find the six polarizations for speeds
different from $1$:

$\displaystyle h_{b}$ $\displaystyle=\frac{1}{2}\left(h_{11}+h_{22}\right),$
$\displaystyle h_{l}$
$\displaystyle=\frac{1}{c_{X}^{2}}h_{00}+\frac{2}{c_{X}}h_{03}+h_{33},$ (442)
$\displaystyle h_{x}$ $\displaystyle=\frac{1}{c_{X}}h_{01}+h_{13},$
$\displaystyle h_{y}$ $\displaystyle=\frac{1}{c_{X}}h_{02}+h_{23},$ (443)
$\displaystyle h_{+}$ $\displaystyle=\frac{1}{2}\left(h_{11}-h_{22}\right),$
$\displaystyle h_{\times}$ $\displaystyle=h_{12}.$ (444)

Here we treat the inverse mode speed $1/c_{X}$ as an operator, defined as
follows:

$\displaystyle\frac{1}{c_{X}}\psi$
$\displaystyle=\psi\sqrt{\frac{\Delta\psi}{\ddot{\psi}}}.$ (445)

In other words, the concrete value that $c_{X}$ takes will depend on the entry
of $h_{\mu\nu}$ that it is acting on.

We now want to use the gauge invariant polarizations to express our metric
perturbations in the synchronous gauge. The polarizations $h_{b}$, $h_{+}$ and
$h_{\times}$ naturally stay the same before and after transformation since we
leave $h_{11}$, $h_{12}$ and $h_{22}$ untouched. For the remaining three
polarizations the transformed components of $h_{\mu\nu}$, marked with a prime,
must be equal to their counterparts from before the transformation:

$\displaystyle\frac{1}{c_{X}}h^{\prime}_{01}+h^{\prime}_{13}$
$\displaystyle=\frac{1}{c_{X}}h_{01}+h_{13},$ (446)
$\displaystyle\frac{1}{c_{X}}h^{\prime}_{02}+h^{\prime}_{23}$
$\displaystyle=\frac{1}{c_{X}}h_{02}+h_{23},$ (447)
$\displaystyle\frac{1}{c_{X}^{2}}h^{\prime}_{00}+\frac{2}{c_{X}}h^{\prime}_{03}+h^{\prime}_{33}$
$\displaystyle=\frac{1}{c_{X}^{2}}h_{00}+\frac{2}{c_{X}}h_{03}+h_{33}.$ (448)

$h_{\alpha 0}$ vanishes after the transformation while $h_{03}$, $h_{13}$ and
$h_{23}$ vanish before:

$\displaystyle h^{\prime}_{13}$ $\displaystyle=\frac{1}{c_{X}}h_{01},$
$\displaystyle h^{\prime}_{23}$ $\displaystyle=\frac{1}{c_{X}}h_{02},$
$\displaystyle h^{\prime}_{33}$
$\displaystyle=\frac{1}{c_{X}^{2}}h_{00}+h_{33}.$ (449)

We know that $h_{01}$ and $h_{02}$ are proportional to $\gamma_{1}$ and
$\gamma_{2}$, so we can replace the $1/c_{X}$ acting on them with $1/c_{V}$.
Similarly, we know that $h_{00}$ is proportional to $F$. In this case, we
replace $c_{X}$ with $c_{S}$.

$\displaystyle h^{\prime}_{13}$ $\displaystyle=\frac{1}{c_{V}}h_{01},$
$\displaystyle h^{\prime}_{23}$ $\displaystyle=\frac{1}{c_{V}}h_{02},$
$\displaystyle h^{\prime}_{33}$
$\displaystyle=\frac{1}{c_{S}^{2}}h_{00}+h_{33}.$ (450)

We apply this transformation to the metric perturbation:

$\displaystyle h_{\mu\nu}$
$\displaystyle=\left(\begin{array}[]{ccccccc}0&&0&&0&&0\\\ \\\
0&&\phi_{11}+\dfrac{F}{2}&&\phi_{12}&&\dfrac{\gamma_{1}}{c_{V}}\\\ \\\
0&&\phi_{12}&&-\phi_{11}+\dfrac{F}{2}&&\dfrac{\gamma_{2}}{c_{V}}\\\ \\\
0&&\dfrac{\gamma_{1}}{c_{V}}&&\dfrac{\gamma_{2}}{c_{V}}&&\alpha
F\end{array}\right).$ (458)

Here we have defined
$\alpha=1/\left(c_{S}^{2}c_{14}\right)-(1+c_{2})/c_{123}$. We can also easily
express the polarizations in terms of our five modes:

$\displaystyle h_{b}$ $\displaystyle=\frac{F}{2},$ $\displaystyle h_{x}$
$\displaystyle=\frac{\gamma_{1}}{c_{V}},$ $\displaystyle h_{+}$
$\displaystyle=\phi_{11},$ (459) $\displaystyle h_{l}$ $\displaystyle=\alpha
F,$ $\displaystyle h_{y}$ $\displaystyle=\frac{\gamma_{2}}{c_{V}},$
$\displaystyle h_{\times}$ $\displaystyle=\phi_{12}.$ (460)

Since $F$ excites a mixture of breathing and longitudinal mode we study it in
a bit more detail. Assuming that only $F$ is excited, the 3D line element
takes the following form:

$\displaystyle ds^{2}$
$\displaystyle=\left(1+\frac{F}{2}\right)\left(dx^{2}+dy^{2}\right)+\left(1+\alpha
F\right)dz^{2}.$ (461)

We can rewrite this in spherical coordinates with two parameters $\chi$ and
$\beta$:

$\displaystyle ds^{2}$ $\displaystyle=\chi^{2}\left(1+\beta\cos
2\theta\right)dr^{2},$ (462) $\displaystyle\chi$
$\displaystyle=\sqrt{\frac{1}{2}\left(2+\frac{F}{2}+\alpha F\right)},$ (463)
$\displaystyle\beta$
$\displaystyle=\frac{F\left(\alpha-1/2\right)}{2+F/2+\alpha F}.$ (464)

Note that while $\chi$ and $\beta$ on their own can be complex and singular
the relevant quantities $\chi^{2}$ and $\chi^{2}\beta$ never are.

We find that the equation for ${\rm d}s^{2}$ can be approximated as an ellipse
equation. Scalar Einstein-aether waves we deform a sphere into an ellipsoid.
In the $xz$ or $yz$ plane, we have a semimajor axis
$a=\chi\sqrt{1+\beta}\,{\rm d}r$ and an area of $\pi\chi^{2}\,{\rm d}r^{2}$.
Unlike the tensor or vector modes, the area of the ellipse cross section can
change depending on $\chi$ since $F$ excites the breathing mode. The total
effect of the scalar wave will depend on the value of $\alpha$, which depends
on the coupling parameters of Einstein-aether and can be positive, negative or
zero.

##  Appendix D Mathematica details

Throughout this paper, results were mostly first calculated by hand and then
double-checked using Wolfram Mathematica and the OGRe (Object-Oriented General
Relativity) package for Mathematica [99]. This package was also used for the
calculations in Sec. IV. For those calculations we kept the high frequency
approximation. This means all derivatives yield zero if they do not act on the
perturbation. Since no products of perturbations appear, we can then simply
treat $\partial_{\alpha}$ as a multiplier for the sake of speeding up our
calculation.

To shorten the result for the full kinetic matrix we used Mathematica’s
Simplify function and switched back and forth between using $\mu_{0}$ or
inserting $\mu_{0}=\sqrt{1+\mu_{1}^{2}+\mu_{2}^{2}+\mu_{3}^{2}}$. Still, the
smallest size we could obtain in Mathematica was 936.752 bytes. Obtaining the
exact mode speeds would entail finding the eigenvalues of the kinetic matrix
and setting them to zero. This is equivalent to solving $\det\hat{K}=0$. The
large kinetic matrix causes the determinant to blow up to 16 megabytes. This
means that finding the exact mode speeds without using more computational
power or different algorithms is unfeasible: all attempts ran for many hours
before it was decided to manually terminate them.

In Sec. IV.3 we applied the functions Simplify and FullSimplify to each order
of the kinetic matrix in order to shorten expressions.

## References

  * Berti  _et al._ [2015] E. Berti _et al._ , Testing General Relativity with Present and Future Astrophysical Observations, Class. Quant. Grav. 32, 243001 (2015), arXiv:1501.07274 [gr-qc] .
  * Clifton  _et al._ [2012] T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, Modified Gravity and Cosmology, Phys. Rept. 513, 1 (2012), arXiv:1106.2476 [astro-ph.CO] .
  * Shankaranarayanan and Johnson [2022] S. Shankaranarayanan and J. P. Johnson, Modified theories of gravity: Why, how and what?, Gen. Rel. Grav. 54, 44 (2022), arXiv:2204.06533 [gr-qc] .
  * Jacobson and Mattingly [2001] T. Jacobson and D. Mattingly, Gravity with a dynamical preferred frame, Phys. Rev. D 64, 024028 (2001), arXiv:gr-qc/0007031 .
  * Mattingly [2005] D. Mattingly, Modern tests of Lorentz invariance, Living Rev. Rel. 8, 5 (2005), arXiv:gr-qc/0502097 .
  * Oost  _et al._ [2018] J. Oost, S. Mukohyama, and A. Wang, Constraints on Einstein-aether theory after GW170817, Phys. Rev. D 97, 124023 (2018), arXiv:1802.04303 [gr-qc] .
  * Sarbach  _et al._ [2019] O. Sarbach, E. Barausse, and J. A. Preciado-López, Well-posed Cauchy formulation for Einstein-æther theory, Class. Quant. Grav. 36, 165007 (2019), arXiv:1902.05130 [gr-qc] .
  * Eling and Jacobson [2004] C. Eling and T. Jacobson, Static postNewtonian equivalence of GR and gravity with a dynamical preferred frame, Phys. Rev. D 69, 064005 (2004), arXiv:gr-qc/0310044 .
  * Eling and Jacobson [2006a] C. Eling and T. Jacobson, Spherical solutions in Einstein-aether theory: Static aether and stars, Class. Quant. Grav. 23, 5625 (2006a), [Erratum: Class.Quant.Grav. 27, 049801 (2010)], arXiv:gr-qc/0603058 .
  * Eling and Jacobson [2006b] C. Eling and T. Jacobson, Black Holes in Einstein-Aether Theory, Class. Quant. Grav. 23, 5643 (2006b), [Erratum: Class.Quant.Grav. 27, 049802 (2010)], arXiv:gr-qc/0604088 .
  * Jacobson [2010] T. Jacobson, Extended Horava gravity and Einstein-aether theory, Phys. Rev. D 81, 101502 (2010), [Erratum: Phys.Rev.D 82, 129901 (2010)], arXiv:1001.4823 [hep-th] .
  * Carroll and Lim [2004] S. M. Carroll and E. A. Lim, Lorentz-violating vector fields slow the universe down, Phys. Rev. D 70, 123525 (2004), arXiv:hep-th/0407149 .
  * Graesser  _et al._ [2005] M. L. Graesser, A. Jenkins, and M. B. Wise, Spontaneous Lorentz violation and the long-range gravitational preferred-frame effect, Phys. Lett. B 613, 5 (2005), arXiv:hep-th/0501223 .
  * Foster and Jacobson [2006] B. Z. Foster and T. Jacobson, Post-Newtonian parameters and constraints on Einstein-aether theory, Phys. Rev. D 73, 064015 (2006), arXiv:gr-qc/0509083 .
  * Yagi  _et al._ [2014] K. Yagi, D. Blas, E. Barausse, and N. Yunes, Constraints on Einstein-Æther theory and Hořava gravity from binary pulsar observations, Phys. Rev. D 89, 084067 (2014), [Erratum: Phys.Rev.D 90, 069902 (2014), Erratum: Phys.Rev.D 90, 069901 (2014)], arXiv:1311.7144 [gr-qc] .
  * Jacobson and Mattingly [2004] T. Jacobson and D. Mattingly, Einstein-Aether waves, Phys. Rev. D 70, 024003 (2004), arXiv:gr-qc/0402005 .
  * Abbott  _et al._ [2021a] R. Abbott _et al._ (LIGO Scientific, Virgo), Tests of general relativity with binary black holes from the second LIGO-Virgo gravitational-wave transient catalog, Phys. Rev. D 103, 122002 (2021a), arXiv:2010.14529 [gr-qc] .
  * Abbott  _et al._ [2021b] R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), Tests of General Relativity with GWTC-3 (2021b), arXiv:2112.06861 [gr-qc] .
  * Takeda  _et al._ [2022] H. Takeda, S. Morisaki, and A. Nishizawa, Search for scalar-tensor mixed polarization modes of gravitational waves, Phys. Rev. D 105, 084019 (2022), arXiv:2105.00253 [gr-qc] .
  * Takeda  _et al._ [2023] H. Takeda, Y. Manita, H. Omiya, and T. Tanaka, Scalar polarization window in gravitational-wave signals, PTEP 2023, 073E01 (2023), arXiv:2304.14430 [gr-qc] .
  * Eling [2006] C. Eling, Energy in the Einstein-aether theory, Phys. Rev. D 73, 084026 (2006), [Erratum: Phys.Rev.D 80, 129905 (2009)], arXiv:gr-qc/0507059 .
  * Elliott  _et al._ [2005] J. W. Elliott, G. D. Moore, and H. Stoica, Constraining the new Aether: Gravitational Cerenkov radiation, JHEP 2005 (08), 066, arXiv:hep-ph/0505211 .
  * Foster [2006] B. Z. Foster, Radiation damping in Einstein-aether theory, Phys. Rev. D 73, 104012 (2006), [Erratum: Phys.Rev.D 75, 129904 (2007)], arXiv:gr-qc/0602004 .
  * Foster [2007] B. Z. Foster, Strong field effects on binary systems in Einstein-aether theory, Phys. Rev. D 76, 084033 (2007), arXiv:0706.0704 [gr-qc] .
  * Ezquiaga and Zumalacárregui [2017] J. M. Ezquiaga and M. Zumalacárregui, Dark Energy After GW170817: Dead Ends and the Road Ahead, Phys. Rev. Lett. 119, 251304 (2017), arXiv:1710.05901 [astro-ph.CO] .
  * Creminelli and Vernizzi [2017] P. Creminelli and F. Vernizzi, Dark Energy after GW170817 and GRB170817A, Phys. Rev. Lett. 119, 251302 (2017), arXiv:1710.05877 [astro-ph.CO] .
  * Baker  _et al._ [2017] T. Baker, E. Bellini, P. G. Ferreira, M. Lagos, J. Noller, and I. Sawicki, Strong constraints on cosmological gravity from GW170817 and GRB 170817A, Phys. Rev. Lett. 119, 251301 (2017), arXiv:1710.06394 [astro-ph.CO] .
  * Sakstein and Jain [2017] J. Sakstein and B. Jain, Implications of the Neutron Star Merger GW170817 for Cosmological Scalar-Tensor Theories, Phys. Rev. Lett. 119, 251303 (2017), arXiv:1710.05893 [astro-ph.CO] .
  * Lombriser and Taylor [2016] L. Lombriser and A. Taylor, Breaking a Dark Degeneracy with Gravitational Waves, JCAP 2016 (03), 031, arXiv:1509.08458 [astro-ph.CO] .
  * Bettoni  _et al._ [2017] D. Bettoni, J. M. Ezquiaga, K. Hinterbichler, and M. Zumalacárregui, Speed of Gravitational Waves and the Fate of Scalar-Tensor Gravity, Phys. Rev. D 95, 084029 (2017), arXiv:1608.01982 [gr-qc] .
  * Abbott  _et al._ [2017a] B. P. Abbott _et al._ (LIGO Scientific, Virgo), GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral, Phys. Rev. Lett. 119, 161101 (2017a), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ [2017b] B. P. Abbott _et al._ (LIGO Scientific, Virgo, Fermi-GBM, INTEGRAL), Gravitational Waves and Gamma-rays from a Binary Neutron Star Merger: GW170817 and GRB 170817A, Astrophys. J. Lett. 848, L13 (2017b), arXiv:1710.05834 [astro-ph.HE] .
  * Abbott  _et al._ [2017c] B. P. Abbott _et al._ (LIGO Scientific, Virgo, Fermi GBM, INTEGRAL, IceCube, AstroSat Cadmium Zinc Telluride Imager Team, IPN, Insight-Hxmt, ANTARES, Swift, AGILE Team, 1M2H Team, Dark Energy Camera GW-EM, DES, DLT40, GRAWITA, Fermi-LAT, ATCA, ASKAP, Las Cumbres Observatory Group, OzGrav, DWF (Deeper Wider Faster Program), AST3, CAASTRO, VINROUGE, MASTER, J-GEM, GROWTH, JAGWAR, CaltechNRAO, TTU-NRAO, NuSTAR, Pan-STARRS, MAXI Team, TZAC Consortium, KU, Nordic Optical Telescope, ePESSTO, GROND, Texas Tech University, SALT Group, TOROS, BOOTES, MWA, CALET, IKI-GW Follow-up, H.E.S.S., LOFAR, LWA, HAWC, Pierre Auger, ALMA, Euro VLBI Team, Pi of Sky, Chandra Team at McGill University, DFN, ATLAS Telescopes, High Time Resolution Universe Survey, RIMAS, RATIR, SKA South Africa/MeerKAT), Multi-messenger Observations of a Binary Neutron Star Merger, Astrophys. J. Lett. 848, L12 (2017c), arXiv:1710.05833 [astro-ph.HE] .
  * Abbott  _et al._ [2019] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Tests of General Relativity with GW170817, Phys. Rev. Lett. 123, 011102 (2019), arXiv:1811.00364 [gr-qc] .
  * Takeda  _et al._ [2021] H. Takeda, S. Morisaki, and A. Nishizawa, Pure polarization test of GW170814 and GW170817 using waveforms consistent with modified theories of gravity, Phys. Rev. D 103, 064037 (2021), arXiv:2010.14538 [gr-qc] .
  * Gong  _et al._ [2018] Y. Gong, S. Hou, D. Liang, and E. Papantonopoulos, Gravitational waves in Einstein-æther and generalized TeVeS theory after GW170817, Phys. Rev. D 97, 084040 (2018), arXiv:1801.03382 [gr-qc] .
  * Gupta  _et al._ [2021] T. Gupta, M. Herrero-Valea, D. Blas, E. Barausse, N. Cornish, K. Yagi, and N. Yunes, New binary pulsar constraints on Einstein-æther theory after GW170817, Class. Quant. Grav. 38, 195003 (2021), arXiv:2104.04596 [gr-qc] .
  * Schumacher  _et al._ [2023] K. Schumacher, S. E. Perkins, A. Shaw, K. Yagi, and N. Yunes, Gravitational wave constraints on Einstein-æther theory with LIGO/Virgo data, Phys. Rev. D 108, 104053 (2023), arXiv:2304.06801 [gr-qc] .
  * Ezquiaga and Zumalacárregui [2018] J. M. Ezquiaga and M. Zumalacárregui, Dark Energy in light of Multi-Messenger Gravitational-Wave astronomy, Front. Astron. Space Sci. 5, 44 (2018), arXiv:1807.09241 [astro-ph.CO] .
  * Reyes  _et al._ [2010] R. Reyes, R. Mandelbaum, U. Seljak, T. Baldauf, J. E. Gunn, L. Lombriser, and R. E. Smith, Confirmation of general relativity on large scales from weak lensing and galaxy velocities, Nature 464, 256 (2010), arXiv:1003.2185 [astro-ph.CO] .
  * Collett  _et al._ [2018] T. E. Collett, L. J. Oldham, R. J. Smith, M. W. Auger, K. B. Westfall, D. Bacon, R. C. Nichol, K. L. Masters, K. Koyama, and R. van den Bosch, A precise extragalactic test of General Relativity, Science 360, 1342 (2018), arXiv:1806.08300 [astro-ph.CO] .
  * Mukherjee  _et al._ [2020] S. Mukherjee, B. D. Wandelt, and J. Silk, Probing the theory of gravity with gravitational lensing of gravitational waves and galaxy surveys, Mon. Not. Roy. Astron. Soc. 494, 1956 (2020), arXiv:1908.08951 [astro-ph.CO] .
  * Goyal  _et al._ [2021] S. Goyal, K. Haris, A. K. Mehta, and P. Ajith, Testing the nature of gravitational-wave polarizations using strongly lensed signals, Phys. Rev. D 103, 024038 (2021), arXiv:2008.07060 [gr-qc] .
  * Magaña Hernandez [2022] I. Magaña Hernandez, Measuring the polarization content of gravitational waves with strongly lensed binary black hole mergers (2022), arXiv:2211.01272 [gr-qc] .
  * Collett and Bacon [2017] T. E. Collett and D. Bacon, Testing the speed of gravitational waves over cosmological distances with strong gravitational lensing, Phys. Rev. Lett. 118, 091101 (2017), arXiv:1602.05882 [astro-ph.HE] .
  * Fan  _et al._ [2017] X.-L. Fan, K. Liao, M. Biesiada, A. Piorkowska-Kurpas, and Z.-H. Zhu, Speed of Gravitational Waves from Strongly Lensed Gravitational Waves and Electromagnetic Signals, Phys. Rev. Lett. 118, 091102 (2017), arXiv:1612.04095 [gr-qc] .
  * Finke  _et al._ [2021] A. Finke, S. Foffa, F. Iacovelli, M. Maggiore, and M. Mancarella, Probing modified gravitational wave propagation with strongly lensed coalescing binaries, Phys. Rev. D 104, 084057 (2021), arXiv:2107.05046 [gr-qc] .
  * Narola  _et al._ [2023] H. Narola, J. Janquart, L. Haegel, K. Haris, O. A. Hannuksela, and C. Van Den Broeck, How well can modified gravitational wave propagation be constrained with strong lensing? (2023), arXiv:2308.01709 [gr-qc] .
  * Takeda and Tanaka [2024] H. Takeda and T. Tanaka, Strong lensing of gravitational waves with modified propagation (2024), arXiv:2404.10809 [gr-qc] .
  * Chung and Li [2021] A. K.-W. Chung and T. G. F. Li, Lensing of gravitational waves as a novel probe of graviton mass, Phys. Rev. D 104, 124060 (2021), arXiv:2106.09630 [gr-qc] .
  * Deka  _et al._ [2024] U. Deka, S. Chakraborty, S. J. Kapadia, M. A. Shaikh, and P. Ajith, Probing black hole charge with gravitational microlensing of gravitational waves (2024), arXiv:2401.06553 [gr-qc] .
  * Ezquiaga and Zumalacárregui [2020] J. M. Ezquiaga and M. Zumalacárregui, Gravitational wave lensing beyond general relativity: birefringence, echoes and shadows, Phys. Rev. D 102, 124048 (2020), arXiv:2009.12187 [gr-qc] .
  * Goyal  _et al._ [2023] S. Goyal, A. Vijaykumar, J. M. Ezquiaga, and M. Zumalacárregui, Probing lens-induced gravitational-wave birefringence as a test of general relativity, Phys. Rev. D 108, 024052 (2023), arXiv:2301.04826 [gr-qc] .
  * Zhao  _et al._ [2020] W. Zhao, T. Zhu, J. Qiao, and A. Wang, Waveform of gravitational waves in the general parity-violating gravities, Phys. Rev. D 101, 024002 (2020), arXiv:1909.10887 [gr-qc] .
  * Wang  _et al._ [2022] Y.-F. Wang, S. M. Brown, L. Shao, and W. Zhao, Tests of gravitational-wave birefringence with the open gravitational-wave catalog, Phys. Rev. D 106, 084005 (2022), arXiv:2109.09718 [astro-ph.HE] .
  * Jenks  _et al._ [2023] L. Jenks, L. Choi, M. Lagos, and N. Yunes, Parametrized parity violation in gravitational wave propagation, Phys. Rev. D 108, 044023 (2023), arXiv:2305.10478 [gr-qc] .
  * Lagos  _et al._ [2024] M. Lagos, L. Jenks, M. Isi, K. Hotokezaka, B. D. Metzger, E. Burns, W. M. Farr, S. Perkins, K. W. K. Wong, and N. Yunes, Birefringence tests of gravity with multi-messenger binaries (2024), arXiv:2402.05316 [gr-qc] .
  * Will [1993] C. M. Will,  _Theory and Experiment in Gravitational Physics_ (Cambridge University Press, 1993).
  * Jacobson [2007] T. Jacobson, Einstein-aether gravity: A Status report, PoS QG-PH, 020 (2007), arXiv:0801.1547 [gr-qc] .
  * Takeda [2023] H. Takeda, Testing theories of gravity in extreme environments with gravitational-wave polarization from compact binary systems (2023), unpublished slides of a presentation for the Cosmology meeting of the Max Planck Institute for Gravitational Physics.
  * Nishizawa  _et al._ [2009] A. Nishizawa, A. Taruya, K. Hayama, S. Kawamura, and M.-a. Sakagami, Probing non-tensorial polarizations of stochastic gravitational-wave backgrounds with ground-based laser interferometers, Phys. Rev. D 79, 082002 (2009), arXiv:0903.0528 [astro-ph.CO] .
  * Ghosh  _et al._ [2021] R. Ghosh, A. K. Mishra, and S. Sarkar, Overcharging extremal black holes, Phys. Rev. D 104, 104043 (2021), arXiv:2106.10667 [gr-qc] .
  * Nordtvedt [1968] K. Nordtvedt, Equivalence Principle for Massive Bodies. 2. Theory, Phys. Rev. 169, 1017 (1968).
  * Nordtvedt [1969] K. Nordtvedt, Equivalence Principle for Massive Bodies Including Rotational Energy and Radiation Pressure, Phys. Rev. 180, 1293 (1969).
  * Will [1971] C. M. Will, Theoretical Frameworks for Testing Relativistic Gravity. 2. Parametrized Post-Newtonian Hydrodynamics, and the Nordtvedt Effect, Astrophys. J. 163, 611 (1971).
  * Will and Nordtvedt [1972] C. M. Will and K. Nordtvedt, Jr., Conservation Laws and Preferred Frames in Relativistic Gravity. I. Preferred-Frame Theories and an Extended PPN Formalism, Astrophys. J. 177, 757 (1972).
  * Muller  _et al._ [1996] J. Muller, K. Nordtvedt, and D. Vokrouhlicky, Improved constraint on the alpha-1 PPN parameter from lunar motion, Phys. Rev. D 54, R5927 (1996).
  * Shao and Wex [2012] L. Shao and N. Wex, New tests of local Lorentz invariance of gravity with small-eccentricity binary pulsars, Class. Quant. Grav. 29, 215018 (2012), arXiv:1209.4503 [gr-qc] .
  * Shao  _et al._ [2013] L. Shao, R. N. Caballero, M. Kramer, N. Wex, D. J. Champion, and A. Jessner, A new limit on local Lorentz invariance violation of gravity from solitary pulsars, Class. Quant. Grav. 30, 165019 (2013), arXiv:1307.2552 [gr-qc] .
  * Muller  _et al._ [1999] U. Muller, C. Schubert, and A. M. E. van de Ven, A Closed formula for the Riemann normal coordinate expansion, Gen. Rel. Grav. 31, 1759 (1999), arXiv:gr-qc/9712092 .
  * Laboratory [2015] L. Laboratory, LIGO R&D (2015), accessed: 2023-04-14.
  * Bland-Hawthorn and Gerhard [2016] J. Bland-Hawthorn and O. Gerhard, The Galaxy in Context: Structural, Kinematic, and Integrated Properties, Annu. Rev. Astron. Astrophys. 54, 529 (2016), arXiv:1602.07702 [astro-ph.GA] .
  * Tambalo  _et al._ [2023a] G. Tambalo, M. Zumalacárregui, L. Dai, and M. H.-Y. Cheung, Gravitational wave lensing as a probe of halo properties and dark matter, Phys. Rev. D 108, 103529 (2023a), arXiv:2212.11960 [astro-ph.CO] .
  * Schneider  _et al._ [1992] P. Schneider, J. Ehlers, and E. E. Falco,  _Gravitational Lenses_ (Springer, 1992).
  * Takahashi and Nakamura [2003] R. Takahashi and T. Nakamura, Wave effects in gravitational lensing of gravitational waves from chirping binaries, Astrophys. J. 595, 1039 (2003), arXiv:astro-ph/0305055 .
  * Tambalo  _et al._ [2023b] G. Tambalo, M. Zumalacárregui, L. Dai, and M. H.-Y. Cheung, Lensing of gravitational waves: Efficient wave-optics methods and validation with symmetric lenses, Phys. Rev. D 108, 043527 (2023b), arXiv:2210.05658 [gr-qc] .
  * Dalang  _et al._ [2021] C. Dalang, P. Fleury, and L. Lombriser, Scalar and tensor gravitational waves, Phys. Rev. D 103, 064075 (2021), arXiv:2009.11827 [gr-qc] .
  * Boran  _et al._ [2018] S. Boran, S. Desai, E. O. Kahya, and R. P. Woodard, GW170817 Falsifies Dark Matter Emulators, Phys. Rev. D 97, 041501 (2018), arXiv:1710.06168 [astro-ph.HE] .