Testing the Speed of Gravity with Black Hole Ringdown

.

Measuring the speed of gravitational waves, c 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 GW different from the speed of light, and indeed it has been shown that binary compact object mergers can place powerful constraints on c 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 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 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 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. ∼ 20 − 2000 Hz.Expressed as an energy scale this corresponds to ∼ 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 GW and dark energy theories that do affect c GW generically come with a cutoff around 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. 1 This high energy completion will naturally enforce c 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. 2 In other words, in theories that do affect c GW at cosmological scales, one therefore naturally expects a frequency-dependent transition back to c GW = c upon approaching the LVK band.With frequencies in the LISA band ∼ 10 −4 − 10 −1 Hz (and the corresponding energies ∼ 10 −19 − 10 −16 eV) being significantly lower, upcoming LISA [27] and TianQin [28] measurements therefore provide a much cleaner probe of c GW in such dark energy related theories.
Existing and upcoming constraints on c GW : Existing constraints on c 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 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 |α T | ≲ O(10 −2 ) for frequencies f ∼ 10 −5 Hz [21].Here we have conveniently expressed bounds on c GW in terms of the dimensionless α T parameter which we will use throughout this paper.Bounds from even lower frequencies f ∼ 10 −18 − 10 −14 Hz come from the cosmic microwave background and large scale structure measurements (see  and references therein) and require |α T | ≲ O (1).Finally there are already a number of c GWrelated 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 whitedwarf binary will be able to constrain |α T | ≲ 10 −12 (in the event of a non-detection of any α T -related effect).
2) For the case when there is a significant frequencydependence for c 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 |α T | ≲ 10 −4 .
3) Also for frequency-dependent c GW within the LISA band, [30] forecasted that a bound |α T | ≲ 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 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 |α T | ≲ 10 −15 (again in the event of a non-detection). 3ooking 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 |α T | ≲ 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. 4Here 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 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 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 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 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 φ , 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] 5 , which is governed by the following action Here we have introduced the shorthands φ µ ≡ ∇ µ φ and φ µν ≡ ∇ ν ∇ µ φ , and the G i are free functions of φ and X, where X ≡ − 1 2 φ µ φ µ .G iX denotes the partial derivative of G i with respect to X.Most relevant for our purposes will be the (Xdependent parts of the) G 4 interactions and the G 5 interactions, since (as we will discuss below) these are the only interactions affecting α T .Also note that, for simplicity, we will often focus on the case where G 4 is φ -independent, so that G 4φ = 0 (where G iφ denotes the partial derivative of G i with respect to φ ) -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 GW on cosmological scales, then this theory only applies up to its cutoff, expected at or below the aforementioned 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 GW and hence α 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 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 GW is effectively frequency-independent within this window.
Rigorously computing analogous predictions in frequency-windows where there is significant frequency-dependence for c 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 α 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 α 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 welldescribed 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 µν = 0 = R) here, in particular Schwarzschild spacetimes.We therefore write the background metric ḡµν as where A, B,C are general functions of the radial coordinate r and dΩ 2 is the line element of the standard 2-sphere.We now consider metric perturbations h around this background, where 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 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ωt .Since perturbations of the scalar φ 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.

A. 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]. 6A natural starting point are therefore Schwarzschild spacetimes with a constant scalar field background profile φ as e.g.investigated by [82][83][84] 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 [82,83,[85][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 6 This 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].
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.

B. Hairy black holes: Background
If φ acquires a non-trivial background profile, this will provide a medium for gravitational waves (i.e.here in particular h odd µν ) to travel through and hence can affect c GW .Probing c 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 ∼ 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 Here δ A i , δC i , δ φ i are functions of r and ε is simply a useful order parameter, since we will work perturbatively up to quadratic order in the (small) hair δ φε has no physical meaning beyond this. 7Note that we will denote quantities which are evaluated on the background (so h µν is set to zero, recall odd scalar perturbations vanish identically) with a bar, so φ denotes the scalar field as evaluated on the background.Quantities where in addition the small (background) scalar hair δ φ i is set to zero are denoted by a hat, so e.g.φ denotes the scalar field as evaluated on the background in the absence of any non-trivial scalar hair.As a consequence we have X = 0, while X here acquires non-zero contributions via the δ φ i .While ( 7) is a very general parametrisation, for our 7 It 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 δ φ 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 odd µν and perturbations in the background hair.We will work up to linear order in h µν (at the level of the equations of motion) and up to quadratic order in the hair.purposes we will be able to work with a highly simplified subset.We are interested in probing the effect of c GW (or equivalently α T ) on the ringdown phase.Since deviations from c GW = c (or equivalently α 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 ε any α T -dependent contribution only depends on δ φ 1 in (7) and not on δ A i , δC i , or δ φ 2 .We collect results showing this explicitly in appendix A, but here we will therefore proceed by working with the much simpler parametrised ansatz Recall that δ φ is a small deviation in the background solution φ .This will allow us to identify the leading order contributions imprinted by a non-zero α T , so is ideally suited for our purposes.We will later discuss to what extent the constraints we will derive on α T may be contaminated/weakened in the presence of non-zero δ A i , δ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], 8i.e.we are working with a so-called 'stealth' solution for the metric.

C. 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] 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π/(2ℓ + 1) coming from angular integration.The expressions for the āi agree with those found by [101,102] and satisfy (7), but not δ A 2 , δ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 GW on top of a parametrised background ansatz.See appendix A for a more in-depth discussion of what happens when δ A 1 = 0 = δC 1 , but δ A 2 , δC 2 are non-zero and fully taken into account.
where the ā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).ε A,B are contributions that vanish on-shell, and 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, satisfying 9 Re-writing the quadratic action in terms of Q in tortoise coordinates r * (defined as dr = Bdr * ), one then finds where the potential is given by Note that, we have liberally used our background ansatz (8) to simplify the āi , etc. in comparison to the more general expressions in [101,102] -we collect those general expressions in appendix B for comparison.

D. 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 9 In 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].
We assume that the time dependence of Q is given by e −iωt , substitute F , G , H and V for our background (8), and finally obtain the modified Regge-Wheeler equation where α T satisfies Note that, given the background we are considering, α T naturally is a function of r as well as of the Schwarzschild mass M.
In general there are further contributions to α T depending on G 5X , but these only enter at cubic order in ε as can be deduced from (11), so do not contribute here. 10V RW is the well-known Regge-Wheeler potential in GR and δV is given by 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φ = 0 = G 4φ φ here, as would e.g.be the case in shift-symmetric theories.We do this to isolate the effect of α T on the ringdown spectrum, but will further discuss how G 4φ ̸ = 0 ̸ = G 4φ φ 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 ω.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], 11 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.

A. Parametrized ringdown
There are a number of techniques one can use to obtain the QNM themselves (see e.g.[104][105][106][107][108][109][110]) but we refer to [65] for an extensive review of those. 12In 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 This can easily be achieved by absorbing the ε 2 ω 2 α T term into δV in (16).Because this term is a small correction to δ V , we can take ω to be the unperturbed frequencies ω 0 of the unmodified Regge-Wheeler equation characteristic of GR, around which we will compute the leading order δ ω corrections below. 13Doing so, we obtain a modified Regge-Wheeler equation in the form of (20) with with δV being given by (19).It is then instructive to express the modification to the potential δ V as an expansion in powers of (2M/r) Once expressed in this form, [111] show that the quasinormal frequencies are determined by the same a j coefficients as follows given that a smallness criterion on the coefficients |a j | ≪ (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 Ifor 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 δ V enters with different (increasing) powers of ∼ 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 ∼ 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 ω 0 , but this does not entail that the sum of all corrections has to be parametrically suppressed.

B. Parametrizing scalar hair and α T
Before proceeding with the QNM computation and applying the above formalism to our (21), we require more information about the functional form of δ φ and α 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 δ φ , in the main text we will follow [100] and focus on a scalar hair profile parametrised as where ϕ c is a constant. 14In appendix C we discuss the more general parametrisation δ φ = ϕ c 2M r 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 Here we plot α T , the deviation of the speed of gravitational waves from that of light defined via α T ≡ (c 2 GW − c 2 )/c 2 , as a function of 2M/r.We show α 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 α T (r) = 0 at spatial infinity and at the horizon, but displays nontrivial behaviour in-between with the overall amplitude and sign determined by A Ti while the radial dependence depends on the choice of i considered.n = 1 scalar hair profile (24), we will discuss how different profiles affect eventual constraints on α T in section IV D. Having parametrised δ φ , we turn our attention to the one remaining function of r affecting δ V , namely α T .To this end it will be useful to separate out the dependence on the scalar hair background profile and other geometric factors from α T in (17) as follows Here the dimensionless G T parameter has been defined to isolate the dependence of α 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 α T can tell us about scalar hair and vice versa. 15Because we have a nontrivial scalar profile, all the G i (and hence also G T ) are functions of r and so to fully specify the r-dependence of α T we finally also expand G T in powers of r along the same lines as discussed for δ φ above.Doing so we can write where each of the G Ti are constant coefficients.Putting everything together, i.e. substituting δ φ (24) and G T (26) into the expression for α T (25), we can finally write 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 α T and satisfy 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 α T , we illustrate the dependence of α T on 2M/r for various choices of the amplitude coefficients A Ti in figure 1.

C. Parametrized QNMs
Having parametrized all the functional freedom encoded within δ V above, it is now straightforward to combine the above expressions.Doing so we can express δ V as Note that we have dropped the order parameter ε 2 at this point.From this we can read off the a-coefficients defined in ( 22) 16  and from (23) we also obtain the following expression for the 16 For a given i the contributions to these coefficients are By using the 'element sign' ∈ 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.
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 I. Real and imaginary components of the e j 'basis' functions for ℓ = 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 ℓ up to ℓ = 10, together with the 'basis' for even-gravitational and even-scalar perturbations, see [111].
quasinormal frequencies where E 1 i has been defined for convenience, and its subscript 1 refers to n = 1.A set of E n i 'basis' functions for general n are provided in appendix C.
Anticipating some of our later discussion, we can already see that, for small ℓ and 2Mω 0 ∼ 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 δ ω, see appendix C and especially figure 4.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 | ≪ (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 T 0 ≪ 1.5. 17Including 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.

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 α T .We would now like to use this to forecast how well future GW experiments will be able to constrain α 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.

A. Fisher forecast setup
We begin by modelling the waveform as where h +,× represent the strain in the two polarisations of the gravitational wave.These are in principle functions of all coordinates, i.e. h +,× (t, r, θ , φ ).However, to distinguish between time and frequency domains, we will only make t (or ν) explicit and take the dependence on r, θ , φ as understood.F +,× 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 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 + ℓm , 18 18 The strain functions h +/× appear with different normalisation factors in the literature depending on the setup in question, e.g. with a factor of 1/2 √ 10π, an extra geometrical 3/4 for LISA or a 1 r factor [103,112,115,116].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 ρ, for which the appropriate detector-specific values will be discussed and used in the next section.and where f ℓm and τ ℓm characterise the real and imaginary parts of ω ℓm in the following way where Q ℓm is the quality factor.{A + ℓm , A × ℓm = A + ℓm N × , φ + ℓm , φ × ℓm } are the amplitudes and phases for the two polarisations.Finally, S ℓm are spheroidal functions carrying angular dependencies.Because modes with different (ℓ, m) do not mix due to the nature of our background, the ℓm indices in {ω ℓm , f ℓm , τ ℓm , Q ℓm , S ℓm , A + ℓm , φ + ℓm , φ × ℓm } will play no role for the time being, so we will obviate them to simplify our notation (and explicitly discuss which ℓm modes are of interest when this becomes relevant below).
Using the above strain functions, we compute the signal-tonoise-ratio (SNR) with the usual where S h (ν) is the noise spectral density characteristic of the detector and h(ν) is the Fourier transform of h(t). 19We now make use of the following set of simplifying assumptions: We also make use of the fact that we can approximate S h (ν) to be constant.For details on (and explicit checks of) these assumptions see [112].Using these assumptions, ( 35) can be re-expressed as20 To derive error estimates we make use of the Fisher Information Matrix, given by where θ a is the set of parameters for our theory and the noiseweighted product (•|•) is defined as Then, we can calculate the parameter errors by inverting the Fisher matrix (which gives the covariance matrix Σ).The error for a parameter a is given by As an initial estimate we will here study the simplified case where all the usual parameters of the waveform are known (A, φ + , ...) 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 ω (i.e.inside f and Q) 21 , 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 constraints 22 where the prime denotes a derivative with respect to A Ti , and for double-parameter constraints where again a prime denotes a derivative with respect to A Ti , and a dot represents a derivative with respect to A T j .This matches analogous expressions in [100] (albeit for non-α Trelated parameters there).Before deriving error estimates on different parameter combinations, let us briefly return to the question of which (ℓ, m) modes are of most interest. 23As discussed above, while the QNM spectrum does not depend on initial conditions, the amplitude of individual modes does.The dominant observable 21 Note that we need two variables to represent the real and imaginary parts of ω and, similarly to previous literature [100] we choose to work with the pair { f , Q} rather than { f , τ}.We note that, if working with the latter, the SNR equation ( 37) would be further simplified to 22 This simple expression also implicitly makes use of the large-Q limit (or equivalently large damping times τ).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]. 23There 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 |Im(ω)|.
contributions, i.e. the modes with the largest amplitudes for astrophysical binary compact object mergers, generically are the ℓ = m modes, more specifically the (2, 2) mode [112,[115][116][117][118][119]. 24Note 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 ℓ = 3) grows as the mass ratio q and angular momentum j of the remnant black hole increases [115,116,120,121] -also see those references for discussions related to the detectability of such modes.Nonetheless, the ℓ = 2 mode still generically dominates in all scenarios and higher ℓ modes decay more quickly, see table II.Note that the damping time τ goes as the inverse of the imaginary component of ω, which increases for higher ℓ 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 ℓ = 2 mode will be additionally enhanced relative to other modes [122][123][124].While it is straightforward to repeat the analysis for other higher ℓ modes, the above rationale truly singles out the ℓ = 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 → 0 recovers the Schwarzschild line element) so one expects that the nonrotating scenario still captures the leading order information in the quasinormal frequencies for sufficiently slow rotation.

B. Constraining A T 0
We begin by considering a minimal setup, where there is only a single relevant A Ti parameter, namely A T 0 .From (31) we then find the QNM shift to be given by where E 1 0 is shown in (31) and we quote it here for reference Substituting in the numerical values for the e j from table I, we obtain In evaluating this, we have also used the ℓ = 2 mode in table II.This now allows us to obtain parametric expressions for the α T -induced deviations in the QNM spectrum.From ( 46) and table II we find the following percentage differences for the real and imaginary parts, respectively Finally, we are also in a position to extract an expression for the accuracy with which an experiment with ringdown SNR ρ will be able to measure A T 0 .Reading off f and Q from ( 46), as defined in (34), and substituting them into the singleparameter error expression (42), we obtain an estimate on its detectability in the same fashion as [100]. 25This gives us26 The numerical value of this error calculation, as well as the analogous ones which will follow in this section, will be translated into α T constraints for specific detectors in section V. There, we will compare our results with other existing and forecasted constraints.

C. Constraining multiple A Ti
Having considered the single-parameter case above, a natural next step is to consider a more complex functional form for α T and hence for the A Ti .Here we consider the case where α T is controlled by two parameters, A T 0 as before and a second parameter A T 1 .Proceeding as before, we then have Reading off expressions for f and Q from equation ( 49) as before, we find the following error estimates from ( 43) 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 [65,118,130], this will of course allow constraining a correspondingly larger parameter space.Equation (50) shows that A T 0 and A T 1 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 T 0 and A T 1 are strongly correlated, as can be seen from the off-diagonal elements of the covariance matrix for these two parameters 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. 27Under standard matrix diagonalization procedures we obtain Σab = S −1 Σ ab S ∼ 10 3 5 0 0 303 .
This transformation amounts to identifying the combinations of A T 0 and A T 1 that yield uncorrelated parameters A TA and A T B , i.e. we have performed the parameter transformation (A T 0 , A T 1 ) → (A TA , A T B ) 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.84AT 0 − 0.54A T 1 and A T B = −0.54AT 0 + 0.84A T 1 . 28Finally, the errors for the new parameters are where we indeed see that A TA can be constrained more tightly than any parameter in the previous basis.

D. Dependence on scalar hair profile
Above we have derived expressions for the precision with which a generic future experiment with SNR ρ 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 27 We thank Sigurd Naess for related discussions. 28Equivalently, the matrix S is built with the eigenvectors of ( 51) Here we show forecasted errors on the strength of the interactions contributing to c GW as parametrised by G Ti (26), where we focus on G T 0 as an example.The corresponding error σ G T 0 is shown as a function of the scalar hair amplitude ϕ c and of the detector SNR as quantified by ρ.We see that σ G T 0 improves as the scalar hair amplitude grows and as ρ increases, as expected from ( 56).More concretely, for an SNR of ρ ∼ 10 x and a scalar hair amplitude of order ϕ c ∼ 10 y , we find σ G T 0 ∼ 10 2−x−2y .From table III, at lower frequencies for LISA we would therefore expect σ G T 0 ∼ 10 −3−2y constraints, while for LVK one would need y ≳ 8 to yield constraints on the underlying interactions there that are competitive with or stronger than existing bounds in this band.
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 Here ϕ 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 ϕ 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 ϕ 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 ϕ c as a free parameter. 29It is interesting, then, to disentangle the effect of ϕ 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 ϕ c the constraint on the G Ti becomes stronger.More explicitly 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 GW (in other words: in the event of a future detection of a c GW ̸ = c from another probe), and II) infer a constraint of the G Ti , given other independent information about the amplitude ϕ 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 T 0 case.The corrections to the quasinormal frequencies are now given by where for n = 2 the basis E n i (C7) becomes leading to the following expression for the error on A T 0 We see that the precision with which A T 0 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 T 0 and A T 1 .
29 ϕ 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].
Detector   (27), as an example.The error on α T , σ α T , is quoted as one order-of-magnitude less than the corresponding error on A T 0 , as observed in figure 1.We stress that the precise mapping of underlying amplitude parameters to α 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 α T are qualitatively similar -see e.g.table IV.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 ∼ O( 103 ) [146,147], while [145] forecast SNRs up to ∼ 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.

V. FORECASTING OBSERVATIONAL CONSTRAINTS
In the previous section we derived parametric expressions for the precision with which the parameters controlling the behaviour of α T and hence c GW can be measured for probes with a general ringdown SNR ρ.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 III. 30  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 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 30 It 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 [115,117,[148][149][150][151][152] for discussions on this.
faithful description of the relevant physics.Since (2) gives rise to a frequency-independent c GW , we are therefore implicitly assuming that at the very least in the frequency-window spanning the ringdown frequencies in question, c 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 GW effectively becomes a constant c GW ̸ = 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 GW = c (0) GW 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 GW = c at high frequencies due to the Lorentz invariant nature of the UV completion.The frequency-dependence in c 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 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 frequencydependence of c 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 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 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 GW and hence α 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 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 [30,154] -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]. 31In 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 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 GW within the band under investigation (i.e. the LISA forecasts assume a frequency-independent c GW in the LISA band and so on).

A. LISA band forecasts
As motivated above and in the introduction, the LISA band is particularly promising in terms of testing for deviations of c GW from the speed of light, given that a frequency-dependent transition of c 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 III 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 ∼ O(10 5 ) [137,145].Mapping this back to α T itself (27), this implies one will be able to detect deviations down to the α T ∼ O(10 −4 ) level from LISA band ringdown alone in the context of the models we consider, 32 i.e.
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 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 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 σ LISA/TianQin α T ∼ 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.

B. 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 GW specific to this frequency band.From measuring the coincidence of the GW170817 signal in GW and optical counterpart observations, one finds that α T ≲ 10 −15 [2][3][4][5][6].When even a very mild frequency-dependence of c GW is present in the LVK band, this bound can be strengthened to α T ≲ 10 −17 [30].Contrast this with bounds from ringdown observations alone, where the A Ti can be constrained at the ∼ 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 III and see [138,139] and [140,141] for ET and CE, respectively.Again mapping this to constraints on α T itself, we therefore ultimately expect 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 O(10 4 − 10 5 ) events with a ringdown SNR of 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 σ ET/CE α T ∼ 10 −3 .In the LVK context it is also interesting to point out that existing (non-ringdown-specific) constraints on c GW from GW waveforms in the LVK band have already seen similar improvements by stacking events.More specifically, when comparing I) constraints on c 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 GW ), this improves these bounds on c GW by around two orders of magnitude. 33t first sight (61) as well as the improved σ ET/CE α T ∼ 10 −3 bound reachable by stacking events are rather weak, albeit complementary, constraints on c 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 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 α T from the beginning while not precluding a more sizeable α 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 D -we show this in figure 2.More specifically, from (56) and in the event of a scalar hair amplitude ϕ c ∼ 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 ϕ c .Reversing the argument, if future observations were to identify a small but non-zero cosmological α 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 Λ.In cosmology this scale is typically chosen to be Λ = Λ 3 ≡ (M Pl H 2 0 ) 1/3 , where this choice ensures those interactions give O(1) contributions to cosmology.However, if a different Λ is chosen, the cosmological α T (and hence G T ) scales as α T ∼ (Λ 3 /Λ) 6 .So raising the interaction scale Λ by just three orders of magnitude suppresses the cosmological α T down to a level of 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 ∼ O(10 −18 ), while future ringdown constraints from ET/CE along the lines investigated here do not yield evidence for a non-zero α T .From (28) this would allow us to derive a bound on the scalar hair ϕ c ≲ O(10 9 ) for frequencies in the LVK band.Note that other complementary bounds on ϕ c may be obtainable e.g. from considering even perturbations or going beyond linear theory.

C. 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 GW in the LISA band by pointing out that a frequency-dependent transition from a nearly constant c GW = c at LVK frequencies to a different low-frequency c 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 α T down to a level of O(10 −15 ) [30,58].In addition, once AEDGE/DECIGO observations are available, direct constraints on c GW from this band will be obtainable in analogy to the LVK/LISA analyses discussed above.Whenever there is significant frequency-dependence for c 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 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 ∼ O(10 −1 ) level, c.f. table III.Mapping this to constraints on α T itself, as before, this implies 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 O(10 3 ) should be observable with AEDGE per year [144], so optimistically an improvement up to σ AEDGE/DECIGO α T ∼ 10 −4 appears feasible eventually.

VI. CONCLUSIONS
In this paper we have investigated how the speed of gravitational waves c 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 GW from the speed of light only affect the QNMs in the presence of a non-trivial scalar hair profile φ = φ (r) in agreement with the results of [100].Any deviations from c 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 GW at the O(10 −4 ) level.For AEDGE/DECIGO the equivalent precision will be 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 GW e.g. from GW170817, the importantly probe different frequency ranges.This is particularly relevant in the context of testing c GW , given large classes of dark energy models naturally give rise to a frequency-dependent transition in c GW below the LVK band.
• With ringdown constraints we are testing the effect of deviations from c GW = c on a different background solution than that relevant for GW propagation constraints on c GW (black hole vs. cosmological space-times).The precise dependence of c 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 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 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 GW for various detectors throughout the O(10 −4 ) − O( 103 ) 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 GW ̸ = 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 GW is approximately constant in specific frequency win-dows/bands when forecasting constraints for those respective bands.A more comprehensive analysis, extending the present work, investigating degeneracies and constraints in higherdimensional 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 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 O(1) level -see e.g.[31,32] and references therein -which are expected to improve to O(10 −1 ) in the near future [38].While we have not mandated a specific sign for any potential deviation of c 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 GW ≥ c for large classes of models [171,172].We close by re-emphasising that we have mostly focused on investigating how well c 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.
where the potential perturbations are given by 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 δ A 1,2 = δC 1,2 = G 4φ = G 4φ φ = 0.However, we see that δ A 1 , δC 1 would contribute at lower order (as well as δ φ 1 , if G 4φ ̸ = 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 ε is very small, such lower-order-in-ε contributions would be expected to dominate over any α 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 α 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.δ A 1 = δ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φ ̸ = 0 ̸ = G 4φ φ .It is interesting to point out that interactions contributing to G 4φ in our background can be removed with a conformal transformation of the metric.This is because we still have X = 0, 34 meaning that all terms in G 4φ are X-independent or, in other words, those interactions are of the form f (φ )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φ = 0 = G 4φ φ .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φ = 0 = G 4φ φ as an additional simplifying assumption.Note that this assumption is trivially satisfied in theories where the scalar φ is endowed with a shift symmetry φ → φ + c.
To conclude this appendix, let us briefly consider the case where one or many of G 4φ , G 4φ φ , δ A 2 , δ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.quadratic action is given by with the coefficients being where again the bar denotes that the āi are evaluated on the background and, to avoid clutter, bars are implied everywhere on the right hand side.ε A,B are contributions that vanish onshell 35 and As shown in [101,102], this action can be rewritten to make the presence of only one degree of freedom explicit where the potential is given by 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 δ φ = ϕ c 2M r n , (C1) 35 Expressions for ε A,B , which are obtained by varying the action with respect to A and B, are fully provided in [102].as a function of 2M/r.i = 0 has been set such that A T 0 is the only non-zero parameter (and fixed to a fiducial value of 1).We see that α T (r) = 0 at spatial infinity and at the horizon, but again observe non-trivial behaviour in the intermediate region.The size of α T (partially controlled by A T 0 ) is considerably enhanced by increasing n.This is mainly due to the factor n 2 accompanying all δ ω (as can be seen from ( C6) and (C7)).
which means α T is now given by where the new n-dependence is plotted in figure 3. Substituting this back into δ V (21) we get δ V = n 2M 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 n i .In diagram 4 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 α T ) is straightforwardly repeated in the same fashion as shown in section IV.We show in table IV 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 IV is that errors decrease for increasing i and n (at least for the single-parameter cases).The last column (i : 0, 1) displays the errors for the two parameters (A TA , A T B ), obtained from (A T 0 , A T 1 ) via the same diagonalization procedure as shown in section IV C.
FIG.1.Here we plot α T , the deviation of the speed of gravitational waves from that of light defined via α T ≡ (c 2 GW − c 2 )/c 2 , as a function of 2M/r.We show α 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 α T (r) = 0 at spatial infinity and at the horizon, but displays nontrivial behaviour in-between with the overall amplitude and sign determined by A Ti while the radial dependence depends on the choice of i considered.

E 1 0 6 +
= (2Mω 0 ) 2 e 4 − (ℓ(ℓ + 1) − 9)e FIG. 2.Here we show forecasted errors on the strength of the interactions contributing to c GW as parametrised by G Ti(26), where we focus on G T 0 as an example.The corresponding error σ G T 0 is shown as a function of the scalar hair amplitude ϕ c and of the detector SNR as quantified by ρ.We see that σ G T 0 improves as the scalar hair amplitude grows and as ρ increases, as expected from(56).More concretely, for an SNR of ρ ∼ 10 x and a scalar hair amplitude of order ϕ c ∼ 10 y , we find σ G T 0 ∼ 10 2−x−2y .From table III, at lower frequencies for LISA we would therefore expect σ G T 0 ∼ 10 −3−2y constraints, while for LVK one would need y ≳ 8 to yield constraints on the underlying interactions there that are competitive with or stronger than existing bounds in this band.

4 FIG. 3 .
FIG.3.Here we show α T (C2) for different choices of n in (C1) as a function of 2M/r.i = 0 has been set such that A T 0 is the only non-zero parameter (and fixed to a fiducial value of 1).We see that α T (r) = 0 at spatial infinity and at the horizon, but again observe non-trivial behaviour in the intermediate region.The size of α T (partially controlled by A T 0 ) is considerably enhanced by increasing n.This is mainly due to the factor n 2 accompanying all δ ω (as can be seen from (C6) and (C7)).
IV. Values for σ A Ti ρ for a scalar profile δ φ = ϕ c 2M r n .

TABLE III .
Achievable order-of-magnitude ringdown SNRs for a single observed event for different GW detectors and the corresponding order-of-magnitude errors on α T .Errors in this table are computed assuming A T 0 is the only amplitude parameter contributing to α T in