# Probing the speed of gravity with LVK, LISA, and joint observations

Ian Harry  Corresponding author<EMAIL_ADDRESS>Institute of Cosmology &
Gravitation, University of Portsmouth, Portsmouth, PO1 3FX, U.K.  Johannes
Noller  Corresponding author<EMAIL_ADDRESS>Institute of
Cosmology & Gravitation, University of Portsmouth, Portsmouth, PO1 3FX, U.K.
DAMTP, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, U.K.

###### Abstract

Theories of dark energy that affect the speed of gravitational waves $c_{\rm
GW}$ on cosmological scales naturally lead to a frequency-dependent transition
of that speed close to the LIGO/Virgo/KAGRA (LVK) band. While observations
such as GW170817 assure us that $c_{\rm GW}$ is extremely close to the speed
of light in the LVK band, a frequency-dependent transition below the LVK band
is a smoking-gun signal for large classes of dynamical dark energy theories.
Here we discuss 1) how the remnants of such a transition can be constrained
with observations in the LVK band, 2) what signatures are associated with such
a transition in the LISA band, and 3) how joint observations in the LVK and
LISA bands allow us to place tight constraints on this transition and the
underlying theories. We find that deviations of $c_{\rm GW}$ can be
constrained down to a level of $\sim 10^{-17}$ in the LVK and LISA bands even
for mild frequency-dependence, much stronger than existing bounds for
frequency-independent $c_{\rm GW}\neq c$. We use the strain data from GW170817
to bound the deviation of $c_{\rm GW}$ to be less than $10^{-17}$ at 100 Hz
and less than $10^{-18}$ at 500 Hz. We also identify a particularly
interesting type of transition in between the LVK and LISA bands and show how
multi-band observations can constrain this further. Finally, we discuss what
these current and forecasted constraints imply for the underlying dark energy
theories.

###### Contents

  1. I Introduction
  2. II From theory to templates
    1. II.1 Dark energy and $c_{\rm GW}\neq c$
    2. II.2 The reach of a theory and frequency-dependent $c_{\rm GW}$
    3. II.3 Templates for $c_{\rm GW}$
  3. III The gravitational wave signal
    1. III.1 Waveform stretching or squeezing
    2. III.2 Arrival times
    3. III.3 Intrinsic source effects
  4. IV Probing $c_{\rm GW}$ in the LVK band
    1. IV.1 Including waveform stretching/squeezing in the gravitational wave signal
    2. IV.2 Parameterizing $\delta c_{\rm GW}(f)$
    3. IV.3 Approximate measurement of observable $c_{\rm GW}$ deviation
    4. IV.4 Bayesian parameter estimation constraints of $c_{\rm GW}$ using GW170817
  5. V Probing $c_{\rm GW}$ in the LISA band
  6. VI Probing $c_{\rm GW}$ with joint LVK/LISA observations
  7. VII Conclusions
  8. A $c_{\rm GW}$ in Horndeski gravity
  9. B More $c_{\rm GW}$ templates
  10. C Intrinsic source evolution with non-constant $c_{\rm GW}$
  11. D Constraints using only LVK or ET observations

##  I Introduction

One of the most fundamental observables in testing our understanding of
gravity is the speed with which gravity propagates. In a Lorentz-invariant
solution, gravitational waves will propagate at the speed of light, but this
need not be the case if Lorentz invariance is spontaneously broken, as it is
in our Universe. Here, just as for any other waves, the speed of gravitational
waves $c_{\rm GW}$ depends on the medium they are propagating through.
Measuring $c_{\rm GW}$ therefore lets us probe the nature of this medium. With
dark energy accounting for close to 70% of the energy budget of our current
Universe Planck:2018vyg , constraining $c_{\rm GW}$ provides an especially
powerful tool in hunting for new degrees of freedom potentially associated
with dark energy or a modified theory of gravity. Such new degrees of freedom,
which non-trivially affect the background evolution and hence the universal
‘medium’, are a generic consequence of departures from modelling dark energy
as a cosmological constant $\Lambda{}$. Tight constraints found on some such
theories Baker:2017hug ; Creminelli:2017sry ; Sakstein:2017xjx ;
Ezquiaga:2017ekz  – also see earlier related work Amendola:2012ky ;
Amendola:2014wma ; Deffayet:2010qz ; Linder:2014fna ; Raveri:2014eea ;
Saltas:2014dha ; Lombriser:2015sxa ; Lombriser:2016yzn ; Jimenez:2015bwa ;
Bettoni:2016mij ; Sawicki:2016klv  – following the joint observations of
GW170817 and GRB 170817A TheLIGOScientific:2017qsa ; 2041-8205-848-2-L14 ;
2041-8205-848-2-L15 ; LIGOScientific:2017zic ; LIGOScientific:2017ync  neatly
illustrate this point.

Current constraints on the speed of gravitational waves, $c_{\rm GW}$,
strongly depend on the frequency range in question and are conveniently
expressed in terms of the parameter

$\displaystyle\alpha_{T}\equiv(c_{\rm GW}^{2}-c^{2})/c^{2},$ (1)

where $c$ is the speed of light. Note that $\alpha_{T}$ therefore is a measure
of the ‘tensor speed excess’, i.e. it is positive for $c_{\rm GW}>c$ and
negative when $c_{\rm GW}<c$. We will remain agnostic here and will not assume
anything about the sign of $\alpha_{T}$ in what follows111It is useful to
recall that for a spontaneously Lorentz-breaking solution, $c_{\rm GW}\leq c$
need not be imposed by causality requirements – in fact the opposite can be
the case, see e.g. deRham:2019ctd ; deRham:2021fpu  and references therein..
Constraints from the cosmic microwave background and large scale structure
(see Noller:2018wyv ; Bellini:2015xja ; Hu:2013twa ; Raveri:2014cka ;
Gleyzes:2015rua ; Kreisch:2017uet ; Zumalacarregui:2016pph ; Alonso:2016suf ;
Arai:2017hxj ; Frusciante:2018jzw ; Reischke:2018ooh ; Mancini:2018qtb ;
Noller:2018eht ; Arai:2019zul ; Brando:2019xbv ; Arjona:2019rfn ;
Raveri:2019mxg ; Perenon:2019dpc ; Frusciante:2019xia ; Arai:2019zul ;
SpurioMancini:2019rxy ; Bonilla:2019mbm ; Baker:2020apq ; Joudaki:2020shz ;
Noller:2020lav ; Noller:2020afd ; Traykova:2021hbr  and references therein),
corresponding to frequencies $f\sim 10^{-18}-10^{-14}\;{\rm Hz}$, require
$|\alpha_{T}|\lesssim{\cal O}(1)$. Binary pulsars, in particular the Hulse-
Taylor binary, impose $|\alpha_{T}|\lesssim{\cal O}(10^{-2})$ for frequencies
$f\sim 10^{-4}\;{\rm Hz}$ Jimenez:2015bwa . Gravitational wave observations of
the afore-mentioned GW170817 event then place an extremely tight bound of
$|\alpha_{T}|\lesssim{\cal O}(10^{-15})$ for frequencies $f\sim
10^{1}-10^{3}\;{\rm Hz}$ TheLIGOScientific:2017qsa ; 2041-8205-848-2-L14 ;
2041-8205-848-2-L15 ; LIGOScientific:2017zic ; LIGOScientific:2017ync  , while
the presence of high energy cosmic rays (and inferred absence of gravitational
Cherenkov radiation) leads to a $\alpha_{T}\gtrsim-{\cal O}(10^{-15})$ bound
at energies of $\sim 10^{10}$ GeV Moore:2001bv , i.e. for very high
frequencies $f\sim 10^{34}\;{\rm Hz}$. Note that the last bound does depend on
the sign of $\alpha_{T}$. Finally, several forecasted bounds exist for future
LISA observations Littenberg:2019mob ; Baker:2022rhh . Since these are
particularly relevant in the context of this paper, we will discuss them in
more detail in section V. The frequency-dependence of the bounds on
$\alpha_{T}$ is particularly important, since the theoretical predictions for
this parameter are also frequency-dependent. More specifically, for large
classes of dark-energy related models that yield a significantly non-zero
$\alpha_{T}$ on cosmological scales, i.e. at very low frequencies, one
naturally expects a transition back towards $\alpha_{T}=0$ close to or
somewhat below $10^{2}\;{\rm Hz}$ deRham:2018red . This means that the tight
GW170817 constraint on $\alpha_{T}$ is perfectly consistent with models that
yield a detectable, non-zero $\alpha_{T}$ at lower frequencies.

Current ground-based observatories, Advanced LIGO LIGOScientific:2014pky ,
Advanced Virgo VIRGO:2014yos  and KAGRA Aso:2013eba  have sensitivity between
approximately $20\;{\rm Hz}$ and $\sim 2000\;{\rm Hz}$. In the coming years,
these observatories will continue to run a series of increasingly sensitive
observing runs, with lower frequency sensitivity reaching no lower than
$10\;{\rm Hz}$ KAGRA:2013rdx . It is expected that these runs will yield
numerous observations of binary neutron star mergers, and hopefully numerous
multimessenger observations KAGRA:2013rdx , with which powerful tests of the
speed of gravity can be performed. However, to perform tests outside of this
frequency range we must wait for the next generation of gravitational-wave
observatories, expected to become operational in the 2030s. The third
generation of ground-based observatories, Einstein Telescope Punturo:2010zza
and Cosmic Explorer Reitze:2019iox , will increase overall sensitivity, and
increase the range of sensitive frequencies, covering from $\sim 1\;{\rm Hz}$
up to $\sim 10^{4}\;{\rm Hz}$. Perhaps of most interest in this context though
are the proposed space-based observatories, LISA LISA:2017pwj  and TianQin
TianQin:2015yph , which will cover a range of frequencies inaccessible from
the ground, covering $\sim 10^{-4}\;{\rm Hz}$ to $\sim 10^{-1}\;{\rm Hz}$.

In this paper we investigate how observations from the LIGO/Virgo/KAGRA
observatories (henceforth LVK), and GW170817 in particular, can be leveraged
to provide additional constraints on a frequency-dependent transition of
$c_{\rm GW}$, extending existing constraints that assume a frequency-
independent $c_{\rm GW}$. We also explore how one can use lower frequency
measurements of gravitational waves in the LISA band to test for the existence
and nature of such a transition and how one might be able to combine
information from both space-based and ground-based observation bands to
provide further constraints.

The paper is structured as follows. In section II we discuss in what sense the
frequency-dependence of $c_{\rm GW}$ is a generic phenomenon in models that do
lead to a non-zero $\alpha_{T}$ on cosmological scales. As a key outcome of
this discussion, we introduce and motivate templates to capture this frequency
dependence. In section III we investigate what the imprint of a frequency-
dependent $c_{\rm GW}$ is on the gravitational wave signal observed,
highlighting the resulting stretching and squeezing of the waveform as the
main associated observable. In the following sections we then compute
constraints and forecasts on the frequency-dependence of $c_{\rm GW}$: In
section IV we present constraints computed for the LVK band, that extend
existing bounds on a constant and hence frequency-independent $c_{\rm GW}$. In
section V we present analogous forecasts for the LISA band, and in section VI
we discuss how joint observations in both bands will allow us to place highly
precise constraints on $c_{\rm GW}$ over the whole range of relevant frequency
and energy scales. We conclude and summarise our findings in VII and collect
additional details in the appendices.

##  II From theory to templates

In this section we will illustrate how a non-luminal $c_{\rm GW}$ can arise in
theories of dark energy (II.1), why this generically leads to a $c_{\rm GW}$
that is frequency-dependent, where (if indeed present) one would expect this
to be especially pronounced close to the frequencies observed by present (LVK)
and near-future (LISA/TianQin/Einstein Telescope/Cosmic Explorer) detectors
(II.2), and how one can use this information to build theoretically well-
motivated templates to be used in constraining and searching for frequency-
dependent $c_{\rm GW}$ (II.3).

###  II.1 Dark energy and $c_{\rm GW}\neq c$

Since GR is the single consistent theory of a massless spin-2 field, testing
for (potentially dark energy-related) deviations away from it generically
amounts to probing the presence of new gravitational degrees of freedom.
Scalar-tensor theories are a minimal deviation from GR in this sense, as they
only introduce a single additional degree of freedom. Accordingly, Horndeski
gravity Horndeski:1974wa ; Deffayet:2011gz 222For the equivalence between the
formulations of Horndeski:1974wa  and Deffayet:2011gz , see Kobayashi:2011nu
., the most general Lorentz-invariant scalar-tensor action that gives rise to
second order equations of motion, has recently been the main workhorse in
testing for deviations from GR. In appendix A we will summarise the key
relevant results for the complete Horndeski theory, but here we will
illustrate the salient points with the following example scalar-tensor
theory333We emphasise that we merely use (2) as an example to illustrate
generic $\alpha_{T}$ features in scalar-tensor theories. While they will not
be important here, we refer the reader curious about some of the interesting
stability/screening/positivity properties of (2) to Kobayashi:2011nu ;
Bellini:2014fua ; Berezhiani:2013dw ; Koyama:2013paa ; Melville:2019wyy  and
references therein.

$\displaystyle{\cal L}=\frac{M_{\rm
Pl}^{2}}{2}R-\frac{1}{2}\nabla_{\mu}\phi\nabla^{\mu}\phi+\frac{M_{\rm
Pl}^{2}}{\Lambda^{3}}g(\phi)G_{\mu\nu}\nabla^{\mu}\nabla^{\nu}\phi.$ (2)

The first two terms are the standard kinetic interactions for the massless
tensor $g_{\mu\nu}$ and the scalar $\phi$, respectively, while the last term
is a higher-order interaction involving both the scalar and tensor.
$G_{\mu\nu}$ is the Einstein tensor and $g(\phi)$ is a dimensionless function
of $\phi/M_{\rm Pl}$. There are two key mass scales: the Planck mass $M_{\rm
Pl}$, and the scale $\Lambda$ associated with second derivatives of the
scalar. In cosmology these are linked to the value of the Hubble constant
today, $H_{0}$, and are conventionally taken to satisfy $\Lambda^{3}=M_{\rm
Pl}H_{0}^{2}$. This choice ensures that all interactions can give ${{\cal
O}(1)}$ contributions to the (cosmological) background evolution, so is
phenomenologically motivated in order to allow the scalar $\phi$ to act as
dark energy.

We can now work out the speed of gravitational waves $c_{\rm GW}$ for the
theory (2). In a scalar-tensor theory as we are considering here,
gravitational waves propagate through a ‘medium’ provided by the presence of
the scalar $\phi$. The interactions of $\phi$ with $g_{\mu\nu}$ encoded in (2)
can affect the refractive index of this medium and hence $c_{\rm GW}$. More
specifically, we find

$\displaystyle\alpha_{T}$
$\displaystyle=\frac{2(\nabla\bar{\phi})^{2}g^{\prime}(\bar{\phi})}{M_{\rm
Pl}\Lambda^{3}-(\nabla\bar{\phi})^{2}g^{\prime}(\bar{\phi})},$ (3)

where a prime denotes a derivative with respect to the argument of the
function (we recall this is $\phi/M_{\rm Pl}$ for $g$ and so both $g$ and
$g^{\prime}$ are dimensionless) and $\bar{\phi}$ is the background solution
for the scalar. This makes it clear that for a Lorentz-invariant background
solution with $\bar{\phi}=0$ we obtain $\alpha_{T}=0$ and hence gravitational
waves travel at the speed of light there, as expected. However, when
considering the cosmological background solution $\bar{\phi}=\bar{\phi}(t)$ of
interest here, this time-dependent scalar evolution spontaneously breaks
Lorentz invariance and provides a non-trivial medium for gravitational waves
to travel through. In this case, $\alpha_{T}$ will generically be non-zero and
we obtain the following dispersion relation

$\displaystyle\omega^{2}=c_{\rm GW}^{2}(t)k^{2},$ (4)

where $c_{\rm GW}^{2}=(1+\alpha_{T})c^{2}$ as above. We emphasise that this
means $\alpha_{T}$ is generically time-dependent, but (so far) frequency-
independent – we will see how frequency-dependence enters below. Regarding the
time-dependence, note that the characteristic timescale for this in theories
of dark energy is a Hubble time ($1/H_{0}\sim 10^{10}$ years). Compared to the
characteristic travel time for the signals seen by gravitational wave
detectors this is approx. the same order of magnitude for a source at a
distance of $\sim$ Gpc (travel time $\sim 3\times 10^{9}$ years) and about two
orders of magnitude larger for closer sources (the travel time is $\sim
10^{8}$ years for a source at $40$ Mpc), assuming the signal travels at close
to the speed of light. While taking into account this time-dependence can
therefore have a noticeable effect for far-away sources, e.g. modulating some
of the relevant observables such as arrival times, this will not change the
qualitative, order-of-magnitude constraints we are focusing on here. We will
therefore ignore any time-dependence here and leave an exploration of this
effect to future work444Note that, while dark energy only becomes a major
driver of the cosmological evolution from $z\sim{\cal O}(1)$ onwards, it only
takes light approx. twice as much time to travel to us from redshift 1 vs.
from redshift 1000, so this difference is immaterial to what we are doing
here..

###  II.2 The reach of a theory and frequency-dependent $c_{\rm GW}$

Let us return to the example theory (2). The higher-order interactions encoded
in the final term mean this theory is predictive at most up to the energy
scale $\Lambda$ and cannot resolve energy scales close to or larger than this
– see e.g. Georgi:1993mps ; Donoghue:1995cz ; Burgess:2007pt  for reviews of
the underlying field theoretic reasoning. This is precisely analogous to what
happens in GR, which also encodes higher-order interactions within the
Einstein-Hilbert term and where predictivity is lost as one approaches the
corresponding energy scale there, namely the Planck scale. For cosmologically
motivated models, where the scalar $\phi$ is linked to dark energy, we recall
that $\Lambda^{3}=M_{\rm Pl}H_{0}^{2}$. Expressed as a frequency this amounts
to

$\displaystyle\Lambda\sim 260\;{\rm Hz}.$ (5)

This is the largest possible energy/frequency scale, a so-called cutoff, where
(2) stops being applicable. Making predictions around or above this scale then
requires knowledge of a more complete high energy theory, a UV completion,
that extends the regime of validity of this theory beyond the frequency/energy
scale $\Lambda$. In the absence of knowledge about such a UV completion, i.e.
in the field’s present situation, we simply do not know how to relate a
measurement at energies $\gtrsim\Lambda$ to a theory such as (2). Note that
this conclusion does not rely on the specific form of (2) and affects all
sectors of the theory involving $\phi$ interactions555More specifically, it is
precisely the interactions that affect $\alpha_{T}$ which depend on $\Lambda$
and hence limit the applicability of the theory to energy scales below
$\Lambda$. In other words, if $\phi$ interactions in a Horndeski theory
contribute to the dark energy evolution today and affect $\alpha_{T}$ on
cosmological scales, this naturally leads to a cutoff (5). So this conclusion
does indeed not rely on the specific example (2). Some (but not all) well-
known dark energy theories that do not affect $\alpha_{T}$ can have a much
larger cutoff, e.g. quintessence Ratra:1987rm ; Wetterich:1987fm ;
Ferreira:1997hj ; Caldwell:1997ii  or k-essence ArmendarizPicon:1999rj ;
Garriga:1999vw ; ArmendarizPicon:2000dh ; ArmendarizPicon:2000ah  theories..
It is worth emphasizing that $\Lambda$ is the largest possible scale up to
which such a theory can be predictive. Where precisely this threshold is
depends on the specific nature of the UV completion. Compare this with the
Fermi theory for weak interactions, where the naive cutoff scale is $\sim$
TeV, yet new physics associated with a UV completion already enters an order
of magnitude below at the scale of the W-boson at $\sim 80$ GeV. The LISA and
lower end of the LVK bands (as well as intermediate bands) are therefore
particularly well-motivated regions to look for the onset of new physics
associated with UV completions of dark energy theories.

How is the above discussion of regimes of validity for scalar-tensor theories
relevant to the main observable we are considering in this paper, $c_{\rm
GW}$? This is discussed in detail in deRham:2018red  and we here summarise the
essential results relevant for this paper. When deriving (3) and (4) we were
implicitly working on cosmological scales, firmly within the regime of
validity of (2). However, as one approaches the cutoff of the theory
$\lesssim\Lambda$, the existence of a UV completion will manifest itself by
additional interactions (highly suppressed on cosmological scales) becoming
important. Precisely what these interactions are is determined by the
(unknown) UV completion, but we can infer some types of candidate additional
interactions from investigating radiative corrections to generic Horndeski
theories Luty:2003vm ; Nicolis:2004qq ; deRham:2010eu ; Burrage:2010cu ;
deRham:2014wfa ; Pirtskhalava:2015nla ; deRham:2012ew ; Goon:2016ihr ;
Saltas:2016nkg ; Noller:2018eht ; Heisenberg:2019wjv ; Heisenberg:2020cyi ;
Goon:2020myi . The example set of such interactions focused on in
deRham:2018red  is

$\displaystyle{\cal L}^{(n)}=M_{\rm
Pl}^{2}G_{\mu\nu}\frac{\Box^{n}}{M^{2n+4}}(\nabla^{\mu}\phi\nabla^{\nu}\phi),$
(6)

where $M\lesssim\Lambda$ is the scale associated with these UV-induced new
operators and $n\geq 2$. Of particular relevance for our context is the
(representative and generic) higher-derivative nature of these interactions.
Including these interactions modifies the dispersion relation (4) and one now
finds deRham:2018red

$\displaystyle\omega^{2}\sim c_{\rm GW}^{2}(0)k^{2}+\sum_{n\geq
2}\frac{c_{n}\Lambda^{6}}{3M^{2n+4}}(c^{2}k^{2}-\omega^{2})^{n-1}(\omega^{2}+{\cal
O}(k^{2})),$ (7)

where $\omega$ is the angular frequency, $k$ the wavenumber, the $c_{n}$ are
some order one coefficients (whose precise value is not important here) and we
have ignored suppressed $H$-dependent corrections, since $H^{2}/\omega^{2}\ll
1$ for the frequencies relevant here. Here we explicitly write $c_{\rm GW}$ as
a function of $k$ and hence denote the asymptotic value of $c_{\rm GW}$ for
low energies/frequencies, i.e. effectively its cosmological value derived in
(3), by $c_{\rm GW}(0)$ – see deRham:2018red  for details on the derivation of
the full dispersion relation for this example as well as other related
examples. Crucially, (7) asymptotes to the cosmological $c_{\rm GW}=c_{\rm
GW}(0)$ for small $k$ (low frequencies) and asymptotes to luminal $c_{\rm
GW}=c$ for large $k$ (high frequencies), where the
$(c^{2}k^{2}-\omega^{2})^{n-1}$ drives this large $k$ approach to luminality.
Here it is important to emphasise the point made by deRham:2018red  that
asymptoting to $c_{\rm GW}=c$ at high energies, as is entailed by (7), is a
generic consequence of any UV completion as long as Lorentz invariance is
restored at scales $\gtrsim M$. So this feature should certainly not be seen
as an accident of the example interactions shown here.

Finally, it is instructive to consider the asymptotic scaling of $c_{\rm GW}$
with $k$ deriving from (7). We focus on a single $n$ interaction in (6) and
denote the values of $k$ and $f$ corresponding to the transition scale
$M\lesssim\Lambda$ by $k_{\star}$ and $f_{\star}$, respectively, in what
follows. Doing so, we find

$\displaystyle k\ll k_{\star}:$ $\displaystyle c_{\rm GW}-c_{\rm GW}(0)$
$\displaystyle\propto k^{2n-2},$ $\displaystyle k\gg k_{\star}:$
$\displaystyle c-c_{\rm GW}$ $\displaystyle\propto k^{-2}.$ (8)

We see that there is a clean power law scaling in both limits, where the
relevant power for small $k$ depends on the precise interaction term present,
while there is a universal $1/k^{2}$ approach to luminal $c_{\rm GW}$ for
large $k$. With an eye on constructing templates below, it is very important
to separate what is a generic feature any template should recover vs. what is
an accident of the example chosen. At scales below the transition scale
$k_{\star}$, a UV completion will indeed generically give rise to higher-
derivative interactions such as (6). In the dispersion relation these will
manifest themselves as higher powers of $k$ suppressed by the scale $M$, so
generically one does indeed expect that a power law scaling is an excellent
approximation for the low energy approach into the $c_{\rm GW}$ transition.
However, once reaching energy scales $\sim M$ (or equivalently $k_{\star}$) no
firm prediction of the transition behaviour can be made without the (unknown)
UV completion. In particular, the transition can in principle happen
arbitrarily quickly and need not even follow a power law. So the universal
$1/k^{2}$ scaling found for the simple example above should not be hard-coded
into any template666As we shall see in more detail later, transitions that do
have a ‘slow’ $1/k^{2}$ scaling at high energies are severely constrained by
LVK measurements such as GW170817 even when the analogue of $f_{\star}$ (and
hence $k_{\star}$) are firmly within the LISA band. But while this is of
course an interesting finding for transitions with this scaling in its own
right, it should not be mistaken for a generic conclusion.. Furthermore, when
going beyond the simple example above and considering a generic situation
where all $c_{n}$ can be non-zero, a different large $k$ scaling can be
obtained for a given $c_{\rm GW}(0)$ by carefully tuning the $c_{n}$. Judging
whether such a tuning is natural again requires detailed knowledge of the UV
completion, but this already provides further reason why no fixed power-law
dependence for large $k$ should be hard-coded into any template777We thank
Scott Melville for related discussions..

###  II.3 Templates for $c_{\rm GW}$

In the sections below we will use current/forecasted observational data from
LVK and LISA bands to place constraints on the frequency-dependence of $c_{\rm
GW}$. Doing so generically requires a template for the waveform and, in
particular, for $c_{\rm GW}$ itself. Given the above considerations and
plethora of potentially contributing interactions, the precise form of $c_{\rm
GW}(f)$ in principle depends on a large number of parameters. This includes
the $c_{n}$ in (7), and crucially parameters associated to the unknown UV
completion that generically becomes important in or close to the LVK/LISA
bands for dark energy theories leading to a non-luminal $c_{\rm GW}$ at very
low frequencies, i.e. for cosmology.

We would therefore like to work with a template for $c_{\rm GW}$ that, while
simplifying the parametric dependence and making an analysis practical, still
captures the salient features of this frequency-dependence. As discussed
above, such a template ought to asymptote to a luminal $c_{\rm GW}$ at large
frequencies (to enable consistency with LVK bounds), asymptote to a constant
$c_{\rm GW}(0)$ at small frequencies (consistent with (3)), while also
allowing different scalings with $k$ throughout the transition. A minimal
template will depend on at least three parameters

  * •

A speed $c_{\rm GW}(0)$, which corresponds to the asymptotic speed of
gravitational waves at low frequencies. Effectively this is the cosmological
$c_{\rm GW}$ that can differ significantly from $c$.

  * •

A frequency scale $f_{\star}$, denoting the ‘central frequency’ of the
transition from $c_{\rm GW}(0)$ to $c$. In terms of the underlying physics
this is set by the scale where new physics associated with the UV completion
becomes important, e.g. $M$ in (6).

  * •

A parameter $\sigma$, which controls how quickly the transition takes place.
This is effectively a measure of the nature of the new interactions present
due to the UV completion.

Note that there can in principle be many more parameters controlling
interactions that determine how quickly the transition takes place, so
introducing a single parameter for this should be seen as a lowest order
approximation.

Using the parameters listed above, we can build the following useful
hyperbolic fitting function

$\displaystyle\Theta_{\pm}(f,\sigma,f_{\star})$
$\displaystyle\equiv\tfrac{1}{2}\pm\tfrac{1}{2}\tanh\left[\sigma\cdot\log\left(f/f_{\star}\right)\right],$
(9)

where $\Theta_{+}$ transitions between 0 and 1 around $f_{\star}$, while
$\Theta_{-}$ transitions from 1 to 0 around the same frequency. $\sigma$
controls how quickly this transition takes place, as advertised. With this
function we can now build a straightforward template for $c_{\rm GW}$ as a
function of frequency $f$, that captures the essential features outlined above

$\displaystyle c_{\rm GW}(f,\sigma,f_{\star})=c_{\rm GW}^{(0)}+(c-c_{\rm
GW}^{(0)})\Theta_{+}.$ (10)

This behaviour is illustrated in Figure 1.

It is again instructive to consider asymptotic scalings with $k$. From (10) we
find

$\displaystyle k\ll k_{\star}:$ $\displaystyle c_{\rm GW}-c_{\rm GW}(0)$
$\displaystyle\propto k^{2\sigma},$ $\displaystyle k\gg k_{\star}:$
$\displaystyle c-c_{\rm GW}$ $\displaystyle\propto k^{-2\sigma}.$ (11)

This clearly shows that the steepness parameter $\sigma$ directly controls the
power law scaling in both asymptotic limits, allowing us to mimic different
UV-completion-induced scalings at high energies. Due to the simplicity of the
template, the single parameter $\sigma$ controls both asymptotes and hence we
have symmetric $k^{2\sigma}$ and $k^{-2\sigma}$ scalings for this template.
There is no fundamental reason why the transition ought to be symmetric in
this way and it is straightforward to refine the template to include
asymmetric transitions upon introducing an additional parameter – see appendix
B888We also emphasise that the most highly UV sensitive (and hence least
robust) features of the transition, e.g. the precise functional form around
$f_{\star}$ in (7) (which are generically different for the small and large
frequency asymptote) are not reproduced by the template by design.. However,
for the constraint analysis in the following sections we will work with the
minimal template (10) and leave a more detailed analysis of more refined
templates including additional corrections for future work.

Figure 1:  The speed of gravitational waves, $c_{\rm GW}$, as a function of
frequency $f$ for the templates discussed in section II. The different curves
correspond to different choices of the ‘steepness’ parameter $\sigma$ that
also controls the asymptotic power law scalings of the $c_{\rm GW}$ – see
(11). $c_{\rm GW}$ neatly interpolates between $c_{\rm GW}(0)$, the speed of
GWs on cosmological scales, and the speed of light $c$ at large frequencies.
$f_{\star}$ physically gets set by the scale $M\lesssim\Lambda$ that marks the
onset of the UV completion for the dark energy theory in question, cf. (6).

##  III The gravitational wave signal

Before proceeding to probing $c_{\rm GW}$ using current and forecasted data,
we will find it useful to understand some key features of the GW signal
itself.

###  III.1 Waveform stretching or squeezing

If any detectable form of frequency-dependence (we will be more quantitative
below) takes place within the observable LVK and/or LISA bands, where the
source is tens to thousands of megaparsecs away from the detectors, the
leading-order effect of a frequency-dependent $c_{\rm GW}$ will be that the
observed signal is stretched (for negative $\alpha_{T}$) or squeezed (for
positive $\alpha_{T}$) as the travel time varies with frequency.

The travel time of a signal between a source and an observer on the Earth can
be expressed as

$\displaystyle t=\frac{D_{L}}{c_{\rm GW}(f)}\equiv\frac{D_{L}}{c(1-\delta
c_{\rm GW}(f))}\approx\frac{D_{L}}{c}+\frac{D_{L}}{c}\delta c_{\rm GW}(f),$
(12)

where ${D_{L}}$ is the luminosity distance to the source. We define here
$\delta c_{\rm GW}(f)\equiv 1-c_{\rm GW}(f)/c$, which is the fractional
variation in the gravitational-wave speed from $c$. This is related to
$\alpha_{T}$ according to $\alpha_{T}=-2\delta c_{\rm GW}+(\delta c_{\rm
GW})^{2}$, the latter term being negligible in the cases that we will
consider. $\frac{{D_{L}}}{c}\delta c_{\rm GW}(f)$ then denotes the frequency
dependent time delay that is incurred if $c_{\rm GW}$ is less than $c$. As the
frequency of radiation emitted by the source increases, waves emitted at
different times (and hence different frequencies) will have a different
associated $\delta c_{\rm GW}(f)$ (and hence travel times). This is what leads
to the stretching/squeezing of the observed signal, where importantly the
distance to the source affects how pronounced this effect is. We note that
this squeezing and stretching affect is the same as the frequency-domain
dephasing introduced in Mirshekari:2011yq , albeit expressed in a different
formalism. In the case of a large distortion in the observed waveform, this
can lead to effects analogous to the ‘inverse chirping’ of scalar waves in
scalar-tensor theories Sperhake:2017itk ; Aurrekoetxea:2022ika . For distances
of 100Mpc a value of only $\delta c_{\rm GW}(f)=1\times 10^{-16}$ is required
to add a shift of 1s to a signal. While 1 second may seem like a small amount,
shifting a signal by 1 second in the tens of seconds it takes a signal to move
through the LVK observing band is a significant shift, as we will explore in
more detail in later sections.

###  III.2 Arrival times

A second case of interest is when there is no detectable frequency-dependence
in the LVK or LISA bands themselves, yet there is a sharp transition in
between the two bands. In that case systems that are first observable in LISA
and then enter the LVK band are an excellent probe of the transition. The
primary signature here is not in the form of the signal itself, but instead
even a small $\delta c_{\rm GW}$ will result in an ‘arrival time’ of the
signal in the LVK band significantly different from what would be predicted
from LISA observations of the signal itself. The idea of multi-band
observations—using information from both LISA and the LVK band to measure the
properties of compact binary mergers—was first demonstrated in Sesana:2016ljz
. In that work it was predicted that LISA observations could measure the
arrival time of the signal in the LVK band to an accuracy of tens of seconds.
This finding was later supported by more detailed followup work Marsat:2020rtl
; Toubiana:2020cqv . However, in an independent study Klein:2022rbf , a more
pessismistic prediction was provided. There it is found that for stellar-mass
black systems, like GW150914, one can predict the arrival time in the LVK band
with a precision of $\sim$ hours from observations in the LISA band. This
discrepancy has not yet been resolved in the literature999We thank the authors
of Toubiana:2020cqv ; Baker:2022eiz  for related discussions.. Nevertheless,
even an uncertainty on the arrival time of the signal in the LVK band
constrained to hours can in principle be used to place very strong constraints
on $\delta c_{\rm GW}$ in the case of multi-band observations, as we explore
in section VI.

###  III.3 Intrinsic source effects

We will explore the effect of the waveform stretching/squeezing and arrival
time further in the following sections, but first we will briefly state that
the emission of gravitational waves is affected when a non-trivial $c_{\rm
GW}$ is present even when one is very near to the source. The full derivation
can be found in appendix C, but we show that the intrinsic evolution of the
source, in terms of the angular frequency ($\omega$), is changed in the
presence of a non-constant $c_{\rm GW}$ according to

$\displaystyle\dot{\omega}(t)=\left(\frac{G_{\rm GW}c}{G_{N}c_{\rm
GW}}\right)\frac{3456^{1/3}}{5}\left(\frac{G_{N}M_{c}}{c^{3}}\right)^{5/3}\omega^{11/3}.$
(13)

Here $G_{\rm GW}$ is the effective gravitational ‘constant’ $G$ seen by tensor
perturbations in the metric and $G_{N}$ is the ‘standard’ gravitational
constant as e.g. entering the Poisson equation or Kepler’s laws101010In other
words, $G_{N}$ is the gravitational constant relevant for the interaction
between (test) masses, so e.g. for the computation of planetary orbits or
torsion balance tests, while (in the context of the Horndeski theories
considered here) gravitational waves feel a different effective gravitational
constant $G_{\rm GW}$ due to the presence of the extra scalar degree of
freedom. For Horndeski theories these satisfy the relationship $G_{N}/G_{\rm
GW}=c_{\rm GW}^{2}/c^{2}$ Jimenez:2015bwa  – see appendix C for details.. This
is the standard general relativity solution, computed to leading order (“0
Post-Newtonian order”), with an additional term containing $c_{\rm GW}$ and
$G_{\rm GW}$. However, as shown in appendix C, this term/correction will be
negligible in comparison to the arrival time delays (in the region we are
exploring), so we will neglect this effect in what follows.

##  IV Probing $c_{\rm GW}$ in the LVK band

The joint observation of GW170817 and GRB 170817A allowed for a simple, but
very powerful, test of the speed of gravity. The peak of gravitational-wave
emission and the gamma-ray burst arrived at the Earth within seconds of each
other. If one assumes that these were emitted from the source at roughly the
same time, then any deviation between their arrival times at the Earth would
be due to differences in the speed of the two messengers. In the absence of
any observable delay, a very strong bound on $c_{\rm GW}$ can be placed.
Indeed, this observation constrained $\delta c_{\rm GW}$ in the 10 - 1000Hz
frequency band to $\lesssim 3\times 10^{-15}$ LIGOScientific:2017zic .

This limit was placed assuming a _constant_ $c_{\rm GW}$. In the scalar-tensor
theories we consider here we assume that $\delta c_{\rm GW}$ is allowed to be
nonzero, but it must be asymptoting to 0 in the $>10\mathrm{Hz}$ frequency
band in which GW170817 was observed by LIGO and Virgo. We will demonstrate
that a much tighter constraint on $\delta c_{\rm GW}$ across the LVK sensitive
band can be placed whenever there is a frequency-dependent $c_{\rm GW}$. To do
this we will first formulate how we will model the gravitational-wave signal
that ground-based observatories would see assuming a variable $c_{\rm GW}$, we
will then define how we will model our specific example, and then use these
models to identify what range of $c_{\rm GW}$ is consistent with the data
observed around GW170817.

###  IV.1 Including waveform stretching/squeezing in the gravitational wave
signal

We first need to define how we will apply a frequency dependent $\delta c_{\rm
GW}$ to the waveforms that we will be using. As discussed in section III.3, we
neglect any intrinsic source evolution effects. However, the variable travel
time delay as the waveform evolves from 10Hz towards 1000Hz and merger is not
negligible, and we must incorporate this effect into our waveform models.

To do this we begin by using the “TaylorF2" post-Newtonian waveform
Buonanno:2009zt  as implemented in lalsuite lalsuite , and apply a frequency-
dependent time shift to it. We first identify time before merger ($\tau$) as a
function of frequency using the leading order general relativity approximation

$\displaystyle
f=\left(\frac{256}{5c^{3}}M_{c}^{5/3}\pi^{8/3}\tau\right)^{-3/8}.$ (14)

We then identify the time delay as a function of frequency according to

$\displaystyle\delta\tau(f)=\frac{{D_{L}}}{c}\delta c_{\rm GW}(f).$ (15)

From these two expressions we can get $\delta\tau(\tau)$—the arrival time
delay as a function of time before merger—which can be used to produce the
timeshifted waveform we might expect to observe under these assumptions. We
apply the frequency dependent time shift by Fourier transforming the
frequency-domain waveform into the time-domain, applying the time shift to all
sample points, and interpolating the waveform back to an equally spaced
timeseries for use in analysis.

We note that our interpolation technique is computationally expensive. The
technique of adding this squeezing and stretching in terms of phase offset in
the frequency domain Mirshekari:2011yq  is more efficient, but would need to
be generalized from the power-law case considered there (although we do use a
power-law approximation in this section, we use the form of equation 10 in
later sections). We also only consider the dominant mode of gravitational-wave
emission in our studies, we refer the reader to Ezquiaga:2022nak  for a
discussion on how higher-order modes could be included.

###  IV.2 Parameterizing $\delta c_{\rm GW}(f)$

We next need a functional form of $\delta c_{\rm GW}(f)$. We use the formalism
we discussed in section II.3 in the limit where $\delta c_{\rm GW}(f)$ is
asymptoting to 0, i.e. in what we previously called the $k\gg k_{\star}$
limit. We model this asymptotic behaviour with two parameters—$\delta c_{\rm
GW}(f_{\mathrm{ref}})$ and $\sigma$—according to

$\displaystyle\delta c_{\rm GW}(f)=\delta c_{\rm
GW}(f_{\mathrm{ref}})\times\frac{f_{\mathrm{ref}}^{2\sigma}}{f^{2\sigma}}.$
(16)

Here e.g. $\delta c_{\rm GW}(30)$ corresponds to the value of $\delta c_{\rm
GW}$ at $30\;{\rm Hz}$, and we assume that $\delta c_{\rm GW}$ is approaching
zero over the band of observation following a power law with slope $2\sigma$.
$f_{\mathrm{ref}}$ is some reference frequency, which we stress is distinct
from $f_{\star}$. Note that this emergent power law behaviour neatly links up
with power-law parametrisations for (modifications of) the dispersion relation
previously employed in the literature Mirshekari:2011yq ;
LIGOScientific:2018dkp ; Baker:2022rhh . However, it is worth emphasizing that
in our context this behaviour only emerges asymptotically and it will
eventually be important to keep track of the full functional form of (10) when
inferring constraints on $\delta c_{\rm GW}(f)$ throughout the entire range of
frequencies considered here. We will do so in the following sections.

###  IV.3 Approximate measurement of observable $c_{\rm GW}$ deviation

Figure 2:  (Top) A prediction of how well we would be able to constrain
$\delta c_{\rm GW}$ at 30Hz, 100Hz and 500Hz as a function of the power law
slope parameter $\sigma$. This power law models how quickly $\delta c_{\rm
GW}(f)$ asymptotes to 0 as the frequency increases. This prediction assumes a
signal with the same masses (1.36$M_{\odot}$ for both bodies) as GW170817, and
the same signal-to-noise ratio (32.4). The Advanced LIGO noise curve and the
Einstein Telescope predicted noise curve lalsuite  are used while computing
these predictions. (Bottom) Bayesian inference exclusion posterior points,
projected down to $\sigma$ and $\delta c_{\rm GW}(30)$. An animated version of
this panel can be found at THIS LINK

We then want to quantify if such a deviation in the waveform is observable
with respect to a waveform where always $c_{\rm GW}=c$. A way to approximate
what values of $\delta c_{\rm GW}(f_{\mathrm{ref}})$ and $\sigma$ are required
to produce a distinguishable effect is to use the "distinguishability
criteria" from Lindblom:2008cm . This states that two waveforms are considered
to be distinguishable if, after subtracting one from the other, the residual
has a signal-to-noise ratio that is larger than one. That is if

$\left\langle h_{2}-h_{1}|h_{2}-h_{1}\right\rangle\geq 1$ (17)

then two waveforms could be distinguished from each other. This uses the
definition of the noise-weighted inner product between two waveforms $h_{1}$
and $h_{2}$, which is commonly used in gravitational-wave data analysis

$\left\langle
h_{1}|h_{2}\right\rangle=4\Re\int_{0}^{\infty}\frac{h^{*}_{1}(f)h_{2}(f)}{S_{n}(f)}df.$
(18)

Then, to find the level of deviation which is observable, we first choose a
fiducial physical system and generate the waveform assuming GR. We then
generate the same waveform, but assuming some parameterized deviation from GR
and measure the "distinguishability criterion". We then vary the values of the
GR-deviation parameters until we get a value of 1. We note that this measure
can provide overly optimistic predictions in the case where there are
correlations between parameters, for example it might be possible to change
other physical parameters to mimic the waveform stretching/squeezing. We will
check the validity of this approximation in the following subsection.

In the top panel of Figure 2, we demonstrate this for a signal with the same
masses and signal-to-noise ratio as GW170817, and using the Advanced LIGO
noise curve. This allows us to predict that such a signal would allow us to
constrain $\delta c_{\rm GW}$ at 30Hz to be no larger than $\sim 10^{-19}$,
for a slope with power law $\sigma=1$. For larger values of $n$ the exclusion
quickly becomes much tighter. We also predict the exclusion of $\delta c_{\rm
GW}$ at 100Hz and 500Hz. Exclusions at these frequencies are considerably
tighter, as would be expected when we model $\delta c_{\rm GW}$ as asymptoting
to 1. Finally, we predict the exclusion that would be possible if GW170817
were observed with the future Einstein Telescope observatory, these exclusions
are around two orders of magnitude tighter at $\sigma=1$ and decay
considerably quicker, reflecting the improved sensitivity of Einstein
Telescope at lower frequencies. This figure is generated assuming a subluminal
$c_{\rm GW}$, however we find that the results are identical to within
numerical precision when using a superluminal $c_{\rm GW}$111111In the next
section we do find slightly different behaviour between subluminal and
superluminal $c_{\rm GW}$. This is because a partial degeneracy between source
masses and the deviation parameters becomes important, and at some points we
see divergence between the two cases..

For full reproducibility of this, and other figures in this work, the code
used to produce this figure can be found by clicking on this link.

###  IV.4 Bayesian parameter estimation constraints of $c_{\rm GW}$ using
GW170817

We next consider what constraints can be placed using the recorded strain data
around the GW170817 observation when a frequency-varying $c_{\rm GW}$ is
assumed. This will allow us to verify the predicted curve in Figure 2 and
check to see if there are correlations between the parameters parameterizing
$\delta c_{\rm GW}$ and other physical parameters of the system. We use
Bayesian inference to place bounds on the $\sigma$ and $\delta c_{\rm GW}(30)$
parameters. We take the publicly released data from GWOSC from around GW170817
LIGOScientific:2019lzm  and follow closely the analysis described in
Nitz:2020oeq . We take the known sky location from the optical counterpart to
GW170817 LIGOScientific:2017ync  and assume a luminosity distance of 40Mpc
Cantiello:2018ffy  when applying the frequency-dependent time shift. We do not
consider spins in the analysis as neutron star spins are expected to be small
Brown:2012qf , which allows us to reduce computational cost and we also do not
include tidal terms in the waveform.

We choose a prior on $\sigma$ that is uniformly distributed between 0.5 and 5,
and a uniform in log prior on $\delta c_{\rm GW}(30)$ between $1\times
10^{-13}$ and $1\times 10^{-20}$. Other prior ranges match those of
Nitz:2020oeq . We use the emcee parallel tempered sampler Foreman-
Mackey:2012any  within PyCBC inference Biwer:2018osg  to generate results. We
run the analysis twice, once allowing for only subluminal values of $c_{\rm
GW}$ and once allowing for only superluminal values of $c_{\rm GW}$. The
results are then combined for visualization. For full reproducability, the
configuration files and code versions can be found in our data release page.

We can see the results of the Bayesian analysis in the bottom panel of Figure
2. We observe that the exclusion plots strongly depend on the frequency
dependence of $c_{\rm GW}$, i.e., on the power law slope $\sigma$ in Equation
11. The vast majority of samples we obtain have a $\delta c_{\rm GW}$ at 30Hz
smaller than $1\times 10^{-16.5}$, although there is one curious spike where
$\sigma\sim 1.5$, where a shift in the masses of the source, coupled with this
particular frequency-dependent time shift, seems to reproduce the GW170817
waveform well. This occurs for both superluminal and subluminal $c_{\rm GW}$,
although the shifted source masses differ in the two cases. We also provide an
animated version of this plot, where we show how the constraints improve as a
function of frequency. We find that all samples have $\delta c_{\rm GW}$
smaller than $1\times 10^{-17}$ at 100Hz, and all samples are smaller than
$1\times 10^{-18}$ at 500Hz. We can also compare the results between the two
panels of Figure 2 to validate the accuracy of the distinguishability
criterion used in the top panel. While we do notice that the top panel is more
optimistic than the bottom panel, it does agree to within 1-2 orders of
magnitude over the range of parameters considered. We therefore consider the
distinguishability criterion sufficient for making qualititative predictions
and we will use this again in future sections.

These constraints are considerably stronger than the limit of $\lesssim
3\times 10^{-15}$ placed with GW170817 LIGOScientific:2017zic . This can be
understood by remembering that this limit was based off of a time delay of
$\sim 5$ seconds. However, a shift of even milliseconds in a gravitational-
wave signal can be observable with matched-filtering and so we are able to
probe much smaller values of $\delta c_{\rm GW}$.

The LVK collaboration papers have also considered a variable gravitational-
wave speed in terms of a “modified dispersion relation” LIGOScientific:2017bnn
; LIGOScientific:2018dkp ; LIGOScientific:2019fpa ; LIGOScientific:2020tif ;
LIGOScientific:2021sio . Bounds in these papers have not been directly quoted
in terms of $\delta c_{\rm GW}$. However, it is possible to translate between
the constraining power on, for example, the graviton mass quoted in these
papers and infer the implied constraining power for $\delta c_{\rm
GW}$121212We thank Nathan Johnson-McDaniel for pointing this out to us and for
discussions about how the LVK results compare to our Figure 2.. The LVK
analyses only compute bounds at specific values of $\sigma$, and do not
consider $\sigma$ values larger than 1. However, the results from analysing
GW170817 data are compatible with our results, in regions where our analyses
overlap LIGOScientific:2018dkp . Tighter bounds on the graviton mass have been
placed on the LVK catalog using 43 confident binary black hole mergers in
LIGOScientific:2021sio . If this data were used to place bounds on $c_{\rm
GW}$ it would result in bounds that are around two orders of magnitude tighter
than those quoted here. This further backs up the main point of this section,
which is that ground based observations of gravitational-wave mergers can
place, and indeed have placed, very tight constraints on $c_{\rm GW}$ in the
$\sim 10-1000$Hz region.

##  V Probing $c_{\rm GW}$ in the LISA band

Figure 3:  The values of $\delta c_{\rm GW}^{(0)}$ that would be
distinguishable from a purely GR signal, as a function of the transition
frequency ($f_{\star}$) and the steepness of the transition ($\alpha$). (Left)
This is computed assuming a single LISA observation of a system of two
supermassive black holes each having a mass of $4.154\times 10^{6}M_{\odot}$
at a distance of 1Gpc. (Middle) Computed assuming the combination of an LVK
observation of a GW170817-like source and a single LISA observation of two
supermassive black holes. (Right) Computed assuming the combination of an ET
observation of a GW170817-like source and a single LISA observation of two
supermassive black holes. In all cases contour lines correspond to $\delta
c_{\rm GW}^{(0)}=10^{-3},10^{-6},10^{-9},10^{-12}$ and $10^{-15}$.

In the previous section we have demonstrated that we can place much tighter
constraints upon deviations from general relativity in the LVK band if we
assume a model where the speed of gravity diverges from the speed of light,
but varies as a function of frequency. Even a very small variation as a
function of frequency can cause a measurable frequency-dependent delay in
arrival time in the observed gravitational-wave signal. Similar tests would
also be possible for observations in the 1mHz - 100mHz with future space-based
observatories, such as LISA LISA:2017pwj  and TianQin TianQin:2015yph . Here
we focus on the LISA observatory. Perhaps the most useful astrophysical source
of gravitational waves for probing deviations from general relativity with
LISA is the merger of two supermassive black holes. Such systems would have
very large signal-to-noise ratio, would have relatively small mass
ratios—where reliability of waveform models is highest—and would cover a broad
range of frequencies in the scale of O(few years) that LISA might observe such
a source Klein:2015hvg ; Katz:2018dgn .

Within the LISA context, bounds on $\delta c_{\rm GW}$ have previously been
considered in two different settings. First, Littenberg:2019mob  investigated
testing a frequency-independent $\delta c_{\rm GW}$ in the LISA band with a
multi-messenger observation akin to GW170817. They forecasted a resulting
bound $|\delta c_{\rm GW}|\lesssim 10^{-12}$ in the event of a non-detection.
As before, here we instead focus on candidate signatures of a frequency-
dependent $\delta c_{\rm GW}$ and (given our findings in the LVK band above)
expect to place significantly stronger bounds on such a signal. Secondly,
while this paper was being completed Baker:2022rhh  also investigated a
frequency-dependent $c_{\rm GW}$ in the LISA band, similarly motivated by the
observation of deRham:2018red  that a frequency-dependent $c_{\rm GW}$ can
lead to interesting phenomenology close to the LISA band and exploring its
effect on waveform models and on (forecasted) parameter constraints. In terms
of $\delta c_{\rm GW}$ constraints, their results are driven by redshift-
induced frequency dependence (e.g. that the frequency of a monochromatic wave
differs at the source and at the observer in an expanding Universe). This
effect also is a generic consequence of a frequency-dependent $c_{\rm GW}$ and
Baker:2022rhh  uses this to place bounds on $|\delta c_{\rm GW}|\lesssim
10^{-4}$. Here we show that these bounds can be improved by several orders of
magnitude, more specifically to $|\delta c_{\rm GW}|\lesssim 10^{-17}$, by
considering the effects of $c_{\rm GW}\neq c$ outlined in section III.

To demonstrate the potential capability of LISA to observe such effects we
again use the distinguishability criterion introduced in the previous section.
Here we consider a system with a pair of black holes both having a mass of
$4.154\times 10^{6}M_{\odot}$ (the mass of the black hole at the centre of the
Milky Way). We situate the system at a distance of 1 Gpc, which we take to be
at a redshift of 0.20. We assume that the sky position and orientation of the
system is such that it will induce the optimal response for the LISA detector
and we assume that we observe the system from an initial gravitational wave
frequency of $5\times 10^{-5}$Hz, corresponding to roughly 1.5 months before
merger, the period in which the majority of the signal power would be accrued.
We model the LISA sensitivity curve using the LISA Sensitivity
code131313https://github.com/eXtremeGravityInstitute/LISA_Sensitivity. Such a
system would be recovered with a signal-to-noise ratio of around 5000 with
LISA. It is possible that the transition frequency would occur within the LISA
band of observation, so rather than modelling $\delta c_{\rm GW}$ as an
asymptoting function, we use the 3-parameter model defined earlier in equation
10. We then, as a function of the transition frequency ($f_{\star}$) and the
steepness of the transition ($\alpha$), predict the values of $\delta c_{\rm
GW}^{(0)}$ that would be distinguishable from a signal modelled assuming GR.

The result of this can be seen in Figure 3. We can see that in the absence of
a detection of a deviation from GR with such an observation, it would be
possible to place very tight bounds over a large range of the parameter space
that we consider. However, if $c_{\rm GW}$ has a steep transition with a
transition frequency above $0.01$ to $0.1$Hz it would produce a signal in LISA
that looks almost identical to a purely GR signal. It is important to remember
that the detectability criterion we apply here was demonstrated to be
optimistic by 1 to 2 orders of magnitude in the LVK band, and we might expect
a similar behaviour here as well. Nevertheless, this will not affect the main
features of our results shown in Figure 3. In this context also note that the
bounds forecasted here are for a single source and, from the analogous LVK
band discussion in the previous section, one may expect these bounds to
improve in a similar manner when considering a catalog of LISA observations.

##  VI Probing $c_{\rm GW}$ with joint LVK/LISA observations

It is interesting to compare how the bounds placed from simulated observations
in the LVK band would combine with analogous bounds in the LISA band. If we
re-compute our approximate bounds on $c_{\rm GW}$ used in Figure 2 in terms of
the 3 parameters used in the previous section, it is easy to combine these
plots and produce a visualization of the area of parameter space that would be
excluded from observations in both bands. This can be seen in Figure 3. We can
see that this combined exclusion includes much tighter constraints at larger
values of the transition frequency, and would place very strong constraints
for $\sigma<2$ over the full range of transition frequencies considered.
However, it would still permit cases where the transition frequency is between
the two bounds (around 0.1 to 1 Hz) with a very steep transition curve. For
completeness, we also show in Figure 3 how the constraints would improve with
the observation of a GW170817-like source with the Einstein Telescope.
Finally, we also include in Appendix D the uncombined exclusions, with only an
LVK or ET observation of a GW170817-like source. Indeed, there we see that LVK
band observations by themselves already strongly constrain $\delta c_{\rm GW}$
throughout the LISA band for $\sigma=1$.

In addition to observations made with only ground-based, or only with space-
based observatories, it will be possible, in the coming decades, to observe a
gravitational-wave source _both_ with ground-based and space-based
observatories. This was first explored in Sesana:2016ljz  where it was
demonstrated that a GW150914-like merger would be observable to LISA in the
years before it merged in the ground-based observation window. Such _multi-
band_ observations offer the ability to perform particularly stringent tests
of general relativity from observations of sources over a wide range of
frequencies. In particular, consider the case that $c_{\rm GW}$ asymptotes to
1 before entering the ground-based observation window, and $c_{\rm GW}$
asymptotes to a constant $c_{\rm GW}(0)\neq c$ in the LISA band of observation
(at levels of accuracy sufficient to evade the detection of any frequency-
dependence in the LVK and LISA bands themselves – see previous sections). Such
a deviation would not be directly observable in either band, but in the gap
between the signal leaving the LISA band and being observed in the LVK band
there would be a shift in the travel time required for the signal to reach an
observer. Relating this to our model fit of equation 10 this would require a
large value of $\sigma$ with $f_{\star}$ in the frequency gap between the LISA
and LVK bands, which is the region that is poorly constrained in Figure 3.

If we consider GW150914 at a distance of 400Mpc, the travel time, in the
Earth’s reference frame, for the signal to reach the Earth is $4\times
10^{16}$ seconds. Then, as with the previous cases, we can see that even very
tiny changes in $c_{\rm GW}$ will produce large variations in the arrival
time. To produce time shifts of $\sim 10000$ seconds—roughly the accuracy at
which an observation in the LISA band would be able to constrain the merger
time as predicted in Klein:2022rbf 141414We remind the reader that there are
inconsistent predictions of the ability to measure the arrival time of such a
signal in the LVK band with LISA, as we discussed in section III.2—would
require $|\delta c_{\rm GW}|=3\times 10^{-13}$. To produce time shifts of
$\sim 10$ seconds—as predicted in Sesana:2016ljz ; Marsat:2020rtl ;
Toubiana:2020cqv —would require $|\delta c_{\rm GW}|=3\times 10^{-16}$. More
generally speaking, if in the future the merger time in the LVK band can be
predicted with an accuracy of $\sim 10^{x}$ seconds from LISA band
observations, $\delta c_{\rm GW}\sim 10^{-14+x-y}$ will be discernible for a
source at $10^{y}$ Mpc.

An interesting observation can be made when considering larger values of
$|\delta c_{\rm GW}|$. In particular the time delay due to $|\delta c_{\rm
GW}|$ changing between the LISA and LVK bands can quickly become very large. A
value of only $|\delta c_{\rm GW}|=8\times 10^{-10}$ is required before the
time delay is larger than a year. In such cases, with a subluminal $c_{\rm
GW}$, we might observe the signal with LISA at _the same time_ as the signal
arriving at ground based observatories. For values of $|\delta c_{\rm GW}|$
that are much larger than this, we would not observe a multi-band signal at
all, and indeed the absence of multi-band signals, when such observations are
expected, could be an indicator of $|\delta c_{\rm GW}|$ being large in the
LISA band. With existing upper bounds on $c_{\rm GW}$ relevant for the LISA
band imposing $|\delta c_{\rm GW}|\lesssim 10^{-2}$ Jimenez:2015bwa , this
means there is in principle a window $|\delta c_{\rm GW}|\sim 10^{-8}-10^{-2}$
for which, given a sufficiently fast transition in between the LVK and LISA
bands, the absence of multi-band signals would be the main signature in near-
future measurements. This would also be an especially interesting target for
proposed future experiments targeting this intermediate region, such as AEDGE
AEDGE:2019nxb . Finally note that, while this draft was in internal review, we
were made aware of ongoing work Baker:2022eiz  that will also explore similar
multiband constraints.

##  VII Conclusions

Summary: In this paper we have explored how frequency-dependent transitions in
the speed of gravitational waves $c_{\rm GW}$ – a generic consequence in large
classes of dark energy theories – can be constrained by current and future
observations in the LIGO/Virgo/KAGRA (LVK) band, in the LISA band and with
joint observations in both bands. Since such dark energy-related transitions
naturally occur close to those two bands, the observations place powerful
constraints on $c_{\rm GW}$. We would like to highlight the following key
results:

  * •

We find that deviations away from $c_{\rm GW}=c$ can be constrained down to a
level of $|\delta c_{\rm GW}|\sim 10^{-17}$ in both the LVK and LISA bands
even for mild frequency-dependence, much stronger than existing bounds for
frequency-independent $c_{\rm GW}$ from GW170817 or indeed analogous
forecasted bounds for multi-messenger observations in the LISA band
Littenberg:2019mob . Our constraints are driven by $c_{\rm GW}$-induced
frequency-dependent time shifts in the observed signals and we have discussed
in detail how they depend on the precise frequency(-ies) associated with the
transition as well as on its functional form. Fig. 3 summarises the main
constraints we find.

  * •

We have identified a class of interesting very fast transitions taking place
around frequencies of $\sim 10^{-2}-10\;{\rm Hz}$, which can almost completely
evade the aforementioned bounds when they proceed sufficiently quickly151515In
terms of a power law scaling, this e.g. means $c_{\rm GW}=c$ is approached at
large frequencies as $f^{-2\sigma}$, where $\sigma\gtrsim 5$.. However, joint
observations of sources visible in both LVK and LISA bands would be able to
either: 1) constrain $|\delta c_{\rm GW}|\lesssim 10^{-15}$, when LVK
observations ‘see’ sources at the times expected from prior LISA observations
of the same sources within expected uncertainties (note that this constraint
can weaken to $|\delta c_{\rm GW}|\lesssim 10^{-12}$ for more pessimistic
forecasted uncertainties – see section VI for details), 2) measure $|\delta
c_{\rm GW}|$ anywhere in the range $|\delta c_{\rm GW}|\sim 10^{-15}-10^{-9}$,
or 3) indicate that $|\delta c_{\rm GW}|\sim 10^{-8}-10^{-2}$, when LVK
observations do not ‘see’ any expected counterparts for sources expected from
prior LISA observations within ${\cal O}$(years).

  * •

As a precursor for deriving these constraints, we have constructed a number of
theoretically well-motivated templates for a fiducial frequency-dependent
$c_{\rm GW}$. Underlying most of our results is a simple three-parameter
ansatz that captures the essential features of candidate transitions (see
section II), but we collect a discussion of how to construct more
sophisticated templates for future analyses in appendix B.

Implications for dark energy: In the immediate aftermath of GW170817 the
consequences of tight bounds on $c_{\rm GW}$ – $|\delta c_{\rm GW}|\lesssim
10^{-15}$ – on dark energy theories were explored in detail (see e.g.
Baker:2017hug ; Creminelli:2017sry ; Sakstein:2017xjx ; Ezquiaga:2017ekz  and
references therein), assuming these bounds can simply be ported to much lower
cosmological frequencies. The resulting conclusions would be significantly
more robust in the presence of an analogous precision measurement at LISA
frequencies, as discussed and forecasted here. This is because we can
consistently describe both cosmological and LISA scales in dark energy
theories that lead to $c_{\rm GW}\neq c$ on cosmological scales, while this is
much more challenging for LVK scales deRham:2018red .

In the context of Horndeski gravity explored in this paper, a tight bound on
$\delta c_{\rm GW}$ from LISA would therefore firmly reduce the set of
surviving scalar-tensor dark energy theories to Baker:2017hug ;
Creminelli:2017sry ; Sakstein:2017xjx ; Ezquiaga:2017ekz  (also see
Amendola:2012ky ; Amendola:2014wma ; Deffayet:2010qz ; Linder:2014fna ;
Raveri:2014eea ; Saltas:2014dha ; Lombriser:2016yzn ; Lombriser:2015sxa ;
Jimenez:2015bwa ; Bettoni:2016mij ; Bettoni:2016mij ; Sawicki:2016klv  for
closely related prior work)

$\displaystyle{\cal L}=G_{2}(\phi,X)+G_{3}(\phi,X)\Box\phi+G_{4}(\phi)R,$ (19)

where the $G_{i}$ are free functions of the scalar $\phi$ and its first
derivative via $X\equiv-\tfrac{1}{2}\nabla_{\mu}\phi\nabla^{\mu}\phi$.
Conversely, within the Horndeski context a detection of a non-zero $\delta
c_{\rm GW}$ would constitute a measurement of significant $G_{4}$ and/or
$G_{5}$ interactions contributing to $c_{\rm GW}$. Note that, by constraining
the templates discussed here and in the event of a detection of non-zero
$\delta c_{\rm GW}$, not only would we obtain information about $\delta c_{\rm
GW}$ itself and hence the value of $c_{\rm GW}$ on cosmological scales, but
also about the $f_{\star}$ and $\sigma$ parameters detailed in section II.
$f_{\star}$ encodes information about the energy/frequency-scale associated
with a UV completion for dark energy, while $\sigma$ encodes information about
the specific novel physics and interactions entailed by such a candidate UV
completion.

Implications for future gravitational wave searches: Throughout most of this
paper our search for and constraints on $c_{\rm GW}$ with separate or
combined, present and or forecasted LVK/LISA observations has been primarily
motivated by large classes of dark energy theories that generically lead to a
frequency-dependent $c_{\rm GW}$ in or close to those bands. What do our
results imply beyond this dark energy-specific context? Ultimately a non-zero
$\delta c_{\rm GW}$ would be evidence for novel gravitational interactions,
new physics, that provide a non-trivial medium for gravitational waves to
travel through. The template(s) we have worked with here do assume 1) an
asymptotically constant $c_{\rm GW}$ at very low frequencies and 2) a
frequency-dependent transition to $c_{\rm GW}=c$ at some frequency
$f_{\star}$, whose form we have motivated by general field-theoretic
arguments. Importantly, this means no further dark energy-specific input
affects the templates themselves, so these can be implemented in future
searches and offer a robust probe of any novel physics affecting $c_{\rm GW}$
in the LVK/LISA bands and (at least approximately) consistent with the above
two assumptions, whether dark energy-related or otherwise.

## Acknowledgments

We especially thank Nathan Johnson-McDaniel and Scott Melville for several
helpful discussions and shared insights. We also thank Tessa Baker, Anson
Chen, Claudia de Rham, Macarena Lagos, Eugene Lim, Mauro Pieroni, Anand
Sengupta and Gianmassimo Tasinato for useful discussions and comments on a
draft. JN is supported by an STFC Ernest Rutherford Fellowship (ST/S004572/1).
IH is supported by STFC grants ST/T000333/1 and ST/V005715/1. In deriving the
results of this paper, we have used: xAct xAct , PyCBC alex_nitz_2022_6583784
, LALSuite lalsuite , emcee Foreman-Mackey:2012any , numpy Harris:2020xlr  and
scipy Virtanen:2019joe .

## Data availability

All data used in this work, and the information necessary to fully reproduce
any results and figures presented here, can be found at our data release page
https://icg-gravwaves.github.io/probing_speed_of_gravity/.

## Open access

This article is licensed under a Creative Commons Attribution 4.0
International License, which permits use, sharing, adaptation, distribution
and reproduction in any medium or format, as long as you give appropriate
credit to the original author(s) and the source, provide a link to the
Creative Commons licence, and indicate if changes were made. The images or
other third party material in this article are included in the article’s
Creative Commons licence, unless indicated otherwise in a credit line to the
material. If material is not included in the article’s Creative Commons
licence and your intended use is not permitted by statutory regulation or
exceeds the permitted use, you will need to obtain permission directly from
the copyright holder. To view a copy of this licence, visit
http://creativecommons.org/licenses/by/4.0/.

##  Appendix A $c_{\rm GW}$ in Horndeski gravity

In section II of the main text we focused on a specific example theory, namely
(2). However, the key features relevant for this paper remain true for the
much more general class of Horndeski theories, the most general Lorentz-
invariant scalar-tensor action that gives rise to second order equations of
motion. It is described by the following action

$\displaystyle S_{H}=\int\mathrm{d}^{4}x\sqrt{-g}\left\\{\sum_{i=2}^{5}{\cal
L}_{i}[\phi,g_{\mu\nu}]\right\\},$ (20)

where we write the scalar-tensor Lagrangians ${\cal L}_{i}$ (for a scalar
$\phi$ and a massless tensor $g_{\mu\nu}$) as

$\displaystyle{\cal L}_{2}$
$\displaystyle=\Lambda_{2}^{4}\,G_{2}~{},\quad\quad\quad\quad{\cal
L}_{3}=\frac{\Lambda_{2}^{4}}{\Lambda_{3}^{3}}G_{3}\cdot[\Phi]\,,$
$\displaystyle{\cal L}_{4}$
$\displaystyle=\frac{\Lambda_{2}^{8}}{\Lambda_{3}^{6}}G_{4}R+\frac{\Lambda_{2}^{4}}{\Lambda_{3}^{6}}~{}G_{4,X}\left([\Phi]^{2}-[\Phi^{2}]\right)\,,$
(21) $\displaystyle{\cal L}_{5}$
$\displaystyle=\frac{\Lambda_{2}^{8}}{\Lambda_{3}^{9}}G_{5}G_{\mu\nu}\Phi^{\mu\nu}-\frac{1}{6}\frac{\Lambda_{2}^{4}}{\Lambda_{3}^{9}}G_{5,X}([\Phi]^{3}-3[\Phi][\Phi^{2}]+2[\Phi^{3}]).$

Here we adopt a dimensionless definition
$X\equiv-\tfrac{1}{2}\nabla_{\mu}\phi\nabla^{\mu}\phi/\Lambda_{2}^{4}$ for
what is essentially the scalar kinetic term, and we have defined
$\Phi^{\mu}_{\;\;\nu}\equiv\nabla^{\mu}\nabla_{\nu}\phi$. The $G_{i}$ are
dimensionless functions of $\phi/M_{\rm Pl}$ and $X$, and $G_{i,\phi}$ and
$G_{i,X}$ denote the partial derivatives of the $G_{i}$ (with respect to these
dimensionless arguments). Square brackets denote the trace, e.g.
$[\Phi^{2}]\equiv\nabla^{\mu}\nabla_{\nu}\phi\nabla^{\nu}\nabla_{\mu}\phi$,
and we have three mass scales: $M_{\rm Pl},\Lambda_{2}$ and $\Lambda$. In
cosmology they are conventionally taken to satisfy $\Lambda_{2}^{2}=M_{\rm
Pl}H_{0}$ and $\Lambda_{3}^{3}=M_{\rm Pl}H_{0}^{2}$, which ensures that all
interactions can give ${{\cal O}(1)}$ contributions to the (cosmological)
background evolution. General Relativity is recovered when
$G_{2}=G_{3}=G_{5}=0$ and $G_{4}=1$. Note that the example theory (2)
effectively amounts to a minimal choice of $G_{2}$ and $G_{4}$, while
$G_{5}=g(\phi)$ and $G_{3}=0$. For the general Horndeski theory(21), one can
then work out the effect on $c_{\rm GW}$ and find

$\displaystyle\alpha_{T}$
$\displaystyle=2\frac{X}{M^{2}}\left[2G_{4,X}-2G_{5,\phi}-\left(\frac{\ddot{\phi}}{H_{0}^{2}}-\frac{\dot{\phi}H}{H_{0}^{2}}\right)G_{5,X}\right]\,.$
(22)

Here
$M^{2}/2=G_{4}-2XG_{4,X}+XG_{5,\phi}-\tfrac{\dot{\phi}H}{H_{0}^{2}}XG_{5,X}$
is the so-called effective Planck mass (albeit dimensionless in the way
written here), $H$ is the Hubble scale that measures the expansion of the
Universe and satisfies $H=\dot{a}/a$ (where a is the scale factor of the
Universe), and $H_{0}$ is the value of the Hubble scale today. Overdots denote
time-derivatives with respect to proper time $t$ and we recall that, on a
cosmological background, the scalar $\phi$ is a function of time only, i.e.
$\phi=\phi(t)$. From this it is clear that the presence of the scalar can
affect $c_{\rm GW}$ through non-trivial $G_{4}$ and $G_{5}$ interactions.
Conversely, if $G_{5}$ vanishes and $G_{4}$ is at most a function of $\phi$
(but not of its derivatives), then $\alpha_{T}=0$ and there is no effect on
$c_{\rm GW}$ and no deviation of this speed from the speed of light.

##  Appendix B More $c_{\rm GW}$ templates

Figure 4:  The speed of gravitational waves, $c_{\rm GW}$, as a function of
frequency $f$ for the extended template (24). The different curves correspond
to different choices of the ‘steepness’ parameters $\sigma$ and $\beta$ that
control the asymptotic power law scalings of $c_{\rm GW}$ – see (25). We have
set $\gamma=6$ for this example. This choice does not affect the asymptotes,
but does determine how quickly the transition from a $\sigma$-dominated
evolution to a $\beta$-dominated one occurs.

The $c_{\rm GW}$ template (10) we discussed in the main text is controlled by
a single steepness parameter $\sigma$ that controls both the small and large
$k$ asymptotes. The simplicity of the template results in a symmetric
$k^{2\sigma}$ vs. $k^{-2\sigma}$ scaling for those asymptotes, but this can
straightforwardly be refined by considering higher order (in the number of
parameters) generalisations of this template, in particular ones that allow
for different power law scalings in the asymptotic small and large frequency
regimes. Here we will discuss how this can be done more explicitly.

Using the $\Theta_{\pm}$ functions defined in (9) we can define

$\displaystyle\hat{\Theta}(\sigma,\beta,\gamma)\equiv\Theta_{+}(\gamma)\Theta_{+}(\beta)+\Theta_{-}(\gamma)\Theta_{+}(\sigma),$
(23)

where it is understood that any $\Theta_{\pm}$ depends on $f$ and $f_{\star}$,
but we no longer write this explicitly. We can now again build a template for
$c_{\rm GW}$ as before, writing

$\displaystyle c_{\rm GW}(f)=c_{\rm GW}^{(0)}+(c-c_{\rm
GW}^{(0)})\hat{\Theta}.$ (24)

Note that we require $\gamma>\\{\sigma,\beta\\}$, but since it will not affect
the large or small frequency asymptotes the precise value of $\gamma$ is
immaterial161616Nevertheless note that it does affect the detailed transition
evolution around $f_{\star}$. Since those details are of lesser importance for
the phenomenology investigated throughout this paper, we leave a more precise
investigation of their effects to future work.. In this sense (24) is really a
four parameter template controlled by $\\{c_{\rm
GW}(0),f_{\star},\sigma,\beta\\}$. Investigating its asymptotic scalings
explicitly, we find

$\displaystyle k\ll k_{\star}:$ $\displaystyle c_{\rm GW}-c_{\rm GW}(0)$
$\displaystyle\propto k^{2\sigma},$ $\displaystyle k\gg k_{\star}:$
$\displaystyle c-c_{\rm GW}$ $\displaystyle\propto k^{-2\beta}.$ (25)

So $\sigma$ and $\beta$ control the relevant power law exponents in the large
and small $k$ regime, respectively, allowing different scalings in those two
regimes. As discussed in the main text, a sufficiently fast (polynomial or
otherwise) scaling as $c_{\rm GW}$ asymptotes back to unity is particularly
important for ensuring consistency with measurements in the LVK band. So
allowing for this parametric freedom for the large $k$ scaling is especially
relevant there, ultimately capturing one specific aspect of the lack of
knowledge we have about the precise nature of any would-be UV completion of
relevance for our analysis.

Finally, let us point out that the above template for asymmetric transitions
is of course not unique. One example of an alternative template (albeit with
less clean asymptotics) is a logistic/Verhulst function of the following type

$\displaystyle\Theta_{V}(\sigma,\beta)\equiv
1-\left(1+(2^{\beta}-1)\left(\frac{f}{f_{\star}}\right)^{\sigma}\right)^{-1/\beta},$
(26)

from which we can build the template

$\displaystyle c_{\rm GW}(f)=c_{\rm GW}^{(0)}+(c-c_{\rm GW}^{(0)})\Theta_{V}.$
(27)

Again we have introduced another parameter, $\beta$, that controls the
asymmetry in the transition, allowing asymmetric scalings for small and large
$k$ in the transition.

##  Appendix C Intrinsic source evolution with non-constant $c_{\rm GW}$

In this appendix we would like to understand how a non-constant $c_{\rm GW}$
affects GW emission and propagation in the source frame, in particular how the
intrinsic evolution of (and power emitted by) a binary system is affected if
$c_{\rm GW}$ deviates from its GR prediction. We will compute this at leading,
0-PN, order, to demonstrate that these effects exist, but are negligible in
our work. Higher order corrections would be needed if these effects were
important. Closely following Jimenez:2015bwa ; Maggiore:1900zz , we can write
the quadratic action for the standard metric tensor perturbations $\gamma$ as

$\displaystyle S_{h}=\frac{1}{64\pi G_{\rm GW}}\int
d^{4}x\left(\dot{\gamma}_{ij}\dot{\gamma}^{ij}-k^{2}c_{\rm
GW}^{2}{\gamma}_{ij}{\gamma}^{ij}\right),$ (28)

where $G_{\rm GW}$ denotes the effective gravitational ‘constant’ $G$ seen by
these tensor perturbations. We can then express the instantaneous power
emitted as

$\frac{dE}{dt}=\frac{r^{2}c^{4}}{32\pi c_{\rm GW}G_{\rm GW}}\int
d\Omega\left\langle\partial_{t}\gamma_{ij}\partial_{t}\gamma_{ij}\right\rangle\,,$
(29)

where we integrate over solid angle $\Omega$ and average over a region of
spacetime much larger than the GW wavelength (denoted by
$\langle\dots\rangle$). In order to use the above we furthermore solve for the
amplitude $\gamma_{ij}$ of the radiated GWs by performing a multipole
expansion and working to lowest (quadrupole) order in velocity. Solving the
linearised Einstein equation and proceeding along the usual Green’s function
solution (see e.g. Maggiore:1900zz ), we recover the expression
Jimenez:2015bwa

$[\gamma_{ij}]_{quad}\ =\ \frac{2G_{\rm
GW}}{rc^{4}}\ddot{Q}_{ij}^{TT}\\!\left(t-\frac{r}{c_{\rm GW}}\right),$ (30)

at leading order in velocity. Here we have defined the usual quadrupole moment
$Q^{ij}\equiv M^{ij}-\tfrac{1}{3}\delta^{ij}M^{k}_{k}$ in terms of the momenta
of $T^{00}/c^{2}$, $M$. Note that, in performing the above computation it is
important that the time it takes for a GW to traverse the source is negligible
compared to the time variation scale of $c_{\rm GW}$, which allows us to
effectively treat $c_{\rm GW}$ as constant in the relevant integrals.

Now we would like to work out the total power emitted, starting from (29). Re-
arranging, we obtain

$\displaystyle\frac{dP}{d\Omega}$ $\displaystyle=\frac{r^{2}c^{4}}{32\pi
c_{\rm GW}G_{\rm
GW}}\left\langle\partial_{t}\gamma_{ij}\partial_{t}\gamma_{ij}\right\rangle\
=\frac{G_{\rm GW}}{8\pi c_{\rm
GW}c^{4}}\left\langle\dddot{Q}_{ij}^{TT}\dddot{Q}_{ij}^{TT}\right\rangle\,.$
(31)

Here $P\equiv dE/dt$ and we have substituted for $\gamma_{ij}$ using (30).
Instead of writing this in terms of its transverse-traceless projection
$Q_{ij}^{TT}$, we would like to express this in terms of the quadrupole moment
of the source $Q_{ij}$ itself. Doing so and integrating over solid angle
again, we find the total (quadrupolar) power emitted (at a given time)

$\displaystyle P_{\rm quad}=\frac{G_{\rm GW}}{5c_{\rm
GW}c^{4}}\left\langle\dddot{Q}_{kl}\dddot{Q}_{mn}\right\rangle\,,$ (32)

again recovering the result from Jimenez:2015bwa . We emphasise that this
expression holds for a time/frequency-varying $c_{\rm GW}$, since it is the
total power emitted at a given time/frequency, but care needs to be taken when
integrating this expression to e.g. obtain the total energy emitted throughout
the evolution.

We now focus on a circular inspiral and explicitly work out the power emitted
in terms of the parameters describing this system. Establishing some notation,
we have

$\displaystyle\omega_{s}^{2}$ $\displaystyle=\frac{G_{N}m}{R^{3}},$
$\displaystyle\omega_{\rm GW}$ $\displaystyle=2\omega_{s},$ $\displaystyle
f_{\rm GW}$ $\displaystyle=\frac{\omega_{\rm GW}}{2\pi},$ $\displaystyle
M_{c}$ $\displaystyle=\mu^{3/5}m^{2/5},$ (33)

where we recall that

$\displaystyle\mu$ $\displaystyle\equiv\frac{m_{1}m_{2}}{m_{1}+m_{2}},$
$\displaystyle m$ $\displaystyle\equiv m_{1}+m_{2},$ (34)

and $m_{1}$ and $m_{2}$ are the masses of the two compact objects constituting
the inspiralling binary in question. Importantly, in (33), note that the
expression for $\omega_{s}$ is effectively Kepler’s law and so the appropriate
gravitational coupling constant here is $G_{N}$, i.e. the gravitational
constant that enters in the Poisson equation and which is generically
different from $G_{\rm GW}$. In this setup we can solve for $\gamma$ and find

$\displaystyle\gamma_{\times}(t,\theta,\phi)=\frac{4G_{\rm
GW}\mu\omega_{s}^{2}R^{2}}{rc^{3}}c_{\rm GW}\cos\theta\sin(2\omega_{s}t_{\rm
ret}+2\phi)$ (35)

where we note the factor of $G_{\rm GW}$ that enters following on from (30)
and we only show the $\times$ polarisation for brevity here – see (3.330) and
(3.331) in Maggiore:1900zz  for an analogous $+$ polarisation expression,
which can be adapted along the same lines. Note that $t_{\rm ret}$ has an
implicit $c_{\rm GW}$ dependence as well. Now using (33) to replace $R$ with
$\omega_{s}$ and explicitly introducing the chirp mass, we can write

$\displaystyle\gamma_{\times}(t,\theta,\phi)$ $\displaystyle=\frac{G_{\rm
GW}c}{G_{N}c_{\rm
GW}}\frac{4}{r}\left(\frac{G_{N}M_{c}}{c^{2}}\right)^{5/3}\left(\frac{\pi
f_{\rm GW}}{c}\right)^{2/3}$
$\displaystyle\cdot\cos\theta\sin(2\omega_{s}t_{\rm ret}+2\phi).$ (36)

The first fraction encodes all dark energy-related modifications we consider
here, while everything that remains is as in standard GR. Finally, we can
follow the same logic as above and compute the resulting total power emitted
at any given point in time

$\displaystyle P=\frac{G_{\rm GW}c}{G_{N}c_{\rm
GW}}\frac{32c^{5}}{5G_{N}}\left(\frac{G_{N}M_{c}\omega_{\rm
GW}}{2c^{3}}\right)^{\frac{10}{3}}.$ (37)

The power lost to gravitational radiation will be equal to loss of orbital
energy, and, assuming Keplerian orbits, this can expressed in terms of $M_{c}$
and $\omega$ as Maggiore:1900zz

$\displaystyle
P=-\frac{dE_{orbit}}{dt}=\frac{2}{3}\left(\frac{G_{N}^{2}M_{c}^{5/3}}{32}\right)^{1/3}\omega^{-1/3}\dot{\omega}.$
(38)

This can be rearranged in terms of $\dot{\omega}$ to give

$\displaystyle\dot{\omega}(t)=\left(\frac{G_{\rm GW}c}{G_{N}c_{\rm
GW}}\right)\frac{3456^{1/3}}{5}\left(\frac{G_{N}M_{c}}{c^{3}}\right)^{5/3}\omega^{11/3}.$
(39)

In the generic case we are considering where $c_{\rm GW}$ and $G_{\rm GW}$
depend on frequency we cannot simplify further without specifying that
dependency. In the simple case that $c_{\rm GW}=1$ and $G_{\rm GW}=G_{N}$ the
extra term on the right hand side vanishes and this collapses to the standard
0-PN GR solution. In the case where $c_{\rm GW}$ and $G_{\rm GW}$ deviate from
the GR values but are constant over the range of observation this can be
integrated analytically and rearranged to find $h_{+}(t)$ and $h_{\times}(t)$.
In this case we find that both the phasing evolution and the overall amplitude
depend on $c_{\rm GW}$ and $G_{\rm GW}$, but do so in a way that is degenerate
with chirp mass. Therefore a signal emitted in such a model would appear to be
consistent with general relativity, but the recovered measurement of the chirp
mass would be offset from the true value – also see Mirshekari:2011yq  for a
related discussion of phase shifts and degeneracies. In the case of a general
solution for $c_{\rm GW}(\omega)$ and $G_{\rm GW}(\omega)$ this could be
integrated numerically to find a model of the emitted waveform.

Figure 5:  The values of $\delta c_{\rm GW}^{(0)}$ that would be
distinguishable from a purely GR signal, as a function of the transition
frequency ($f_{\star}$) and the steepness of the transition ($\alpha$). (Left)
This is computed assuming a single LISA observation of a system of two
supermassive black holes each having a mass of $4.154\times 10^{6}M_{\odot}$
at a distance of 1Gpc. (Middle) Computed assuming an LVK observation of a
GW170817-like source. (Right) Computed assuming an ET observation of a
GW170817-like source. In all cases contour lines correspond to $\delta c_{\rm
GW}^{(0)}=10^{-3},10^{-6},10^{-9},10^{-12}$ and $10^{-15}$.

In the above we have kept $G_{\rm GW}$ and $G_{N}$, i.e. two different
gravitational coupling constants for gravitational waves and for matter,
generic. However, within the context of (21) and inside a screened regime as
is relevant for the emission of gravitational waves considered here, these are
found to satisfy $G_{N}/G_{\rm GW}=c_{\rm GW}^{2}/c^{2}$ Jimenez:2015bwa
171717This relationship may break down once approaching the cutoff of the
cosmological theory close to $\Lambda$, so care needs to be taken if this is
to be modelled in detail.. So here any modification of (39) away from its GR
limit will be controlled by a single, dimensionless parameter: $c_{\rm GW}/c$.
However, bounds obtained from this effect will be significantly weaker than
the ones discussed in the main text. Compare the $|\delta c_{\rm GW}|\lesssim
10^{-2}$ bound obtained by Jimenez:2015bwa  from the Hulse-Taylor binary using
analogous reasoning with the $|\delta c_{\rm GW}|\lesssim 10^{-17}$ and better
constraints discussed in the main text. We therefore do not pursue bounds from
this effect further here.

##  Appendix D Constraints using only LVK or ET observations

In Figure 3 we showed the potential constraints that could be placed on
$c_{\rm GW}$ when combining an observation of a super-massive black hole
merger with LISA with a single observation of a GW170817-like source with
LVK/ET. For completeness, in Figure 5, we also show what the constraints would
look like if _only_ considering a single GW170817-like source observed with
either LVK or ET, alongside the constraints from LISA alone as comparison.

## References

  * (1) Planck Collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641 (2020) A6, arXiv:1807.06209 [astro-ph.CO].
  * (2) 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 (2017) no. 25, 251301, arXiv:1710.06394 [astro-ph.CO].
  * (3) P. Creminelli and F. Vernizzi, “Dark Energy after GW170817 and GRB170817A,” Phys. Rev. Lett. 119 (2017) no. 25, 251302, arXiv:1710.05877 [astro-ph.CO].
  * (4) J. Sakstein and B. Jain, “Implications of the Neutron Star Merger GW170817 for Cosmological Scalar-Tensor Theories,” Phys. Rev. Lett. 119 (2017) no. 25, 251303, arXiv:1710.05893 [astro-ph.CO].
  * (5) J. M. Ezquiaga and M. Zumalacárregui, “Dark Energy After GW170817: Dead Ends and the Road Ahead,” Phys. Rev. Lett. 119 (2017) no. 25, 251304, arXiv:1710.05901 [astro-ph.CO].
  * (6) L. Amendola, M. Kunz, M. Motta, I. D. Saltas, and I. Sawicki, “Observables and unobservables in dark energy cosmologies,” Phys. Rev. D87 (2013) no. 2, 023501, arXiv:1210.0439 [astro-ph.CO].
  * (7) L. Amendola, G. Ballesteros, and V. Pettorino, “Effects of modified gravity on B-mode polarization,” Phys. Rev. D90 (2014) 043009, arXiv:1405.7004 [astro-ph.CO].
  * (8) C. Deffayet, O. Pujolas, I. Sawicki, and A. Vikman, “Imperfect Dark Energy from Kinetic Gravity Braiding,” JCAP 1010 (2010) 026, arXiv:1008.0048 [hep-th].
  * (9) E. V. Linder, “Are Scalar and Tensor Deviations Related in Modified Gravity?,” Phys. Rev. D90 (2014) no. 8, 083536, arXiv:1407.8184 [astro-ph.CO].
  * (10) M. Raveri, C. Baccigalupi, A. Silvestri, and S.-Y. Zhou, “Measuring the speed of cosmological gravitational waves,” Phys. Rev. D91 (2015) no. 6, 061501, arXiv:1405.7974 [astro-ph.CO].
  * (11) I. D. Saltas, I. Sawicki, L. Amendola, and M. Kunz, “Anisotropic Stress as a Signature of Nonstandard Propagation of Gravitational Waves,” Phys. Rev. Lett. 113 (2014) no. 19, 191101, arXiv:1406.7139 [astro-ph.CO].
  * (12) L. Lombriser and A. Taylor, “Breaking a Dark Degeneracy with Gravitational Waves,” JCAP 1603 (2016) no. 03, 031, arXiv:1509.08458 [astro-ph.CO].
  * (13) L. Lombriser and N. A. Lima, “Challenges to Self-Acceleration in Modified Gravity from Gravitational Waves and Large-Scale Structure,” Phys. Lett. B765 (2017) 382–385, arXiv:1602.07670 [astro-ph.CO].
  * (14) J. Beltran Jimenez, F. Piazza, and H. Velten, “Evading the Vainshtein Mechanism with Anomalous Gravitational Wave Speed: Constraints on Modified Gravity from Binary Pulsars,” Phys. Rev. Lett. 116 (2016) no. 6, 061101, arXiv:1507.05047 [gr-qc].
  * (15) D. Bettoni, J. M. Ezquiaga, K. Hinterbichler, and M. Zumalacárregui, “Speed of Gravitational Waves and the Fate of Scalar-Tensor Gravity,” Phys. Rev. D95 (2017) no. 8, 084029, arXiv:1608.01982 [gr-qc].
  * (16) I. Sawicki, I. D. Saltas, M. Motta, L. Amendola, and M. Kunz, “Nonstandard gravitational waves imply gravitational slip: On the difficulty of partially hiding new gravitational degrees of freedom,” Phys. Rev. D95 (2017) no. 8, 083520, arXiv:1612.02002 [astro-ph.CO].
  * (17) LIGO Scientific, Virgo Collaboration, B. Abbott et al., “GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral,” Phys. Rev. Lett. 119 (2017) no. 16, 161101, arXiv:1710.05832 [gr-qc].
  * (18) A. Goldstein et al., “An Ordinary Short Gamma-Ray Burst with Extraordinary Implications: Fermi-GBM Detection of GRB 170817A,” Astrophys. J. Lett. 848 (2017) no. 2, L14, arXiv:1710.05446 [astro-ph.HE].
  * (19) V. Savchenko et al., “INTEGRAL Detection of the First Prompt Gamma-Ray Signal Coincident with the Gravitational-wave Event GW170817,” Astrophys. J. Lett. 848 (2017) no. 2, L15, arXiv:1710.05449 [astro-ph.HE].
  * (20) LIGO Scientific, Virgo, Fermi-GBM, INTEGRAL Collaboration, B. P. Abbott et al., “Gravitational Waves and Gamma-rays from a Binary Neutron Star Merger: GW170817 and GRB 170817A,” Astrophys. J. Lett. 848 (2017) no. 2, L13, arXiv:1710.05834 [astro-ph.HE].
  * (21) B. P. Abbott et al., “Multi-messenger Observations of a Binary Neutron Star Merger,” Astrophys. J. Lett. 848 (2017) no. 2, L12, arXiv:1710.05833 [astro-ph.HE].
  * (22) C. de Rham and A. J. Tolley, “Speed of gravity,” Phys. Rev. D 101 (2020) no. 6, 063518, arXiv:1909.00881 [hep-th].
  * (23) C. de Rham, S. Melville, and J. Noller, “Positivity bounds on dark energy: when matter matters,” JCAP 08 (2021) 018, arXiv:2103.06855 [astro-ph.CO].
  * (24) J. Noller and A. Nicola, “Cosmological parameter constraints for Horndeski scalar-tensor gravity,” Phys. Rev. D 99 (2019) no. 10, 103502, arXiv:1811.12928 [astro-ph.CO].
  * (25) E. Bellini, A. J. Cuesta, R. Jimenez, and L. Verde, “Constraints on deviations from $\Lambda$CDM within Horndeski gravity,” JCAP 02 (2016) 053, arXiv:1509.07816 [astro-ph.CO]. [Erratum: JCAP 06, E01 (2016)].
  * (26) B. Hu, M. Raveri, N. Frusciante, and A. Silvestri, “Effective Field Theory of Cosmic Acceleration: an implementation in CAMB,” Phys. Rev. D 89 (2014) no. 10, 103530, arXiv:1312.5742 [astro-ph.CO].
  * (27) M. Raveri, B. Hu, N. Frusciante, and A. Silvestri, “Effective Field Theory of Cosmic Acceleration: constraining dark energy with CMB data,” Phys. Rev. D90 (2014) no. 4, 043513, arXiv:1405.1022 [astro-ph.CO].
  * (28) J. Gleyzes, D. Langlois, M. Mancarella, and F. Vernizzi, “Effective Theory of Dark Energy at Redshift Survey Scales,” JCAP 1602 (2016) no. 02, 056, arXiv:1509.02191 [astro-ph.CO].
  * (29) C. D. Kreisch and E. Komatsu, “Cosmological Constraints on Horndeski Gravity in Light of GW170817,” arXiv:1712.02710 [astro-ph.CO].
  * (30) M. Zumalacárregui, E. Bellini, I. Sawicki, J. Lesgourgues, and P. G. Ferreira, “hi_class: Horndeski in the Cosmic Linear Anisotropy Solving System,” JCAP 1708 (2017) no. 08, 019, arXiv:1605.06102 [astro-ph.CO].
  * (31) D. Alonso, E. Bellini, P. G. Ferreira, and M. Zumalacárregui, “Observational future of cosmological scalar-tensor theories,” Phys. Rev. D95 (2017) no. 6, 063502, arXiv:1610.09290 [astro-ph.CO].
  * (32) S. Arai and A. Nishizawa, “Generalized framework for testing gravity with gravitational-wave propagation. II. Constraints on Horndeski theory,” Phys. Rev. D97 (2018) no. 10, 104038, arXiv:1711.03776 [gr-qc].
  * (33) N. Frusciante, S. Peirone, S. Casas, and N. A. Lima, “The road ahead of Horndeski: cosmology of surviving scalar-tensor theories,” arXiv:1810.10521 [astro-ph.CO].
  * (34) R. Reischke, A. S. Mancini, B. M. Schafer, and P. M. Merkel, “Investigating scalar-tensor-gravity with statistics of the cosmic large-scale structure,” arXiv:1804.02441 [astro-ph.CO].
  * (35) A. Spurio Mancini, R. Reischke, V. Pettorino, B. M. Schafer, and M. Zumalacárregui, “Testing (modified) gravity with 3D and tomographic cosmic shear,” Mon. Not. Roy. Astron. Soc. 480 (2018) 3725, arXiv:1801.04251 [astro-ph.CO].
  * (36) J. Noller and A. Nicola, “Radiative stability and observational constraints on dark energy and modified gravity,” Phys. Rev. D 102 (2020) no. 10, 104045, arXiv:1811.03082 [astro-ph.CO].
  * (37) S. Arai, P. Karmakar, and A. Nishizawa, “Cosmological evolution of viable models in the generalized scalar-tensor theory,” Phys. Rev. D 102 (2020) no. 2, 024003, arXiv:1912.01768 [gr-qc].
  * (38) G. Brando, F. T. Falciano, E. V. Linder, and H. E. S. Velten, “Modified gravity away from a $\Lambda$CDM background,” JCAP 11 (2019) 018, arXiv:1904.12903 [astro-ph.CO].
  * (39) R. Arjona, W. Cardona, and S. Nesseris, “Designing Horndeski and the effective fluid approach,” Phys. Rev. D 100 (2019) no. 6, 063526, arXiv:1904.06294 [astro-ph.CO].
  * (40) M. Raveri, “Reconstructing Gravity on Cosmological Scales,” Phys. Rev. D 101 (2020) no. 8, 083524, arXiv:1902.01366 [astro-ph.CO].
  * (41) L. Perenon, J. Bel, R. Maartens, and A. de la Cruz-Dombriz, “Optimising growth of structure constraints on modified gravity,” JCAP 06 (2019) 020, arXiv:1901.11063 [astro-ph.CO].
  * (42) N. Frusciante and L. Perenon, “Effective field theory of dark energy: A review,” Phys. Rept. 857 (2020) 1–63, arXiv:1907.03150 [astro-ph.CO].
  * (43) A. Spurio Mancini, F. Köhlinger, B. Joachimi, V. Pettorino, B. M. Schäfer, R. Reischke, E. van Uitert, S. Brieden, M. Archidiacono, and J. Lesgourgues, “KiDS + GAMA: constraints on horndeski gravity from combined large-scale structure probes,” Mon. Not. Roy. Astron. Soc. 490 (2019) no. 2, 2155–2177, arXiv:1901.03686 [astro-ph.CO].
  * (44) A. Bonilla, R. D’Agostino, R. C. Nunes, and J. C. N. de Araujo, “Forecasts on the speed of gravitational waves at high $z$,” JCAP 03 (2020) 015, arXiv:1910.05631 [gr-qc].
  * (45) T. Baker and I. Harrison, “Constraining Scalar-Tensor Modified Gravity with Gravitational Waves and Large Scale Structure Surveys,” JCAP 01 (2021) 068, arXiv:2007.13791 [astro-ph.CO].
  * (46) S. Joudaki, P. G. Ferreira, N. A. Lima, and H. A. Winther, “Testing Gravity on Cosmic Scales: A Case Study of Jordan-Brans-Dicke Theory,” arXiv:2010.15278 [astro-ph.CO].
  * (47) J. Noller, L. Santoni, E. Trincherini, and L. G. Trombetta, “Scalar-tensor cosmologies without screening,” JCAP 01 (2021) 045, arXiv:2008.08649 [astro-ph.CO].
  * (48) J. Noller, “Cosmological constraints on dark energy in light of gravitational wave bounds,” Phys. Rev. D 101 (2020) no. 6, 063524, arXiv:2001.05469 [astro-ph.CO].
  * (49) D. Traykova, E. Bellini, P. G. Ferreira, C. García-García, J. Noller, and M. Zumalacárregui, “Theoretical priors in scalar-tensor cosmologies: Shift-symmetric Horndeski models,” Phys. Rev. D 104 (2021) no. 8, 083502, arXiv:2103.11195 [astro-ph.CO].
  * (50) G. D. Moore and A. E. Nelson, “Lower bound on the propagation speed of gravity from gravitational Cherenkov radiation,” JHEP 09 (2001) 023, arXiv:hep-ph/0106220.
  * (51) T. B. Littenberg and N. J. Cornish, “Prospects for Gravitational Wave Measurement of ZTF J1539+5027,” Astrophys. J. Lett. 881 (2019) no. 2, L43, arXiv:1908.00678 [astro-ph.IM].
  * (52) T. Baker et al., “Measuring the propagation speed of gravitational waves with LISA,” arXiv:2203.00566 [gr-qc].
  * (53) C. de Rham and S. Melville, “Gravitational Rainbows: LIGO and Dark Energy at its Cutoff,” Phys. Rev. Lett. 121 (2018) no. 22, 221101, arXiv:1806.09417 [hep-th].
  * (54) LIGO Scientific Collaboration, J. Aasi et al., “Advanced LIGO,” Class. Quant. Grav. 32 (2015) 074001, arXiv:1411.4547 [gr-qc].
  * (55) VIRGO Collaboration, F. Acernese et al., “Advanced Virgo: a second-generation interferometric gravitational wave detector,” Class. Quant. Grav. 32 (2015) no. 2, 024001, arXiv:1408.3978 [gr-qc].
  * (56) KAGRA Collaboration, Y. Aso, Y. Michimura, K. Somiya, M. Ando, O. Miyakawa, T. Sekiguchi, D. Tatsumi, and H. Yamamoto, “Interferometer design of the KAGRA gravitational wave detector,” Phys. Rev. D 88 (2013) no. 4, 043007, arXiv:1306.6747 [gr-qc].
  * (57) KAGRA, LIGO Scientific, Virgo Collaboration, B. P. Abbott et al., “Prospects for observing and localizing gravitational-wave transients with Advanced LIGO, Advanced Virgo and KAGRA,” Living Rev. Rel. 21 (2018) no. 1, 3, arXiv:1304.0670 [gr-qc].
  * (58) M. Punturo et al., “The third generation of gravitational wave observatories and their science reach,” Class. Quant. Grav. 27 (2010) 084007.
  * (59) D. Reitze et al., “Cosmic Explorer: The U.S. Contribution to Gravitational-Wave Astronomy beyond LIGO,” Bull. Am. Astron. Soc. 51 (2019) no. 7, 035, arXiv:1907.04833 [astro-ph.IM].
  * (60) LISA Collaboration, P. Amaro-Seoane et al., “Laser Interferometer Space Antenna,” arXiv:1702.00786 [astro-ph.IM].
  * (61) TianQin Collaboration, J. Luo et al., “TianQin: a space-borne gravitational wave detector,” Class. Quant. Grav. 33 (2016) no. 3, 035010, arXiv:1512.02076 [astro-ph.IM].
  * (62) G. W. Horndeski, “Second-order scalar-tensor field equations in a four-dimensional space,” Int. J. Theor. Phys. 10 (1974) 363–384.
  * (63) C. Deffayet, X. Gao, D. A. Steer, and G. Zahariade, “From k-essence to generalised Galileons,” Phys. Rev. D84 (2011) 064039, arXiv:1103.3260 [hep-th].
  * (64) T. Kobayashi, M. Yamaguchi, and J. Yokoyama, “Generalized G-inflation: Inflation with the most general second-order field equations,” Prog. Theor. Phys. 126 (2011) 511–529, arXiv:1105.5723 [hep-th].
  * (65) E. Bellini and I. Sawicki, “Maximal freedom at minimum cost: linear large-scale structure in general modifications of gravity,” JCAP 1407 (2014) 050, arXiv:1404.3713 [astro-ph.CO].
  * (66) L. Berezhiani, G. Chkareuli, and G. Gabadadze, “Restricted Galileons,” Phys. Rev. D88 (2013) 124020, arXiv:1302.0549 [hep-th].
  * (67) K. Koyama, G. Niz, and G. Tasinato, “Effective theory for the Vainshtein mechanism from the Horndeski action,” Phys. Rev. D 88 (2013) 021502, arXiv:1305.0279 [hep-th].
  * (68) S. Melville and J. Noller, “Positivity in the Sky: Constraining dark energy and modified gravity from the UV,” Phys. Rev. D 101 (2020) no. 2, 021502, arXiv:1904.05874 [astro-ph.CO].
  * (69) H. Georgi, “Effective field theory,” Ann. Rev. Nucl. Part. Sci. 43 (1993) 209–252.
  * (70) J. F. Donoghue, “Introduction to the effective field theory description of gravity,” in Advanced School on Effective Theories.  6, 1995.  arXiv:gr-qc/9512024.
  * (71) C. P. Burgess, “Introduction to Effective Field Theory,” Ann. Rev. Nucl. Part. Sci. 57 (2007) 329–362, arXiv:hep-th/0701053.
  * (72) B. Ratra and P. J. E. Peebles, “Cosmological Consequences of a Rolling Homogeneous Scalar Field,” Phys. Rev. D 37 (1988) 3406.
  * (73) C. Wetterich, “Cosmology and the Fate of Dilatation Symmetry,” Nucl. Phys. B 302 (1988) 668–696, arXiv:1711.03844 [hep-th].
  * (74) P. G. Ferreira and M. Joyce, “Cosmology with a primordial scaling field,” Phys. Rev. D 58 (1998) 023503, arXiv:astro-ph/9711102.
  * (75) R. R. Caldwell, R. Dave, and P. J. Steinhardt, “Cosmological imprint of an energy component with general equation of state,” Phys. Rev. Lett. 80 (1998) 1582–1585, arXiv:astro-ph/9708069.
  * (76) C. Armendariz-Picon, T. Damour, and V. F. Mukhanov, “k - inflation,” Phys. Lett. B458 (1999) 209–218, arXiv:hep-th/9904075 [hep-th].
  * (77) J. Garriga and V. F. Mukhanov, “Perturbations in k-inflation,” Phys. Lett. B458 (1999) 219–225, arXiv:hep-th/9904176 [hep-th].
  * (78) C. Armendariz-Picon, V. F. Mukhanov, and P. J. Steinhardt, “A Dynamical solution to the problem of a small cosmological constant and late time cosmic acceleration,” Phys. Rev. Lett. 85 (2000) 4438–4441, arXiv:astro-ph/0004134 [astro-ph].
  * (79) C. Armendariz-Picon, V. F. Mukhanov, and P. J. Steinhardt, “Essentials of k essence,” Phys. Rev. D63 (2001) 103510, arXiv:astro-ph/0006373 [astro-ph].
  * (80) M. A. Luty, M. Porrati, and R. Rattazzi, “Strong interactions and stability in the DGP model,” JHEP 09 (2003) 029, arXiv:hep-th/0303116 [hep-th].
  * (81) A. Nicolis and R. Rattazzi, “Classical and quantum consistency of the DGP model,” JHEP 06 (2004) 059, arXiv:hep-th/0404159 [hep-th].
  * (82) C. de Rham and A. J. Tolley, “DBI and the Galileon reunited,” JCAP 05 (2010) 015, arXiv:1003.5917 [hep-th].
  * (83) C. Burrage, C. de Rham, D. Seery, and A. J. Tolley, “Galileon inflation,” JCAP 01 (2011) 014, arXiv:1009.2497 [hep-th].
  * (84) C. de Rham and R. H. Ribeiro, “Riding on irrelevant operators,” JCAP 1411 (2014) no. 11, 016, arXiv:1405.5213 [hep-th].
  * (85) D. Pirtskhalava, L. Santoni, E. Trincherini, and F. Vernizzi, “Weakly Broken Galileon Symmetry,” JCAP 1509 (2015) no. 09, 007, arXiv:1505.00007 [hep-th].
  * (86) C. de Rham, G. Gabadadze, L. Heisenberg, and D. Pirtskhalava, “Nonrenormalization and naturalness in a class of scalar-tensor theories,” Phys. Rev. D 87 (2013) no. 8, 085017, arXiv:1212.4128 [hep-th].
  * (87) G. Goon, K. Hinterbichler, A. Joyce, and M. Trodden, “Aspects of Galileon Non-Renormalization,” JHEP 11 (2016) 100, arXiv:1606.02295 [hep-th].
  * (88) I. D. Saltas and V. Vitagliano, “Covariantly Quantum Galileon,” Phys. Rev. D 95 (2017) no. 10, 105002, arXiv:1611.07984 [hep-th].
  * (89) L. Heisenberg and C. F. Steinwachs, “Geometrized quantum Galileons,” JCAP 02 (2020) 031, arXiv:1909.07111 [hep-th].
  * (90) L. Heisenberg, J. Noller, and J. Zosso, “Horndeski under the quantum loupe,” JCAP 10 (2020) 010, arXiv:2004.11655 [hep-th].
  * (91) G. Goon, S. Melville, and J. Noller, “Quantum corrections to generic branes: DBI, NLSM, and more,” JHEP 01 (2021) 159, arXiv:2010.05913 [hep-th].
  * (92) S. Mirshekari, N. Yunes, and C. M. Will, “Constraining Generic Lorentz Violation and the Speed of the Graviton with Gravitational Waves,” Phys. Rev. D 85 (2012) 024041, arXiv:1110.2720 [gr-qc].
  * (93) U. Sperhake, C. J. Moore, R. Rosca, M. Agathos, D. Gerosa, and C. D. Ott, “Long-lived inverse chirp signals from core collapse in massive scalar-tensor gravity,” Phys. Rev. Lett. 119 (2017) no. 20, 201103, arXiv:1708.03651 [gr-qc].
  * (94) J. C. Aurrekoetxea, P. G. Ferreira, K. Clough, E. A. Lim, and O. J. Tattersall, “Where is the ringdown? Reconstructing quasinormal modes from dispersive waves,” arXiv:2205.15878 [gr-qc].
  * (95) A. Sesana, “Prospects for Multiband Gravitational-Wave Astronomy after GW150914,” Phys. Rev. Lett. 116 (2016) no. 23, 231102, arXiv:1602.06951 [gr-qc].
  * (96) S. Marsat, J. G. Baker, and T. Dal Canton, “Exploring the Bayesian parameter estimation of binary black holes with LISA,” Phys. Rev. D 103 (2021) no. 8, 083011, arXiv:2003.00357 [gr-qc].
  * (97) A. Toubiana, S. Marsat, S. Babak, J. Baker, and T. Dal Canton, “Parameter estimation of stellar-mass black hole binaries with LISA,” Phys. Rev. D 102 (2020) 124037, arXiv:2007.08544 [gr-qc].
  * (98) A. Klein et al., “The last three years: multiband gravitational-wave observations of stellar-mass binary black holes,” arXiv:2204.03423 [astro-ph.HE].
  * (99) T. Baker, E. Barausse, A. Chen, C. de Rham, M. Pieroni, and G. Tasinato, “Testing gravitational wave propagation with multiband detections,” arXiv:2209.14398 [gr-qc].
  * (100) A. Buonanno, B. Iyer, E. Ochsner, Y. Pan, and B. S. Sathyaprakash, “Comparison of post-Newtonian templates for compact binary inspiral signals in gravitational-wave detectors,” Phys. Rev. D 80 (2009) 084043, arXiv:0907.0700 [gr-qc].
  * (101) LIGO Scientific Collaboration, “LIGO Algorithm Library - LALSuite.” Free software (gpl), 2018.
  * (102) J. M. Ezquiaga, W. Hu, M. Lagos, M.-X. Lin, and F. Xu, “Modified gravitational wave propagation with higher modes and its degeneracies with lensing,” arXiv:2203.13252 [gr-qc].
  * (103) LIGO Scientific, Virgo Collaboration, B. P. Abbott et al., “Tests of General Relativity with GW170817,” Phys. Rev. Lett. 123 (2019) no. 1, 011102, arXiv:1811.00364 [gr-qc].
  * (104) L. Lindblom, B. J. Owen, and D. A. Brown, “Model Waveform Accuracy Standards for Gravitational Wave Data Analysis,” Phys. Rev. D 78 (2008) 124020, arXiv:0809.3844 [gr-qc].
  * (105) R. Abbott et al., “Open data from the first and second observing runs of Advanced LIGO and Advanced Virgo,” SoftwareX 13 (2021) 100658, arXiv:1912.11716 [gr-qc].
  * (106) A. H. Nitz, T. Dent, G. S. Davies, S. Kumar, C. D. Capano, I. Harry, S. Mozzon, L. Nuttall, A. Lundgren, and M. Tápai, “2-OGC: Open Gravitational-wave Catalog of binary mergers from analysis of public Advanced LIGO and Virgo data,” Astrophys. J. 891 (2020) 123, arXiv:1910.05331 [astro-ph.HE].
  * (107) M. Cantiello et al., “A Precise Distance to the Host Galaxy of the Binary Neutron Star Merger GW170817 Using Surface Brightness Fluctuations,” Astrophys. J. Lett. 854 (2018) no. 2, L31, arXiv:1801.06080 [astro-ph.GA].
  * (108) D. A. Brown, I. Harry, A. Lundgren, and A. H. Nitz, “Detecting binary neutron star systems with spin in advanced gravitational-wave detectors,” Phys. Rev. D 86 (2012) 084017, arXiv:1207.6406 [gr-qc].
  * (109) D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, “emcee: The MCMC Hammer,” Publ. Astron. Soc. Pac. 125 (2013) 306–312, arXiv:1202.3665 [astro-ph.IM].
  * (110) C. M. Biwer, C. D. Capano, S. De, M. Cabero, D. A. Brown, A. H. Nitz, and V. Raymond, “PyCBC Inference: A Python-based parameter estimation toolkit for compact binary coalescence signals,” Publ. Astron. Soc. Pac. 131 (2019) no. 996, 024503, arXiv:1807.10312 [astro-ph.IM].
  * (111) LIGO Scientific, VIRGO Collaboration, B. P. Abbott et al., “GW170104: Observation of a 50-Solar-Mass Binary Black Hole Coalescence at Redshift 0.2,” Phys. Rev. Lett. 118 (2017) no. 22, 221101, arXiv:1706.01812 [gr-qc]. [Erratum: Phys.Rev.Lett. 121, 129901 (2018)].
  * (112) LIGO Scientific, Virgo Collaboration, B. P. Abbott et al., “Tests of General Relativity with the Binary Black Hole Signals from the LIGO-Virgo Catalog GWTC-1,” Phys. Rev. D 100 (2019) no. 10, 104036, arXiv:1903.04467 [gr-qc].
  * (113) LIGO Scientific, Virgo Collaboration, R. Abbott et al., “Tests of general relativity with binary black holes from the second LIGO-Virgo gravitational-wave transient catalog,” Phys. Rev. D 103 (2021) no. 12, 122002, arXiv:2010.14529 [gr-qc].
  * (114) LIGO Scientific, VIRGO, KAGRA Collaboration, R. Abbott et al., “Tests of General Relativity with GWTC-3,” arXiv:2112.06861 [gr-qc].
  * (115) A. Klein et al., “Science with the space-based interferometer eLISA: Supermassive black hole binaries,” Phys. Rev. D 93 (2016) no. 2, 024003, arXiv:1511.05581 [gr-qc].
  * (116) M. L. Katz and S. L. Larson, “Evaluating Black Hole Detectability with LISA,” Mon. Not. Roy. Astron. Soc. 483 (2019) no. 3, 3108–3118, arXiv:1807.02511 [gr-qc].
  * (117) AEDGE Collaboration, Y. A. El-Neaj et al., “AEDGE: Atomic Experiment for Dark Matter and Gravity Exploration in Space,” EPJ Quant. Technol. 7 (2020) 6, arXiv:1908.00802 [gr-qc].
  * (118) J. M. Martín-García, “xAct 2002-2014,” http://www.xact.es/ .
  * (119) A. Nitz et al., “PyCBC: v2.0.3 release of PyCBC.” Zenodo, May, 2022\.  https://doi.org/10.5281/zenodo.6583784.
  * (120) C. R. Harris et al., “Array programming with NumPy,” Nature 585 (2020) no. 7825, 357–362, arXiv:2006.10256 [cs.MS].
  * (121) P. Virtanen et al., “SciPy 1.0–Fundamental Algorithms for Scientific Computing in Python,” Nature Meth. 17 (2020) 261, arXiv:1907.10121 [cs.MS].
  * (122) M. Maggiore, Gravitational Waves. Vol. 1: Theory and Experiments.  Oxford Master Series in Physics. Oxford University Press, 2007.