# Testing the Speed of Gravity with Black Hole Ringdown

Sergi Sirera Lahoz  Institute of Cosmology & Gravitation, University of
Portsmouth, Portsmouth, PO1 3FX, U.K.  Johannes Noller  Institute of Cosmology
& Gravitation, University of Portsmouth, Portsmouth, PO1 3FX, U.K.  DAMTP,
University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, U.K.

###### Abstract

We investigate how the speed of gravitational waves, $c_{\rm GW}$, can be
tested by upcoming black hole ringdown observations. We do so in the context
of hairy black hole solutions, where the hair is associated with a new scalar
degree of freedom, forecasting that LISA and TianQin will be able to constrain
deviations of $c_{\rm GW}$ from the speed of light at the ${\cal O}(10^{-4})$
level from a single supermassive black hole merger. We discuss how these
constraints depend on the nature of the scalar hair, what different aspects of
the underlying physics they are sensitive to in comparison with constraints
derived from gravitational wave propagation effects, which observable systems
will place the most stringent bounds, and that constraints are expected to
improve by up to two orders of magnitude with multiple observations. This is
especially interesting for dark energy-related theories, where existing bounds
from GW170817 need not apply at lower frequencies and where upcoming bounds
from lower-frequency missions will therefore be especially powerful. As such,
we also forecast analogous bounds for the intermediate-frequency AEDGE and
DECIGO missions. Finally, we discuss and forecast analogous black hole
ringdown constraints at higher frequencies (so from LVK, the Einstein
Telescope and Cosmic Explorer) and in what circumstances they can yield new
information on top of existing constraints on $c_{\rm GW}$. All calculations
performed in this paper are reproducible via a companion Mathematica notebook
[1].

###### Contents

  1. I Introduction
  2. II Black hole perturbation theory
    1. II.1 Schwarzschild black holes without hair
    2. II.2 Hairy black holes: Background
    3. II.3 Hairy black holes: Quadratic action
    4. II.4 Modified Regge-Wheeler equation
  3. III Quasinormal modes
    1. III.1 Parametrized ringdown
    2. III.2 Parametrizing scalar hair and $\alpha_{T}$
    3. III.3 Parametrized QNMs
  4. IV Parametrized constraints
    1. IV.1 Fisher forecast setup
    2. IV.2 Constraining $A_{T0}$
    3. IV.3 Constraining multiple $A_{Ti}$
    4. IV.4 Dependence on scalar hair profile
  5. V Forecasting observational constraints
    1. V.1 LISA band forecasts
      1. V.2 LVK band forecasts
        1. V.3 Intermediate band forecasts
          1. VI Conclusions
            1. A Parametrised hair
            2. B Coefficients in quadratic action
            3. C General scalar profile

##  I Introduction

Measuring the speed of gravitational waves, $c_{\rm GW}$, places strong
constraints on the ‘medium’ gravitational waves are propagating through and
hence on the particle content of the Universe. In the strong gravity regime,
binary compact object mergers – e.g. binary black hole (BBH) or binary neutron
star (BNS) mergers – are one of the cleanest probes of this particle content.
Here interactions associated with novel particles can leave an imprint in the
inspiral, merger and ringdown phases. These systems can therefore act as a
particle detector, identifying or constraining the new physics that would be a
consequence of such particles and associated ‘fifth forces’. One of the
smoking gun signals for the presence of such new physics is a $c_{\rm GW}$
different from the speed of light, and indeed it has been shown that binary
compact object mergers can place powerful constraints on $c_{\rm GW}$ [2, 3,
4, 5, 6] and hence on the presence and potential dynamics of new degrees of
freedom – see e.g. [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21,
22, 23, 24] and references therein. These previous constraints on $c_{\rm GW}$
from binary compact object mergers have mostly focused on propagation effects
(see e.g. [7, 8, 9, 10] and references therein) or emission effects during the
inspiral phase of such systems (see e.g. [21]). In this paper we instead
investigate what bounds can be derived on $c_{\rm GW}$ from the ringdown phase
alone. This phase is particularly amenable to being understood perturbatively
and hence promises an especially clean analytic understanding. While strong
constraints on $c_{\rm GW}$ exist, most notably from the binary neutron star
merger GW170817 [2, 3, 4, 5, 6], it is important to keep in mind that these
are for frequencies in the LIGO-Virgo-KAGRA (LVK) band, i.e. $\sim 20-2000$
Hz. Expressed as an energy scale this corresponds to $\sim 10^{-14}-10^{-12}$
eV. This range of values is important, because dark energy-related physics is
one of the primary targets that can be constrained with measurements of
$c_{\rm GW}$ and dark energy theories that do affect $c_{\rm GW}$ generically
come with a cutoff around ${\cal O}(10^{2})$ Hz [25]. This means that, for
such theories, an (unknown) high energy completion of the fiducial new dark
energy physics ought to take over as one approaches this cutoff, i.e. close to
or somewhat below the LVK band.111Note that the cutoff is the largest possible
energy/frequency scale, where the high energy completion can take over, but
this can already take place at significantly lower energies/frequencies.
Theoretically predicting the precise scale would require detailed knowledge
about such a fiducial (currently unknown) high energy completion. This high
energy completion will naturally enforce $c_{\rm GW}=c$ at high energies if it
permits Lorentz invariant solutions, so LVK measurements such as GW170817 may
simply confirm that feature of the high energy completion instead of probing
the original (low energy) dark energy physics itself.222The same is true for
bounds from (the absence of) gravitational Cherenkov radiation [26], which
place a lower bound on $c_{\rm GW}$ at energy scales of order $\sim 10^{10}$
GeV, i.e. far above the energy scales probed by gravitational wave detectors.
In other words, in theories that do affect $c_{\rm GW}$ at cosmological
scales, one therefore naturally expects a frequency-dependent transition back
to $c_{\rm GW}=c$ upon approaching the LVK band. With frequencies in the LISA
band $\sim 10^{-4}-10^{-1}$ Hz (and the corresponding energies $\sim
10^{-19}-10^{-16}$ eV) being significantly lower, upcoming LISA [27] and
TianQin [28] measurements therefore provide a much cleaner probe of $c_{\rm
GW}$ in such dark energy related theories.

Existing and upcoming constraints on $c_{\rm GW}$:  Existing constraints on
$c_{\rm GW}$ from lower frequencies, i.e. frequencies below the LVK band, are
comparatively weak, so it is particularly interesting to forecast constraints
from and for LISA. Here it is worth emphasising that a frequency-dependent
$c_{\rm GW}$, so e.g. different speeds in the LVK and LISA bands, is a generic
consequence of the afore-mentioned dark energy theories – see [25, 29, 30] for
more detailed discussions on this point. The existing relevant constraints
closest to (in fact, just below) the LISA band are from binary pulsars, in
particular from the Hulse-Taylor binary, and place a bound of
$|\alpha_{T}|\lesssim{\cal O}(10^{-2})$ for frequencies $f\sim 10^{-5}\;{\rm
Hz}$ [21]. Here we have conveniently expressed bounds on $c_{\rm GW}$ in terms
of the dimensionless $\alpha_{T}$ parameter

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

which we will use throughout this paper. Bounds from even lower frequencies
$f\sim 10^{-18}-10^{-14}\;{\rm Hz}$ come from the cosmic microwave background
and large scale structure measurements (see [31, 32, 33, 34, 35, 36, 37, 38,
39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 44, 50, 51, 52, 53, 54, 55, 56]
and references therein) and require $|\alpha_{T}|\lesssim{\cal O}(1)$. Finally
there are already a number of $c_{\rm GW}$-related forecasts for upcoming
measurements in the LISA band:

  * 1)

[57] forecasted that a multi-messenger observation in the LISA band using
observations of an eclipsing white-dwarf binary will be able to constrain
$|\alpha_{T}|\lesssim 10^{-12}$ (in the event of a non-detection of any
$\alpha_{T}$-related effect).

  * 2)

For the case when there is a significant frequency-dependence for $c_{\rm GW}$
already within the LISA band, [29] used redshift-induced frequency dependence
imprinted on waveforms to be observed in the LISA band (i.e. without the need
for an optical counterpart) to forecast a constraint of $|\alpha_{T}|\lesssim
10^{-4}$.

  * 3)

Also for frequency-dependent $c_{\rm GW}$ within the LISA band, [30]
forecasted that a bound $|\alpha_{T}|\lesssim 10^{-17}$ can be placed by using
the fact that waveforms to be observed by LISA will be
squeezed/stretched/scrambled due to the different speeds with which different
frequencies will propagate (for frequency-dependent $c_{\rm GW}$ and again
without the need for an optical counterpart).

  * 4)

Finally, if there is no detectable frequency-dependence in both the LISA or
LVK bands individually, but a transition in between, [30, 58] showed that
multiband observations using systems such as GW150914 – that are first
observable in the LISA band and later enter the LVK band [59] – will constrain
$|\alpha_{T}|\lesssim 10^{-15}$ (again in the event of a non-detection).333
Several of the forecasted constraints listed here were computed considering
single waveforms. For upcoming future detectors signal overlap will likely be
a regular occurrence, so understanding to what extent this impacts the above
bounds will be an interesting issue to explore going forward.

Looking forward to upcoming LISA observations this leaves us with the
following situation when looking for the strongest possible upcoming bounds.
If there is any significant frequency-dependence in the LISA band, a strong
$|\alpha_{T}|\lesssim 10^{-17}$ bound will very quickly be established once a
single sufficiently loud SMBH (super massive black hole) merger has been
observed. No optical counterpart or multi-band observation is required for
this. If no such frequency-dependence is present, multi-messenger events and
multi-band observations will eventually place bounds at the $10^{-12}$ level
and $10^{-15}$ level, respectively.444Note that the galactic eclipsing white
dwarf binary considered in [57] is a known system which is expected to be
clearly observable in LISA, whereas detection rates for multi-messenger events
more akin to GW170817 (i.e. compact object mergers with a clearly identifiable
optical counterpart that pinpoints the merger itself) are highly uncertain
[60]. Multi-band observations as discussed in [59, 30, 58] will take several
years to constrain $c_{\rm GW}$, given the signal has to ‘migrate’ from the
LISA to the LVK frequency band for the constraint to arise. Here we will show
that additional bounds at the $10^{-4}$ level can be derived from the ringdown
phase of an observed SMBH merger. These bounds are more model-dependent (we
will detail how below), but will effectively be obtainable as soon as LISA
goes online (given an expected ${\cal O}(10-100)$ observable SMBH mergers per
year [61, 62, 63, 64, 65, 66, 67] ). In the LISA context these bounds are
therefore most relevant in the event that no significant frequency-dependence
is detectable within the LISA band itself, e.g. when $c_{\rm GW}$ quickly
asymptotes to a constant value for high and low frequencies and its frequency-
dependence and transition between those asymptotes is effectively localised to
a narrow band between the frequencies accessible by LISA and LVK. We will
further discuss this setup – as pointed out above, this is the same basic
setup as explored in [30, 58] in the context of multi-band observations –
below, as well as how our analysis is affected when frequency-dependence leaks
into the frequency band under investigation. As we will show, the ringdown
bounds on $c_{\rm GW}$ discussed here are also eventually expected to tighten
by up to two orders of magnitude when stacking observations of multiple
events. We will also highlight that such bounds are not just a complementary
and independent constraint on $c_{\rm GW}$, but the fact that they are derived
for a different background space-time compared with constraints from
gravitational wave propagation (black hole vs. cosmological space-times) also
allows us to extract novel insights about the underlying physics.

Scalar-tensor theories: We will focus on theories where the fiducial new
physics is minimal in the sense that it is described by a single scalar degree
of freedom $\phi$, so that we are dealing with a scalar-tensor theory. The
most general such theory which results in second-order equations of motion is
commonly know as Horndeski scalar-tensor theory [68, 69]555For the equivalence
between the formulations of [68] and [69], see [70]., which is governed by the
following action

$\displaystyle S=\int$ $\displaystyle
d^{4}x\sqrt{-g}\Big{[}G_{2}+G_{3}\Box\phi+G_{4}R+$ $\displaystyle
G_{4X}\left[(\square\phi)^{2}-\phi^{\mu\nu}\phi_{\mu\nu}\right]+G_{5}G_{\mu\nu}\phi^{\mu\nu}-$
$\displaystyle\frac{1}{6}G_{5X}\left[(\square\phi)^{3}-3\phi^{\mu\nu}\phi_{\mu\nu}\square\phi+2\phi_{\mu\nu}\phi^{\mu\sigma}\phi_{\sigma}^{\nu}\right]\Big{]}.$
(2)

Here we have introduced the shorthands $\phi_{\mu}\equiv\nabla_{\mu}\phi$ and
$\phi_{\mu\nu}\equiv\nabla_{\nu}\nabla_{\mu}\phi$, and the $G_{i}$ are free
functions of $\phi$ and $X$, where $X\equiv-\frac{1}{2}\phi_{\mu}\phi^{\mu}$.
$G_{iX}$ denotes the partial derivative of $G_{i}$ with respect to $X$. Most
relevant for our purposes will be the ($X$-dependent parts of the) $G_{4}$
interactions and the $G_{5}$ interactions, since (as we will discuss below)
these are the only interactions affecting $\alpha_{T}$. Also note that, for
simplicity, we will often focus on the case where $G_{4}$ is
$\phi$-independent, so that $G_{4\phi}=0$ (where $G_{i\phi}$ denotes the
partial derivative of $G_{i}$ with respect to $\phi$) – we discuss what this
assumption entails in more detail in appendix A.

It is important to highlight that, just like General Relativity (GR), the
Horndeski scalar-tensor theory (2) is an effective field theory (EFT), so has
a limited range of validity. When (2) is taken to be a fiducial dark energy
theory that does affect $c_{\rm GW}$ on cosmological scales, then this theory
only applies up to its cutoff, expected at or below the aforementioned ${\cal
O}(10^{2})$ Hz. This is precisely analogous to the way in which GR is at most
a valid description of gravitational phenomena up to the Planck scale. These
observations have an important practical implication when computing BBH merger
observables as we do here: $c_{\rm GW}$ and hence $\alpha_{T}$ derived from
(2) are in fact frequency-independent as a consequence of the structure of (2)
imposed by the requirement of 2nd order equations of motion. The frequency-
dependence of $c_{\rm GW}$ alluded to above only enters as a consequence of
the unknown UV (high energy) completion of (2), in other words once we are
about to leave the regime of validity of (2). Throughout most of this paper we
will compute and analyse ringdown predictions derived from (2), so we are
implicitly assuming that we are operating within a frequency-window where I)
(2) is firmly within its regime of validity, and hence II) $c_{\rm GW}$ is
effectively frequency-independent within this window. Rigorously computing
analogous predictions in frequency-windows where there is significant
frequency-dependence for $c_{\rm GW}$ would require incorporating at least
some of the effects of the UV completion and hence supplementing/replacing (2)
with the relevant interactions. We will point out the implications of this
assumption in more detail below as well as when one can extrapolate to more
general scenarios.

Outline: With the above setup in place, this paper is organised as follows. We
collect and discuss the relevant results from black hole perturbation theory
in section II, both for ‘bald’ and ‘hairy’ black hole solutions. We extract
the observable quasinormal spectrum from the relevant solutions in section
III, discussing issues related to the parametrisation of $\alpha_{T}$ in the
process. Parametrized constraints are then presented in section IV, where we
analytically compute the precision with which upcoming ringdown observations
will be able to probe $\alpha_{T}$ for a generic observation. We discuss
correlations between different constrainable parameters and how constraints
depend on the underlying interactions. Forecasted constraints for a range of
specific missions and experiments are then discussed in section V. We conclude
in section VI and collect further relevant details in the appendices.

##  II Black hole perturbation theory

Since the ringdown phase of BBH mergers can be well-described perturbatively,
we first ought to discuss the relevant setup in black hole perturbation
theory. We will consider static and spherically symmetric background solutions
that are Ricci-flat ($R_{\mu\nu}=0=R$) here, in particular Schwarzschild
spacetimes. We therefore write the background metric $\bar{g}_{\mu\nu}$ as

$ds^{2}=\bar{g}_{\mu\nu}dx^{\mu}dx^{\nu}=-A(r)dt^{2}+\frac{1}{B(r)}dr^{2}+C(r)d\Omega^{2},$
(3)

where $A,B,C$ are general functions of the radial coordinate $r$ and
$d\Omega^{2}$ is the line element of the standard 2-sphere. We now consider
metric perturbations $h$ around this background, where

$\displaystyle g_{\mu\nu}$ $\displaystyle=\bar{g}_{\mu\nu}+h_{\mu\nu}.$ (4)

Around the static and spherically symmetric backgrounds considered here such
perturbations can be decomposed into odd and even parity perturbations (under
rotations), which decouple from one another at linear order  (i.e. they evolve
independently from one another and can therefore be treated separately). We
will work to leading (linear) order in this paper, but note that this
decoupling does not hold at higher orders – see [71, 72, 73, 74, 75, 76] for
details on the behaviour of higher order modes. In this paper we will
exclusively focus on odd perturbations, which can be written as

$h_{\mu\nu}^{\rm odd}=\left(\begin{array}[]{cccc}0&0&0&h_{0}\\\ 0&0&0&h_{1}\\\
0&0&0&0\\\
h_{0}&h_{1}&0&0\end{array}\right)\sin{\theta}\partial_{\theta}Y_{\ell m}$ (5)

where we have used the Regge-Wheeler gauge [77] and, since we assume a static
background metric, we have set $m=0$ without loss of generality. $h_{0}$ and
$h_{1}$ are functions of $(r,t)$, where the $t$-dependence will be taken to be
of the form $e^{-i\omega t}$. Since perturbations of the scalar $\phi$ are
even under parity transformations and we focus on parity odd modes, we will
therefore only be concerned with metric perturbations. These perturbations are
however affected by the background solution they are propagating on, so odd
metric perturbations will nevertheless be sensitive to the new physics encoded
by the (background solution of the) fiducial scalar degree of freedom we are
probing here.

###  II.1 Schwarzschild black holes without hair

No-hair theorems guaranteeing a trivial scalar field profile exist for a wide
range of scalar-tensor theories [78, 79, 80, 81].666This was first shown for
stationary black holes in minimally coupled Brans-Dicke theories [78], and
subsequently extended to a more general class of scalar-tensor theories
including self-interactions of the scalar [79], to spherically symmetric
static black holes in Galilean-invariant theories [80], and for slowly
rotating black holes in more general shift-symmetric theories [81]. A natural
starting point are therefore Schwarzschild spacetimes with a constant scalar
field background profile $\bar{\phi}$ as e.g. investigated by [82, 83, 84]

$\displaystyle ds^{2}$
$\displaystyle=-\left(1-\tfrac{2M}{r}\right)dt^{2}+\frac{1}{\left(1-\tfrac{2M}{r}\right)}dr^{2}+d\Omega^{2},$
$\displaystyle\bar{\phi}$ $\displaystyle={\rm constant}.$ (6)

Note that, when we mention the ‘background’ or ‘background solution’ going
forward, we refer to both the metric and scalar background solutions, as e.g.
provided in (6). Around the background (6) odd metric perturbations trivially
behave just as in GR, since they are unaffected by the even sector (where
scalar perturbations do induce non-trivial effects) and also do not feel any
effects from the scalar background solution (since this is trivial in the
present no-hair setup). So in order to explore potentially observable effects
induced by the scalar, one ought to either investigate different background
solutions or consider even perturbations. For detailed discussions of the
second option we refer to [85, 82, 83, 86, 87, 88, 89, 90, 91, 92, 93] for
work in the context of Horndeski gravity, and to [86, 94, 95, 96, 97, 98] for
work in the context of other theories (scalar-tensor or otherwise). However,
here we will proceed along the first route, considering the dynamics of odd
perturbations around different background solutions.  We leave an
investigation of how the speed of gravity impacts quasi-normal modes in the
even sector in the presence of a non-trivial background solution (i.e.
combining the two options discussed above) for future work.

###  II.2 Hairy black holes: Background

If $\phi$ acquires a non-trivial background profile, this will provide a
medium for gravitational waves (i.e. here in particular $h^{\rm
odd}_{\mu\nu}$) to travel through and hence can affect $c_{\rm GW}$. Probing
$c_{\rm GW}$ therefore constitutes a powerful test for departures from GR in
such cases, as neatly illustrated in the aforementioned cosmological context.
For the black hole solutions we focus on here, a well-known scalar-tensor
theory example that can have scalar hair are scalar-Gauss-Bonnet (sGB)
theories [81, 99]. In the context of Horndeski theories these are described by
an action that (in addition to a standard kinetic terms) contains a $G_{5}$
interaction where $G_{5}\sim\ln|X|$ [70]. However, instead of focusing on a
specific hairy solution, we will here follow the approach of [100] and
parametrise the scalar-induced hair in a perturbative fashion, but otherwise
remain agnostic about the precise nature of the hair. More specifically, we
will consider a no-hair Schwarzschild black hole solution at lowest order and
introduce small hairy deviations away from this. These can in principle
manifest themselves both in the background solution for the metric as well as
in the scalar profile, so [100] proposed the following parametrised ansatz

$\displaystyle A(r)$ $\displaystyle=B(r)=1-\frac{2M}{r}+\epsilon\delta
A_{1}+\epsilon^{2}\delta A_{2}+{\cal O}(\epsilon^{3}),$ $\displaystyle C(r)$
$\displaystyle=(1+\epsilon\delta C_{1}+\epsilon^{2}\delta C_{2})r^{2}+{\cal
O}(\epsilon^{3})$ $\displaystyle\bar{\phi}$
$\displaystyle=\widehat{\phi}+\epsilon\delta\phi_{1}+\epsilon^{2}\delta\phi_{2}+{\cal
O}(\epsilon^{3}).$ (7)

Here $\delta A_{i},\delta C_{i},\delta\phi_{i}$ are functions of $r$ and
$\epsilon$ is simply a useful order parameter, since we will work
perturbatively up to quadratic order in the (small) hair $\delta\phi$ –
$\epsilon$ has no physical meaning beyond this.777It is worth emphasising an
important subtlety here. As mentioned above, since we are focusing on the odd
parity sector of perturbations, there are no scalar perturbations contributing
in our setup. The $\delta\phi_{i}$ in (7) therefore describe small deviations
in the background solution for the scalar away from a lowest order constant
scalar profile. There are therefore implicitly two perturbative hierarchies at
play here. Metric perturbations $h_{\mu\nu}^{\rm odd}$ and perturbations in
the background hair. We will work up to linear order in $h_{\mu\nu}$ (at the
level of the equations of motion) and up to quadratic order in the hair. Note
that we will denote quantities which are evaluated on the background (so
$h_{\mu\nu}$ is set to zero, recall odd scalar perturbations vanish
identically) with a bar, so $\bar{\phi}$ denotes the scalar field as evaluated
on the background. Quantities where in addition the small (background) scalar
hair $\delta\phi_{i}$ is set to zero are denoted by a hat, so e.g.
$\widehat{\phi}$ denotes the scalar field as evaluated on the background in
the absence of any non-trivial scalar hair. As a consequence we have
$\widehat{X}=0$, while $\bar{X}$ here acquires non-zero contributions via the
$\delta\phi_{i}$. While (7) is a very general parametrisation, for our
purposes we will be able to work with a highly simplified subset. We are
interested in probing the effect of $c_{\rm GW}$ (or equivalently
$\alpha_{T}$) on the ringdown phase. Since deviations from $c_{\rm GW}=c$ (or
equivalently $\alpha_{T}=0$) arise due to a non-trivial scalar field profile
acting as a medium for gravitational waves passing through, it is unsurprising
that at lowest order in $\epsilon$ any $\alpha_{T}$-dependent contribution
only depends on $\delta\phi_{1}$ in (7) and not on $\delta A_{i},\delta
C_{i}$, or $\delta\phi_{2}$. We collect results showing this explicitly in
appendix A, but here we will therefore proceed by working with the much
simpler parametrised ansatz

$\displaystyle ds^{2}$
$\displaystyle=-\left(1-\tfrac{2M}{r}\right)dt^{2}+\frac{1}{\left(1-\frac{2M}{r}\right)}dr^{2}+d\Omega^{2},$
$\displaystyle\bar{\phi}$ $\displaystyle=\widehat{\phi}+\epsilon\delta\phi.$
(8)

Recall that $\delta\phi$ is a small deviation in the background solution
$\bar{\phi}$. This will allow us to identify the leading order contributions
imprinted by a non-zero $\alpha_{T}$, so is ideally suited for our purposes.
We will later discuss to what extent the constraints we will derive on
$\alpha_{T}$ may be contaminated/weakened in the presence of non-zero $\delta
A_{i},\delta C_{i}$, but for now proceed with (8) as a proof of principle.
However, do note that our simplified ansatz is (partially) motivated by sGB-
like hair. There, when working perturbatively in a small sGB coupling, at
leading order only the scalar background acquires a non-trivial contribution,
while the metric remains Schwarzschild [81, 99],888This motivation is only
‘partial’, since 1) it would correspond to setting $\delta A_{1}=0=\delta
C_{1}$ in (7), but not $\delta A_{2},\delta C_{2}$ and 2) because the guiding
principle here is not to explore the consequences of any specific theory, but
rather to explore the consequences of a non-trivial $c_{\rm GW}$ on top of a
parametrised background ansatz. See appendix A for a more in-depth discussion
of what happens when $\delta A_{1}=0=\delta C_{1}$, but $\delta A_{2},\delta
C_{2}$ are non-zero and fully taken into account. i.e. we are working with a
so-called ‘stealth’ solution for the metric.

###  II.3 Hairy black holes: Quadratic action

In order to extract the ringdown signal we need to compute the behaviour of
(odd parity) perturbations on top of the background (8). Working out the
quadratically perturbed action, substituting the components of (5) as well as
our background solution (8), integrating over the angular coordinates and
performing several integrations by parts, we recover the action [101, 102]

$\displaystyle S^{(2)}=\int
dtdr\left[\bar{a}_{1}h_{0}^{2}+\bar{a}_{2}h_{1}^{2}+\bar{a}_{3}\left(\dot{h}_{1}^{2}+h_{0}^{\prime
2}-2\dot{h}_{1}h_{0}^{\prime}+\frac{4}{r}\dot{h}_{1}h_{0}\right)\right],$ (9)

where a dot and a prime denote derivatives with respect to $t$ and $r$,
respectively, and we have dropped an overall multiplicative factor of
$2\pi/(2\ell+1)$ coming from angular integration. The expressions for the
$\bar{a}_{i}$ agree with those found by [101, 102] and satisfy

$\displaystyle\bar{a}_{1}=\frac{\ell(\ell+1)}{2r^{2}}\left[\left(r\mathcal{H}\right)^{\prime}+\frac{(\ell-1)(\ell+2)\mathcal{F}}{2B}+\frac{r^{2}}{B}\epsilon_{A}\right],$
$\displaystyle\bar{a}_{2}=-\frac{\ell(\ell+1)}{2}B\left[\frac{(\ell-1)(\ell+2)\mathcal{G}}{2r^{2}}+\epsilon_{B}\right],$
$\displaystyle\bar{a}_{3}=\frac{\ell(\ell+1)}{4}\mathcal{H}.$ (10)

where the $\bar{a}_{i}$ are to be evaluated on the background (8) (to avoid
clutter bars are implied, but not written explicitly, for all expressions on
the right hand side). $\epsilon_{A,B}$ are contributions that vanish on-shell,
and

$\displaystyle\mathcal{F}$
$\displaystyle=2\left(G_{4}+\frac{1}{2}B\phi^{\prime}X^{\prime}G_{5X}-XG_{5\phi}\right),$
$\displaystyle\mathcal{G}$
$\displaystyle=2\left[G_{4}-2XG_{4X}+X\left(\frac{B^{\prime}}{2}\phi^{\prime}G_{5X}+G_{5\phi}\right)\right],$
$\displaystyle\mathcal{H}$
$\displaystyle=2\left[G_{4}-2XG_{4X}+X\left(\frac{B}{r}\phi^{\prime}G_{5X}+G_{5\phi}\right)\right].$
(11)

The quadratic action (9) contains two fields $(h_{0},h_{1})$, but describes
only one dynamical degree of freedom. [101, 102] show how the action can be
rewritten to make this manifest. To do so an auxiliary field $q$ is defined
and then re-defined into a field $Q$, satisfying999In short, the one-field
quadratic action is obtained by introducing an auxiliary field $q$ into (9)
while leaving the dynamics invariant, varying with respect to $h_{0,1}$ to
obtain $h_{0,1}$ in terms of $q$, substituting them back into the action and
finally performing a field redefinition $q(Q)$. This is explicitly shown in
the accompanying notebook [1].

$\displaystyle h_{0}$
$\displaystyle=-\frac{(r^{2}\bar{a}_{3}q)^{\prime}}{r^{2}\bar{a}_{1}-2(r\bar{a}_{3})^{\prime}},$
$\displaystyle h_{1}$ $\displaystyle=\frac{\bar{a}_{3}}{\bar{a}_{2}}\dot{q},$
$\displaystyle q$ $\displaystyle=\frac{\sqrt{\mathcal{F}}}{r\mathcal{H}}Q.$
(12)

Re-writing the quadratic action in terms of $Q$ in tortoise coordinates
$r_{*}$ (defined as $dr=Bdr_{*}$), one then finds

$\displaystyle S^{(2)}=\frac{\ell(\ell+1)}{4(\ell-1)(\ell+2)}\int
dtdr_{*}\left[\frac{\mathcal{F}}{\mathcal{G}}\dot{Q}^{2}-\left(\frac{dQ}{dr_{*}}\right)^{2}-V(r)Q^{2}\right],$
(13)

where the potential is given by

$\displaystyle V$
$\displaystyle=(\ell+2)(\ell-1)\frac{B}{r^{2}}\frac{\mathcal{F}}{\mathcal{H}}-\frac{r^{3}}{2}\left(\frac{B^{2}}{r^{4}}\right)^{\prime}-\frac{r^{4}\mathcal{F}^{2}}{4\mathcal{F}^{\prime}}\left(\frac{B^{2}\mathcal{F}^{\prime
2}}{r^{4}\mathcal{F}^{3}}\right)^{\prime}.$ (14)

Note that, we have liberally used our background ansatz (8) to simplify the
$\bar{a}_{i}$, etc. in comparison to the more general expressions in [101,
102] – we collect those general expressions in appendix B for comparison.

###  II.4 Modified Regge-Wheeler equation

In order to obtain the analogue of the Regge-Wheeler equation, we now vary the
action with respect to the field $Q$ and find

$\frac{\partial^{2}Q}{\partial
r_{*}^{2}}-\frac{\mathcal{F}}{\mathcal{G}}\frac{\partial^{2}Q}{\partial
t^{2}}-VQ=0.$ (15)

We assume that the time dependence of Q is given by $e^{-i\omega t}$,
substitute $\mathcal{F}$, $\mathcal{G}$, $\mathcal{H}$ and $V$ for our
background (8), and finally obtain the modified Regge-Wheeler equation

$\frac{d^{2}Q}{dr_{*}^{2}}+\Biggl{[}\omega^{2}(1+\epsilon^{2}\alpha_{T})-B(V_{RW}+\epsilon^{2}\delta
V)\Biggr{]}Q=0,$ (16)

where $\alpha_{T}$ satisfies

$\alpha_{T}=-B\frac{G_{4X}-G_{5\phi}}{G_{4}}\delta\phi^{\prime 2}.$ (17)

Note that, given the background we are considering, $\alpha_{T}$ naturally is
a function of $r$ as well as of the Schwarzschild mass $M$. In general there
are further contributions to $\alpha_{T}$ depending on $G_{5X}$, but these
only enter at cubic order in $\epsilon$ as can be deduced from (11), so do not
contribute here.101010This also implies that (17) can be simplified further by
integrating the contributing $G_{5}$ terms by parts in the original theory,
but as we will work with a general $\alpha_{T}$ here, this does not affect our
subsequent expressions. $V_{RW}$ is the well-known Regge-Wheeler potential in
GR

$V_{RW}=\frac{\ell(\ell+1)}{r^{2}}-\frac{6M}{r^{3}},$ (18)

and $\delta V$ is given by

$\begin{split}\delta
V=\alpha_{T}&\Biggl{[}\frac{M(2r-5M)}{r^{3}(r-2M)}+\frac{(\ell+2)(\ell-1)}{r^{2}}\\\
&-\frac{r-2M}{2r}\left(\left(\frac{\delta\phi^{\prime\prime}}{\delta\phi^{\prime}}\right)^{2}-\frac{\delta\phi^{\prime\prime\prime}}{\delta\phi^{\prime}}\right)+\frac{r-5M}{r^{2}}\frac{\delta\phi^{\prime\prime}}{\delta\phi^{\prime}}\Biggr{]}.\end{split}$
(19)

While it has been re-arranged into a more concise form here, this as well as
the above expressions in this subsection agree with the corresponding results
given in [100], when specialised to our ansatz (8). Note that we have
implicitly assumed that $G_{4\phi}=0=G_{4\phi\phi}$ here, as would e.g. be the
case in shift-symmetric theories. We do this to isolate the effect of
$\alpha_{T}$ on the ringdown spectrum, but will further discuss how
$G_{4\phi}\neq 0\neq G_{4\phi\phi}$ would affect our results in appendix A.

##  III Quasinormal modes

Having derived and collected the relevant results from black hole perturbation
theory in the previous section, we are now in a position to extract the key
observable in the black hole ringdown context: the quasinormal modes (QNMs).
As before, we will be focusing on the perturbations of odd modes and the
modified Regge-Wheeler equation (16) governing them. This equation can now be
solved to obtain the frequencies of the associated quasinormal modes $\omega$.
Unlike normal modes, these frequencies are complex numbers, where the real
part represents the physical oscillation frequency and the imaginary part
represents the exponential damping due to dissipation in the system. The QNM
spectrum only depends on the properties of the final black hole (mass, angular
momentum, charge) as well as on the structure of the underlying theory.
Detecting and measuring this spectrum is hence a powerful way to constrain the
presence of novel degrees of freedom and interactions, as well as to generally
test the Kerr hypothesis [103],111111The hypothesis states that the spacetime
around a black hole after gravitational collapse is well described by the Kerr
metric and therefore contains no hair. a scheme that has received the name of
black hole spectroscopy. Note that, while the QNM spectrum does not depend on
initial conditions, the amplitude of individual modes does and this will be
relevant for us in section IV.

###  III.1 Parametrized ringdown

There are a number of techniques one can use to obtain the QNM themselves (see
e.g. [104, 105, 106, 107, 107, 108, 109, 110]) but we refer to [65] for an
extensive review of those.121212The rationale for obtaining numerical
solutions for QNM is schematically the following. One imposes boundary
conditions on the horizon and spatial infinity (corresponding to
$r_{*}\rightarrow\pm\infty$) such that on the horizon wavepackets are moving
inwards and at infinity wavepackets are moving outwards. The imposition of
boundary conditions will then select ‘quantised’ values of $\omega$ (poles in
the Green’s function) which correspond to the QNM. In this paper, we will make
use of the parametrized ringdown formalism [111], the relevant key aspects of
which we will now summarise. In order to apply this formalism, our modified
Regge-Wheeler equation has to be re-cast into the following form

$B\frac{d}{dr}\left(B\frac{dQ}{dr}\right)+\Bigl{[}\omega^{2}-B(V_{RW}+\tilde{\delta
V)}\Bigr{]}Q=0.$ (20)

This can easily be achieved by absorbing the
$\epsilon^{2}\omega^{2}\alpha_{T}$ term into $\delta V$ in (16). Because this
term is a small correction to $\tilde{\delta V}$, we can take $\omega$ to be
the unperturbed frequencies $\omega_{0}$ of the unmodified Regge-Wheeler
equation characteristic of GR, around which we will compute the leading order
$\delta\omega$ corrections below.131313As we will clarify more explicitly in
the next section, there are multiple modes encoded within $\omega_{0}$. Doing
so, we obtain a modified Regge-Wheeler equation in the form of (20) with

$\tilde{\delta V}=\epsilon^{2}\left(\delta
V-\frac{1}{B}\omega_{0}^{2}\alpha_{T}\right),$ (21)

with $\delta V$ being given by (19). It is then instructive to express the
modification to the potential $\tilde{\delta V}$ as an expansion in powers of
$(2M/r)$

$\tilde{\delta
V}=\frac{1}{(2M)^{2}}\sum_{j=0}^{\infty}a_{j}\left(\frac{2M}{r}\right)^{j}.$
(22)

Once expressed in this form, [111] show that the quasinormal frequencies are
determined by the same $a_{j}$ coefficients as follows

$\omega=\omega_{0}+\delta\omega\equiv\omega_{0}+\sum_{j=0}^{\infty}a_{j}e_{j},$
(23)

given that a smallness criterion on the coefficients
$|a_{j}|\ll(1+1/j)^{j}(j+1)$ is satisfied. The $e_{j}$ are a complex ‘basis’
and we summarise the low-order $e_{j}$ most relevant here in table 1 – for
more details and an explicit computation of this basis see [111].

At this point we can already appreciate an important subtlety from the
structure of equations (22) and (23). Each coefficient $a_{j}$ contributing to
$\tilde{\delta V}$ enters with different (increasing) powers of $\sim 1/r$
(22). While this does mean that those contributions to the potential are
suppressed in the far distance limit, i.e. far away from the horizon, it does
not entail that these contributions are providing a sub-dominant contribution
to the frequency spectrum for the QNMs. Indeed, from (23) we explicitly see
that this $\sim 1/r^{j}$ suppression does not play a role in determining the
QNM frequencies. The $j$-th correction enters as $a_{j}e_{j}$ and, while the
$e_{j}$’s tend to slowly decrease in size as $j$ increases, there is no
parametric suppression of higher $j$ contributions. Also note that situations
where (some) higher $j$ contributions dominate over lower $j$ contributions do
arise rather generically – we will see explicit examples below. Finally, note
that the smallness criterion mentioned above guarantees that the $j$-th
contribution to the QNM frequencies is a small correction to $\omega_{0}$, but
this does not entail that the sum of all corrections has to be parametrically
suppressed.

###  III.2 Parametrizing scalar hair and $\alpha_{T}$

Figure 1: Here we plot $\alpha_{T}$, the deviation of the speed of
gravitational waves from that of light defined via $\alpha_{T}\equiv(c_{\rm
GW}^{2}-c^{2})/c^{2}$, as a function of $2M/r$. We show $\alpha_{T}$ for
different example choices of $i$ in (27), where only the amplitude $A_{Ti}$
corresponding to this specific $i$ is non-zero (and fixed to a fiducial value
of $1$) for each plotted curve. Note that $2M/r=0$ therefore corresponds to
spatial infinity whereas $2M/r=1$ corresponds to the Schwarzschild radius, so
we are interested in this range of values/distances. One can clearly see that
$\alpha_{T}(r)=0$ at spatial infinity and at the horizon, but displays non-
trivial behaviour in-between with the overall amplitude and sign determined by
$A_{Ti}$ while the radial dependence depends on the choice of $i$ considered.

Before proceeding with the QNM computation and applying the above formalism to
our (21), we require more information about the functional form of
$\delta\phi$ and $\alpha_{T}$. In close analogy to the above discussion, it is
natural to to view these functions as an expansion in powers of $(2M/r)$ as
well. Starting with the scalar hair function $\delta\phi$, in the main text we
will follow [100] and focus on a scalar hair profile parametrised as

$\delta\phi=\varphi_{c}\left(\frac{2M}{r}\right),$ (24)

where $\varphi_{c}$ is a constant.141414Note that such a profile is indeed
recovered in scalar Gauss-Bonnet theories [81, 99] (in the long distance limit
and when working perturbatively in the Gauss-Bonnet coupling). In appendix C
we discuss the more general parametrisation
$\delta\phi=\varphi_{c}\left(\frac{2M}{r}\right)^{n}$ (where one remains
agnostic of the leading order in $1/r$ at which the hair enters) as well as
superpositions of different $r$-dependencies in the scalar hair profile. We
leave an investigation of even more general (non-power-law) parametrisations
for future investigation. Note that, while in this section we will focus on
the $n=1$ scalar hair profile (24), we will discuss how different profiles
affect eventual constraints on $\alpha_{T}$ in section IV.4. Having
parametrised $\delta\phi$, we turn our attention to the one remaining function
of $r$ affecting $\tilde{\delta V}$, namely $\alpha_{T}$. To this end it will
be useful to separate out the dependence on the scalar hair background profile
and other geometric factors from $\alpha_{T}$ in (17) as follows

$\displaystyle\alpha_{T}$ $\displaystyle=-B(2M)^{2}G_{T}\delta\phi^{\prime
2},$ $\displaystyle G_{T}$
$\displaystyle\equiv\frac{1}{(2M)^{2}}\frac{G_{4X}-G_{5\phi}}{G_{4}}.$ (25)

Here the dimensionless $G_{T}$ parameter has been defined to isolate the
dependence of $\alpha_{T}$ on the Lagrangian $G_{i}$ functions, as opposed to
the $r$-dependence following directly from the scalar profile or via the
dependence on the Schwarzschild function $B(r)$. We will find this separation
especially useful later on when investigating what constraints on $\alpha_{T}$
can tell us about scalar hair and vice versa.151515Note that we have used the
Schwarzschild mass $M$ as a mass scale to define a dimensionless $G_{T}$ here,
but in principle this mass scale is arbitrary. Because we have a non-trivial
scalar profile, all the $G_{i}$ (and hence also $G_{T}$) are functions of $r$
and so to fully specify the $r$-dependence of $\alpha_{T}$ we finally also
expand $G_{T}$ in powers of $r$ along the same lines as discussed for
$\delta\phi$ above. Doing so we can write

$G_{T}=\sum_{i}G_{Ti}\left(\frac{2M}{r}\right)^{i},$ (26)

where each of the $G_{Ti}$ are constant coefficients. Putting everything
together, i.e. substituting $\delta\phi$ (24) and $G_{T}$ (26) into the
expression for $\alpha_{T}$ (25), we can finally write

$\displaystyle\alpha_{T}$
$\displaystyle=-\sum_{i=0}^{\infty}G_{Ti}\varphi_{c}^{2}\left(1-\frac{2M}{r}\right)\left(\frac{2M}{r}\right)^{i+4},$
$\displaystyle=-\sum_{i=0}^{\infty}A_{Ti}\left(1-\frac{2M}{r}\right)\left(\frac{2M}{r}\right)^{i+4}.$
(27)

In the final line we have implicitly defined a final shorthand as part of our
notational setup. The dimensionless amplitude parameters $A_{Ti}$ neatly
encapsulate the coefficients controlling $\alpha_{T}$ and satisfy

$\displaystyle A_{Ti}\equiv G_{Ti}\varphi_{c}^{2}.$ (28)

As one may expect, these are also the effective constant parameters that, as
we will find below, QNM observations will constrain observationally. To
provide some intuition on the relationship between $A_{Ti}$ and $\alpha_{T}$,
we illustrate the dependence of $\alpha_{T}$ on $2M/r$ for various choices of
the amplitude coefficients $A_{Ti}$ in figure 1.

###  III.3 Parametrized QNMs

Having parametrized all the functional freedom encoded within $\tilde{\delta
V}$ above, it is now straightforward to combine the above expressions. Doing
so we can express $\tilde{\delta V}$ as

$\displaystyle\tilde{\delta
V}=\frac{1}{(2M)^{2}}\sum_{i=0}^{\infty}A_{Ti}\Bigg{[}$
$\displaystyle\left(\frac{2M}{r}\right)^{4+i}(2M\omega_{0})^{2}$
$\displaystyle+$
$\displaystyle\left(\frac{2M}{r}\right)^{6+i}\left(-\ell(\ell+1)+9\right)$
$\displaystyle+$
$\displaystyle\left(\frac{2M}{r}\right)^{7+i}\left(\ell(\ell+1)-20\right)$
$\displaystyle+$
$\displaystyle\left(\frac{2M}{r}\right)^{8+i}\frac{45}{4}\Bigg{]},$ (29)

Note that we have dropped the order parameter $\epsilon^{2}$ at this point.
From this we can read off the $a$-coefficients defined in (22)161616For a
given $i$ the contributions to these coefficients are  $\displaystyle a_{4+i}\
\ \rotatebox[origin={c}]{-180.0}{$\in$}\ \ A_{Ti}(2M\omega_{0})^{2},$
$\displaystyle a_{6+i}\ \ \rotatebox[origin={c}]{-180.0}{$\in$}\ \
A_{Ti}(-\ell(\ell+1)+9),$ $\displaystyle a_{7+i}\ \
\rotatebox[origin={c}]{-180.0}{$\in$}\ \ A_{Ti}(\ell(\ell+1)-20),$
$\displaystyle a_{8+i}\ \ \rotatebox[origin={c}]{-180.0}{$\in$}\ \
A_{Ti}\frac{45}{4}.$ (30) By using the ‘element sign’ $\in$ we stress that a
given $a_{j}$ can be built from contributions from different $i$’s. For
instance, $a_{6}$ obtains contributions from $i=0$ and $i=2$. and from (23) we
also obtain the following expression for the quasinormal frequencies

$\displaystyle\delta\omega$ $\displaystyle=\sum_{i=0}^{\infty}A_{Ti}\cdot
E_{i}^{1},$
$\displaystyle=\sum_{i=0}^{\infty}A_{Ti}\Big{[}(2M\omega_{0})^{2}e_{4+i}-(\ell(\ell+1)-9)e_{6+i}$
$\displaystyle\hskip
42.67912pt+(\ell(\ell+1)-20)e_{7+i}+\frac{45}{4}e_{8+i}\Big{]},$ (31)

where $E_{i}^{1}$ has been defined for convenience, and its subscript $1$
refers to $n=1$. A set of $E_{i}^{n}$ ‘basis’ functions for general $n$ are
provided in appendix C.

Anticipating some of our later discussion, we can already see that, for small
$\ell$ and $2M\omega_{0}\sim{\cal O}(1)$, the higher $e_{j}$ terms are
enhanced relative to smaller $j$ terms (for a given $A_{Ti}$). For additional
details on how different contributions enter into $\delta\omega$, see appendix
C and especially figure V. There we also explicitly discuss how this picture
changes for different scalar field profiles. Also note that the consistency
criterion we alluded to above, $|a_{j}|\ll(1+1/j)^{j}(j+1)$, places an
implicit bound on the amplitudes $A_{Ti}$ for the perturbative treatment we
have outlined to be valid. As an example, for the case where only $i=0$ terms
contribute this bound requires that $A_{T0}\ll 1.5$.171717Note that in
practice this bound is set by the $a_{7}$ coefficient, which places a stronger
bound than the other $a_{j}$. Including higher order $i$’s will generate
similar joint constraints on different $A_{Ti}$. As we will see, observational
bounds will constrain the $A_{Ti}$ at the $10^{-1}$ level or stronger. Given
we are measuring deviations away from the standard GR expectation $A_{Ti}=0$,
we therefore expect these consistency bounds to be satisfied in all relevant
scenarios here.

| Re($2Me_{j}$) | Im($2Me_{j}$)
---|---|---
$j=4$ | 0.03668 | -0.00044
$j=5$ | 0.02404 | 0.00273
$j=6$ | 0.01634 | 0.00484
$j=7$ | 0.01136 | 0.00601
$j=8$ | 0.00795 | 0.00654
Table 1: Real and imaginary components of the $e_{j}$ ‘basis’ functions for
$\ell=2$, taken from [111]. Note that we start with $j=4$ as this is our
lowest-order non-zero $a_{j}$. For the full collection up to $j=50$ for each
$\ell$ up to $\ell=10$, together with the ‘basis’ for even-gravitational and
even-scalar perturbations, see [111].

##  IV Parametrized constraints

In the previous section we derived an analytic expression for the QNM
frequencies, assuming these to be close to the corresponding GR frequencies
for a Schwarzschild black hole with the small perturbations encoding
information about interactions in the underlying scalar-tensor theory, in
particular about $\alpha_{T}$. We would now like to use this to forecast how
well future GW experiments will be able to constrain $\alpha_{T}$ using
ringdown. More specifically, we perform a Fisher forecast to estimate the
error in the $A_{Ti}$ parameters (28). In this section we therefore derive
general expressions for the resulting constraints and discuss their overall
features, following this up in the next section by forecasting and discussing
constraints for specific upcoming experiments. Note that, throughout this
section, we ubiquitously use the techniques developed in [112] for our
analysis.

###  IV.1 Fisher forecast setup

| Re($2M\omega_{0}$) | Im($2M\omega_{0}$)
---|---|---
$\ell=2$ | 0.7474 | -0.1779
$\ell=3$ | 1.1989 | -0.1854
$\ell=4$ | 1.6184 | -0.1883
$\ell=5$ | 2.0246 | -0.1897
Table 2: Real and imaginary components of the quasinormal frequencies
$\omega_{0}$ of a Schwarzschild BH in GR for $\ell=2,3,4,5$. Quasinormal data
is provided online [113, 114].

We begin by modelling the waveform as

$h=h_{+}F_{+}+h_{\times}F_{\times},$ (32)

where $h_{+,\times}$ represent the strain in the two polarisations of the
gravitational wave. These are in principle functions of all coordinates, i.e.
$h_{+,\times}(t,r,\theta,\phi)$. However, to distinguish between time and
frequency domains, we will only make $t$ (or $\nu$) explicit and take the
dependence on $r,\theta,\phi$ as understood. $F_{+,\times}$ are functions
encoding the geometry of the problem (i.e. they depend on the angles
specifying the orientation of the source with respect to the detector). The
strain functions for the ringdown are given by

$\displaystyle h_{+}(t)$ $\displaystyle=\sum_{\ell m}A_{\ell
m}^{+}e^{-\frac{\pi tf_{\ell m}}{Q_{\ell m}}}S_{\ell m}\cos{(\phi^{+}_{\ell
m}+2\pi tf_{\ell m})},$ $\displaystyle h_{\times}(t)$
$\displaystyle=\sum_{\ell m}A_{\ell m}^{+}N_{\times}e^{-\frac{\pi tf_{\ell
m}}{Q_{\ell m}}}S_{\ell m}\sin{(\phi^{\times}_{\ell m}+2\pi tf_{\ell m})},$
(33)

where these are the strain functions as emitted by the source and we will
implicitly assume that these trivially propagate to the detector here and will
briefly discuss what this assumption entails and when propagation effects can
be relevant in the next section. In (33) we have absorbed any overall constant
normalization factors into the amplitude parameters $A_{\ell m}^{+}$,181818The
strain functions $h_{+/\times}$ appear with different normalisation factors in
the literature depending on the setup in question, e.g. with a factor of
$1/2\sqrt{10\pi}$, an extra geometrical $\sqrt{3/4}$ for LISA or a
$\frac{1}{r}$ factor [112, 115, 116, 103]. We choose to remain general and
absorb all such factors into the amplitudes $A^{+}$. This does not affect the
calculations presented in this section, as these factors only enter trough the
signal-to-noise-ratio $\rho$, for which the appropriate detector-specific
values will be discussed and used in the next section.  and where $f_{\ell m}$
and $\tau_{\ell m}$ characterise the real and imaginary parts of $\omega_{\ell
m}$ in the following way

$\displaystyle\omega_{\ell m}=2\pi f_{\ell m}+\frac{i}{\tau_{\ell m}},$
$\displaystyle Q_{\ell m}=\pi f_{\ell m}\tau_{\ell m},$ (34)

where $Q_{\ell m}$ is the quality factor. $\\{A^{+}_{\ell m},A^{\times}_{\ell
m}=A^{+}_{\ell m}N_{\times},\phi^{+}_{\ell m},\phi^{\times}_{\ell m}\\}$ are
the amplitudes and phases for the two polarisations. Finally, $S_{\ell m}$ are
spheroidal functions carrying angular dependencies. Because modes with
different $(\ell,m)$ do not mix due to the nature of our background, the $\ell
m$ indices in $\\{\omega_{\ell m},f_{\ell m},\tau_{\ell m},Q_{\ell m},S_{\ell
m},A^{+}_{\ell m},\phi^{+}_{\ell m},\phi^{\times}_{\ell m}\\}$ will play no
role for the time being, so we will obviate them to simplify our notation (and
explicitly discuss which $\ell m$ modes are of interest when this becomes
relevant below).

Using the above strain functions, we compute the signal-to-noise-ratio (SNR)
with the usual

$\rho^{2}=(h|h)=4\int_{0}^{\infty}d\nu\frac{\tilde{h}(\nu)^{*}\tilde{h}(\nu)}{S_{h}(\nu)},$
(35)

where $S_{h}(\nu)$ is the noise spectral density characteristic of the
detector and $\tilde{h}(\nu)$ is the Fourier transform of $h(t)$.191919Note
that in (35) we use $\nu$ rather than $f$ for the frequency domain
representation (or Fourier transform) of the time coordinate. This is to
distinguish it from the real component of the quasinormal modes $f_{\ell m}$
as defined in (34), especially since we will be omitting the $\ell m$ indices.
We now make use of the following set of simplifying assumptions: $\langle
F_{+}\rangle=\langle F_{\times}\rangle=1/5$, $\langle
F_{+}F_{\times}\rangle=0$, $\langle|S|^{2}\rangle=1/4\pi$, $A^{+}=A$. We also
make use of the fact that we can approximate $S_{h}(\nu)$ to be constant. For
details on (and explicit checks of) these assumptions see [112]. Using these
assumptions, (35) can be re-expressed as202020Note that to obtain this simple
expression for $\rho^{2}$ we have further assumed that
$\phi^{+}=\phi^{\times}$ and $N_{\times}=1$ (i.e. $A^{\times}=A^{+}=A$) [112].
However, we stress that this is not a necessary assumption to recover the
expression for the single-parameter error (42) laid out in the following
subsections. We do find it necessary to recover the double-parameter error
expressions (43). Without this assumption, we have
$\rho^{2}=\frac{QA^{2}((1+N_{\times}^{2})(1+4Q^{2})+\cos{2\phi^{+}}-N_{\times}^{2}\cos{2\phi^{\times}})}{(1+N_{\times}^{2})(1+4Q^{2})\pi
fS_{h}}.$ (36)

$\rho^{2}=\frac{QA^{2}}{\pi fS_{h}}.$ (37)

To derive error estimates we make use of the Fisher Information Matrix, given
by

$\Gamma_{ab}=\left(\frac{\delta h}{\delta\theta^{a}}\Big{|}\frac{\delta
h}{\delta\theta^{b}}\right),$ (38)

where $\theta^{a}$ is the set of parameters for our theory and the noise-
weighted product $(\cdot|\cdot)$ is defined as

$(h_{1}|h_{2})=2\int_{0}^{\infty}d\nu\frac{\tilde{h_{1}}^{*}\tilde{h_{2}}+\tilde{h_{2}}^{*}\tilde{h_{1}}}{S_{h}(\nu)}.$
(39)

Then, we can calculate the parameter errors by inverting the Fisher matrix
(which gives the covariance matrix $\Sigma$). The error for a parameter $a$ is
given by

$\sigma_{a}=\sqrt{\Sigma_{aa}}=\sqrt{\Gamma_{aa}^{-1}}.$ (40)

As an initial estimate we will here study the simplified case where all the
usual parameters of the waveform are known $(A,\phi^{+},...)$ and our only
free parameters are the $A_{Ti}$’s (28).  We leave forecasting full joint
constraints to future work. This simplified setup means the estimates which we
will compute below effectively are upper bounds on the precision one can
expect. For a setup as considered here, where the only waveform parameters we
want to constrain are those appearing inside the quasinormal frequencies
$\omega$ (i.e. inside $f$ and $Q$)212121Note that we need two variables to
represent the real and imaginary parts of $\omega$ and, similarly to previous
literature [100] we choose to work with the pair $\\{f,Q\\}$ rather than
$\\{f,\tau\\}$. We note that, if working with the latter, the SNR equation
(37) would be further simplified to  $\rho^{2}=\frac{\tau A^{2}}{S_{h}}.$ (41)
,  general expressions for the errors can be analytically derived. These only
depend on the number of parameters one wants to constrain. In this paper we
will constrain up to two $A_{Ti}$ together so we provide here the expression
for single-parameter constraints222222This simple expression also implicitly
makes use of the large-$Q$ limit (or equivalently large damping times $\tau$).
More terms appear at the $Q^{-4}$ order. For details on the validity of this
approximation we again refer to [112] and the full ‘unapproximated’
expressions are available in the companion notebook [1].

$\sigma_{A_{Ti}}^{2}\rho^{2}=\frac{1}{2}\left(\frac{f}{Qf^{\prime}}\right)^{2},$
(42)

where the prime denotes a derivative with respect to $A_{Ti}$, and for double-
parameter constraints

$\displaystyle\sigma_{A_{Ti}}^{2}\rho^{2}=\frac{\dot{f}^{2}}{2}\frac{(2Q)^{2}+(1-\frac{f\dot{Q}}{\dot{f}Q})^{2}}{(\dot{Q}f^{\prime}-\dot{f}Q^{\prime})^{2}},$
$\displaystyle\sigma_{A_{Tj}}^{2}\rho^{2}=\frac{f^{\prime
2}}{2}\frac{(2Q)^{2}+(1-\frac{fQ^{\prime}}{f^{\prime}Q})^{2}}{(\dot{Q}f^{\prime}-\dot{f}Q^{\prime})^{2}},$
(43)

where again a prime denotes a derivative with respect to $A_{Ti}$, and a dot
represents a derivative with respect to $A_{Tj}$. This matches analogous
expressions in [100] (albeit for non-$\alpha_{T}$-related parameters there).

Before deriving error estimates on different parameter combinations, let us
briefly return to the question of which $(\ell,m)$ modes are of most
interest.232323There is a third index characterising the QNM spectrum, the
overtone number $n$. Here we only focus on the ‘fundamental mode’ $n=0$. Modes
with higher $n$’s (i.e. overtones) are more suppressed by virtue of having
increasing values of $|\text{Im}(\omega)|$. As discussed above, while the QNM
spectrum does not depend on initial conditions, the amplitude of individual
modes does. The dominant observable contributions, i.e. the modes with the
largest amplitudes for astrophysical binary compact object mergers,
generically are the $\ell=m$ modes, more specifically the $(2,2)$ mode [112,
115, 117, 116, 118, 119].242424Note that the discussion in section II only
applies for $\ell\geq 1$ modes and the dipole perturbation $\ell=1$ requires
special treatment, as the Regge-Wheeler gauge used in section II does not
fully fix all gauge degrees of freedom for this mode, see [101] for details.
However, the contribution from the $\ell=1$ mode can be shown to be negligible
for the background solutions we are probing [101], so this is of no concern
here. Note that, for a non-rotating black hole solution as we are focusing on
here, the equations of motion are independent of $m$ [86]. So while $m=0$ is
typically fixed in such setups for simplicity, as we have done here, the
results derived apply for any $m$. The relative amplitude of subdominant modes
(in particular $\ell=3$) grows as the mass ratio $q$ and angular momentum $j$
of the remnant black hole increases [116, 115, 120, 121] –  also see those
references for discussions related to the detectability of such modes.
Nonetheless, the $\ell=2$ mode still generically dominates in all scenarios
and higher $\ell$ modes decay more quickly, see table 2. Note that the damping
time $\tau$ goes as the inverse of the imaginary component of $\omega$, which
increases for higher $\ell$ modes. So in addition to generically possessing a
smaller amplitude, these modes also decay faster. Finally, also notice that,
for binary systems that have orbited each other for a sufficiently long time
for orbits to have approximately circularised, the $\ell=2$ mode will be
additionally enhanced relative to other modes [122, 123, 124]. While it is
straightforward to repeat the analysis for other higher $\ell$ modes,  the
above rationale truly singles out the $\ell=2$ mode as the observationally
most relevant. We will therefore focus on this mode in what follows. Having
said this, a multiple mode analysis will of course be a powerful tool to probe
higher dimensional parameter spaces using ringdown alone tests in the future.
While, as we have seen, quasinormal modes are independent of $m$ for static
black holes, all astrophysical black holes do in fact rotate. For those,
$(2,2)$ is truly the dominant mode, and hence we will focus on this one to
perform the Fisher forecast and leave a more detailed study of constraints in
the presence of several detected modes for the future. Extending the
quasinormal mode calculations in sections II and III to rotating black holes
is an interesting way forward for which some machinery already exists, at
least for slowly-rotating black holes (see e.g. [125, 126, 127, 128, 129]).
However, such metrics and the Schwarzschild metric are smoothly connected
(i.e. taking the limit of zero rotation $j\rightarrow 0$ recovers the
Schwarzschild line element) so one expects that the non-rotating scenario
still captures the leading order information in the quasinormal frequencies
for sufficiently slow rotation.

###  IV.2 Constraining $A_{T0}$

We begin by considering a minimal setup, where there is only a single relevant
$A_{Ti}$ parameter, namely $A_{T0}$. From (31) we then find the QNM shift to
be given by

$\displaystyle\delta\omega=A_{T0}\cdot E_{0}^{1}.$ (44)

where $E_{0}^{1}$ is shown in (31) and we quote it here for reference

$\displaystyle E_{0}^{1}=\Big{[}$
$\displaystyle(2M\omega_{0})^{2}e_{4}-(\ell(\ell+1)-9)e_{6}$
$\displaystyle+(\ell(\ell+1)-20)e_{7}+\frac{45}{4}e_{8}\Big{]}.$ (45)

Substituting in the numerical values for the $e_{j}$ from table 1, we obtain

$M\delta\omega=-[0.00070+0.00306i]\cdot A_{T0}.$ (46)

In evaluating this, we have also used the $\ell=2$ mode in table 2. This now
allows us to obtain parametric expressions for the $\alpha_{T}$-induced
deviations in the QNM spectrum. From (46) and table 2 we find the following
percentage differences for the real and imaginary parts, respectively

$\displaystyle\frac{\delta\omega_{R}}{\omega_{0R}}\approx-0.19\cdot A_{T0}\%,$
$\displaystyle\frac{\delta\omega_{I}}{\omega_{0I}}\approx 3.44\cdot A_{T0}\%.$
(47)

Finally, we are also in a position to extract an expression for the accuracy
with which an experiment with ringdown SNR $\rho$ will be able to measure
$A_{T0}$. Reading off $f$ and $Q$ from (46), as defined in (34), and
substituting them into the single-parameter error expression (42), we obtain
an estimate on its detectability in the same fashion as [100].252525Note that,
in evaluating the final expression, we set $A_{T0}$ to zero. This should
simply be understood as capturing the leading order contributions to the error
– depending on the actual value of $A_{T0}$ the precise error can differ by
$\lesssim{\cal O}$(10%). This gives us262626More precise results are provided
in [1]. Ultimately, we will only be interested in the robust order-of-
magnitude constraints here, so e.g. in table V we will approximate
$\sigma_{A_{T0}}\rho\approx{\cal O}(10^{2})$.

$\sigma_{A_{T0}}\rho\approx 181.$ (48)

The numerical value of this error calculation, as well as the analogous ones
which will follow in this section, will be translated into $\alpha_{T}$
constraints for specific detectors in section V. There, we will compare our
results with other existing and forecasted constraints.

###  IV.3 Constraining multiple $A_{Ti}$

Having considered the single-parameter case above, a natural next step is to
consider a more complex functional form for $\alpha_{T}$ and hence for the
$A_{Ti}$. Here we consider the case where $\alpha_{T}$ is controlled by two
parameters, $A_{T0}$ as before and a second parameter $A_{T1}$. Proceeding as
before, we then have

$\displaystyle\delta\omega=A_{T0}\cdot E_{0}^{1}+A_{T1}\cdot E_{1}^{1}.$ (49)

Reading off expressions for $f$ and $Q$ from equation (49) as before, we find
the following error estimates from (43)

$\displaystyle\sigma_{A_{T0}}\rho$ $\displaystyle\approx 302,$
$\displaystyle\sigma_{A_{T1}}\rho$ $\displaystyle\approx 465.$ (50)

One may wonder how we can constrain two parameters with the measurement of a
single mode. To this end note that the measurement of a single mode carries
information about the oscillation frequency of that mode as well as for the
associated damping time and these independent pieces of information allow
constraining two parameters here. Once future observations are capable of
measuring multiple modes [118, 65, 130], this will of course allow
constraining a correspondingly larger parameter space.

Equation (50) shows that $A_{T0}$ and $A_{T1}$ can be constrained to a similar
order of precision. This re-iterates that terms at higher order in a $1/r$
expansion are not parametrically suppressed in their contribution to the QNM
frequency spectrum, so a $1/r$ expansion is not an ideal basis in terms of
observational constraints. Indeed, upon closer inspection, we find that
constraints on $A_{T0}$ and $A_{T1}$ are strongly correlated, as can be seen
from the off-diagonal elements of the covariance matrix for these two
parameters

$\Sigma_{ab}\sim 10^{4}\begin{pmatrix}9&-16\\\ -16&22\end{pmatrix},$ (51)

We can therefore diagonalize the covariance matrix to obtain the eigenmodes
that will be constrained by the data, i.e. a more optimal basis from a
detectability point of view.272727We thank Sigurd Naess for related
discussions. Under standard matrix diagonalization procedures we obtain

$\tilde{\Sigma}_{ab}=S^{-1}\Sigma_{ab}S\sim 10^{3}\begin{pmatrix}5&0\\\
0&303\end{pmatrix}.$ (52)

This transformation amounts to identifying the combinations of $A_{T0}$ and
$A_{T1}$ that yield uncorrelated parameters $A_{TA}$ and $A_{TB}$, i.e. we
have performed the parameter transformation
$(A_{T0},A_{T1})\rightarrow(A_{TA},A_{TB})$ such that the covariance matrix of
the latter is the one given by equation (52). More explicitly, the relevant
eigenmodes here are $A_{TA}=-0.84A_{T0}-0.54A_{T1}$ and
$A_{TB}=-0.54A_{T0}+0.84A_{T1}$.282828Equivalently, the matrix $S$ is built
with the eigenvectors of (51)
$S=\begin{pmatrix}\vec{s_{1}}&\vec{s_{2}}\end{pmatrix}\sim\begin{pmatrix}-0.84&-0.54\\\
-0.54&0.84\end{pmatrix}.$ (53) Finally, the errors for the new parameters are

$\displaystyle\sigma_{A_{TA}}\rho$ $\displaystyle\approx 68,$
$\displaystyle\sigma_{A_{TB}}\rho$ $\displaystyle\approx 550,$ (54)

where we indeed see that $A_{TA}$ can be constrained more tightly than any
parameter in the previous basis.

###  IV.4 Dependence on scalar hair profile

Figure 2: Here we show forecasted errors on the strength of the interactions
contributing to $c_{\rm GW}$ as parametrised by $G_{Ti}$ (26), where we focus
on $G_{T0}$ as an example. The corresponding error $\sigma_{G_{T0}}$ is shown
as a function of the scalar hair amplitude $\varphi_{c}$ and of the detector
SNR as quantified by $\rho$. We see that $\sigma_{G_{T0}}$ improves as the
scalar hair amplitude grows and as $\rho$ increases, as expected from (56).
More concretely, for an SNR of $\rho\sim 10^{x}$ and a scalar hair amplitude
of order $\varphi_{c}\sim 10^{y}$, we find $\sigma_{G_{T0}}\sim 10^{2-x-2y}$.
From table V, at lower frequencies for LISA we would therefore expect
$\sigma_{G_{T0}}\sim 10^{-3-2y}$ constraints, while for LVK one would need
$y\gtrsim 8$ to yield constraints on the underlying interactions there that
are competitive with or stronger than existing bounds in this band.

Above we have derived expressions for the precision with which a generic
future experiment with SNR $\rho$ will be able to constrain the relevant
parameter combinations affecting QNM frequencies, namely the $A_{Ti}$. Here we
would like to investigate to what extent the specific form of the scalar hair
profile affects this. As we will argue, in certain cases this argument can
then also be inverted to place constraints on the scalar hair itself. Recall
that we are parametrising the scalar hair profile as

$\delta\phi=\varphi_{c}\left(\frac{2M}{r}\right)^{n}.$ (55)

Here $\varphi_{c}$ effectively captures the scalar field amplitude, while $n$
carries information about the radial dependence of this profile. Until now we
have set $n=1$.

Amplitude: The QNM frequencies derived above are functions of the $A_{Ti}$,
which we recall depend both on the scalar amplitude $\varphi_{c}$ as well as
on the $G_{Ti}$ (i.e. the interactions in the underlying theory) via (28).
This has an immediate important consequence, namely that a detection of the
specific QNM shifts discussed here implies both a detection of scalar hair and
of non-trivial $G_{4}$ and/or $G_{5}$ interactions contributing to the
$G_{Ti}$ – cf. (25). The scalar amplitude $\varphi_{c}$, analogously to the
amplitudes of QNMs, will depend on the ‘initial conditions’ for the ringdown
phase. It is worth emphasising that, at present, it is not yet well understood
how the non-linear merger stage affects this amplitude in scalar-tensor
theories of interest, so we will leave $\varphi_{c}$ as a free parameter.
292929$\varphi_{c}$ may be significantly enhanced or suppressed during the
non-linear merger stage, so in the absence of comprehensive numerical (merger)
simulations for the theories in question, even an order of magnitude estimate
appears premature. Note that, in cases where the scalar hair does affect the
black hole geometry (so unlike the ‘stealth’ solutions (8) we consider here),
this effect on the geometry can be used to place additional constraints on the
nature and amplitude of the hair e.g. along the lines presented in [131, 132,
133, 134, 135]. It is interesting, then, to disentangle the effect of
$\varphi_{c}$ and of $G_{Ti}$ on the constrained $A_{Ti}$ parameter(s). This
is shown in Figure 2. As can clearly be seen, and indeed as expected from
(28), in the presence of a larger scalar hair amplitude $\varphi_{c}$ the
constraint on the $G_{Ti}$ becomes stronger. More explicitly

$\displaystyle\sigma_{G_{Ti}}=\sigma_{A_{Ti}}\varphi_{c}^{-2}$ (56)

Interestingly, this implies that one can I) infer a constraint on the scalar
hair amplitude from measurements of the QNMs, given another non-trivial bound
on the $G_{Ti}$ e.g. from non-ringdown related constraints on $c_{\rm GW}$ (in
other words: in the event of a future detection of a $c_{\rm GW}\neq c$ from
another probe), and II) infer a constraint of the $G_{Ti}$, given other
independent information about the amplitude $\varphi_{c}$ (in other words: in
the event of a complementary detection of scalar hair).

Radial dependence: Having considered the effect of the scalar hair amplitude
above, we now investigate how our analysis is affected when the functional
form of the scalar hair, i.e. its $r$-dependence and hence $n$ in (55),
changes. We have considered $n=1$ above and here we repeat the Fisher analysis
for $n=2$ as a complementary example – for further details on the effect of
general $n$ see appendix C. For concreteness, we again consider the single-
parameter $A_{T0}$ case. The corrections to the quasinormal frequencies are
now given by

$\displaystyle\delta\omega=A_{T0}\cdot E_{0}^{2},$ (57)

where for $n=2$ the basis $E_{i}^{n}$ (76) becomes

$\displaystyle E_{i}^{2}=4\Bigg{[}$
$\displaystyle(2M\omega_{0})^{2}e_{6+i}+\left(-\ell(\ell+1)+\frac{31}{2}\right)e_{8+i}$
$\displaystyle+\left(\ell(\ell+1)-\frac{69}{2}\right)e_{9+i}+\frac{76}{4}e_{10+i}\Bigg{]},$
(58)

leading to the following expression for the error on $A_{T0}$

$\sigma_{A_{T0}}\rho\approx 132.$ (59)

We see that the precision with which $A_{T0}$ can be measured has improved for
the $n=2$ case compared to $n=1$. We find this tightening of constraints with
increasing $n$ to be generic and discuss it more in appendix C, also for the
two parameter case with $A_{T0}$ and $A_{T1}$.

##  V Forecasting observational constraints

[1pt] Detector(s) | Ringdown SNR ($\rho$) |  Error on $\alpha_{T}$
---|---|---
[1pt][1pt] LVK |  $10$ [136, 137, 138] |  $1$
[1pt] ET / CE |  $10^{2}$ [139, 138, 140, 141] |  $10^{-1}$
[1pt][1pt] pre-DECIGO |  $10^{2}$ [142] |  $10^{-1}$
[1pt] DECIGO / AEDGE |  $10^{3}$ [143, 144]* |  $10^{-2}$
[1pt] [1pt] LISA |  $10^{5}$ [137, 145] |  $10^{-4}$
[1pt] TianQin |  $10^{5}$ [145] |  $10^{-4}$
[1pt] AMIGO |  $10^{5}$ [130] |  $10^{-4}$
[1pt] |  |

Table 3: Achievable order-of-magnitude ringdown SNRs for a single observed
event for different GW detectors and the corresponding order-of-magnitude
errors on $\alpha_{T}$. Errors in this table are computed assuming $A_{T0}$ is
the only amplitude parameter contributing to $\alpha_{T}$ in (27), as an
example. The error on $\alpha_{T}$, $\sigma_{\alpha_{T}}$, is quoted as one
order-of-magnitude less than the corresponding error on $A_{T0}$, as observed
in figure 1. We stress that the precise mapping of underlying amplitude
parameters to $\alpha_{T}$ mildly depends on the precise functional form of
the scalar hair and the underlying interactions, but note that errors on other
$A_{Ti}$ parameters and hence $\alpha_{T}$ are qualitatively similar – see
e.g. table 4. A star (*) denotes that the quoted forecasted SNR is not
ringdown-specific. For ET/CE we have quoted the ringdown-specific ET forecast
[138], in the current absence (to our knowledge) of an analogous forecast for
CE. For LISA, we note that the quoted SNR is significantly larger than typical
event SNRs in the LISA Mock Data Challenge which go up to $\sim{\cal
O}(10^{3})$ [146, 147], while [145] forecast SNRs up to $\sim{\cal O}(10^{5})$
for (sufficiently nearby and massive) events. This also illustrates that there
is still significant variance in the forecasted SNRs relevant for the missions
considered here.

In the previous section we derived parametric expressions for the precision
with which the parameters controlling the behaviour of $\alpha_{T}$ and hence
$c_{\rm GW}$ can be measured for probes with a general ringdown SNR $\rho$. In
this section we now summarise and briefly discuss what this concretely implies
for a range of current and future missions, spanning the frequency range from
$10^{-4}$ Hz to $10^{3}$ Hz. The main results are collected in table
V.303030It is worth pointing out that the ringdown SNRs quoted implicitly
depend on when precisely the transition from merger to ringdown phase is
assumed to take place. While discussing this in detail is beyond the scope of
this paper, we point the interested reader to [117, 148, 115, 149, 150, 151,
152] for discussions on this.

Before discussing forecasted constraints in detail for the respective missions
and frequency bands, this is a good point to recall our introductory
discussion in section I about how and where the frequency-dependence in
$c_{\rm GW}$ may be localised and how this is tied to the regime (i.e.
frequency range) where the underlying theoretical framework is valid. Rather
obviously, any prediction derived from our starting point – the Horndeski
scalar-tensor action (2) – is only trustworthy when (2) is a faithful
description of the relevant physics. Since (2) gives rise to a frequency-
independent $c_{\rm GW}$, we are therefore implicitly assuming that at the
very least in the frequency-window spanning the ringdown frequencies in
question, $c_{\rm GW}$ is constant as a function of frequency to high
accuracy. A natural scenario to consider would therefore be the one alluded to
in the introduction: $c_{\rm GW}$ effectively becomes a constant $c_{\rm
GW}^{(0)}\neq c$ at low frequencies where (2) applies and may indeed be
intimately linked to dark energy phenomenology on cosmological scales. Now we
consider the ringdown following a SMBH merger observable in the LISA band and
effectively have $c_{\rm GW}=c_{\rm GW}^{(0)}$ there. We can therefore
straightforwardly use (2) to compute this ringdown signal. In this scenario we
also assume (2) stops being an accurate description of the relevant
gravitational physics between the LISA and LVK bands and its unknown UV (high
energy) completion takes over there, resulting in a transition back to $c_{\rm
GW}=c$ at high frequencies due to the Lorentz invariant nature of the UV
completion. The frequency-dependence in $c_{\rm GW}$, induced by the UV
completion is sharply localised in frequency-space between the LISA and LVK
bands and so fully consistent with existing bounds on $c_{\rm GW}$ from the
LVK band. Now this scenario – as explored in detail in the context of
forecasting upcoming multi-band constraints in [30, 58] – is only illustrative
and the frequency-dependence of $c_{\rm GW}$ and the regime of validity of (2)
can easily be altered depending on the UV completion and if the connection to
dark energy is loosened or severed completely. We refer to [25, 29, 30, 58]
and [153] for more detailed discussions of those two points, respectively, and
note that in this paper this is especially relevant in the context of
forecasts for frequency-bands above (i.e. at higher frequencies than) the LISA
band. We will come back to this point below.

What would it take to extrapolate/extend the results from the above sections
to cases where $c_{\rm GW}$ is frequency-dependent in the frequency-window
associated with ringdown signals of interest? On the theoretical side, we
already pointed out that this would involve supplementing/replacing (2) with
the interactions inducing the frequency-dependence of $c_{\rm GW}$, which
requires knowledge of (or assumptions about) the UV completion of (2). The
resulting action could then be used to repeat the analysis for this frequency-
window. It is worth highlighting that the results of sections III and IV only
know about the Horndeski scalar-tensor action by assuming the corresponding
modified form of the Regge-Wheeler equation (16). So any UV completion that
does not modify this form other than inducing a frequency-dependent $c_{\rm
GW}$ and hence $\alpha_{T}$ is covered by the analysis in sections III and IV.
We leave an exploration of how UV completions might otherwise affect the
modified Regge-Wheeler equation and how this affects the subsequent analysis
for future work. On the observational side, a frequency-dependent $c_{\rm GW}$
would introduce another challenge in the ringdown analysis. Since different
parts (i.e. frequencies) of the waveform then travel at different speeds, the
received signal at the detector will be stretched/squeezed/scrambled with
respect to the signal emitted at the source [154, 30] – also see formally
related discussions in [155, 156]. Specifically in the ringdown context, this
can make identifying the correct frequencies more challenging and this
therefore requires a dedicated analysis [154].313131We thank Josu Aurrekoetxea
and Pedro Ferreira for related discussions. In practice this means that the
strain functions (33) accurately describe the signal at emission but will be
altered via non-trivial dispersion effects by the time they reach the
detector, so this needs to be taken into account to correctly forecast
constraints when a frequency-dependent $c_{\rm GW}$ affects the frequency-
window associated with the signal under investigation. We will leave such a
dedicated analysis to future work and (as also motivated by the theoretical
considerations above) in this section forecast constraints for different
frequency bands, assuming an effectively frequency-independent $c_{\rm GW}$
within the band under investigation (i.e. the LISA forecasts assume a
frequency-independent $c_{\rm GW}$ in the LISA band and so on).

###  V.1 LISA band forecasts

As motivated above and in the introduction, the LISA band is particularly
promising in terms of testing for deviations of $c_{\rm GW}$ from the speed of
light, given that a frequency-dependent transition of $c_{\rm GW}$ just or
somewhat below the LVK band is a natural prediction in a range of candidate
dark energy models. In terms of the amplitude parameters $A_{Ti}$, we see from
table V that one expects the leading order such parameter to be constrained at
the $10^{-3}$ level with future LISA/TianQin observations that are forecasted
to yield a ringdown SNR of $\sim{\cal O}(10^{5})$ [137, 145]. Mapping this
back to $\alpha_{T}$ itself (27), this implies one will be able to detect
deviations down to the $\alpha_{T}\sim{\cal O}(10^{-4})$ level from LISA band
ringdown alone in the context of the models we consider,323232Note that,
$\alpha_{T}$ is generically about one order of magnitude smaller than the
dominant $A_{Ti}$ – see figure 1. Also, as should be obvious from (27), this
mapping is mildly dependent on the $i$ coefficient. i.e.

$\displaystyle\sigma_{\alpha_{T}}^{\rm LISA/TianQin}\sim 10^{-4}.$ (60)

Note that present forecasts for far-future missions such as AMIGO predict the
same order-of-magnitude ringdown SNR, so this would not qualitatively alter
constraints on $c_{\rm GW}$ in comparison with those expected from
LISA/TianQin for a single event.

It is worth emphasising that the main bounds discussed here are forecasted for
a single ringdown observation with the SNR achievable by the relevant
detector. It is reasonable to expect that qualitatively improved constraints
will be obtained when combining multiple observations. Indeed, for
sufficiently large $N$ (where $N$ is the number of detected events) the
measurement precision for QNMs is expected to improve as $N^{-1/2}$ [157,
158],  assuming $N$ identical events. For the LISA band, expected event rates
are somewhat uncertain, but most estimates lie in the ${\cal{O}}(10-100)$ per
year range for SMBH mergers – see e.g. [61, 62, 63, 64, 65, 66, 67, 159]. An
improvement of up to two orders of magnitude on the above constraints
therefore seems achievable after several years of operation, so that one may
hope to ultimately reach a precision of close to $\sigma_{\alpha_{T}}^{\rm
LISA/TianQin}\sim 10^{-6}$.  We again emphasise that the $N^{-1/2}$ scaling
discussed here assumes an idealised case with $N$ identical events with the
large SNRs considered here and so the above should be taken as an optimistic
bound (e.g. many events to be detected will be at higher redshifts and have
correspondingly reduced SNRs), also depending on the precise SMBH merger and
detection rates as discussed above.

###  V.2 LVK band forecasts

Having summarised results for the LISA band above, let us consider the LVK
band. The situation here is qualitatively different, given that there already
are tight constraints on $c_{\rm GW}$ specific to this frequency band. From
measuring the coincidence of the GW170817 signal in GW and optical counterpart
observations, one finds that $\alpha_{T}\lesssim 10^{-15}$ [2, 3, 4, 5, 6].
When even a very mild frequency-dependence of $c_{\rm GW}$ is present in the
LVK band, this bound can be strengthened to $\alpha_{T}\lesssim 10^{-17}$
[30]. Contrast this with bounds from ringdown observations alone, where the
$A_{Ti}$ can be constrained at the $\sim{\cal O}(10)$ level with LVK
observations with an improvement by approximately an additional order of
magnitude to be expected from the future Einstein Telescope(ET)/Cosmic
Explorer(CE) missions – cf. table V and see [139, 138] and [140, 141] for ET
and CE, respectively. Again mapping this to constraints on $\alpha_{T}$
itself, we therefore ultimately expect

$\displaystyle\sigma_{\alpha_{T}}^{\rm ET/CE}\sim 10^{-1}.$ (61)

once ET/CE are collecting data in the future.

The bound (61) given above is again for a single event with the SNR achievable
by ET/CE. Taking into account the $N^{-1/2}$ improvement of the measurement
precision for $N$ detected events discussed above, we can again extrapolate
how this precision might be improved over time. For ET
${\cal{O}}(10^{4}-10^{5})$ events with a ringdown SNR of ${\cal O}(10)$ are
expected per year [103]. One may therefore reasonably expect that constraints
can eventually be improved by about two orders of magnitude to
$\sigma_{\alpha_{T}}^{\rm ET/CE}\sim 10^{-3}$. In the LVK context it is also
interesting to point out that existing (non-ringdown-specific) constraints on
$c_{\rm GW}$ from GW waveforms in the LVK band have already seen similar
improvements by stacking events. More specifically, when comparing I)
constraints on $c_{\rm GW}$ from GW170817 data alone (i.e. without using an
optical counterpart) [12] with II) constraints obtained using a LVK catalog of
43 confident binary black hole mergers (used to obtain bounds on the graviton
mass in [160], but straightforwardly re-interpretable to place bounds on
$c_{\rm GW}$), this improves these bounds on $c_{\rm GW}$ by around two orders
of magnitude.333333This improvement, while still partially driven by the
larger number $N$ of observations included, is also partially due to other
events having higher individual SNRs than GW170817.

At first sight (61) as well as the improved $\sigma_{\alpha_{T}}^{\rm
ET/CE}\sim 10^{-3}$ bound reachable by stacking events are rather weak, albeit
complementary, constraints on $c_{\rm GW}$ when compared with the existing GW
propagation bounds from GW170817 discussed above. Also note that, for the
purposes of this subsection and as discussed in detail above, we are assuming
that (2) is a valid description of the underlying physics in (at least part
of) the LVK band. As discussed, in dark energy-related theories within (2)
where $c_{\rm GW}$ receives order one corrections on cosmological scales one
would not expect this to be the case.  One can remedy this (i.e. ‘return’ the
LVK band to within the regime of validity of (2)) in two different ways.
First, by severing the connection to cosmology/dark energy and looking at the
constraints derived here in their own right. Or second, by suppressing the
cosmological $\alpha_{T}$ from the beginning while not precluding a more
sizeable $\alpha_{T}$ around black hole space-times. We will briefly recap a
specific scenario related to the second case below. However, a more general
related point is the following: The fact that the constraints derived in this
paper are computed for a different background solution than cosmological
background GW-propagation constraints derived e.g. from GW170817 means that
they nevertheless contain some interesting new information on the scalar hair
profile and the underlying interactions encoded in $G_{T}$ along the lines
discussed in section IV.4 – we show this in figure 2. More specifically, from
(56) and in the event of a scalar hair amplitude $\varphi_{c}\sim{\cal
O}(10^{8})$, the constraint on the underlying interactions will be as strong
as constraints on the same interactions from GW170817 and even stronger for a
larger amplitude $\varphi_{c}$. Reversing the argument, if future observations
were to identify a small but non-zero cosmological $\alpha_{T}$, this would
allow placing a bound on the scalar hair amplitude from the ringdown
constraints investigated here. For concreteness, consider the following setup:
The higher derivative scalar interactions in the $G_{3,4,5}$ terms in (2) come
with an implicit mass scale $\Lambda$. In cosmology this scale is typically
chosen to be $\Lambda=\Lambda_{3}\equiv(M_{\rm Pl}H_{0}^{2})^{1/3}$, where
this choice ensures those interactions give ${\cal O}(1)$ contributions to
cosmology. However, if a different $\Lambda$ is chosen, the cosmological
$\alpha_{T}$ (and hence $G_{T}$) scales as
$\alpha_{T}\sim(\Lambda_{3}/\Lambda)^{6}$. So raising the interaction scale
$\Lambda$ by just three orders of magnitude suppresses the cosmological
$\alpha_{T}$ down to a level of ${\cal O}(10^{-18})$, comfortably consistent
with bounds from GW170817. This setup also allows the full LVK band to be
within the regime of validity of the physics described by (2) – see [153] for
further details on this scenario. Now suppose that a future constraint indeed
establishes $G_{T}\sim{\cal O}(10^{-18})$, while future ringdown constraints
from ET/CE along the lines investigated here do not yield evidence for a non-
zero $\alpha_{T}$. From (28) this would allow us to derive a bound on the
scalar hair $\varphi_{c}\lesssim{\cal O}(10^{9})$ for frequencies in the LVK
band. Note that other complementary bounds on $\varphi_{c}$ may be obtainable
e.g. from considering even perturbations or going beyond linear theory.

###  V.3 Intermediate band forecasts

With LISA and LVK forecasts discussed above, the intermediate frequency band
stands out as a third region of interest. Here the upcoming AEDGE [144] and
DECIGO [161, 162] experiments will detect and investigate GWs in the future.
In the introduction we motivated probing $c_{\rm GW}$ in the LISA band by
pointing out that a frequency-dependent transition from a nearly constant
$c_{\rm GW}=c$ at LVK frequencies to a different low-frequency $c_{\rm GW}$
naturally occurs just or somewhat below the LVK band in large classes of dark
energy theories. This motivation of course equally applies to the frequencies
probed by AEDGE/DECIGO. Candidate transitions in this intermediate band may
‘leak out’ into the LISA and/or LVK bands, in which case the considerations
outlined above for those bands already promise tight constraints. However,
another interesting class of transitions are those investigated by [30, 58],
where the transition is effectively completely contained within the
intermediate frequency band and no detectable frequency-dependence leaks out
into the LISA and/or LVK bands. In such a case multiband observations using
systems such as GW150914 that are first observable in the LISA band and
eventually enter the LVK band can be used to obtain an integrated constraint
on any features residing at intermediate frequencies and indeed will be able
to constrain $\alpha_{T}$ down to a level of ${\cal O}(10^{-15})$ [30, 58]. In
addition, once AEDGE/DECIGO observations are available, direct constraints on
$c_{\rm GW}$ from this band will be obtainable in analogy to the LVK/LISA
analyses discussed above. Whenever there is significant frequency-dependence
for $c_{\rm GW}$ in band, a complementary ringdown-specific analysis faces
similar theoretical challenges as discussed for the LVK band above, as well as
the observational modelling challenges mentioned earlier in this section. So,
as before, the bounds forecasted in this subsection will be for the case where
(2) applies within (at least part of) the AEDGE/DECIGO band and hence $c_{\rm
GW}$ is frequency-independent in this band to high accuracy. With these
assumptions and from ringdown alone, we find that AEDGE/DECIGO will be able to
constrain the $A_{Ti}$ at the $\sim{\cal O}(10^{-1})$ level, c.f. table V.
Mapping this to constraints on $\alpha_{T}$ itself, as before, this implies

$\displaystyle\sigma_{\alpha_{T}}^{\rm AEDGE/DECIGO}\sim 10^{-2}.$ (62)

While weaker in magnitude than the integrated multiband constraints discussed
above, these bounds are complementary in the same sense as discussed in the
LVK section above. Note that one may again expect this bound to be improved
significantly when stacking multiple observed events: Several dozen
intermediate mass black hole (IMBH) mergers with an SNR ${\cal O}(10^{3})$
should be observable with AEDGE per year [144], so optimistically an
improvement up to $\sigma_{\alpha_{T}}^{\rm AEDGE/DECIGO}\sim 10^{-4}$ appears
feasible eventually.

##  VI Conclusions

In this paper we have investigated how the speed of gravitational waves
$c_{\rm GW}$ can be probed using black hole ringdown observations. Focusing on
scalar-tensor theories of the Horndeski type and on odd parity quasinormal
modes (QNMs), our key findings are as follows:

  * •

In the context of non-rotating black holes where the metric background
solution is given by Schwarzschild, we find that deviations of $c_{\rm GW}$
from the speed of light only affect the QNMs in the presence of a non-trivial
scalar hair profile $\phi=\phi(r)$ in agreement with the results of [100]. Any
deviations from $c_{\rm GW}=c$ are then proportional to the square of the
amplitude of the scalar hair.

  * •

For a single event, ringdown observations from LISA and TianQin will be able
to constrain $c_{\rm GW}$ at the ${\cal O}(10^{-4})$ level. For AEDGE/DECIGO
the equivalent precision will be ${\cal O}(10^{-2})$. When stacking
observations over several years, both constraints may be improved by up to two
orders of magnitude, depending on precise event rates. While those constraints
are weaker than existing constraints on $c_{\rm GW}$ e.g. from GW170817, the
importantly probe different frequency ranges. This is particularly relevant in
the context of testing $c_{\rm GW}$, given large classes of dark energy models
naturally give rise to a frequency-dependent transition in $c_{\rm GW}$ below
the LVK band.

  * •

With ringdown constraints we are testing the effect of deviations from $c_{\rm
GW}=c$ on a different background solution than that relevant for GW
propagation constraints on $c_{\rm GW}$ (black hole vs. cosmological space-
times).  The precise dependence of $c_{\rm GW}$ on interactions in the
underlying theory is different for these two backgrounds. We have highlighted
examples where, in the presence of a sufficiently large scalar hair profile,
ringdown observations can provide novel constraints on those interactions.
Likewise, given complementary information on those underlying interactions, we
have shown how ringdown observations can constrain the nature of scalar hair.
We stress that therefore even LVK band ringdown measurements, where we find
that ${\cal O}(10^{-1})$ level will be obtainable from the Einstein
Telecope/Cosmic Explorer for a single event, can yield valuable information
complementary to existing constraints on $c_{\rm GW}$.

Overall we have therefore derived forecasts for the precision with which
ringdown observations will be able to constrain the speed of gravitational
waves $c_{\rm GW}$ for various detectors throughout the ${\cal
O}(10^{-4})-{\cal O}(10^{3})$ Hz frequency range. Our study has been idealised
in the sense that we have assumed I) the ‘usual’ binary black hole merger
parameters (masses, amplitudes, phases) to be known and focused on the effect
of novel parameters associated with $c_{\rm GW}\neq c$, II) focused on a
specific background solution for the black hole geometry and scalar hair
profile, and III) by working with Horndeski scalar tensor theories, we have
implicitly assumed that $c_{\rm GW}$ is approximately constant in specific
frequency windows/bands when forecasting constraints for those respective
bands. A more comprehensive analysis, extending the present work,
investigating degeneracies and constraints in higher-dimensional parameter
spaces as well as a wider range of hairy black hole solutions and underlying
theoretical setups, will therefore be an interesting next step. As has been
mentioned, another promising route to make our setup more physically realistic
is to extend the quasinormal mode calculations to rotating black hole
solutions [125, 126, 127, 128, 129]. Another step in this direction would be
to include in the analysis surrounding matter fields that dynamically interact
with the black hole in a way that also affects the emitted quasinormal modes –
see e.g. [163, 164, 165, 166, 167, 168, 169, 170]. It is also worth
emphasising that there are several complementary probes of $c_{\rm GW}$ in
addition to the gravitational wave probes discussed throughout this paper and
corresponding to energy/frequency scales outside of the range considered here.
These include constraints from cosmological large scale structure, currently
at the ${\cal O}(1)$ level – see e.g. [32, 31] and references therein – which
are expected to improve to ${\cal O}(10^{-1})$ in the near future [38]. While
we have not mandated a specific sign for any potential deviation of $c_{\rm
GW}$ away from $c$, theoretical bounds from requiring causality, locality and
unitarity at high energies can further yield information on these deviations
at the (comparatively) low energies probed by gravitational waves and
cosmology, noticeably mandating $c_{\rm GW}\geq c$ for large classes of models
[171, 172]. We close by re-emphasising that we have mostly focused on
investigating how well $c_{\rm GW}$ can be tested by ringdown observations,
for a single detected odd-parity quasinormal mode. As more sources and modes
are detected in the future and the theoretical machinery to analyse them is
further developed, we fully expect further tightened constraints to become
obtainable.

## Acknowledgments

We thank Josu Aurrekoetxea, David Bacon, Emanuele Berti, Pedro Ferreira, Ian
Harry, Kazuya Koyama, Andrew Lundgren, Sigurd Naess, Thomas Sotiriou and
Oliver Tattersall for useful discussions. SSL is supported by an STFC
studentship. JN is supported by an STFC Ernest Rutherford Fellowship
(ST/S004572/1). In deriving the results of this paper, we have used xAct [173]
as well as the xAct notebook repository by Reginald Bernardo [174, 175, 90].

##  Appendix A Parametrised hair

In this appendix we quote the results of [100], which consider the full
parametrised hair ansatz (7). The modified Regge-Wheeler equation becomes

$\displaystyle\left[\frac{d^{2}}{dr_{\ast}^{2}}+\omega^{2}\left(1+\epsilon^{2}\alpha_{T}(r)\right)-A(r)\left(\frac{\ell(\ell+1)}{r^{2}}-\frac{6M}{r^{3}}+\epsilon\delta
V_{1}+\epsilon^{2}\delta V_{2}\right)\right]Q=0,$ (63)

where the potential perturbations are given by

$\displaystyle\delta V_{1}=$ $\displaystyle\;\frac{1}{2r^{2}}\left[4\delta
A_{1}-2r\delta A_{1}^{\prime}-2(\ell+2)(\ell-1)\delta C_{1}+2(r-3M)\delta
C_{1}^{\prime}-r(r-2M)\delta
C_{1}^{\prime\prime}-\frac{G_{4\phi}}{G_{4}}\left(r\left(r-2M\right)\delta\phi_{1}^{\prime\prime}-2(r-3M)\delta\phi_{1}^{\prime}\right)\right]$
(64a) $\displaystyle\delta V_{2}=$
$\displaystyle\;\frac{1}{4r^{2}}\left[8\delta A_{2}-4r\delta
A_{2}^{\prime}+4(\ell+2)(\ell-1)\left(\delta C_{1}^{2}-\delta
C_{2}\right)+3r(r-2M)\delta C_{1}^{\prime 2}+4(r-3M)\delta
C_{2}^{\prime}-2r(r-2M)\delta C_{2}^{\prime\prime}+4r\delta A_{1}\delta
C_{1}^{\prime}\right.$ $\displaystyle\left.-2r^{2}\left(\delta
A_{1}^{\prime}\delta C_{1}^{\prime}+\delta A_{1}\delta
C_{1}^{\prime\prime}\right)-4(r-3M)\delta C_{1}\delta
C_{1}^{\prime\prime}+2r(r-2M)\delta C_{1}\delta C_{1}^{\prime\prime}\right]$
$\displaystyle-\frac{1}{2r^{2}}\frac{G_{4\phi}}{G_{4}}\left[-2(r-3M)\delta\phi_{2}^{\prime}+r\left(r\delta
A_{1}^{\prime}\delta\phi_{1}^{\prime}-\delta
A_{1}\left(2\delta\phi_{1}^{\prime}-r\delta\phi_{1}^{\prime\prime}\right)+(r-2M)\left(\delta\phi_{2}^{\prime\prime}-\delta
C_{1}^{\prime}\delta\phi_{1}^{\prime}\right)\right)\right]$
$\displaystyle+\frac{1}{4r^{2}}\left(\frac{G_{4\phi}}{G_{4}}\right)^{2}\left[3r(r-2M)\delta\phi_{1}^{\prime
2}+2\delta\phi_{1}\left(r(r-2M)\delta\phi_{1}^{\prime\prime}-2(r-3M)\delta\phi_{1}^{\prime}\right)\right]$
$\displaystyle-\frac{1}{2r^{2}}\frac{G_{4\phi\phi}}{G_{4}}\left[r(r-2M)\delta\phi_{1}^{\prime
2}+\delta\phi_{1}\left(r(r-2M)\delta\phi_{1}^{\prime\prime}-2(r-3M)\delta\phi_{1}^{\prime}\right)\right]$
$\displaystyle-\frac{\alpha_{T}(r)}{2r^{3}}\left[-5M+Mr(r-2M)^{-1}-2r(\ell+2)(\ell-1)+r^{2}(r-2M)\left(\frac{\delta\phi_{1}^{\prime\prime}}{\delta\phi_{1}^{\prime}}\right)^{2}+r\left(r(r-2M)\frac{\delta\phi_{1}^{\prime\prime\prime}}{\delta\phi_{1}^{\prime}}-2(r-5M)\frac{\delta\phi_{1}^{\prime\prime}}{\delta\phi_{1}^{\prime}}\right)\right].$
(64b)

There are several observations we can make from these expressions. Let us
first point out that the modified Regge-Wheeler equation in our main text (16)
with (19) can be recovered by setting $\delta A_{1,2}=\delta
C_{1,2}=G_{4\phi}=G_{4\phi\phi}=0$. However, we see that $\delta A_{1},\delta
C_{1}$ would contribute at lower order (as well as $\delta\phi_{1}$, if
$G_{4\phi}\neq 0$). If present, such terms can therefore significantly
contribute to the ringdown signal. Indeed, if the fiducial scalar hair is
highly suppressed, i.e. when $\epsilon$ is very small, such lower-order-
in-$\epsilon$ contributions would be expected to dominate over any
$\alpha_{T}$-induced contributions. The motivation behind our simple setup
then is not to fully explore all parametric effects and degeneracies in a
comprehensive parameter space, but rather to isolate and investigate
observable signatures of $\alpha_{T}$. Having said this, it has been shown
that for known scalar-tensor theories that have hairy black holes, namely
scalar-Gauss-Bonnet theory, the metric at leading perturbative order remains
Schwarzschild (i.e. $\delta A_{1}=\delta C_{1}=0$) [99], suggesting that
several aspects of our simple setup are concretely realised in relevant
theories.

In addition, one could consider cases where $G_{4\phi}\neq 0\neq
G_{4\phi\phi}$. It is interesting to point out that interactions contributing
to $G_{4\phi}$ in our background can be removed with a conformal
transformation of the metric. This is because we still have
$\widehat{X}=0$,343434Recall that $\widehat{X}$ refers to the kinetic term $X$
evaluated on the ‘reference background’ of $\widehat{\phi}$ (i.e.
$\delta\phi=0$). meaning that all terms in $G_{4\phi}$ are $X$-independent or,
in other words, those interactions are of the form $f(\phi)R$ in Jordan frame,
which is well known to be convertible to the usual Einstein-Hilbert term by a
conformal transformation. In the resulting Einstein frame representaton, one
then naturally finds $G_{4\phi}=0=G_{4\phi\phi}$. This would indeed make our
assumptions more general and, because the Horndeski group is closed under
conformal transformations and we are not including matter fields, our
calculations would follow exactly in the same way. The price to pay for
working with the metric in the Einstein frame is that observations in GW
detectors are coupled to matter and therefore measure the metric in Jordan
frame. Hence, forecasts should ultimately be made for gravitational waves as
observed in the detector/Jordan frame. We therefore abstain from removing any
interactions via a conformal transformation here and explicitly highlight
setting $G_{4\phi}=0=G_{4\phi\phi}$ as an additional simplifying assumption.
Note that this assumption is trivially satisfied in theories where the scalar
$\phi$ is endowed with a shift symmetry $\phi\to\phi+c$.

To conclude this appendix, let us briefly consider the case where one or many
of $G_{4\phi},G_{4\phi\phi},\delta A_{2},\delta C_{2}$ are non-vanishing. In
this case an analogous analysis to the one performed in this paper can still
be carried out, where the functional form of any non-vanishing such functions
would need to be specified as above. The additional parameters introduced in
this way mean a higher parameter space would then have to be constrained,
presumably degrading constraints on individual parameters. Breaking
degeneracies in such a higher-dimensional parameter space would likely require
the measurements of multiple QNMs. We leave such a more comprehensive
exploration to future work.

##  Appendix B Coefficients in quadratic action

In this appendix we quote the coefficients in the quadratic action for a
general $A$, $B$ and $C$, as derived in [101, 102]. The quadratic action is
given by

$\displaystyle S^{(2)}=\int
dtdr\left[\bar{a}_{1}h_{0}^{2}+\bar{a}_{2}h_{1}^{2}+\bar{a}_{3}(\dot{h}_{1}^{2}+h_{0}^{\prime
2}-2\dot{h}_{1}h_{0}^{\prime}+2\frac{C^{\prime}}{C}\dot{h}_{1}h_{0})\right]$
(65)

with the coefficients being

$\displaystyle\bar{a}_{1}=\frac{\ell(\ell+1)}{4C}\left[\left(C^{\prime}\sqrt{\frac{B}{A}}\mathcal{H}\right)^{\prime}+\frac{(\ell-1)(\ell+2)}{\sqrt{AB}}\mathcal{F}+\frac{2C}{\sqrt{AB}}\epsilon_{A}\right],$
$\displaystyle\bar{a}_{2}=-\frac{\ell(\ell+1)}{2}\sqrt{AB}\left[\frac{(\ell-1)(\ell+2)}{2C}\mathcal{G}+\epsilon_{B}\right],$
$\displaystyle\bar{a}_{3}=\frac{\ell(\ell+1)}{4}\sqrt{\frac{B}{A}}\mathcal{H},$
(66)

where again the bar denotes that the $\bar{a}_{i}$ are evaluated on the
background and, to avoid clutter, bars are implied everywhere on the right
hand side. $\epsilon_{A,B}$ are contributions that vanish on-
shell353535Expressions for $\epsilon_{A,B}$, which are obtained by varying the
action with respect to $A$ and $B$, are fully provided in [102]. and

$\displaystyle\mathcal{F}$
$\displaystyle=2\left(G_{4}+\frac{1}{2}B\phi^{\prime}X^{\prime}G_{5X}-XG_{5\phi}\right),$
$\displaystyle\mathcal{G}$
$\displaystyle=2\left[G_{4}-2XG_{4X}+X\left(\frac{A^{\prime}}{2A}B\phi^{\prime}G_{5X}+G_{5\phi}\right)\right],$
$\displaystyle\mathcal{H}$
$\displaystyle=2\left[G_{4}-2XG_{4X}+X\left(\frac{C^{\prime}}{2C}B\phi^{\prime}G_{5X}+G_{5\phi}\right)\right].$
(67)

As shown in [101, 102], this action can be rewritten to make the presence of
only one degree of freedom explicit

$\displaystyle S^{(2)}=\frac{\ell(\ell+1)}{4(\ell-1)(\ell+2)}\int
dtdr_{*}\left[\frac{\mathcal{F}}{\mathcal{G}}\dot{Q}^{2}-\left(\frac{dQ}{dr^{*}}\right)^{2}-V(r)Q^{2}\right],$

where the potential is given by

$\displaystyle V$
$\displaystyle=\ell(\ell+1)\frac{A}{C}\frac{\mathcal{F}}{\mathcal{H}}-\frac{C^{\prime
2}}{4C^{\prime\prime}}\left(\frac{ABC^{\prime 2}}{C^{3}}\right)^{\prime}$
$\displaystyle\;\;\;\;-\frac{C^{2}\mathcal{F}^{2}}{4\mathcal{F}^{\prime}}\left(\frac{AB\mathcal{F}^{\prime
2}}{C^{2}\mathcal{F}^{3}}\right)^{\prime}-\frac{2A\mathcal{F}}{C\mathcal{H}}.$
(69)

Our quadratic action (9) and its coefficients as well as the potential (14)
can be recovered from the expressions above by specifying our background:
$A=B$, $C=r^{2}$.

##  Appendix C General scalar profile

Here we repeat the derivation of quasinormal corrections but for a more
general scalar profile, given by

$\delta\phi=\varphi_{c}\left(\frac{2M}{r}\right)^{n},$ (70)

which means $\alpha_{T}$ is now given by

$\displaystyle\alpha_{T}=-\sum_{i=0}^{\infty}A_{Ti}\cdot
n^{2}\left(1-\frac{2M}{r}\right)\left(\frac{2M}{r}\right)^{2n+i+2},$ (71)

Figure 3: Here we show $\alpha_{T}$ (71) for different choices of $n$ in (70)
as a function of $2M/r$. $i=0$ has been set such that $A_{T0}$ is the only
non-zero parameter (and fixed to a fiducial value of $1$). We see that
$\alpha_{T}(r)=0$ at spatial infinity and at the horizon, but again observe
non-trivial behaviour in the intermediate region. The size of $\alpha_{T}$
(partially controlled by $A_{T0}$) is considerably enhanced by increasing $n$.
This is mainly due to the factor $n^{2}$ accompanying all $\delta\omega$ (as
can be seen from (75) and (76)).

where the new $n$-dependence is plotted in figure V. Substituting this back
into $\tilde{\delta V}$ (21) we get

$\begin{split}\tilde{\delta
V}=\left(\frac{n}{2M}\right)^{2}&\sum_{i=0}^{\infty}A_{Ti}\Bigg{[}\left(\frac{2M}{r}\right)^{2n+i+2}(2M\omega_{0})^{2}\\\
+&\left(\frac{2M}{r}\right)^{2n+i+4}\left(-\ell(\ell+1)+\frac{9}{2}+\frac{7}{2}n+n^{2}\right)\\\
+&\left(\frac{2M}{r}\right)^{2n+i+5}\left(\ell(\ell+1)-\left(\frac{19}{2}+\frac{17}{2}n+2n^{2}\right)\right)\\\
+&\left(\frac{2M}{r}\right)^{2n+i+6}\left(\frac{21}{4}+5n+n^{2}\right)\Bigg{]}.\end{split}$
(72)

Note that setting $n=1$ recovers the expression (29). Again, this can be
written as

$\tilde{\delta
V}=\frac{1}{(2M)^{2}}\sum_{j=0}^{\infty}a_{j}\left(\frac{2M}{r}\right)^{j}$
(73)

with the only non-zero a-parameters for a given $i$ contributing as

$\displaystyle a_{2n+i+2}\ \ \rotatebox[origin={c}]{-180.0}{$\in$}\ \
A_{Ti}\cdot n^{2}(2M\omega_{0})^{2}$ $\displaystyle a_{2n+i+4}\ \
\rotatebox[origin={c}]{-180.0}{$\in$}\ \ A_{Ti}\cdot
n^{2}\left(-\ell(\ell+1)+\frac{9}{2}+\frac{7}{2}n+n^{2}\right)$ $\displaystyle
a_{2n+i+5}\ \ \rotatebox[origin={c}]{-180.0}{$\in$}\ \ A_{Ti}\cdot
n^{2}\left(\ell(\ell+1)-\left(\frac{19}{2}+\frac{17}{2}n+2n^{2}\right)\right)$
$\displaystyle a_{2n+i+6}\ \ \rotatebox[origin={c}]{-180.0}{$\in$}\ \
A_{Ti}\cdot n^{2}\left(\frac{21}{4}+5n+n^{2}\right).$ (74)

From these we find the following quasinormal frequency corrections

$\delta\omega=\sum_{j=0}^{\infty}a_{j}e_{j}=\sum_{i=0}^{\infty}A_{Ti}E_{i}^{n},$
(75)

where we have defined the following new basis for convenience

$\displaystyle E_{i}^{n}=n^{2}\Bigg{[}$
$\displaystyle(2M\omega_{0})^{2}e_{2n+i+2}$
$\displaystyle+\left(-\ell(\ell+1)+\frac{9}{2}+\frac{7}{2}n+n^{2}\right)e_{2n+i+4}$
$\displaystyle+\left(\ell(\ell+1)-\left(\frac{19}{2}+\frac{17}{2}n+2n^{2}\right)\right)e_{2n+i+5}$
$\displaystyle+\left(\frac{21}{4}+5n+n^{2}\right)e_{2n+i+6}\Bigg{]}.$ (76)

The basis for $n=1$ (31) and $n=2$ (58) used in the main text can be
straightforwardly recovered from this. Note that here, because we have chosen
to remain agnostic about $n$, we have ended up with two indices, i.e. $n$ and
$i$ that need to be chosen in order to obtain numerical results. This is shown
explicitly in the super and subscripts of the newly defined basis $E_{i}^{n}$.
In diagram V we display the corrections coming from different choices of
$(i,n)$, and make some observations about their structure.

The step from these analytically calculated corrections to the quasinormal
modes to the errors on different $A_{Ti}$ parameters (and hence corrections on
$\alpha_{T}$) is straightforwardly repeated in the same fashion as shown in
section IV. We show in table 4 the results for a few more illustrative cases,
but stress that such an analysis can easily be repeated for any combination
and superposition of $(i,n)$ by adapting the companion notebook provided in
[1]. The general trend we find that can already be appreciated in table 4 is
that errors decrease for increasing $i$ and $n$ (at least for the single-
parameter cases).

| $i=0$ | $i=1$ | $i=\\{0,1\\}$
---|---|---|---
$n=1$ | 181 | 131 | (68, 550)
$n=2$ | 132 | 72 | (27, 250)
$n=\\{1,2\\}$ | 153 | 93 | (40, 344)
Table 4: Values for $\sigma_{A_{Ti}}\rho$ for a scalar profile
$\delta\phi=\varphi_{c}\left(\frac{2M}{r}\right)^{n}$. The last column
($i:0,1$) displays the errors for the two parameters $(A_{TA},A_{TB})$,
obtained from $(A_{T0},A_{T1})$ via the same diagonalization procedure as
shown in section IV.3.

Figure 4: In this diagram we display the contributions to $\alpha_{T}$ (71)
coming from different $n$ values. These are ordered such that all
contributions to a specific $A_{Ti}$ parameter appear in one diagonal, which
in principle is extendible ad infinitum. Contributions from different $n$ also
follow a pattern, which can be appreciated here by focusing on the colour of
the balls. The vertical axis $e_{j}$ corresponds to the paramtrized basis from
[111], which tells us the power of $(1/r)^{j}$ at which each contribution
appears in the potential.

## References

  * [1] S. Sirera Lahoz, “ringdown-calculations.”  https://github.com/sergisl/ringdown-calculations.
  * [2] 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].
  * [3] 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].
  * [4] 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].
  * [5] 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].
  * [6] 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].
  * [7] 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].
  * [8] 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].
  * [9] 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].
  * [10] 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].
  * [11] S. Boran, S. Desai, E. O. Kahya, and R. P. Woodard, “GW170817 Falsifies Dark Matter Emulators,” Phys. Rev. D 97 (2018) no. 4, 041501, arXiv:1710.06168 [astro-ph.HE].
  * [12] 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].
  * [13] 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].
  * [14] 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].
  * [15] 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].
  * [16] 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].
  * [17] 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].
  * [18] 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].
  * [19] L. Lombriser and A. Taylor, “Breaking a Dark Degeneracy with Gravitational Waves,” JCAP 1603 (2016) no. 03, 031, arXiv:1509.08458 [astro-ph.CO].
  * [20] 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].
  * [21] 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].
  * [22] 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].
  * [23] 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].
  * [24] N. Cornish, D. Blas, and G. Nardini, “Bounding the speed of gravity with gravitational wave observations,” arXiv:1707.06101 [gr-qc].
  * [25] 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].
  * [26] 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.
  * [27] LISA Collaboration, P. Amaro-Seoane et al., “Laser Interferometer Space Antenna,” arXiv:1702.00786 [astro-ph.IM].
  * [28] 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].
  * [29] T. Baker et al., “Measuring the propagation speed of gravitational waves with LISA,” arXiv:2203.00566 [gr-qc].
  * [30] I. Harry and J. Noller, “Probing the speed of gravity with LVK, LISA, and joint observations,” arXiv:2207.10096 [gr-qc].
  * [31] 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].
  * [32] 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)].
  * [33] 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].
  * [34] 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].
  * [35] 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].
  * [36] C. D. Kreisch and E. Komatsu, “Cosmological Constraints on Horndeski Gravity in Light of GW170817,” arXiv:1712.02710 [astro-ph.CO].
  * [37] 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].
  * [38] 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].
  * [39] 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].
  * [40] 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].
  * [41] 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].
  * [42] 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].
  * [43] 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].
  * [44] 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].
  * [45] 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].
  * [46] 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].
  * [47] M. Raveri, “Reconstructing Gravity on Cosmological Scales,” Phys. Rev. D 101 (2020) no. 8, 083524, arXiv:1902.01366 [astro-ph.CO].
  * [48] 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].
  * [49] 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].
  * [50] 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].
  * [51] 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].
  * [52] 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].
  * [53] 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].
  * [54] 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].
  * [55] 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].
  * [56] 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].
  * [57] 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].
  * [58] 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].
  * [59] A. Sesana, “Prospects for Multiband Gravitational-Wave Astronomy after GW150914,” Phys. Rev. Lett. 116 (2016) no. 23, 231102, arXiv:1602.06951 [gr-qc].
  * [60] P. Amaro-Seoane et al., “Astrophysics with the Laser Interferometer Space Antenna,” arXiv:2203.06016 [gr-qc].
  * [61] E. Berti, “Lisa observations of massive black hole mergers: event rates and issues in waveform modelling,” Class. Quant. Grav. 23 (2006) S785–S798, arXiv:astro-ph/0602470.
  * [62] A. Sesana, F. Haardt, P. Madau, and M. Volonteri, “The gravitational wave signal from massive black hole binaries and its contribution to the LISA data stream,” Astrophys. J. 623 (2005) 23–30, arXiv:astro-ph/0409255.
  * [63] K. J. Rhook and J. S. B. Wyithe, “Realistic event rates for detection of supermassive black hole coalescence by LISA,” Mon. Not. Roy. Astron. Soc. 361 (2005) 1145–1152, arXiv:astro-ph/0503210.
  * [64] T. Tanaka and Z. Haiman, “The Assembly of Supermassive Black Holes at High Redshifts,” Astrophys. J. 696 (2009) 1798–1822, arXiv:0807.4702 [astro-ph].
  * [65] E. Berti, V. Cardoso, and A. O. Starinets, “Quasinormal modes of black holes and black branes,” Class. Quant. Grav. 26 (2009) 163001, arXiv:0905.2975 [gr-qc].
  * [66] eLISA Collaboration, P. A. Seoane et al., “The Gravitational Universe,” arXiv:1305.5720 [astro-ph.CO].
  * [67] M. Bonetti, A. Sesana, F. Haardt, E. Barausse, and M. Colpi, “Post-Newtonian evolution of massive black hole triplets in galactic nuclei – IV. Implications for LISA,” Mon. Not. Roy. Astron. Soc. 486 (2019) no. 3, 4044–4060, arXiv:1812.01011 [astro-ph.GA].
  * [68] G. W. Horndeski, “Second-order scalar-tensor field equations in a four-dimensional space,” Int. J. Theor. Phys. 10 (1974) 363–384.
  * [69] 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].
  * [70] 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].
  * [71] D. Brizuela, J. M. Martin-Garcia, and G. A. Mena Marugan, “Second and higher-order perturbations of a spherical spacetime,” Phys. Rev. D 74 (2006) 044039, arXiv:gr-qc/0607025.
  * [72] D. Brizuela, J. M. Martin-Garcia, and G. A. M. Marugan, “High-order gauge-invariant perturbations of a spherical spacetime,” Phys. Rev. D 76 (2007) 024004, arXiv:gr-qc/0703069.
  * [73] D. Brizuela, J. M. Martin-Garcia, and M. Tiglio, “A Complete gauge-invariant formalism for arbitrary second-order perturbations of a Schwarzschild black hole,” Phys. Rev. D 80 (2009) 024021, arXiv:0903.1134 [gr-qc].
  * [74] M. Lagos and L. Hui, “Generation and propagation of nonlinear quasi-normal modes of a Schwarzschild black hole,” arXiv:2208.07379 [gr-qc].
  * [75] K. Mitman et al., “Nonlinearities in black hole ringdowns,” arXiv:2208.07380 [gr-qc].
  * [76] M. H.-Y. Cheung et al., “Nonlinear effects in black hole ringdown,” arXiv:2208.07374 [gr-qc].
  * [77] T. Regge and J. A. Wheeler, “Stability of a Schwarzschild singularity,” Phys. Rev. 108 (1957) 1063–1069.
  * [78] S. Hawking, “Black holes in the brans-dicke,” Communications in Mathematical Physics 25 (1972) no. 2, 167–171.
  * [79] T. P. Sotiriou and V. Faraoni, “Black holes in scalar-tensor gravity,” Phys. Rev. Lett. 108 (2012) 081103. https://link.aps.org/doi/10.1103/PhysRevLett.108.081103.
  * [80] L. Hui and A. Nicolis, “No-hair theorem for the galileon,” Phys. Rev. Lett. 110 (2013) 241104. https://link.aps.org/doi/10.1103/PhysRevLett.110.241104.
  * [81] T. P. Sotiriou and S.-Y. Zhou, “Black hole hair in generalized scalar-tensor gravity,” Phys. Rev. Lett. 112 (2014) 251102, arXiv:1312.3622 [gr-qc].
  * [82] O. J. Tattersall and P. G. Ferreira, “Quasinormal modes of black holes in Horndeski gravity,” Phys. Rev. D 97 (2018) no. 10, 104047, arXiv:1804.08950 [gr-qc].
  * [83] O. J. Tattersall and P. G. Ferreira, “Forecasts for Low Spin Black Hole Spectroscopy in Horndeski Gravity,” Phys. Rev. D 99 (2019) no. 10, 104082, arXiv:1904.05112 [gr-qc].
  * [84] H. Motohashi and M. Minamitsuji, “General Relativity solutions in modified gravity,” Phys. Lett. B 781 (2018) 728–734, arXiv:1804.01731 [gr-qc].
  * [85] T. Kobayashi, H. Motohashi, and T. Suyama, “Black hole perturbation in the most general scalar-tensor theory with second-order field equations II: the even-parity sector,” Phys. Rev. D 89 (2014) no. 8, 084042, arXiv:1402.6740 [gr-qc].
  * [86] O. J. Tattersall, P. G. Ferreira, and M. Lagos, “General theories of linear gravitational perturbations to a Schwarzschild Black Hole,” Phys. Rev. D 97 (2018) no. 4, 044021, arXiv:1711.01992 [gr-qc].
  * [87] G. Franciolini, L. Hui, R. Penco, L. Santoni, and E. Trincherini, “Effective Field Theory of Black Hole Quasinormal Modes in Scalar-Tensor Theories,” JHEP 02 (2019) 127, arXiv:1810.07706 [hep-th].
  * [88] S. Datta and S. Bose, “Quasi-normal Modes of Static Spherically Symmetric Black Holes in $f(R)$ Theory,” Eur. Phys. J. C 80 (2020) no. 1, 14, arXiv:1904.01519 [gr-qc].
  * [89] J. Khoury, M. Trodden, and S. S. C. Wong, “Existence and instability of hairy black holes in shift-symmetric Horndeski theories,” JCAP 11 (2020) 044, arXiv:2007.01320 [astro-ph.CO].
  * [90] R. C. Bernardo and I. Vega, “Stealth black hole perturbations in kinetic gravity braiding,” J. Math. Phys. 62 (2021) no. 7, 072501, arXiv:2007.06006 [gr-qc].
  * [91] D. Langlois, K. Noui, and H. Roussille, “Linear perturbations of Einstein-Gauss-Bonnet black holes,” JCAP 09 (2022) 019, arXiv:2204.04107 [gr-qc].
  * [92] M. Minamitsuji, K. Takahashi, and S. Tsujikawa, “Linear stability of black holes in shift-symmetric Horndeski theories with a time-independent scalar field,” Phys. Rev. D 105 (2022) no. 10, 104001, arXiv:2201.09687 [gr-qc].
  * [93] L. Hui, A. Podo, L. Santoni, and E. Trincherini, “Effective Field Theory for the perturbations of a slowly rotating black hole,” JHEP 12 (2021) 183, arXiv:2111.02072 [hep-th].
  * [94] C. de Rham and J. Zhang, “Perturbations of stealth black holes in degenerate higher-order scalar-tensor theories,” Phys. Rev. D 100 (2019) no. 12, 124023, arXiv:1907.00699 [hep-th].
  * [95] D. Langlois, K. Noui, and H. Roussille, “Black hole perturbations in modified gravity,” Phys. Rev. D 104 (2021) no. 12, 124044, arXiv:2103.14750 [gr-qc].
  * [96] K. Glampedakis and H. O. Silva, “Eikonal quasinormal modes of black holes beyond General Relativity,” Phys. Rev. D 100 (2019) no. 4, 044040, arXiv:1906.05455 [gr-qc].
  * [97] C.-Y. Chen and S. Park, “Black hole quasinormal modes and isospectrality in Deser-Woodard nonlocal gravity,” Phys. Rev. D 103 (2021) no. 6, 064029, arXiv:2101.06600 [gr-qc].
  * [98] C.-Y. Chen, M. Bouhmadi-López, and P. Chen, “Lessons from black hole quasinormal modes in modified gravity,” Eur. Phys. J. Plus 136 (2021) no. 2, 253, arXiv:2103.01249 [gr-qc].
  * [99] T. P. Sotiriou and S.-Y. Zhou, “Black hole hair in generalized scalar-tensor gravity: An explicit example,” Phys. Rev. D 90 (2014) 124063, arXiv:1408.1698 [gr-qc].
  * [100] O. J. Tattersall, “Quasi-Normal Modes of Hairy Scalar Tensor Black Holes: Odd Parity,” Class. Quant. Grav. 37 (2020) no. 11, 115007, arXiv:1911.07593 [gr-qc].
  * [101] T. Kobayashi, H. Motohashi, and T. Suyama, “Black hole perturbation in the most general scalar-tensor theory with second-order field equations I: the odd-parity sector,” Phys. Rev. D 85 (2012) 084025, arXiv:1202.4893 [gr-qc]. [Erratum: Phys.Rev.D 96, 109903 (2017)].
  * [102] A. Ganguly, R. Gannouji, M. Gonzalez-Espinoza, and C. Pizarro-Moya, “Black hole stability under odd-parity perturbations in Horndeski gravity,” Class. Quant. Grav. 35 (2018) no. 14, 145008, arXiv:1710.07669 [gr-qc].
  * [103] E. Berti, A. Sesana, E. Barausse, V. Cardoso, and K. Belczynski, “Spectroscopy of Kerr black holes with Earth- and space-based interferometers,” Phys. Rev. Lett. 117 (2016) no. 10, 101102, arXiv:1605.09286 [gr-qc].
  * [104] V. Ferrari and B. Mashhoon, “New approach to the quasinormal modes of a black hole,” Phys. Rev. D 30 (1984) 295–304.
  * [105] C. J. Goebel, “Comments on the “vibrations” of a Black Hole.,” ApJ 172 (1972) L95.
  * [106] B. F. Schutz and C. M. Will, “Black hole normal modes - A semianalytic approach,” ApJ 291 (1985) L33–L36.
  * [107] L. Motl and A. Neitzke, “Asymptotic black hole quasinormal frequencies,” Adv. Theor. Math. Phys. 7 (2003) no. 2, 307–330, arXiv:hep-th/0301173.
  * [108] G. T. Horowitz and V. E. Hubeny, “Quasinormal modes of AdS black holes and the approach to thermal equilibrium,” Phys. Rev. D 62 (2000) 024027, arXiv:hep-th/9909056.