# A parametrized quasi-normal mode framework for modified Teukolsky equations

Pablo A. Cano<EMAIL_ADDRESS>Departament de Física Quàntica i
Astrofísica, Institut de Ciències del Cosmos
Universitat de Barcelona, Martí i Franquès 1, E-08028 Barcelona, Spain
Lodovico Capuano<EMAIL_ADDRESS>SISSA, Via Bonomea 265, 34136 Trieste,
Italy and INFN Sezione di Trieste  IFPU - Institute for Fundamental Physics of
the Universe, Via Beirut 2, 34014 Trieste, Italy  Nicola Franchini
<EMAIL_ADDRESS>Université Paris Cité, CNRS, Astroparticule et
Cosmologie, F-75013 Paris, France  CNRS-UCB International Research Laboratory,
Centre Pierre Binétruy, IRL2007, CPB-IN2P3, Berkeley, CA 94720, USA  Simon
Maenaut<EMAIL_ADDRESS>Institute for Theoretical Physics, KU
Leuven. Celestijnenlaan 200D, B-3001 Leuven, Belgium  Sebastian H. Völkel
<EMAIL_ADDRESS>Max Planck Institute for Gravitational Physics
(Albert Einstein Institute), D-14476 Potsdam, Germany

###### Abstract

Modifications to general relativity lead to effects in the spectrum of quasi-
normal modes of black holes. In this paper, we develop a parametrized
formalism to describe deviations from general relativity in the Teukolsky
equation, which governs linear perturbations of spinning black holes. We do
this by introducing a correction to the effective potential of the Teukolsky
equation in the form of a $1/r$ expansion controlled by free parameters. The
method assumes that a small deviation in the effective potential induces a
small modification in the spectrum of modes and in the angular separation
constants. We isolate and compute the universal linear contribution to the
quasi-normal mode frequencies and separation constants in a set of
coefficients, and test them against known examples in the literature (massive
scalar field, Dudley-Finley equation and higher-derivative gravity). We make
the coefficients publicly available for relevant overtone, angular momentum
and azimuthal numbers of modes and different values of the black hole spin.

##  I Introduction

Gravitational wave astronomy successfully observed more than one hundred
binary black hole (BH) mergers [1]. Such an advancement in the field allows
one to make precision tests of general relativity (GR). In particular, BH
spectroscopy, i.e., the identification of the infinite tower of quasi-normal
modes (QNMs) of which the signal is composed in the linear post-merger phase
known as the ringdown [2, 3, 4], has finally been applied to tens of events
[5, 6, 7]. The detection of two modes simultaneously, despite being
controversial for the first three observing runs [8, 9, 10, 11, 12, 13, 14,
15, 16], is expected to be effective for the current O4a run and future ones.

Detecting a second mode is a crucial ingredient for tests of GR. Due to no-
hair theorems, the QNM spectrum of a Kerr BH depends uniquely on its mass and
spin. If only one mode is detected, this can always be fitted to a QNM
frequency of a Kerr BH with a certain mass and angular momentum. However, the
measurement of any additional modes provides a consistency test of the Kerr
QNM spectrum and hence would allow us to spot deviations from GR. This clear
identification of possible beyond-GR effects makes BH spectroscopy one of the
most promising ways to test GR.

Currently, ringdown tests employ blind deviations from GR in the frequencies
[7], or agnostic deviations constructed assuming small-coupling and slow-spin
parametrization [17, 18, 19]. On the other hand, theory-specific tests are
limited to a handful of cases [20]. This is because the computation of QNMs
for rotating solutions beyond-GR is incomplete. The main difficulties arise
from: absence of analytic background solution, non-separability of the
perturbation equations, additional fields coupled to the metric, different
boundary conditions [4]. It turns out that for a vast class of theories, the
first problem can be solved performing a double simultaneous expansion in the
spin and in the coupling constant of the theory [21, 22, 23, 24, 25]. Then,
one can choose whether to study perturbations giving priority to slow-spin or
small coupling. The former has the advantage of being possible for metric
perturbations, for which couplings between different fields are tractable, at
the cost of not being able to predict precisely QNMs at high spins (which are
relevant for astrophysical purposes) [26, 27, 28, 29, 30, 31, 32].

On the other hand, by assuming small-coupling for the perturbations, one can
work out a modified Teukolsky equation, which is, in principle, reliable at
any spin [33, 34, 35]. The disadvantage comes from the construction itself of
the Teukolsky equation, which is based on curvature perturbations, and once
one has some perturbations of perturbations, as in the case outlined here,
metric reconstruction becomes necessary. This feature strongly hinders one
from going beyond the first order in the coupling expansion.

The necessity of having reliable QNMs at high spins and the fact that
observations seem to narrow down the size of deviations from GR, make the
modified Teukolsky framework preferable for the study of QNMs beyond GR. The
scope of this paper is to develop a general formalism for the quick
computation of QNMs in theories which have a perturbative departure from GR.
The framework is based on an assumption similar to one developed in spherical
symmetry [36, 37, 38, 39]. The advantage of this formalism is that it can be
used also for the inverse problem, i.e., if a modification to GR is detected,
one wants to be able to reconstruct the potential, the metric, or even the
action from which the deviation originated from [38, 40, 41]. In general, this
framework opens the path to the development of a theory-informed description
of QNM agnostic deviations.

The structure of the paper is as follows: we first introduce the modified
Teukolsky equation and the formalism in section II; then we show how to
numerically compute the coefficients with the continued fraction method and
their regime of validity in section III; the formalism has already a few
notable applications, which we use as a check for our computations in section
IV; finally, we outline our conclusions in section V. Throughout the paper we
work with mostly plus signature for the metric, and geometric units
$G_{\text{\tiny{N}}}=c=1$. We also assume that the background, unmodified
metric is a Kerr BH of mass $M=1/2$ and spin $a$. It is worth noting that due
to the choice of units, the spin parameter $a$ ranges between $0$ and $1/2$.

##  II Parametrized formalism

###  II.1 The linear coefficients

Let us start from the radial Teukolsky equation for a spin $s$ field as [42]

$\frac{1}{\Delta^{s}R(r)}\frac{\mathrm{d}}{\mathrm{d}r}\left[\Delta^{s+1}R^{\prime}(r)\right]+V(r)=0\,,$
(1)

where the effective potential reads

$V(r)=2\mathrm{i}s\frac{\mathrm{d}K}{\mathrm{d}r}-\lambda_{\ell
m}+\frac{1}{\Delta}\left(K^{2}-\mathrm{i}sK\frac{\mathrm{d}\Delta}{\mathrm{d}r}\right)\,,$
(2)

and we defined

$\displaystyle\Delta$ $\displaystyle=r^{2}-r+a^{2}\,,\qquad
K=(r^{2}+a^{2})\omega-am\,,$ (3) $\displaystyle\lambda_{\ell m}$
$\displaystyle=B_{\ell m}+a^{2}\omega^{2}-2am\omega\,.$ (4)

It is worth mentioning that the zeros of the function $\Delta$ determine the
location of the BH inner and outer horizons, given by

$r_{\pm}=\frac{1\pm\sqrt{1-4a^{2}}}{2}\,.$ (5)

On the other hand, the angular equation reads [42]

$\begin{split}\frac{1}{S(y)}&\frac{\mathrm{d}}{\mathrm{d}y}\left[(1-y^{2})S^{\prime}(y)\right]+a^{2}\omega^{2}y^{2}\\\
&-2sa\omega y+B_{\ell m}+s-\frac{(m+sy)^{2}}{1-y^{2}}=0\,,\end{split}$ (6)

where $y=\cos\theta$.

We now assume that for a modified theory of gravity whose modifications are
small with respect to GR the equation governing radial perturbations is the
Teukolsky radial equation plus a correction to the potential linear in the
coupling constants

$\frac{1}{\Delta^{s}R(r)}\frac{\mathrm{d}}{\mathrm{d}r}\left[\Delta^{s+1}R^{\prime}(r)\right]+V(r)+\delta
V(r)=0\,,$ (7)

and we assume that the modification is expanded in powers of $r$

$\delta
V(r)=\frac{1}{\Delta}\sum_{k=-K}^{4}\alpha^{(k)}\left(\frac{r}{r_{+}}\right)^{k}\,,$
(8)

where $K$ is the largest negative coefficient of the power series and
$\alpha^{(k)}$ are dimensionful coefficients we assume to be small. This
assumption is justified by the recent developments in obtaining modified
Teukolsky equations by assuming small coupling corrections to GR [33, 34, 35].

On the other hand, we can also assume that the angular equation remains
unchanged. This is due to the fact that the spheroidal harmonics are a
complete basis of the 2-sphere angular variables, and it can be shown that if
one forces an angular expansion of the Weyl scalars in spin-weighted
spheroidal harmonics, all the mixing terms would enter at second order in the
coupling constants [35, 43]. Nevertheless, a modification on the QNM
frequencies will induce a modification to the separation constant $B_{\ell
m}$, hence, we also need to include equation (6) in our analysis.

If we assume the couplings to be small, we are allowed to perform a Taylor
expansion of QNMs and separation constants around their GR values [36, 37].
Hence, we can write

$\begin{split}\omega_{n\ell m}&\simeq\omega_{n\ell
m}^{0}+\sum_{k}d_{\omega,n\ell m}^{(k)}\alpha^{(k)}\,,\\\ B_{\ell
m}(a\omega)&\simeq B_{\ell m}^{0}(a\omega)+\sum_{k}d_{B,\ell
m}^{(k)}\alpha^{(k)}\,.\end{split}$ (9)

In the following steps, we omit the indices $s,n,\ell,m,a,\omega$ for clarity.
The linear coefficients $d_{\omega}$ and $d_{B}$ can be identified as the
derivatives of $\omega$ and $B$ with respect to the single coupling $\alpha$.
To compute them, we perform the following steps. In the GR limit, one finds
$\omega$ and $B_{\ell m}$ as simultaneous roots of two functions constructed
from the radial and the angular equation

$\displaystyle\mathcal{L}_{r}\left[\omega,B_{\ell
m}\right]=0\,,\qquad\mathcal{L}_{\theta}\left[\omega,B_{\ell m}\right]=0\,.$
(10)

The exact form of these functions depends on the chosen numerical method. For
non-zero modifications, we can perform a Taylor expansion of the two functions
around $\alpha=0$

$\left.\mathcal{L}_{j}\right|_{\text{\tiny{GR}}}+\alpha\left.\frac{\mathrm{d}\mathcal{L}_{j}}{\mathrm{d}\alpha}\right|_{\text{\tiny{GR}}}+\mathcal{O}(\alpha)^{2}=0\,,$
(11)

where $j=[r,\theta]$ and we evaluate the derivative around their GR value
$\omega=\omega^{0}$, $B=B^{0}$, $\alpha=0$. By requiring that equations (11)
is satisfied at each order in $\alpha$, and expanding the derivative by chain
rule, we obtain

$\begin{split}\left.\frac{\partial\mathcal{L}_{r}}{\partial\alpha}+\frac{\partial\mathcal{L}_{r}}{\partial\omega}\,d_{\omega}+\frac{\partial\mathcal{L}_{r}}{\partial
B}\,d_{B}\right|_{\text{\tiny{GR}}}&=0\,,\\\
\left.\frac{\partial\mathcal{L}_{\theta}}{\partial\omega}\,d_{\omega}+\frac{\partial\mathcal{L}_{\theta}}{\partial
B}\,d_{B}\right|_{\text{\tiny{GR}}}&=0\,,\end{split}$ (12)

where we identified $d_{\omega}$ and $d_{B}$ from their definition in equation
(9). By solving the conditions above for $d_{\omega}$ and $d_{B}$, we get

$\begin{split}&d_{\omega}=\left.-\frac{\partial\mathcal{L}_{r}}{\partial\alpha}\,\frac{\partial\mathcal{L}_{\theta}}{\partial
B}\,\left(\frac{\partial\mathcal{L}_{r}}{\partial\omega}\frac{\partial\mathcal{L}_{\theta}}{\partial
B}-\frac{\partial\mathcal{L}_{r}}{\partial
B}\frac{\partial\mathcal{L}_{\theta}}{\partial\omega}\right)^{-1}\right|_{\text{\tiny{GR}}}\,,\\\
&d_{B}=\left.\frac{\partial\mathcal{L}_{r}}{\partial\alpha}\,\frac{\partial\mathcal{L}_{\theta}}{\partial\omega}\,\left(\frac{\partial\mathcal{L}_{r}}{\partial\omega}\frac{\partial\mathcal{L}_{\theta}}{\partial
B}-\frac{\partial\mathcal{L}_{r}}{\partial
B}\frac{\partial\mathcal{L}_{\theta}}{\partial\omega}\right)^{-1}\right|_{\text{\tiny{GR}}}\,.\end{split}$
(13)

In section III we show how to numerically define the functions
$\mathcal{L}_{j}$ with Leaver’s continued fraction method.

###  II.2 Maximum number of independent coefficients

In reference [44], Kimura realised that in the case of spherically symmetric
perturbations, there is always an ambiguity in defining the modified
potential, upon a free reparametrization of the field. The same reasoning can
be applied to the Teukolsky equation as well. If we perform the following
transformation in equation (7)

$R(r)\rightarrow\left[1+\varepsilon X(r)\right]R(r)+\varepsilon\Delta
Y(r)R^{\prime}(r)\,,$ (14)

assuming that $\varepsilon\ll 1$, then the equation that $R(r)$ solves is

$\frac{1}{\Delta^{s}R}\frac{\mathrm{d}}{\mathrm{d}r}\left[\Delta^{s+1}R^{\prime}\right]+V+\delta
V+\delta\overline{V}+\delta W\frac{R^{\prime}}{R}=0\,.$ (15)

By imposing $\delta W=0$ we uniquely obtain the free function $X(r)$ as

$X(r)=c+\frac{s}{2}\Delta^{\prime}Y-\frac{1}{2}\Delta Y^{\prime}\,,$ (16)

which yields the "ambiguous" potential in the form

$\begin{split}\delta\overline{V}=\varepsilon\Delta\Bigg{[}&Y^{\prime}\left(\frac{\left(s^{2}-1\right)\left(\Delta^{\prime}\right)^{2}}{2\Delta}+2s-2V-1\right)\\\
&+Y\left(\frac{s(s+1)\Delta^{\prime}}{\Delta}-V^{\prime}-\frac{V\Delta^{\prime}}{\Delta}\right)\\\
&-\frac{1}{2}\Delta
Y^{(3)}-\frac{3\Delta^{\prime}Y^{\prime\prime}}{2}\Bigg{]}\,.\end{split}$ (17)

Now, upon suitable choice of the function $Y$, we can express
$\delta\overline{V}$ in the $r$ basis. It turns out that the ansatz
$Y=Y_{j}=y_{j}(r_{+}/r)^{j}$ yields

$\begin{split}\delta\overline{V}=\frac{\varepsilon
y_{j}}{\Delta}\sum_{k=-3}^{5}r_{+}^{k}\overline{A}_{j}^{(k)}\left(\frac{r}{r_{+}}\right)^{k-j}\,,\end{split}$
(18)

which implies that $j\geq 1$ since the maximum power of $r$ in $V$ is $r^{4}$,
and the full expression of $\overline{A}_{j}^{(k)}$ can be found in appendix
A. In general, one can take a linear combination of the free functions $Y_{j}$
and still get a potential that is equivalent to the starting one. Each term of
this linear combination contains the free parameter $y_{j}$, which can be used
to set to $0$ one of the terms $\alpha^{(k)}$ in equation (8). This reasoning
allows us to fix the negative limit in the power expansion to be $K=3$.

It is possible that by choosing a different ansatz for $Y$ one could further
reduce the number of coefficients in the equation. In fact, for the case of
study of higher derivative gravity that we treat in another publication [45]
the number of independent coefficients reduces to four (being $k=[-2,0,1,2]$)
— see also [46]. Although we could not prove this is a general feature of
arbitrary modifications of the Teukolsky equation, we suspect that it was
possible in that case thanks to the expansion in the spin assumed for every
coefficient. Indeed, we believe that the ansatz for $Y$ that would reduce the
potential to the lowest number of terms would be, perhaps, a rational function
involving powers of $a$ and $r$. To date, we could not find such reduction.

##  III Computation of the coefficients: the continued fraction method

###  III.1 Continued fractions for the Teukolsky equation

We start here by recalling the Leaver method to compute the frequencies and
the separation constant for a Kerr spacetime. The first step to find a
continued fraction expansion is to assume an ansatz for the wavefunctions. Let
us start from the radial equation, where we assume the following ansatz [47]

$R(r)=f^{-\mathrm{i}\sigma-s}(r-r_{-})^{p-1-2s}\mathrm{e}^{qr}\sum_{n=0}^{N}R_{n}f^{n}\,,$
(19)

where $r_{\pm}$ are the zeros of $\Delta$,

$f=\frac{r-r_{+}}{r-r_{-}}\,,$ (20)

and we defined $p=q=\mathrm{i}\omega$ and

$\displaystyle\sigma=\sigma_{\text{\tiny{GR}}}\equiv\frac{r_{+}(\omega-\omega_{c})}{r_{+}-r_{-}}\,,\qquad\omega_{c}=\frac{am}{r_{+}}\,,$
(21) $\displaystyle
r_{\pm}=\frac{1}{2}(1\pm\beta)\,,\qquad\beta=\sqrt{1-4a^{2}}\,.$ (22)

With these definitions, the equation (1) takes the form

$\sum_{n=0}^{N}R_{n}\left(\frac{\alpha^{r}_{n-1}}{f}+\beta^{r}_{n}+\gamma^{r}_{n+1}f\right)f^{n}=0\,,$
(23)

where the coefficients are

$\displaystyle\alpha^{r}_{n}=$
$\displaystyle\,(n+1)\left(n+1-s-2\mathrm{i}\sigma\right)$ (24)
$\displaystyle\beta^{r}_{n}=$
$\displaystyle\,2n(2\mathrm{i}\sigma+p+q\beta-1)-2n^{2}-1-s-B_{\ell m}$
$\displaystyle+q\left(a^{2}q+\beta+s\right)-2i\sigma(p+\beta q+i\omega-1)$
$\displaystyle-p(\beta q+q+s-1)$ (25) $\displaystyle\gamma^{r}_{n}=$
$\displaystyle\,\left(n-p-\mathrm{i}\omega\right)\left(n+s-p-2\mathrm{i}\sigma+\mathrm{i}\omega\right)\,.$
(26)

The equation is satisfied when each term proportional to a power of $f$
vanishes

$\displaystyle\beta^{r}_{0}R_{0}+\alpha^{r}_{0}R_{1}$ $\displaystyle=0\,,$
(27)
$\displaystyle\gamma^{r}_{n}R_{n-1}+\beta^{r}_{n}R_{n}+\alpha^{r}_{n}R_{n+1}$
$\displaystyle=0\quad\text{for }n\geq 1\,.$ (28)

The path for the angular equation is similar. We define an ansatz to be finite
at the regular singular points $y=\pm 1$ [47]

$S(y)=(1+y)^{k_{1}}(1-y)^{k_{2}}\mathrm{e}^{a\omega
y}\sum_{n=0}^{N}S_{n}(1+y)^{n}\,,$ (29)

where $k_{1}=|m-s|/2$ and $k_{2}=|m+s|/2$. We can obtain a similar recurrence
relation by inserting this ansatz into equation (6), which, with an analogous
reshuffling, reads

$\displaystyle\beta^{\theta}_{0}S_{0}+\alpha^{\theta}_{0}S_{1}$
$\displaystyle=0\,,$ (30)
$\displaystyle\gamma^{\theta}_{n}S_{n-1}+\beta^{\theta}_{n}S_{n}+\alpha^{\theta}_{n}S_{n+1}$
$\displaystyle=0\quad\text{for }n\geq 1\,,$ (31)

where the coefficients are

$\displaystyle\alpha^{\theta}_{n}=$
$\displaystyle\,-2(n+1)\left(n+1+2k_{1}\right)\,,$ (32)
$\displaystyle\beta^{\theta}_{n}=$
$\displaystyle\,n(n-1)+2n(k_{1}+k_{2}+1-2a\omega)$
$\displaystyle-2a\omega(2k_{1}+s+1)+(k_{1}+k_{2})(k_{1}+k_{2}+1)$
$\displaystyle-a^{2}\omega^{2}-s(s+1)-B_{\ell m}\,,$ (33)
$\displaystyle\gamma^{\theta}_{n}=$
$\displaystyle\,2a\omega(n+k_{1}+k_{2}+s)\,.$ (34)

To invert the relation we can define the ladder operators which have the
following property $R_{n+1}=-\Lambda^{r}_{n}R_{n}$ and
$S_{n+1}=-\Lambda^{\theta}_{n}S_{n}$ as (the superscript $r/\theta$ is omitted
for clarity)

$\Lambda_{n}=\frac{\gamma_{n+1}}{\beta_{n+1}-\alpha_{n+1}\Lambda_{n+1}}\,.$
(35)

By initializing $\Lambda_{N}$ according to the Nollert expansion (explained in
detail in appendix B), the equations one needs to solve simultaneously to
obtain the eigenfrequency $\omega$ and the separation constant $B_{\ell m}$
are

$\displaystyle\mathcal{L}_{r}=\Lambda_{1}^{r}\alpha^{r}_{0}-\beta^{r}_{0}$
$\displaystyle=0\,,$ (36)
$\displaystyle\mathcal{L}_{\theta}=\Lambda_{1}^{\theta}\alpha^{\theta}_{0}-\beta^{\theta}_{0}$
$\displaystyle=0\,,$ (37)

which are nothing but (27) and (30).

###  III.2 Continued fraction beyond Teukolsky

We now turn our attention to the modified Teukolsky equation. It is always
possible to bring equation (8) into the following form111In section II.2 we
showed that $K=3$, but the following analysis works, in principle, for any
value of $K$, hence we keep it unspecified.

$\begin{split}\delta
V(r)=&\,\frac{A^{(0)}}{\Delta}+\frac{A^{(1)}}{r_{+}(r-r_{-})}+\frac{1}{r_{+}^{2}}\sum_{k=0}^{2}\widetilde{\alpha}^{(k)}\left(\frac{r}{r_{+}}\right)^{k}\\\
&+\frac{1}{\Delta}\sum_{k=1}^{K}\alpha^{(-k)}\left(\frac{r_{+}}{r}\right)^{k}\,,\end{split}$
(38)

where $A^{(0)}$, $A^{(1)}$ and $\widetilde{\alpha}^{(k)}$ are constants that
can be obtained from the constants $\alpha^{(k)}$ appearing in equation (8),
as explained in appendix C. First of all, we notice that the terms multiplied
by $1/\Delta$ modify the behaviour of the equation at the horizon. In order to
take into account of these additional terms, we need to modify the definition
of the exponent $\sigma$ appearing in the ansatz (19). By requesting that the
solution is regular at the horizon, we must replace the value of $\sigma$ into

$\sigma=\frac{\mathrm{i}s}{2}+\sqrt{\left(\sigma_{\text{\tiny{GR}}}-\frac{\mathrm{i}s}{2}\right)^{2}+\frac{1}{\beta^{2}}\sum_{k=-K}^{4}\alpha^{(k)}}\,,$
(39)

where we took the positive sign of the square root in order to obtain the
correct GR limit. On the other hand, the terms $\widetilde{\alpha}^{(1)}$ and
$\widetilde{\alpha}^{(2)}$ modify the behaviour at infinity of the equation.
This leads to a modification of the values of $p$ and $q$ into

$\displaystyle
q=\pm\sqrt{-\frac{\widetilde{\alpha}^{(2)}}{r_{+}^{4}}-\omega^{2}}\,,$ (40)
$\displaystyle
p=-\frac{r_{+}\widetilde{\alpha}^{(1)}+\widetilde{\alpha}^{(2)}-2r_{+}^{4}\left(qs-\mathrm{i}s\omega-\omega^{2}\right)}{2qr_{+}^{4}}\,,$
(41)

where the sign of $q$ is chosen such that $\text{Re}(q)>0$. This asymptotic
behaviour is the reason why we truncate the series in equation (8) at $k=4$.
By repeating the steps done for the GR case, we obtain a modified version of
equation (23)

$\begin{split}\sum_{n=0}^{N}&R_{n}\Bigg{[}\frac{\alpha^{\text{\tiny{bg}}}_{n-1}}{f}+\beta^{\text{\tiny{bg}}}_{n}+\gamma_{n+1}^{\text{\tiny{bg}}}f\\\
&+\frac{1}{\beta^{2}}\frac{(1-f)^{2}}{f}\sum_{k=1}^{K}\alpha^{(-k)}\left(\frac{1-f}{1-\eta
f}\right)^{k}\Bigg{]}f^{n}=0\,,\end{split}$ (42)

where $\eta=r_{-}/r_{+}$ and

$\displaystyle\alpha_{n}^{\text{\tiny{bg}}}=$
$\displaystyle\,\alpha^{r}_{n}-\frac{1}{\beta^{2}}\sum_{k=1}^{K}\alpha^{(-k)}\,,$
(43) $\displaystyle\beta_{n}^{\text{\tiny{bg}}}=$
$\displaystyle\,\beta^{r}_{n}-2\left(\sigma_{\text{\tiny{GR}}}-\sigma\right)\left(\sigma_{\text{\tiny{GR}}}+\sigma-\omega\right)\,,$
$\displaystyle-2\frac{A^{(0)}}{\beta^{2}}+\frac{A^{(1)}}{r_{+}\beta}+\frac{\widetilde{\alpha}^{(0)}}{r_{+}^{2}}$
(44) $\displaystyle\gamma_{n}^{\text{\tiny{bg}}}=$
$\displaystyle\,\gamma^{r}_{n}+2\mathrm{i}\left(\sigma_{\text{\tiny{GR}}}-\sigma\right)(s+\mathrm{i}\omega)$
$\displaystyle-\frac{1}{\beta^{2}}\sum_{k=1}^{K}\alpha^{(-k)}-\frac{A^{(1)}}{r_{+}\beta}\,,$
(45)

and we used the fact that

$\Delta=\beta^{2}\frac{f}{(1-f)^{2}}\,,\qquad r-r_{-}=\frac{\beta}{1-f}\,.$
(46)

If we fix a single modification $k$, we can get rid of the rational behaviour
in $f$ by multiplying the equation by $(1-\eta f)^{k}$, obtaining the
following expression

$\begin{split}\sum_{n=0}^{N}&R_{n}\Bigg{[}\left(\frac{\alpha^{\text{\tiny{bg}}}_{n-1}}{f}+\beta^{\text{\tiny{bg}}}_{n}+\gamma^{\text{\tiny{bg}}}_{n+1}f\right)\left(1-\eta
f\right)^{k}\\\
&+\frac{\alpha^{(-k)}}{\beta^{2}}\frac{(1-f)^{k+2}}{f}\Bigg{]}f^{n}=0\,.\end{split}$
(47)

We can figure out the coefficient relation (equivalent to that of equation
(28)), which at a given $n$ takes the form

$\sum_{j=-1}^{k+1}\left(\widetilde{\gamma}_{n,j-1}+\widetilde{\beta}_{n,j}+\widetilde{\alpha}_{n,j+1}\right)R_{n-j}=0\,.$
(48)

The coefficients appearing in the relation are given by

$\displaystyle\widetilde{\alpha}_{n,j}$
$\displaystyle=\left(-\eta\right)^{j}\binom{k}{j}\alpha^{\text{\tiny{bg}}}_{n-j}+(-1)^{j}\binom{k+2}{j}\frac{\alpha^{(-k)}}{\beta^{2}}\,,$
(49) $\displaystyle\widetilde{\beta}_{n,j}$
$\displaystyle=(-\eta)^{j}\binom{k}{j}\beta^{\text{\tiny{bg}}}_{n-j}\,,$ (50)
$\displaystyle\widetilde{\gamma}_{n,j}$
$\displaystyle=\left(-\eta\right)^{j}\binom{k}{j}\gamma^{\text{\tiny{bg}}}_{n-j}\,.$
(51)

We notice that, from the definition of the binomial,
$\widetilde{\alpha}_{n,j}$ is non-vanishing for $0\leq j\leq k+2$, while
$\widetilde{\beta}_{n,j}$ and $\widetilde{\gamma}_{n,j}$ are non-zero for
$0\leq j\leq k$. Now that we have a $k+3$ terms relation, we can perform a
Gaussian elimination to reduce it to a three-terms relation (details can be
found in the appendix of [38]). Once the three-terms relation is found, one
can re-initialize the ladder operator $\Lambda^{r}_{n}$ and obtain the
modified frequency and separation constant from equations (36)–(37).

###  III.3 Numerical computation of the coefficients

In the previous two sections we explained how to obtain the functions
$\mathcal{L}_{r}\left(\omega,B,\alpha\right)$ and
$\mathcal{L}_{\theta}\left(\omega,B\right)$. To compute the coefficients
$d_{\omega}$ and $d_{B}$ as given in equation (13), we evaluate the
derivatives numerically with a 4-points centered stencil. For each pair of
coefficients, we initialize the ladder operators $\Lambda_{N}$ to some
arbitrary low integer $N$, and then increase it by one until the simultaneous
relative change in $d_{\omega}$ and $d_{B}$ is smaller than a given tolerance
(which we chose to be $10^{-7}$). We computed numerically all the coefficients
for the following values $s=-2$, $n=[0,2]$, $\ell=[2,4]$, $m=[-\ell,\ell]$,
$k=[-3,4]$ in a uniform grid in $a=[0,0.495]$ with spacing $\delta a=0.005$.
The full list of coefficients is available in a public git folder [XXX].

In figure 1 we show the results from this computation for the real and
imaginary parts of the $d_{\omega}^{(k)}$ coefficients for $s=-2$, $n=0$,
$\ell=2$, $m=[-2,2]$ for values of $k=[-3,4]$ and of the spin $a$ comprised
between $0$ and $0.45$, as well as the real and the imaginary part of
$d_{B}^{(k)}$ for the same $n,\ell,m$ and $k=[-1,2]$.

To directly apply our formalism to further studies, e.g., ringdown analysis of
non-linear computations or data analysis, we also provide a python code and a
jupyter notebook with some examples [48]. It allows one to compute the QNMs
and the separation constants as function of $n,\ell,m,a$ and $\alpha^{(k)}$
and can thus, in principle, be efficiently integrated in commonly used code
infrastructure. The code also allows one to access some of the earlier results
for the parametrized QNM framework for modifications to the Regge-Wheeler and
Zerilli potentials, for which coefficients beyond the fundamental mode have
been computed in reference [38]. The GR values for the QNMs have been taken
from reference [3, 49]. For more details about how the code is structured and
how it can be used, we refer to the provided tutorial.

In principle, one should be able to compare the coefficients for $a=0$ with
those computed in [36, 38]. However, we stress that for $a=0$, equation (7)
reduces to the non-spinning limit of the Bardeen-Press equation [50], whereas
the formalism of [36, 38] was developed for the Regge-Wheeler and the Zerilli
equation. The transformation between the Bardeen-Press potential and the
Regge-Wheeler/Zerilli potentials was obtained by Chandrasekhar [51], but
generalizing this to the case of the modified potential with generic
$\alpha^{(k)}$ couplings is non-trivial.

Figure 1: In the top panel we show the real and the imaginary part of
$d_{\omega}$ for $n=0$, $\ell=2$, $m=[-2,2]$, $k=[-3,4]$ and values of the
spin from $a=0$ to $a=0.45$, and each point is on a step of $\delta a=0.005$.
The inset focuses around the coefficients with $k=[-3,0]$. In the bottom panel
we show the real and the imaginary part of $d_{B}$ for $n=0$, $\ell=2$,
$m=[-2,2]$, $k=[-1,2]$ and same values of the spin. The inset focuses around
the coefficients with $m\leq 0$ and $k=[-1,0]$. In both plots the black dot
signals the coefficient value for $a=0$.

###  III.4 Linearized regime of validity

The framework we developed is motivated by the assumption that any
modification of gravity produces only slight deviations from GR in
astrophysical observables. In this section, we expand on the regime of
validity of the formalism, by providing a quantitative assessment of the
accuracy of such approach. It is worth noting that we can only asses the error
made by restricting to linear corrections to the frequencies, as defined in
equation (9) and not taking into account higher-order corrections to the
potential, which are beyond the scope of this paper.

Figure 2: We plot the threshold values $\alpha^{(k)}_{*}$ against the spin for
different values of $k$. With $\alpha^{(k)}_{*}$, we identify the limit value
of a coupling at which an error of $1\%$ the a linear and a non-linear
approximation is obtained. The plot uses as non-linear estimates the full
continued fraction results. The error is evaluated for the $n=0$, $\ell=m=2$
mode and different values of $k$ and of the spin $a$. The black dashed lines
correspond to the estimate (53)

First of all, we give a heuristic motivation on the maximum size of the
coefficients, by requesting that the perturbation equation is not strongly
modified at the boundaries of our dominion and that
$\omega\simeq\mathcal{O}(1)$. At $r\to\infty$, we have seen from equations
(40) and (41) that the only coefficients modifying the asymptotic structure of
the potential are $\alpha^{(3)}$ and $\alpha^{(4)}$. With some simple algebra,
we can infer

$\alpha^{(3)}\lesssim 1\,,\qquad\alpha^{(4)}\lesssim 1\,.$ (52)

On the other hand, the modifications in the potential affect the near-horizon
expansion as in equation (39). The condition is such that the sum of
coefficients must behave as

$\sum_{k}\alpha^{(k)}\lesssim\left|2\beta^{2}\sigma_{\text{\tiny{GR}}}\left(\sigma_{\text{\tiny{GR}}}-\frac{\mathrm{i}s}{2}\right)\right|\,.$
(53)

It is worth noting that when the superradiant condition $\omega=\omega_{c}$ is
activated, one has $\sigma_{\text{\tiny{GR}}}=0$, and we expect that the
formalism is valid only if the sum of the $\alpha^{(k)}$ is approximately 0.

In general, however, each power of $k$ affects in a different way the
effective potential. In order to have a more quantitative estimate of the
allowed regime of validity, we perform two separate analysis. First, we
compare the QNM frequencies computed with the linear approximation against
those obtained with a full continued-fraction method discussed above. We
estimated the error on the frequencies as

$\Delta_{\omega}\equiv\sqrt{\left(\frac{\Delta\omega_{R}}{\omega_{R}}\right)^{2}+\left(\frac{\Delta\omega_{I}}{\omega_{I}}\right)^{2}}\,,$
(54)

where $\omega_{R,I}$ are respectively the real and imaginary parts of the QNM.
We computed the error for several real positive and negative values of the
couplings $\alpha^{(k)}$ and extracted the threshold values $\alpha^{(k)}_{*}$
at which the error reaches $1\%$. The results are represented in figure 2 for
the mode $n=0$, $\ell=m=2$, for selected values of the spins between $a=0$ and
$a=0.45$ and for $k=[-2,4]$.

It can be seen that there is a complicated dependence of the thresholds on the
type of modification we introduce in the potential. However, a main
qualitative feature can be read off, i.e. that, for any modification, the
threshold tends to get smaller for higher spins. The physical interpretation
of that, is that, for a given beyond-GR effect in the modified Teukolsky
equation, rotation tends to exacerbate the deviation of the linear
approximation with respect to the true values of the QNMs. Bearing this caveat
in mind, we will still show in the next section that the linear approximation
provides very good results in a couple of known models of perturbation of
rotating BHs with deviations from Kerr, also for high spin.

Since the computation of QNMs with the continued fraction method is not
immediate nor straightforward to implement, we want to provide a quick
estimate for the errors of the single-$k$ contributions. In this respect, we
compute the diagonal quadratic corrections, as explained in appendix D. We
checked the estimate $\alpha_{*}^{(k)}$ by computing the error
$\Delta_{\omega}$ assuming that the non-linear frequencies are obtained
including quadratic coefficients. By a qualitative comparison, the quadratic
estimate works well to capture the error except for $k=3$, and partially for
$k=4$. Even though it is not as precise as the full non-linear comparison, the
quadratic coefficients can be used as a quick way to understand what threshold
value to take for the couplings.

Lastly, we want to stress that the thresholds that we provided in this
section, are referred to the contribution of a single modification. Hence, it
could be that, depending on the values of the coefficients, the combination of
multiple $k$ would need larger or smaller threshold values. This means that
for a theory-specific case, the bounds on $\alpha^{(k)}$ might differ from
what we inferred in this section, and need to be addressed case-by-case.

##  IV Applications

###  IV.1 Massive scalar perturbations

The first example that we provide to test our formalism is for the computation
of the QNMs of a massive scalar field, a case extensively studied in the
literature [52, 53, 54]. The radial and angular perturbation equations for a
massive scalar field ($s=0$) with mass $\mu$ are

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}r}$ $\displaystyle\left[\Delta
R^{\prime}(r)\right]+\left(\frac{K^{2}}{\Delta}-\lambda_{\ell
m}-\mu^{2}r^{2}\right)R(r)=0\,,$ (55)
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}y}$
$\displaystyle\left[\left(1-y^{2}\right)S^{\prime}(y)\right]$
$\displaystyle+\left[a^{2}\left(\omega^{2}-\mu^{2}\right)y^{2}+B_{\ell
m}-\frac{m^{2}}{1-y^{2}}\right]S(y)=0\,.$ (56)

First of all, we bring the angular equation into the form of equation (6) by
transforming $\omega\to\omega+\frac{\mu^{2}}{2\omega}$. Then, by assuming
$\mu\ll 1$, the radial equation is automatically brought in the form of (7),
with the only non-zero $\alpha^{(k)}$ being

$\alpha^{(1)}=\mu^{2}ar_{+}\left(a-\frac{m}{\omega_{0}}\right)\,,\qquad\alpha^{(3)}=\mu^{2}r_{+}^{3}\,,$
(57)

where $\omega_{0}$ is the unperturbed Kerr frequency. The effect of the mass
on the frequency at linear order in $\mu^{2}$ is given by

$\omega_{\text{\tiny{L}}}=\omega_{0}+\frac{\mu^{2}}{2\omega_{0}}+d_{(1)}\alpha^{(1)}+d_{(3)}\alpha^{(3)}\,.$
(58)

In figure 3 we show the difference
$\delta\omega=|\omega_{\text{\tiny{L}}}-\omega_{\text{\tiny{NL}}}|$ between
the linear results in (58) and the nonlinear QNMs $\omega_{\text{\tiny{NL}}}$
computed in [54] for $\ell=m=2$ modes.

Figure 3:  Relative difference for the real part (solid line) and imaginary
part (dashed line) of the fundamental $\ell=m=2$ mode for a massive scalar
perturbation computed with the linear approximation against the non-linear
results of [54].

###  IV.2 The Dudley-Finley equation

Figure 4: Difference between the real (solid lines) and imaginary (dashed
lines) part of the Dudley-Finley QNMs computed either with the linear
perturbative approach or the full continued fraction method. We show results
for $n=0$ (left panels), $n=1$ (central panels) and $n=2$ (right panels),
$\ell=2$, $m=[-2,2]$ (bottom to top panels) for different values of the spin
and of the electric charge. Note that the three curves have different
endpoints, as for a given $Q$ the maximum value of $a$ is
$a_{\text{\tiny{max}}}=\frac{1}{2}\sqrt{1-4Q^{2}}$

As a second example, we would like to test our formalism against gravitational
perturbation of a Kerr-Newman (KN) BH in the limit of small charge, since the
QNMs for a generic electric charge $Q$ have been computed numerically in [55]
and fits are available in [20]. Unfortunately, the KN perturbation equation is
not explicitly separable, not even in the limit of small charge [56] in which
at least electromagnetic and gravitational perturbations decouple. One could
apply the algorithm of [33, 34, 35] to obtain a modified Teukolsky operator
for the KN solution, but it goes beyond the scopes of this paper. For the sake
of testing the method, we can restrict ourselves to the Dudley-Finley (DF)
equation, a proxy equation for the perturbations of the Kerr-Newman metric
[57, 58, 59]. The DF equation is obtained by taking equation (1) and
performing the following substitution

$\Delta\to\Delta+Q^{2}\,.$ (59)

We can then rescale the equation into the form (7) by assuming that $Q\ll 1/2$
and by defining a new spin parameter

$\widebar{a}=a+\frac{Q^{2}}{2a}\,,$ (60)

such that $\widebar{a}^{2}\simeq a^{2}+Q^{2}$. Retaining only the terms
quadratic in $Q$, we can see that the only non-zero $\alpha^{(k)}$ terms that
contribute to the equation (8) are

$\displaystyle\alpha^{(0)}$
$\displaystyle=Q^{2}\left[\frac{\mathrm{i}s}{2\widebar{a}}(m-2\widebar{a}\widebar{\omega}_{0})-(m-\widebar{a}\widebar{\omega}_{0})^{2}\right]\,,$
(61) $\displaystyle\alpha^{(1)}$
$\displaystyle=Q^{2}r_{+}\left[\widebar{\omega}_{0}\frac{m-\widebar{a}\omega_{0}}{\widebar{a}}-\frac{\mathrm{i}s}{\widebar{a}}(m-2\widebar{a}\widebar{\omega}_{0})\right]\,,$
(62) $\displaystyle\alpha^{(2)}$
$\displaystyle=-Q^{2}r_{+}^{2}\widebar{\omega}_{0}^{2}\,,$ (63)

where $\widebar{\omega}_{0}$ is the Kerr frequency evaluated at spin
$\widebar{a}$. The DF linear frequencies at spin $a$ are obtained by

$\omega_{\text{\tiny{L}}}=\widebar{\omega}_{0}+\sum_{k=0}^{2}\alpha^{(k)}d_{(k)}\,.$
(64)

In Figure 4 we show the real and imaginary part of the absolute difference
$\delta\omega=|\omega_{\text{\tiny{L}}}-\omega_{\text{\tiny{NL}}}|$ between
the linear results in (64) and the nonlinear QNMs $\omega_{\text{\tiny{NL}}}$
computed via the Leaver method in [59], for various spins and different values
of the electric charge. We make the comparison with $\ell=2$ modes, with all
values of $m=[-2,2]$, for the fundamental and first two overtones. The plot
clearly shows that the discrepancy between the linearized QNMs and the full
non-linear results scales with the charge, and the approximation remains valid
for all the different values of $(n,\ell,m)$ surveyed.

Finally, we comment on the fact that the errors grow for small values of the
spin. This is due to the fact that in order to bring the equation in the form
of (7), we performed the transformation (60), which brings a term $1/a$ to the
denominator when $m\neq 0$. In other words, this transformation is valid as
long as $|Q|\ll|a|$. Nevertheless, the smallness of the universal coefficients
$d_{(k)}$ is such that the combination in frequency (64) is finite and
faithful to the non-linear value.

###  IV.3 Higher derivative gravity

Now we want to check the prediction of QNMs in higher derivative gravity using
the parametrized method against the results presented in [46]. In a companion
paper [45], focused on the analysis of QNMs in higher-derivative gravity, we
show how to reduce the radial perturbation equation to the form of equation
(7), with the only non-vanishing values of $\alpha^{(k)}$ being
$k=k^{\text{\tiny{HD}}}=[-2,0,1,2]$

$\delta V^{\pm}=\lambda\sum_{k\in
k^{\text{\tiny{HD}}}}\alpha_{\pm}^{(k)}\left(\frac{r}{r_{+}}\right)^{k}\,,$
(65)

where the $\pm$ refers to the polarization of the perturbation and we
collected out $\lambda$, the normalized coupling constant of the
theory.222cfr. equation (30) of [45], for which
$\alpha^{(k)}=A^{(k)}r_{+}^{k}$, and the coupling constant has been previously
factorized out. Here we use $\lambda$ to refer to the coupling constant, to
avoid misunderstanding with the $\alpha_{\mathrm{q}}$ used in [45]. From this,
we can compute the frequencies deviations, normalized by the coupling constant
$\lambda$

$\delta\omega^{\pm}=\frac{\omega^{\pm}-\omega^{\text{\tiny{Kerr}}}}{\lambda}=\sum_{k\in
k^{\text{\tiny{HD}}}}\alpha^{(k)}_{\pm}d_{(k)}\,,$ (66)

for each parity, and each realization of the theory. In figure 5, we compare
our results against the fits $\delta\omega^{\mathrm{fit}}$ given in [46]. We
truncate the analysis at spin $a=0.35$, since the fits are valid only up to
this value. The plot shows remarkable agreement between the corrections
computed with two different methods, strengthening the validity of the
parametrized formalism. In figure 5 we limited to show polar and axial,
$\ell=m=2$ values for the even-parity cubic theory, labelled as
$\delta\omega_{\mathrm{even}}^{\pm}$, as well as polar, $\ell=m=2$ values for
the odd-parity cubic theory, labelled as $\delta\omega_{\mathrm{odd}}^{+}$.
Details on the definition of these modes can be found in [46] and in the
companion paper where we perform an extensive study of QNMs of rotating BHs in
higher derivative gravity [45].

Figure 5: Real and imaginary part of the absolute difference between the fits
of [46] and the frequencies computed with the parametrized formalism for the
$\pm$ polarizations of even cubic modes and the $+$ polarization of the odd
cubic mode

##  V Conclusions

In this paper we have shown how to connect small deviations parametrized by
powers of the radial coordinate $r$ in the Teukolsky equation to small
deviations in the eigenfrequencies and in the separation constants of modified
Kerr BHs. We proved that for each value of $n,\ell,m$ there are up to nine
independent coefficients in the radial parametrization, but for specific cases
they could be less. We presented a robust method to compute the coefficients
that control the linear corrections to the QNM frequencies and separation
constants through a generalization of Leaver’s continued fractions, and used
it to compute them for $n=[0,2]$, $\ell=[2,4]$, $m=[-\ell,\ell]$ and
$k=[-3,4]$ in a range of spins between $0$ and $0.495$.333We recall that in
our conventions extremality corresponds to $a=1/2$ These results are available
online in a public git repository [48], together with a python code and a
jupyter notebook to compute the QNMs and separation constants. There we also
provide a tutorial demonstrating how to use the code, which can in principle
be applied to compute QNMs with arbitrary $n$, $l$, $m$, $k$ and angular
momentum besides those we computed explicitly here.

We checked the quality of the predictions against three cases known in the
literature: perturbations of a massive scalar field around Kerr, the Dudley-
Finley equation and the QNMs of BHs in higher-derivative gravity. For all the
three cases the frequencies predicted by the formalism show great agreement
with respect to the results in the literature.

These coefficients will be particularly useful for the computation of QNMs of
rotating BHs in alternative theories of gravity for which a modified Teukolsky
equation is obtained, in the spirit of the method developed in [33, 34, 35].
So far, this method has been successfully applied to higher-derivative gravity
[46], but other theories like scalar-Gauss-Bonnet gravity and dynamical-Chern-
Simons gravity [60] are good candidates for this computation. In order to
study those cases, it would be interesting to generalize our parametrized
formalism by including couplings between the Teukolsky equation and a scalar
field, analogous to the analysis of [37] in the case of static BHs.

The modified Teukolsky approach of [33, 34, 35] is currently limited by the
fact that metric reconstruction is only available for GR, meaning that it is
not yet possible to extend the method beyond first order in the coupling. For
this reason, we did not explore further the quadratic coefficients as it was
done in the non-rotating case [37, 38], and we limited the computation of the
diagonal ones just to have a quick estimate of the error of the method itself.

Let us also remark that a different approach to study beyond-GR QNMs based on
spectral methods has recently been introduced and successfully demonstrated
for a wide range of spins in Refs. [61, 62, 63, 64, 65]. These methods have
the advantage of being more flexible, but on the other hand, they have a much
higher computational complexity and cost than the standard perturbative
approaches. In this regard, the results of our method applied to specific
theories may be useful to validate the spectral approaches.

The most intriguing open problem from our analysis is whether one can find a
better way to exploit the potential ambiguity, as done in section II.2. With
the choice we made, we could reduce the number of independent coefficients
$\alpha^{(k)}$ to 8. For the case of higher derivative gravity we have been
able to numerically reduce the number of free coefficients to just 4 [45]. One
may wonder whether higher derivative gravity has a special structure of the
equations, or if there is a fundamental transformation of the potential that
could diminish the number of free parameters.

Such discussion is relevant especially if one wants to use this formalism in a
theory-agnostic setup, e.g., to perform ringdown tests of GR or the inverse
problem. Since the coefficients $d^{(k)}$ have a spin dependence, it would be
interesting to map them to ParSpec [17]. Another useful mapping would be with
the WKB deviation coefficients, as done in [40]. Finally, in the upcoming
analysis, it would be interesting to compare the detectability of beyond-
Teukolsky effects against that of second order QNMs [66, 67, 68].

###### Acknowledgements.

We thank Emanuele Berti for fruitful discussions. We also thank Adrian Ka-Wai
Chung, Pratik Wagle and Nicolas Yunes for useful comments. The work of P. A.
C. received the support of a fellowship from “la Caixa” Foundation (ID
100010434) with code LCF/BQ/PI23/11970032. S. M. acknowledges support from the
Flemish inter-university project IBOF/21/084. S. H. V. acknowledges funding
from the Deutsche Forschungsgemeinschaft (DFG): Project No. 386119226. The
work of L. C. was supported by the European Union’s H2020 ERC Consolidator
Grant “GRavity from Astrophysical to Microscopic Scales” (Grant No.
GRAMS-815673), the PRIN 2022 grant “GUVIRP - Gravity tests in the UltraViolet
and InfraRed with Pulsar timing”, and the EU Horizon 2020 Research and
Innovation Programme under the Marie Sklodowska-Curie Grant Agreement No.
101007855.

##  Appendix A Ambiguity of the potential modifications

Here we list the explicit form of the coefficients of equation (18) when one
transforms the radial Teukolsky function as (14). These values hold for
$s=-2$:

$\displaystyle\overline{A}_{j}^{(-3)}=$
$\displaystyle\,\frac{a^{6}}{2}j(j+1)(j+2)\,,\qquad\overline{A}_{j}^{(-2)}=\,-\frac{3a^{4}}{2}j(j+1)^{2}\,,$
(67) $\displaystyle\overline{A}_{j}^{(-1)}=$
$\displaystyle\,\frac{a^{4}}{2}j\left[3j(j+1)-4B+4m^{2}-8i\omega+10\right]+4ia^{3}jm+\frac{3}{2}a^{2}j\left(j^{2}+j-1\right)\,,$
(68) $\displaystyle\overline{A}_{j}^{(0)}=$
$\displaystyle\,a^{3}(2j-1)\left[a\omega(\omega-4i)-2m(\omega+2i)\right]$
$\displaystyle-a^{2}\left[3j^{3}-j\left(4B-2m^{2}+4i\omega-1\right)+B+2\right]-4iajm-\frac{1}{2}j\left(j^{2}-4\right)\,,$
(69) $\displaystyle\overline{A}_{j}^{(1)}=$
$\displaystyle\,2a^{4}(j-1)\omega^{2}+\frac{1}{2}a^{2}\Big{[}3j^{2}(j-1)+4B+2\omega(\omega-16i)+8+j\left(-8B+4m^{2}-4\omega(\omega-8i)+2\right)\Big{]}$
$\displaystyle+2am\left[2j(\omega+3i)-\omega-2i\right]+\frac{1}{2}j\left[3(j-1)j-4B-11\right]+B+2\,,$
(70) $\displaystyle\overline{A}_{j}^{(2)}=$
$\displaystyle\,a^{2}\omega\left[\omega-16i(j-1)\right]-2a(2j-1)m(\omega+2i)+\frac{1}{2}j(-3(j-2)j+8B-24i\omega+13)-3(B-4i\omega+2)\,,$
(71) $\displaystyle\overline{A}_{j}^{(3)}=$
$\displaystyle\,2a^{2}(2j-3)\omega^{2}+\frac{1}{2}(j-1)\left[(j-2)j-4(B+2)\right]+4i(5j-6)\omega\,,$
(72) $\displaystyle\overline{A}_{j}^{(4)}=$
$\displaystyle\,2\omega\left[2(\omega+3i)-j(\omega+4i)\right]\,,\qquad\overline{A}_{j}^{(5)}=\,2(j-2)\omega^{2}\,.$
(73)

##  Appendix B Nollert’s improvements of continued fraction

The procedure to numerically solve Teukolsky equation through Leaver’s method
requires in practice an initialization for the radial ladder operator
$\Lambda^{r}_{n}$. Such quantity can be expanded for large initialization
number $N$ as

$\Lambda^{r}_{N}=\sum_{j=0}^{J}C_{j}N^{-j/2}+\mathcal{O}\left(N\right)^{-\left(J+1\right)/2}\,.$
(74)

One can initialize the ladder operator just retaining the first term
$C_{0}=-1$. However, this approximation requires in general a very high
initial value for $N$ (which means long computational time) and appears to be
insufficient for frequencies with large imaginary part (namely higher
overtones). In [69], it was shown that adding further corrections to the
initial $\Lambda^{r}_{N}$ improves the accuracy of the method, also allowing
to capture higher overtones. The $(k+3)$-terms recurrence relation of equation
(48) for $n=N$ can be expressed as

$\sum_{j=-1}^{k+1}M_{N,j}\,R_{N-j}=0\,,$ (75)

where we defined

$M_{N,j}\equiv\widetilde{\gamma}_{N,j-1}+\widetilde{\beta}_{N,j}+\widetilde{\alpha}_{N,j+1}\,.$
(76)

Dividing it by $R_{N-k-1}$ and using the definition of ladder operator
$\Lambda^{r}_{N}=-a_{N+1}/a_{N}$ one obtains the equation

$\sum^{k+1}_{j=-1}(-1)^{j}\,M_{N,j}\,\prod_{i=j+1}^{k+2}\Lambda^{r}_{N-i}=0\,.$
(77)

Plugging the definition (74) into the above formula, one can solve for the
coefficients order by order. By the definition of $\widetilde{\gamma}_{N,j}$,
$\widetilde{\beta}_{N,j}$ and $\widetilde{\alpha}_{N,j}$, they scale with $N$
as

$\displaystyle\widetilde{\alpha}_{N,j}$
$\displaystyle=N^{2}+\widetilde{\alpha}_{1}N+\widetilde{\alpha}_{0}$ (78)
$\displaystyle\widetilde{\beta}_{N,j}$
$\displaystyle=-2N^{2}+\widetilde{\beta}_{1}N+\widetilde{\beta}_{0}$ (79)
$\displaystyle\widetilde{\gamma}_{N,j}$
$\displaystyle=N^{2}+\widetilde{\gamma}_{1}N+\widetilde{\gamma}_{0}$ (80)

If we fix $C_{0}=-1$, then the other coefficients up to $J=5$ can be written
as

$\displaystyle C_{1}=$
$\displaystyle\,\pm\sqrt{-\widetilde{\alpha}_{1}-\widetilde{\beta}_{1}-\widetilde{\gamma}_{1}}\,,$
(81) $\displaystyle C_{2}=$
$\displaystyle\,\widetilde{\alpha}_{1}+\frac{\widetilde{\beta}_{1}}{2}-\frac{1}{4}\,,$
(82) $\displaystyle C_{3}=$
$\displaystyle\,\frac{C_{2}^{2}}{2C_{1}}-\frac{C_{2}}{4C_{1}}-\frac{\widetilde{\alpha}_{0}+\widetilde{\beta}_{0}+\widetilde{\gamma}_{0}}{2C_{1}}-\frac{1+2\widetilde{\alpha}_{1}}{4}C_{1}\,,$
(83) $\displaystyle C_{4}=$
$\displaystyle\,\widetilde{\alpha}_{0}-\frac{4\widetilde{\alpha}_{1}-8\widetilde{\beta}_{0}+1}{16}-\frac{1+4\widetilde{\alpha}_{1}}{4}C_{2}-\frac{C_{3}}{2C_{1}}\,,$
(84) $\displaystyle C_{5}=$
$\displaystyle\,\frac{2\widetilde{\alpha}_{1}+3}{4C_{1}}C_{2}^{2}-\frac{8\widetilde{\alpha}_{0}+4\widetilde{\alpha}_{1}+3}{16}C_{1}-\frac{C_{3}^{2}}{2C_{1}}$
$\displaystyle-\frac{4\widetilde{\alpha}_{1}+3}{16C_{1}}C_{2}+C_{3}\left(\frac{C_{2}}{2C_{1}^{2}}-\widetilde{\alpha}_{1}-1\right)$
$\displaystyle+\left(\frac{C_{2}}{C_{1}}-\frac{3}{4C_{1}}\right)C_{4}\,,$ (85)

and the sign of $C_{1}$ is chosen such that $\mathrm{Re}\left(C_{1}\right)>0$.

We can do the same expansion for the angular part, by expanding

$\Lambda^{\theta}_{N}=\sum_{j=0}^{J}D_{j}N^{-j/2}+\mathcal{O}\left(N\right)^{-\left(J+1\right)/2}\,.$
(86)

From equations (32)–(34) we can schematically say that

$\displaystyle\alpha^{\theta}_{N}$
$\displaystyle=-2N^{2}+\alpha_{1}N+\alpha_{0}$ (87)
$\displaystyle\beta^{\theta}_{N}$ $\displaystyle=N^{2}+\beta_{1}N+\beta_{0}$
(88) $\displaystyle\gamma^{\theta}_{N}$ $\displaystyle=\gamma_{1}N+\gamma_{0}$
(89)

By solving perturbatively in $1/N$ the relation (35), we obtain the following
expression for the coefficients $D_{j}$ up to $J=4$

$\displaystyle D_{0}=$ $\displaystyle\,0$ (90) $\displaystyle D_{1}=$
$\displaystyle\,\gamma_{1}$ (91) $\displaystyle D_{2}=$
$\displaystyle\,\gamma_{0}-\gamma_{1}\left(1+\beta_{1}+2\gamma_{1}\right)$
(92) $\displaystyle D_{3}=$
$\displaystyle\,\gamma_{1}^{2}(\alpha_{1}-2)-\gamma_{1}(1+\beta_{0}+\beta_{1})$
$\displaystyle-D_{2}(\beta_{1}+4\gamma_{1}+2)$ (93) $\displaystyle D_{4}=$
$\displaystyle\,\alpha_{0}\gamma_{1}^{2}+\frac{2D_{2}^{2}\left(\beta_{1}+3\gamma_{1}+2\right)}{\gamma_{1}}$
$\displaystyle+D_{2}\left(\beta_{0}+\beta_{1}+2\gamma_{1}+1\right)$
$\displaystyle+D_{3}\left(\frac{2D_{2}}{\gamma_{1}}-\beta_{1}-4\gamma_{1}-2\right)$
(94)

##  Appendix C Splitting of the potential

In this section of the appendix we show how to transform the potential (8)
into the potential (38). We start by splitting equation (8) into

$\displaystyle\delta V(r)$
$\displaystyle=\frac{1}{\Delta}\sum_{k=-K}^{4}\alpha^{(k)}\left(\frac{r}{r_{+}}\right)^{k}$
$\displaystyle=\frac{1}{\Delta}\sum_{k=0}^{4}\alpha^{(k)}\left(\frac{r}{r_{+}}\right)^{k}+\frac{1}{\Delta}\sum_{k=1}^{K}\alpha^{(-k)}\left(\frac{r_{+}}{r}\right)^{k}$
(95)

For $k\geq 1$, the first generic term in $k$ of the sum can be rewritten as

$\displaystyle\frac{\alpha^{(k)}}{\Delta}\left(\frac{r}{r_{+}}\right)^{k}$
$\displaystyle=\frac{\alpha^{(k)}}{\Delta}\left(\frac{r^{k}-r_{+}^{k}}{r_{+}^{k}}+1\right)$
$\displaystyle=\alpha^{(k)}\left[\frac{1}{\Delta}+\frac{1}{r_{+}(r-r_{-})}\sum_{j=0}^{k-1}\left(\frac{r}{r_{+}}\right)^{j}\right]$
(96)

For $k\geq 2$ we can further simplify this term as

$\displaystyle\frac{\alpha^{(k)}}{\Delta}\left(\frac{r}{r_{+}}\right)^{k}$
$\displaystyle=\alpha^{(k)}\left[\frac{1}{\Delta}+\frac{1}{r_{+}(r-r_{-})}\sum_{j=0}^{k-1}\frac{r^{j}-r_{-}^{j}+r_{-}^{j}}{r_{+}^{j}}\right]$
$\displaystyle=\alpha^{(k)}\left[\frac{1}{\Delta}+\frac{1}{r_{+}(r-r_{-})}\sum_{j=0}^{k-1}\left(\frac{r_{-}}{r_{+}}\right)^{j}+\frac{1}{r_{-}r_{+}}\sum_{j=1}^{k-1}\left(\frac{r_{-}}{r_{+}}\right)^{j}\sum_{n=0}^{j-1}\left(\frac{r}{r_{-}}\right)^{n}\right]$
$\displaystyle=\alpha^{(k)}\left[\frac{1}{\Delta}+\frac{1}{r_{+}(r-r_{-})}\sum_{j=0}^{k-1}\left(\frac{r_{-}}{r_{+}}\right)^{j}+\frac{1}{r_{+}^{2}}\sum_{j=0}^{k-2}\left(\frac{r}{r_{-}}\right)^{j}\sum_{n=j}^{k-2}\left(\frac{r_{-}}{r_{+}}\right)^{n}\right]$
(97)

where we obtained the last line by expanding the series and collecting the
terms in $r$ to the same power. Summing over all the non-negative values of
$k$ yields

$\displaystyle\frac{1}{\Delta}\sum_{k=0}^{4}\alpha^{(k)}\left(\frac{r}{r_{+}}\right)^{k}=\frac{1}{\Delta}\sum_{k=0}^{4}\alpha^{(k)}+\frac{1}{r_{+}(r-r_{-})}\sum_{k=1}^{4}\alpha^{(k)}\sum_{j=0}^{k-1}\left(\frac{r_{-}}{r_{+}}\right)^{j}+\frac{1}{r_{+}^{2}}\sum_{k=0}^{2}\left(\frac{r}{r_{+}}\right)^{k}\sum_{j=k}^{2}\alpha^{(j+2)}\sum_{n=0}^{j}\left(\frac{r_{-}}{r_{+}}\right)^{n}$
(98)

Now, we can perform a mapping between the coefficients $\alpha^{(k)}$ of
equation (8) and the coefficients $A^{(0)}$, $A^{(1)}$ and
$\widetilde{\alpha}^{(k)}$ introduced in equations (38). From a direct
comparison we have

$\displaystyle A^{(0)}$ $\displaystyle=\sum_{k=0}^{4}\alpha^{(k)}$ (99)
$\displaystyle A^{(1)}$
$\displaystyle=\sum_{k=1}^{4}\alpha^{(k)}\sum_{j=0}^{k-1}\left(\frac{r_{-}}{r_{+}}\right)^{j}$
(100) $\displaystyle\widetilde{\alpha}^{(k)}$
$\displaystyle=\sum_{j=k}^{2}\alpha^{(j+2)}\sum_{n=0}^{j}\left(\frac{r_{-}}{r_{+}}\right)^{n}$
(101)

##  Appendix D Diagonal quadratic coefficients

We show here how to compute the quadratic diagonal coefficients, defined from
the next-to-leading-order expansion

$\begin{split}\omega&\simeq\omega^{0}+\sum_{k}d_{\omega}^{(k)}\alpha^{(k)}+\frac{1}{2}e_{\omega}^{(k)}{\alpha^{(k)}}^{2}\,,\\\
B&\simeq
B^{0}+\sum_{k}d_{B}^{(k)}\alpha^{(k)}+\frac{1}{2}e_{B}^{(k)}{\alpha^{(k)}}^{2}\,.\end{split}$
(102)

By extending the Taylor expansion (11) to the second order in $\alpha$, we
obtain

$\left.\mathcal{L}_{j}\right|_{\text{\tiny{GR}}}+\alpha\left.\frac{\mathrm{d}\mathcal{L}_{j}}{\mathrm{d}\alpha}\right|_{\text{\tiny{GR}}}+\frac{\alpha^{2}}{2}\left.\frac{\mathrm{d}^{2}\mathcal{L}_{j}}{\mathrm{d}\alpha^{2}}\right|_{\text{\tiny{GR}}}+\mathcal{O}(\alpha)^{3}=0\,.$
(103)

By expanding with the chain rule the total derivative
$\mathrm{d}^{2}\mathcal{L}_{j}/\mathrm{d}\alpha^{2}|_{\text{\tiny{GR}}}$, one
can read off the quadratic coefficients as

$\begin{split}e_{\omega}=&\left(\frac{\partial\mathcal{L}_{r}}{\partial
B}\frac{\partial\mathcal{L}_{\theta}}{\partial\omega}-\frac{\partial\mathcal{L}_{r}}{\partial\omega}\frac{\partial\mathcal{L}_{\theta}}{\partial
B}\right)^{-1}\Bigg{[}\frac{\partial^{2}\mathcal{L}_{r}}{\partial\alpha^{2}}\frac{\partial\mathcal{L}_{\theta}}{\partial
B}+2d_{\omega}\,\frac{\partial^{2}\mathcal{L}_{r}}{\partial\alpha\partial\omega}\,\frac{\partial\mathcal{L}_{\theta}}{\partial
B}+2d_{B}\,\frac{\partial^{2}\mathcal{L}_{r}}{\partial\alpha\partial
B}\,\frac{\partial\mathcal{L}_{\theta}}{\partial B}\\\
&-d_{\omega}^{2}\,\left(\frac{\partial^{2}\mathcal{L}_{\theta}}{\partial\omega^{2}}\frac{\partial\mathcal{L}_{r}}{\partial
B}-\frac{\partial^{2}\mathcal{L}_{r}}{\partial\omega^{2}}\frac{\partial\mathcal{L}_{\theta}}{\partial
B}\right)-d_{B}^{2}\,\left(\frac{\partial^{2}\mathcal{L}_{\theta}}{\partial
B^{2}}\frac{\partial\mathcal{L}_{r}}{\partial
B}-\frac{\partial^{2}\mathcal{L}_{r}}{\partial
B^{2}}\frac{\partial\mathcal{L}_{\theta}}{\partial
B}\right)+2d_{\omega}\,d_{B}\,\left(\frac{\partial^{2}\mathcal{L}_{r}}{\partial\omega\partial
B}\frac{\partial\mathcal{L}_{\theta}}{\partial
B}-\frac{\partial^{2}\mathcal{L}_{\theta}}{\partial\omega\partial
B}\frac{\partial\mathcal{L}_{r}}{\partial B}\right)\Bigg{]}\,,\\\
e_{B}=&-\left(\frac{\partial\mathcal{L}_{r}}{\partial
B}\frac{\partial\mathcal{L}_{\theta}}{\partial\omega}-\frac{\partial\mathcal{L}_{r}}{\partial\omega}\frac{\partial\mathcal{L}_{\theta}}{\partial
B}\right)^{-1}\Bigg{[}\frac{\partial^{2}\mathcal{L}_{r}}{\partial\alpha^{2}}\frac{\partial\mathcal{L}_{\theta}}{\partial\omega}+2d_{\omega}\,\frac{\partial^{2}\mathcal{L}_{r}}{\partial\alpha\partial\omega}\,\frac{\partial\mathcal{L}_{\theta}}{\partial\omega}+2d_{B}\,\frac{\partial^{2}\mathcal{L}_{r}}{\partial\alpha\partial
B}\,\frac{\partial\mathcal{L}_{\theta}}{\partial\omega}\\\
&-d_{\omega}^{2}\,\left(\frac{\partial^{2}\mathcal{L}_{r}}{\partial\omega^{2}}\frac{\partial\mathcal{L}_{\theta}}{\partial\omega}-\frac{\partial^{2}\mathcal{L}_{\theta}}{\partial\omega^{2}}\frac{\partial\mathcal{L}_{r}}{\partial\omega}\right)-d_{B}^{2}\,\left(\frac{\partial^{2}\mathcal{L}_{r}}{\partial
B^{2}}\frac{\partial\mathcal{L}_{\theta}}{\partial\omega}-\frac{\partial^{2}\mathcal{L}_{\theta}}{\partial
B^{2}}\frac{\partial\mathcal{L}_{r}}{\partial\omega}\right)+2d_{\omega}\,d_{B}\,\left(\frac{\partial^{2}\mathcal{L}_{r}}{\partial\omega\partial
B}\frac{\partial\mathcal{L}_{\theta}}{\partial\omega}-\frac{\partial^{2}\mathcal{L}_{\theta}}{\partial\omega\partial
B}\frac{\partial\mathcal{L}_{r}}{\partial\omega}\right)\Bigg{]}\,,\end{split}$
(104)

## References

  * Abbott  _et al._ [2023] R. Abbott _et al._ (KAGRA, VIRGO, LIGO Scientific), GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo during the Second Part of the Third Observing Run, Phys. Rev. X 13, 041039 (2023), arXiv:2111.03606 [gr-qc] .
  * Kokkotas and Schmidt [1999] K. D. Kokkotas and B. G. Schmidt, Quasinormal modes of stars and black holes, Living Rev. Rel. 2, 2 (1999), arXiv:gr-qc/9909058 .
  * Berti  _et al._ [2009] E. Berti, V. Cardoso, and A. O. Starinets, Quasinormal modes of black holes and black branes, Class. Quant. Grav. 26, 163001 (2009), arXiv:0905.2975 [gr-qc] .
  * Franchini and Völkel [2023a] N. Franchini and S. H. Völkel, Testing General Relativity with Black Hole Quasi-Normal Modes (2023a), arXiv:2305.01696 [gr-qc] .
  * Abbott  _et al._ [2016] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Tests of general relativity with GW150914, Phys. Rev. Lett. 116, 221101 (2016), [Erratum: Phys.Rev.Lett. 121, 129902 (2018)], arXiv:1602.03841 [gr-qc] .
  * Abbott  _et al._ [2021a] R. Abbott _et al._ (LIGO Scientific, Virgo), Tests of general relativity with binary black holes from the second LIGO-Virgo gravitational-wave transient catalog, Phys. Rev. D 103, 122002 (2021a), arXiv:2010.14529 [gr-qc] .
  * Abbott  _et al._ [2021b] R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), Tests of General Relativity with GWTC-3 (2021b), arXiv:2112.06861 [gr-qc] .
  * Carullo  _et al._ [2019] G. Carullo, W. Del Pozzo, and J. Veitch, Observational Black Hole Spectroscopy: A time-domain multimode analysis of GW150914, Phys. Rev. D 99, 123029 (2019), [Erratum: Phys.Rev.D 100, 089903 (2019)], arXiv:1902.07527 [gr-qc] .
  * Isi  _et al._ [2019] M. Isi, M. Giesler, W. M. Farr, M. A. Scheel, and S. A. Teukolsky, Testing the no-hair theorem with GW150914, Phys. Rev. Lett. 123, 111102 (2019), arXiv:1905.00869 [gr-qc] .
  * Cotesta  _et al._ [2022] R. Cotesta, G. Carullo, E. Berti, and V. Cardoso, Analysis of Ringdown Overtones in GW150914, Phys. Rev. Lett. 129, 111102 (2022), arXiv:2201.00822 [gr-qc] .
  * Finch and Moore [2022] E. Finch and C. J. Moore, Searching for a ringdown overtone in GW150914, Phys. Rev. D 106, 043005 (2022), arXiv:2205.07809 [gr-qc] .
  * Isi and Farr [2022] M. Isi and W. M. Farr, Revisiting the ringdown of GW150914 (2022), arXiv:2202.02941 [gr-qc] .
  * Carullo  _et al._ [2023] G. Carullo, R. Cotesta, E. Berti, and V. Cardoso, Reply to Comment on ”Analysis of Ringdown Overtones in GW150914”, Phys. Rev. Lett. 131, 169002 (2023), arXiv:2310.20625 [gr-qc] .
  * Crisostomi  _et al._ [2023] M. Crisostomi, K. Dey, E. Barausse, and R. Trotta, Neural posterior estimation with guaranteed exact coverage: The ringdown of GW150914, Phys. Rev. D 108, 044029 (2023), arXiv:2305.18528 [gr-qc] .
  * Pacilio  _et al._ [2024] C. Pacilio, S. Bhagwat, and R. Cotesta, Simulation-based inference of black hole ringdowns in the time domain (2024), arXiv:2404.11373 [gr-qc] .
  * Gennari  _et al._ [2024] V. Gennari, G. Carullo, and W. Del Pozzo, Searching for ringdown higher modes with a numerical relativity-informed post-merger model, Eur. Phys. J. C 84, 233 (2024), arXiv:2312.12515 [gr-qc] .
  * Maselli  _et al._ [2020] A. Maselli, P. Pani, L. Gualtieri, and E. Berti, Parametrized ringdown spin expansion coefficients: a data-analysis framework for black-hole spectroscopy with multiple events, Phys. Rev. D 101, 024043 (2020), arXiv:1910.12893 [gr-qc] .
  * Carullo [2021] G. Carullo, Enhancing modified gravity detection from gravitational-wave observations using the parametrized ringdown spin expansion coeffcients formalism, Phys. Rev. D 103, 124043 (2021), arXiv:2102.05939 [gr-qc] .
  * Maselli  _et al._ [2024] A. Maselli, S. Yi, L. Pierini, V. Vellucci, L. Reali, L. Gualtieri, and E. Berti, Black hole spectroscopy beyond Kerr: Agnostic and theory-based tests with next-generation interferometers, Phys. Rev. D 109, 064060 (2024), arXiv:2311.14803 [gr-qc] .
  * Carullo  _et al._ [2022] G. Carullo, D. Laghi, N. K. Johnson-McDaniel, W. Del Pozzo, O. J. C. Dias, M. Godazgar, and J. E. Santos, Constraints on Kerr-Newman black holes from merger-ringdown gravitational-wave observations, Phys. Rev. D 105, 062009 (2022), arXiv:2109.13961 [gr-qc] .
  * Pani and Cardoso [2009] P. Pani and V. Cardoso, Are black holes in alternative theories serious astrophysical candidates? The Case for Einstein-Dilaton-Gauss-Bonnet black holes, Phys. Rev. D 79, 084031 (2009), arXiv:0902.1569 [gr-qc] .
  * Yunes and Pretorius [2009] N. Yunes and F. Pretorius, Dynamical Chern-Simons Modified Gravity. I. Spinning Black Holes in the Slow-Rotation Approximation, Phys. Rev. D 79, 084043 (2009), arXiv:0902.4669 [gr-qc] .
  * Pani  _et al._ [2011] P. Pani, C. F. B. Macedo, L. C. B. Crispino, and V. Cardoso, Slowly rotating black holes in alternative theories of gravity, Phys. Rev. D 84, 087501 (2011), arXiv:1109.3996 [gr-qc] .
  * Maselli  _et al._ [2015] A. Maselli, P. Pani, L. Gualtieri, and V. Ferrari, Rotating black holes in Einstein-Dilaton-Gauss-Bonnet gravity with finite coupling, Phys. Rev. D 92, 083014 (2015), arXiv:1507.00680 [gr-qc] .
  * Cano and Ruipérez [2019] P. A. Cano and A. Ruipérez, Leading higher-derivative corrections to Kerr geometry, JHEP 05, 189, [Erratum: JHEP 03, 187 (2020)], arXiv:1901.01315 [gr-qc] .
  * Pani [2013] P. Pani, Advanced Methods in Black-Hole Perturbation Theory, Int. J. Mod. Phys. A 28, 1340018 (2013), arXiv:1305.6759 [gr-qc] .
  * Franchini [2023] N. Franchini, Slow rotation black hole perturbation theory, Phys. Rev. D 108, 044079 (2023), arXiv:2305.19313 [gr-qc] .
  * Cano  _et al._ [2022] P. A. Cano, K. Fransen, T. Hertog, and S. Maenaut, Gravitational ringing of rotating black holes in higher-derivative gravity, Phys. Rev. D 105, 024064 (2022), arXiv:2110.11378 [gr-qc] .
  * Pierini and Gualtieri [2021] L. Pierini and L. Gualtieri, Quasi-normal modes of rotating black holes in Einstein-dilaton Gauss-Bonnet gravity: the first order in rotation, Phys. Rev. D 103, 124017 (2021), arXiv:2103.09870 [gr-qc] .
  * Pierini and Gualtieri [2022] L. Pierini and L. Gualtieri, Quasinormal modes of rotating black holes in Einstein-dilaton Gauss-Bonnet gravity: The second order in rotation, Phys. Rev. D 106, 104009 (2022), arXiv:2207.11267 [gr-qc] .
  * Wagle  _et al._ [2022] P. Wagle, N. Yunes, and H. O. Silva, Quasinormal modes of slowly-rotating black holes in dynamical Chern-Simons gravity, Phys. Rev. D 105, 124003 (2022), arXiv:2103.09913 [gr-qc] .
  * Srivastava  _et al._ [2021] M. Srivastava, Y. Chen, and S. Shankaranarayanan, Analytical computation of quasinormal modes of slowly rotating black holes in dynamical Chern-Simons gravity, Phys. Rev. D 104, 064034 (2021), arXiv:2106.06209 [gr-qc] .
  * Li  _et al._ [2023] D. Li, P. Wagle, Y. Chen, and N. Yunes, Perturbations of Spinning Black Holes beyond General Relativity: Modified Teukolsky Equation, Phys. Rev. X 13, 021029 (2023), arXiv:2206.10652 [gr-qc] .
  * Hussain and Zimmerman [2022] A. Hussain and A. Zimmerman, Approach to computing spectral shifts for black holes beyond Kerr, Phys. Rev. D 106, 104018 (2022), arXiv:2206.10653 [gr-qc] .
  * Cano  _et al._ [2023a] P. A. Cano, K. Fransen, T. Hertog, and S. Maenaut, Universal Teukolsky equations and black hole perturbations in higher-derivative gravity, Phys. Rev. D 108, 024040 (2023a), arXiv:2304.02663 [gr-qc] .
  * Cardoso  _et al._ [2019] V. Cardoso, M. Kimura, A. Maselli, E. Berti, C. F. B. Macedo, and R. McManus, Parametrized black hole quasinormal ringdown: Decoupled equations for nonrotating black holes, Phys. Rev. D 99, 104077 (2019), arXiv:1901.01265 [gr-qc] .
  * McManus  _et al._ [2019] R. McManus, E. Berti, C. F. B. Macedo, M. Kimura, A. Maselli, and V. Cardoso, Parametrized black hole quasinormal ringdown. II. Coupled equations and quadratic corrections for nonrotating black holes, Phys. Rev. D 100, 044061 (2019), arXiv:1906.05155 [gr-qc] .
  * Völkel  _et al._ [2022a] S. H. Völkel, N. Franchini, and E. Barausse, Theory-agnostic reconstruction of potential and couplings from quasinormal modes, Phys. Rev. D 105, 084046 (2022a), arXiv:2202.08655 [gr-qc] .
  * Hirano  _et al._ [2024] S. Hirano, M. Kimura, M. Yamaguchi, and J. Zhang, Parametrized black hole quasinormal ringdown formalism for higher overtones, Phys. Rev. D 110, 024015 (2024), arXiv:2404.09672 [gr-qc] .
  * Völkel  _et al._ [2022b] S. H. Völkel, N. Franchini, E. Barausse, and E. Berti, Constraining modifications of black hole perturbation potentials near the light ring with quasinormal modes, Phys. Rev. D 106, 124036 (2022b), arXiv:2209.10564 [gr-qc] .
  * Franchini and Völkel [2023b] N. Franchini and S. H. Völkel, Parametrized quasinormal mode framework for non-Schwarzschild metrics, Phys. Rev. D 107, 124063 (2023b), arXiv:2210.14020 [gr-qc] .
  * Teukolsky [1973] S. A. Teukolsky, Perturbations of a rotating black hole. 1. Fundamental equations for gravitational electromagnetic and neutrino field perturbations, Astrophys. J. 185, 635 (1973).
  * Ghosh  _et al._ [2023] R. Ghosh, N. Franchini, S. H. Völkel, and E. Barausse, Quasinormal modes of nonseparable perturbation equations: The scalar non-Kerr case, Phys. Rev. D 108, 024038 (2023), arXiv:2303.00088 [gr-qc] .
  * Kimura [2020] M. Kimura, Note on the parametrized black hole quasinormal ringdown formalism, Phys. Rev. D 101, 064031 (2020), arXiv:2001.09613 [gr-qc] .
  * Cano  _et al._ [2024] P. A. Cano, L. Capuano, N. Franchini, S. Maenaut, and S. H. Volkel, Black hole quasinormal modes in higher-derivative gravity: rotation, overtones and harmonics, (2024).
  * Cano  _et al._ [2023b] P. A. Cano, K. Fransen, T. Hertog, and S. Maenaut, Quasinormal modes of rotating black holes in higher-derivative gravity, Phys. Rev. D 108, 124032 (2023b), arXiv:2307.07431 [gr-qc] .
  * Leaver [1985] E. W. Leaver, An Analytic representation for the quasi normal modes of Kerr black holes, Proc. Roy. Soc. Lond. A 402, 285 (1985).
  * github [2024] github, parametrized_qnm_framework (2024).
  * Berti  _et al._ [2006] E. Berti, V. Cardoso, and C. M. Will, On gravitational-wave spectroscopy of massive black holes with the space interferometer LISA, Phys. Rev. D 73, 064030 (2006), arXiv:gr-qc/0512160 .
  * Bardeen and Press [1973] J. M. Bardeen and W. H. Press, Radiation fields in the schwarzschild background, J. Math. Phys. 14, 7 (1973).
  * Chandrasekhar [1975] S. Chandrasekhar, On the equations governing the perturbations of the Schwarzschild black hole, Proc. Roy. Soc. Lond. A 343, 289 (1975).
  * Zouros and Eardley [1979] T. J. M. Zouros and D. M. Eardley, Instabilities of massive scalar perturbations of a rotating black hole, Annals Phys. 118, 139 (1979).
  * Detweiler [1980] S. L. Detweiler, Klein-Gordon equation and rotating black holes, Phys. Rev. D 22, 2323 (1980).
  * Dolan [2007] S. R. Dolan, Instability of the massive Klein-Gordon field on the Kerr spacetime, Phys. Rev. D 76, 084001 (2007), arXiv:0705.2880 [gr-qc] .
  * Dias  _et al._ [2022] O. J. C. Dias, M. Godazgar, J. E. Santos, G. Carullo, W. Del Pozzo, and D. Laghi, Eigenvalue repulsions in the quasinormal spectra of the Kerr-Newman black hole, Phys. Rev. D 105, 084044 (2022), arXiv:2109.13949 [gr-qc] .
  * Chandrasekhar [1984] S. Chandrasekhar, The Mathematical Theory of Black Holes, Fundam. Theor. Phys. 9, 5 (1984).
  * Dudley and Finley [1977] A. L. Dudley and J. D. Finley, Separation of Wave Equations for Perturbations of General Type-D Space-Times, Phys. Rev. Lett. 38, 1505 (1977).
  * Dudley and Finley [1979] A. L. Dudley and J. D. Finley, III, Covariant Perturbed Wave Equations in Arbitrary Type $D$ Backgrounds, J. Math. Phys. 20, 311 (1979).
  * Berti and Kokkotas [2005] E. Berti and K. D. Kokkotas, Quasinormal modes of Kerr-Newman black holes: Coupling of electromagnetic and gravitational perturbations, Phys. Rev. D 71, 124008 (2005), arXiv:gr-qc/0502065 .
  * Wagle  _et al._ [2024] P. Wagle, D. Li, Y. Chen, and N. Yunes, Perturbations of spinning black holes in dynamical Chern-Simons gravity: Slow rotation equations, Phys. Rev. D 109, 104029 (2024), arXiv:2311.07706 [gr-qc] .
  * Chung  _et al._ [2023] A. K.-W. Chung, P. Wagle, and N. Yunes, Spectral method for the gravitational perturbations of black holes: Schwarzschild background case, Phys. Rev. D 107, 124032 (2023), arXiv:2302.11624 [gr-qc] .
  * Chung  _et al._ [2024] A. K.-W. Chung, P. Wagle, and N. Yunes, Spectral method for metric perturbations of black holes: Kerr background case in general relativity, Phys. Rev. D 109, 044072 (2024), arXiv:2312.08435 [gr-qc] .
  * Blázquez-Salcedo  _et al._ [2024] J. L. Blázquez-Salcedo, F. S. Khoo, J. Kunz, and L. M. González-Romero, Quasinormal modes of Kerr black holes using a spectral decomposition of the metric perturbations, Phys. Rev. D 109, 064028 (2024), arXiv:2312.10754 [gr-qc] .
  * Chung and Yunes [2024a] A. K.-W. Chung and N. Yunes, Ringing out General Relativity: Quasi-normal mode frequencies for black holes of any spin in modified gravity, (2024a), arXiv:2405.12280 [gr-qc] .
  * Chung and Yunes [2024b] A. K.-W. Chung and N. Yunes, Quasi-normal mode frequencies and gravitational perturbations of black holes with any subextremal spin in modified gravity through METRICS: the scalar-Gauss-Bonnet gravity case, (2024b), arXiv:2406.11986 [gr-qc] .
  * Mitman  _et al._ [2023] K. Mitman _et al._ , Nonlinearities in Black Hole Ringdowns, Phys. Rev. Lett. 130, 081402 (2023), arXiv:2208.07380 [gr-qc] .
  * Cheung  _et al._ [2023] M. H.-Y. Cheung _et al._ , Nonlinear Effects in Black Hole Ringdown, Phys. Rev. Lett. 130, 081401 (2023), arXiv:2208.07374 [gr-qc] .
  * Yi  _et al._ [2024] S. Yi, A. Kuntz, E. Barausse, E. Berti, M. H.-Y. Cheung, K. Kritos, and A. Maselli, Nonlinear quasinormal mode detectability with next-generation gravitational wave detectors, Phys. Rev. D 109, 124029 (2024), arXiv:2403.09767 [gr-qc] .
  * Nollert [1993] H. P. Nollert, Quasinormal modes of Schwarzschild black holes: The determination of quasinormal frequencies with very large imaginary parts, Phys. Rev. D 47, 5253 (1993).