# Perturbative quasinormal mode frequencies

Yasuyuki Hatsuda  Department of Physics, Rikkyo University, Toshima, Tokyo
171-8501, Japan  Masashi Kimura  Department of Informatics and Electronics,
Daiichi Institute of Technology, Tokyo 110-0005, Japan  Department of Physics,
Rikkyo University, Toshima, Tokyo 171-8501, Japan

###### Abstract

We often encounter a situation that black hole solutions can be regarded as
continuous deformations of simpler ones, or modify general relativity by
continuous parameters. We develop a general framework to compute high-order
perturbative corrections to quasinormal mode frequencies in such deformed
problems. Our method has many applications, and allows to compute numerical
values of the high-order corrections very accurately. For several examples, we
perform this computation explicitly, and discuss analytic properties of the
quasinormal mode frequencies for deformation parameters.

††preprint: RUP-23-13

###### Contents

  1. I Introduction
  2. II General framework
  3. III Technical remark: the Bender-Wu approach
    1. III.1 Leading order solution
    2. III.2 First order correction
    3. III.3 On higher order corrections
  4. IV Examples
    1. IV.1 A toy model: the Rosen-Morse potential
    2. IV.2 Massive scalar perturbations
    3. IV.3 Slowly rotating black holes
    4. IV.4 Almost asymptotically flat black holes
    5. IV.5 Reissner-Nordström black holes
      1. IV.5.1 Almost chargeless limit
      2. IV.5.2 Almost extremal limit
      3. IV.5.3 An interpolating function
    6. IV.6 Parameterized black hole QNMs
    7. IV.7 Series expansion method
  5. V Outlook
  6. A Recursion relations among coefficients in parameterized QNM approach
    1. A.1 Parameterized QNM approach
    2. A.2 Ambiguity of effective potential
    3. A.3 Recursion relations for odd parity case
      1. A.3.1 Recursion relations from the Regge-Wheeler potential
      2. A.3.2 Improved recursion relation for $e_{j,k}$
    4. A.4 Reduction of the effective potential

##  I Introduction

Perturbation theory is one of the most powerful tools in physics. We have a
typical situation that a system cannot be solved analytically but its special
limit can be. Perturbation around the special limit provides us a good
approximation method and more importantly a clue to get global information on
the total system by combining with the analytic continuation in complex
analysis or asymptotic analysis. The application range of perturbation theory
is extremely wide. It is important to clarify what we can learn about from
perturbation theory.

In this work, we propose a systematic way to compute high-order perturbative
corrections to quasinormal mode (QNM) frequencies of black holes. QNMs are
solutions to linearized field equations, which satisfy purely ingoing
(outgoing) boundary conditions at the horizon (infinity), around a background
black hole spacetime. It is known that QNMs are related to the late time
behavior of the field dynamics around black holes Nakamura:1987zz ;
Kokkotas:1999bd ; Nollert:1999ji ; berti2009 ; konoplya2011 ; Ferrari:2007dd ;
Leaver:1986gd ; Andersson:1996cm ; Andersson:1995zk ; Nollert:1992ifk ;
Berti:2006wq . In many cases, one can regard some parameters of black hole
solutions as smooth deformation parameters of simpler black holes. We apply
perturbation theory for such deformation parameters.111One should not confuse
it with so-called black hole perturbation theory. We just treat deformation
parameters of black holes as formal perturbation parameters. Similar
situations also happen if one considers possibilities of effective field
theories or modified gravity theories beyond general relativity. Since such
modification parameters are expected to be small, it is natural to expand
physical quantities perturbatively. It is desirable to develop a general
framework widely applicable for such cases.

There are two obstacles to achieve it. One is that we cannot solve the QNM
spectral problem analytically even in spherically symmetric black holes.
Therefore we have only numerical or semi-analytic eigenvalues and
eigenfunctions in this simplest case. The other point is more serious. In the
QNM problem, a set of the eigenfunctions is not complete in the usual sense.
This means that we cannot apply the well-known formula in quantum mechanics to
the computation of perturbative corrections to the QNM spectrum. There is
already an extended formula to compute perturbative corrections to QNM
frequencies leung1997 ; leung1998 ; leung1999 . However it is not clear for us
how to use this formulation for our interested examples systematically and
practically. For this reason, we revisit the similar problem in this work, and
propose another way to get high-order perturbative corrections to QNM
frequencies.

A possible resolution for this problem is simply to use numerical
fittings.222There is another resolution. One can analytically continue
eigenvalue problems from the real line to the complex domain. This is well-
known as the complex scaling method for resonance problems in quantum
mechanics. However it is hard to predict high-order corrections accurately in
this way. Recently, a smart way to compute perturbative quadratic order
corrections to the QNM frequencies was proposed in cardoso2019 ; mcmanus2019 .
We are strongly motivated by these works. We extend them to more general
setups. Our approach is based on the first principle of perturbation theory.
We do not use any numerical fittings to determine the perturbative
coefficients though we need numerical solutions at each order in perturbation.
Our method is quite general and applicable to various situations with smoothly
continuous deformations. In fact, we give, for the first time, the high-
precision perturbative expansion of the QNM frequencies around the extremal
Reissner-Nordström black holes. Combining a recently proposed method
hatsuda2020 , our approach allows to compute numerical values of high-order
corrections very accurately. Once we get the high-order perturbative data, we
can discuss analytic properties (convergence, singularity, analytic
continuation, non-perturbative effect etc.) of the QNM frequencies in
principle.

The organization is as follows. In Section II, we start by explaining a
general framework of our formulation. We illustrate our basic idea to compute
high-order perturbative corrections systematically. In Section III, we present
a technical way to perform the idea in Section II explicitly. In Section IV,
we show various examples in which our method works well. We particularly use
the method proposed in hatsuda2020 , but this is not only the possibility. For
instance, we give another way in Sec. IV.7. In Section V, we consider possible
future directions. In Appendix A, we give some remarks on the so-called
parameterized black hole quasinormal mode approach. These are useful by
combining the results in the main text.

##  II General framework

We first illustrate our idea. In this section, we set up a problem, and
explain a conceptual way to obtain perturbative series of QNM frequencies
systematically. We will show a technical method to achieve it in the next
section. We expect that the problem proposed in this section is solved in many
other ways developed in numerical computations of QNMs, such as Leaver’s
continued fraction method leaver1985 ; konoplya2011 , the direct integration
method Chandrasekhar:1975zza ; konoplya2011  or the pseudospectral method
Jansen:2017oag .

We consider a perturbative deformation of a black hole in a certain theory. We
would like to know perturbative corrections to quasinormal mode frequencies
for a small deformation parameter. Our starting point is the following
(radial) master equation:333Our idea is not restricted to this form. To make
an explanation simpler, we assume it in this paper.

$\displaystyle\left(\frac{d^{2}}{dx^{2}}+\omega^{2}-V(x)\right)\Phi(x)=0,$ (1)

where $x$ is the tortoise coordinate. This is related to the radial variable
$r$ as

$\displaystyle\frac{dx}{dr}=\frac{1}{f(r)}.$ (2)

where $f(r)$ is a function that has a zero at the event horizen $r=r_{H}$.
Explicit forms of $f(r)$, of course, depend on problems. The QNM boundary
condition is then given by444 If the field has a mass term $\mu^{2}$, $\omega$
should be changed into $\sqrt{\omega^{2}-\mu^{2}}$ at $x\to\infty$ in Eq. (3).
In this case, we can still apply the same method. See subsection IV.2.

$\displaystyle\Phi(x)\sim e^{\pm i\omega x}\quad(x\to\pm\infty).$ (3)

We assume that all the quantities in the master equation have smooth
perturbative expansions in a parameter $\alpha$:

$\displaystyle
V(x)=\sum_{k=0}^{\infty}\alpha^{k}V_{k}(x),\qquad\omega^{2}=\sum_{k=0}^{\infty}\alpha^{k}\mathcal{E}_{k},\qquad\Phi(x)=\sum_{k=0}^{\infty}\alpha^{k}\Phi_{k}(x).$
(4)

Typically, the parameter $\alpha$ appears as a deformation parameter of a
black hole or of a modified theory. At this stage, we do not ask its physical
origin for generality. In general, the function $f(r)$ may also depends on
$\alpha$. This dependence causes subtlety on our perturbative treatment. We
will discuss this issue later.

Expanding $\omega$ as a series of $\alpha$,

$\displaystyle\omega=\sum_{k=0}^{\infty}\alpha^{k}\omega_{k},$ (5)

the QNM boundary condition in Eq. (3) can be written as

$\displaystyle\Phi$ $\displaystyle\sim e^{\pm i\omega_{0}x}e^{\pm
i(\alpha\omega_{1}+\alpha^{2}\omega_{2}+\cdots)x}$ $\displaystyle=e^{\pm
i\omega_{0}x}(1+\alpha P^{\pm}_{1}+\alpha^{2}P^{\pm}_{2}+\cdots),$ (6)

where $P^{\pm}_{1},P^{\pm}_{2},\cdots$ are polynomials of $x$. This implies
that the QNM boundary condition for $\Phi_{k}$ is

$\displaystyle\Phi_{k}\sim e^{\pm i\omega_{0}x}\quad(x\to\pm\infty).$ (7)

We solve the master equation perturbatively in $\alpha$. We start with the
zeroth order, at which the eigen-equation is

$\displaystyle\left(\frac{d^{2}}{dx^{2}}+\mathcal{E}_{0}-V_{0}(x)\right)\Phi_{0}(x)=0,\qquad\mathcal{E}_{0}=\omega_{0}^{2}.$
(8)

Note that $\mathcal{E}_{0}$ denotes the zeroth order eigenvalue in
perturbation of $\alpha$, not the fundamental mode eigenvalue. Typically, the
zeroth order equation is the master equation for spherically symmetric black
holes, but our formalism is not restricted to this specific situation. At each
order, we solve the differential equation by requiring proper boundary
conditions, and then get the perturbative corrections to the eigenvalues.

We first solve the zeroth order equation (8) by imposing the ordinary QNM
boundary condition:

$\displaystyle\Phi_{0}(x)\sim e^{\pm i\omega_{0}x}\quad(x\to\pm\infty).$ (9)

This is natural in perturbation theory. There are many techniques to solve (8)
numerically. To go to the next order, we need both the eigenvalue and the
eigenfunction at a fixed energy level (or an overtone number in the QNM
terminology).555This point is quite different from the textbook-like method in
quantum mechanics, in which one needs all the eigenvalues and the
eigenfunctions at the lowest order. Though, in this work, we will use a method
recently proposed in hatsuda2020 , we stress that our idea should work for
many other techniques.

Once we obtain the eigenvalue and the eigenfunction at the zeroth order, we
can proceed to the first order equation. The equation we should solve is

$\displaystyle\left(\frac{d^{2}}{dx^{2}}+\mathcal{E}_{0}-V_{0}(x)\right)\Phi_{1}(x)=(V_{1}(x)-\mathcal{E}_{1})\Phi_{0}(x).$
(10)

We regard this equation as the inhomogeneous differential equation for
$\Phi_{1}(x)$ with the unknown constant $\mathcal{E}_{1}$, while $\Phi_{0}(x)$
and $\mathcal{E}_{0}$ are known. For the function $\Phi_{1}(x)$, we impose the
same QNM boundary condition for $\Phi_{0}(x)$:

$\displaystyle\Phi_{1}(x)\sim e^{\pm i\omega_{0}x}\quad(x\to\pm\infty),$ (11)

as explained in Eq. (7). This is also natural from the perturbative point of
view. As shown in the next section, this inhomogeneous equation is also solved
by the same method as the zeroth order equation. Therefore, we get
$\mathcal{E}_{1}$ and $\Phi_{1}(x)$ at least numerically.

The computations at higher orders are similar. We regard the $k$-th order
equation

$\displaystyle\left(\frac{d^{2}}{dx^{2}}+\mathcal{E}_{0}-V_{0}(x)\right)\Phi_{k}(x)=\sum_{\ell=1}^{k}(V_{\ell}(x)-\mathcal{E}_{\ell})\Phi_{k-\ell}(x).$
(12)

as the inhomogeneous equation for $\mathcal{E}_{k}$ and $\Phi_{k}(x)$ with the
known $\mathcal{E}_{j}$ and $\Phi_{j}(x)$ ($0\leq j\leq k-1$). We solve it
under the boundary condition in Eq. (7). We repeat this computation as many
times as possible.

If the function $f(r)$ depends on the perturbative parameter $\alpha$, there
is a subtle point. In this case, we also expand $f(r)$ in $\alpha$. This gives
a perturbative relation between $r$ and $x$ via the relation (2).
Schematically, we have

$\displaystyle x=x(r,\alpha)=\sum_{k=0}^{\infty}\alpha^{k}x_{k}(r),$ (13)

where $x_{k}(r)$ are functions of $r$. On the other hand, we can inverse this
relation by

$\displaystyle r=r(x,\alpha)=\sum_{k=0}^{\infty}\alpha^{k}r_{k}(x).$ (14)

There is an ambiguity which variable, $r$ or $x$, is fundamental in the
perturbative expansion. In this paper, we regard $x$ as a fundamental
variable, and use (14) to eliminate $r$ to expand the potential
perturbatively. This is because boundary conditions in terms of $x$ seem to be
more natural.

There is a caveat when we apply our framework to a specific system and
calculate the QNM frequencies by numerical calculations. Our framework is
introduced based on the form of the master equation in Eq. (1) which is
written by the tortoise coordinate $x$. However, in many cases, it is
difficult to explicitly write the tortoise coordinate $x$ as a function of $r$
and also the master equation as a function of $x$. This implies that imposing
the boundary condition at each order $\Phi_{k}\sim e^{\pm i\omega_{0}x}$ is
not a trivial task in a concrete example. In that case, the technique to
rewrite the master equation used in cardoso2019  might be useful. When the
function $f$ has a zero at $r=r_{H}$, and it is close to $1-r_{H}/r$, we can
write $f$ as

$\displaystyle f=\left(1-\frac{r_{H}}{r}\right)Z(r;\alpha),$ (15)

where $Z(r;\alpha)$ is a function of $r$ which contains the small parameter
$\alpha$. Choosing $r_{H}$ and $\alpha$ as the fundamental parameters, we can
write the master equation in the form

$\displaystyle\left(1-\frac{r_{H}}{r}\right)\frac{d}{dr}\left(\left(1-\frac{r_{H}}{r}\right)\frac{d\phi}{dr}\right)+(\omega^{2}-\tilde{V})\phi=0,$
(16)

where $\phi=\sqrt{Z}\Phi$ and $\tilde{V}$ is the effective potential which
depends on $\alpha$ cardoso2019 . Regarding this equation as the basic master
equation, we can easily apply our framework to this system because the
tortoise coordinate in this system is explicitly written as
$r+r_{H}\ln(1-r_{H}/r)$. We should note that we do not need to care about this
point as far as we use the Bender-Wu approach introduced in the next section
because the calculation is carried out around potential peak region.

Finally note that our formulation is easily extended to multi-parameter
perturbations. If one wants to consider a two-parameter perturbation:

$\displaystyle
V(x;\alpha,\beta)=V_{0}(x)+\sum_{k=1}^{\infty}(\alpha^{k}V_{k}^{\alpha}(x)+\beta^{k}V_{k}^{\beta}(x)),$
(17)

then the square of the frequency should receive the following perturbative
corrections cardoso2019 ; mcmanus2019 :

$\displaystyle\omega^{2}$
$\displaystyle=\mathcal{E}_{0}+\alpha\mathcal{E}_{1}^{(1,0)}+\beta\mathcal{E}_{1}^{(0,1)}+\alpha^{2}\mathcal{E}_{2}^{(2,0)}+\alpha\beta\mathcal{E}_{2}^{(1,1)}+\beta^{2}\mathcal{E}_{2}^{(0,2)}+\cdots$
(18)
$\displaystyle=\mathcal{E}_{0}+\sum_{k=1}^{\infty}\sum_{\ell=0}^{k}\alpha^{\ell}\beta^{k-\ell}\mathcal{E}_{k}^{(\ell,k-\ell)}.$

To fix the coefficients $\mathcal{E}_{k}^{(\ell,k-\ell)}$, we can choose
various combinations of $(\alpha,\beta)$. For instance, to fix the second
order corrections $\mathcal{E}_{2}^{(2,0)}$, $\mathcal{E}_{2}^{(1,1)}$ and
$\mathcal{E}_{2}^{(0,2)}$, it is sufficient to consider three particular
slices: $(\alpha,\beta)\to(\alpha,0),(\alpha,\alpha),(0,\alpha)$, in which the
problem is reduced to the one-parameter problem. We will return to this issue
in Section IV.

##  III Technical remark: the Bender-Wu approach

In the previous section, we proposed a general idea to compute the
perturbative corrections $\mathcal{E}_{k}$ systematically. The main problem is
of course how we solve the differential equation (12) for our interested QNM
problems. In this section, we see that this is done by the so-called Bender-Wu
approach bender1969  that is recently extended to the QNM computation in
hatsuda2020 ; eniceicu2020 , based on blome1984 ; ferrari1984 ; Ferrari:1984zz
. The main advantage of this approach is that it is widely applicable to many
models, as in the WKB approach mashhoon1983 ; schutz1985 . The Bender-Wu
approach itself also highly depends on perturbation theory. Since we need
eigenfunctions as well as eigenvalues, we review the Bender-Wu approach for
our problem. We follow the notation in sulejmanpasic2018  as much as possible.

###  III.1 Leading order solution

Let us solve the zeroth order equation (8). We first introduce a formal
parameter $g$ by hand,

$\displaystyle\left(-g^{4}\frac{d^{2}}{dx^{2}}+\mathcal{E}_{0}-V_{0}(x)\right)\Phi_{0}(x)=0,\qquad\mathcal{E}_{0}=\omega_{0}^{2}.$
(19)

It is clear to see that $g^{2}$ plays the role of a Planck parameter. Setting
$g=e^{\pi i/4}$, the original equation (8) is reproduced. 666Note that there
is another possibility: $g=e^{-\pi i/4}$. This ambiguity reflects the fact
that the QNM frequencies have two branches for the real part hatsuda2020 . The
basic idea is the following. We first consider the eigenvalue problem for
$g\in\mathbb{R}$. In this case, we have the Schrödinger-type equation with the
inverted potential $-V_{0}(x)$, which admit bound states, and we can apply the
standard perturbative method in quantum mechanics near the minimum of
$-V_{0}(x)$. The important observation in hatsuda2020  is that the boundary
conditions for bound states and QNMs are simply related by the analytic
continuation of $g$. This implies that if we know the bound state energy for
$g\in\mathbb{R}$, we can obtain the QNM eigenvalue by the analytic
continuation $g=e^{\pi i/4}$.

Let $\bar{x}$ be the value of $x$ at which $-V_{0}(x)$ takes the minimal
value. We expand the inverted potential $-V_{0}(x)$ around $x=\bar{x}$:

$\displaystyle-V_{0}(x)=V_{00}+\sum_{j=2}^{\infty}V_{0j}(x-\bar{x})^{j}.$ (20)

We introduce a new variable by $x-\bar{x}=gq$. This change means that as $g$
decreases, we zoom in on the neighborhood of the minimum at $x=\bar{x}$. Then
(19) leads to

$\displaystyle\left(-\frac{1}{2}\frac{d^{2}}{dq^{2}}+\frac{1}{2}\Omega^{2}q^{2}+v_{0}(q)-\epsilon_{0}\right)\psi_{0}(q)=0,$
(21)

where $\Omega:=\sqrt{V_{02}}$,
$\epsilon_{0}:=-(\mathcal{E}_{0}+V_{00})/(2g^{2})$ and

$\displaystyle
v_{0}(q)=\frac{1}{2g^{2}}\sum_{j=3}^{\infty}V_{0j}(gq)^{j}=\sum_{j=1}^{\infty}g^{j}v_{0j}q^{j+2},\quad
v_{0j}:=\frac{V_{0,j+2}}{2}.$ (22)

We denoted $\psi_{0}(q)=\Phi_{0}(\bar{x}+gq)$ to avoid confusion. In this
picture, the Planck constant is unity, and $g$ now plays the role of a
coupling constant in the potential.

We solve (21) perturbatively in $g$ order by order. At the leading order, we
can regard it as the harmonic oscillator with frequency $\Omega$. To eliminate
the exponential factor of the eigenfunction, we rescale
$\psi_{0}(q)=e^{-\Omega x^{2}/2}u_{0}(q)$:

$\displaystyle-\frac{1}{2}u_{0}^{\prime\prime}(q)+\Omega
qu_{0}^{\prime}(q)+\left(\frac{\Omega}{2}+v_{0}(q)-\epsilon_{0}\right)u_{0}(q)=0.$
(23)

We have the following expansions:

$\displaystyle
u_{0}(q)=\sum_{n=0}^{\infty}g^{n}u_{0n}(q),\qquad\epsilon_{0}=\sum_{n=0}^{\infty}g^{n}\epsilon_{0n}.$
(24)

Plugging these expansions into (21), we get

$\displaystyle-\frac{1}{2}u_{0n}^{\prime\prime}+\Omega
qu_{0n}^{\prime}+\frac{\Omega}{2}u_{0n}+\sum_{j=1}^{n}v_{0j}q^{j+2}u_{0,n-j}-\sum_{j=0}^{n}\epsilon_{0j}u_{0,n-j}=0.$
(25)

Let us focus on the ground state for simplicity. The ground state corresponds
to the lowest (or fundamental) overtone mode in the QNM problem. For $n=0$, we
have the trivial solution $u_{00}(q)=1$ and $\epsilon_{00}=\Omega/2$. Using
it, we get

$\displaystyle-\frac{1}{2}u_{0n}^{\prime\prime}+\Omega
qu_{0n}^{\prime}+\sum_{j=1}^{n}(v_{0j}q^{j+2}-\epsilon_{0j})u_{0,n-j}=0,\quad
n\geq 1.$ (26)

The very important fact is that $u_{0n}$ is a polynomial of $q$ whose degree
is $3n$ bender1969 ; sulejmanpasic2018 :

$\displaystyle u_{0n}=\sum_{m=1}^{3n}A_{0n}^{m}q^{m},\quad n\geq 1.$ (27)

As shown in bender1969 , the differential equation (26) determines all the
coefficients $A_{0n}^{m}$ and $\epsilon_{0n}$ recursively. This is what the
Mathematica program in sulejmanpasic2018  is doing. One has to keep in mind
that the above result is valid only for the ground state. For the excited
states, we need to modify it slightly. See sulejmanpasic2018  for these cases.

We finally want to set $g=e^{\pi i/4}$ in the perturbative series. However, in
general, the formal power series (24) are not convergent for any $g\neq 0$.
The substitution of $g=e^{\pi i/4}$ merely gives a meaningless answer. To
avoid it, one needs to truncate all the high-order corrections beyond a
certain optimal order or to use summation methods. Note that the former turns
out to be equivalent to the WKB series in the literature mashhoon1983 ;
schutz1985 . We use the latter, called the Borel summation method, to decode a
meaningful result for finite $g$ from formal divergent series.777An
alternative way is to use Padé approximants matyjasek2017 ; konoplya2019 ;
matyjasek2019 . The conclusion in hatsuda2020  is that the Borel summation of
(24) correctly reproduces the QNM frequencies. We emphasize that the above
method allows us to construct not only the eigenvalue $\mathcal{E}_{0}$ but
also the eigenfunction $\psi_{0}(q)$. In summary, for the ground state, we
have

$\displaystyle\mathcal{E}_{0}$
$\displaystyle=-V_{00}-2g^{2}\sum_{n=0}^{\infty}g^{n}\epsilon_{0n},$ (28)
$\displaystyle\psi_{0}(q)$ $\displaystyle=e^{-\Omega
q^{2}/2}\sum_{n=0}^{\infty}g^{n}u_{0n}(q),\qquad
u_{0n}(q)=\sum_{m=1}^{3n}A_{0n}^{m}q^{m},$

where $\epsilon_{00}=\Omega/2$ and $u_{00}(q)=1$.

###  III.2 First order correction

Let us proceed to the first order correction. We need to solve

$\displaystyle\left(-g^{4}\frac{d^{2}}{dx^{2}}+\mathcal{E}_{0}-V_{0}(x)\right)\Phi_{1}(x)=(V_{1}(x)-\mathcal{E}_{1})\Phi_{0}(x).$
(29)

Note that we already know the zeroth order eigenfunction $\Phi_{0}(x)$ and
eigenvalue $\mathcal{E}_{0}$ in the previous subsection. As in the computation
above, we can rewrite it as

$\displaystyle\left(-\frac{1}{2}\frac{d^{2}}{dq^{2}}+\frac{1}{2}\Omega^{2}q^{2}+v_{0}(q)-\epsilon_{0}\right)\psi_{1}(q)=\frac{V_{1}(x)-\mathcal{E}_{1}}{2g^{2}}\psi_{0}(q),$
(30)

We also expand $-V_{1}(x)$ around $x=\bar{x}$ as

$\displaystyle-V_{1}(x)=\sum_{j=0}^{\infty}V_{1j}(gq)^{j}.$ (31)

Note that $x=\bar{x}$ does not extremize $V_{1}(x)$ in general. As mentioned
in the previous subsection, we have to impose the same boundary conditions for
$\psi_{0}(q)$ and $\psi_{1}(q)$. Therefore we set $\psi_{1}(q)=e^{-\Omega
q^{2}/2}u_{1}(q)$ as well as $\psi_{0}(q)=e^{-\Omega q^{2}/2}u_{0}(q)$, and
get

$\displaystyle-\frac{1}{2}u_{1}^{\prime\prime}+\Omega
qu_{1}^{\prime}+\left(\frac{\Omega}{2}+v_{0}-\epsilon_{0}\right)u_{1}+\left(\frac{V_{11}}{2g}q+v_{1}-\epsilon_{1}\right)u_{0}=0,$
(32)

where

$\displaystyle\epsilon_{1}$
$\displaystyle:=-\frac{\mathcal{E}_{1}+V_{10}}{2g^{2}},$ (33) $\displaystyle
v_{1}(q)$
$\displaystyle:=\frac{1}{2g^{2}}\sum_{j=2}^{\infty}V_{1j}(gq)^{j}=\sum_{j=0}^{\infty}g^{j}v_{1j}q^{j+2},\qquad
v_{1j}=\frac{V_{1,j+2}}{2}.$

We use the zeroth order perturbative solution (24). From the consistency at
the orders $1/g^{2}$ and $1/g$, we should take

$\displaystyle u_{1}(q)=-\frac{V_{11}}{2\Omega
g}q+\sum_{n=0}^{\infty}g^{n}u_{1n}(q),\qquad\epsilon_{1}=\sum_{n=0}^{\infty}g^{n}\epsilon_{1n}.$
(34)

It is observed that for the ground state, $u_{1n}(q)$ is a polynomial of
degree $3n+4$. After putting an ansatz for the polynomial $u_{1n}(q)$, we can
determine all the coefficients of $u_{1n}(q)$ and $\epsilon_{1n}$ from the
perturbative equations. The remaining computation is the same as the zeroth
order one. By performing the Borel summation of $\epsilon_{1}$, we obtain the
first correction $\mathcal{E}_{1}$.

###  III.3 On higher order corrections

The computations for higher orders are straightforward. At the $k$-th order,
we have

$\displaystyle\left(-g^{4}\frac{d^{2}}{dx^{2}}+\mathcal{E}_{0}-V_{0}(x)\right)\Phi_{k}(x)=\sum_{\ell=1}^{k}(V_{\ell}(x)-\mathcal{E}_{\ell})\Phi_{k-\ell}(x).$
(35)

It leads to

$\displaystyle-\frac{1}{2}u_{k}^{\prime\prime}+\Omega
qu_{k}^{\prime}+\left(\frac{\Omega}{2}+v_{0}-\epsilon_{0}\right)u_{k}+\sum_{\ell=1}^{k}\left(\frac{V_{\ell
1}}{2g}q+v_{\ell}-\epsilon_{\ell}\right)u_{k-\ell}=0,$ (36)

where $\Phi_{k}(x)=e^{-\Omega q^{2}/2}u_{k}(q)$ and

$\displaystyle\epsilon_{\ell}:=-\frac{\mathcal{E}_{\ell}+V_{\ell
0}}{2g^{2}},\qquad v_{\ell}(q):=\frac{1}{2g^{2}}\sum_{j=2}^{\infty}V_{\ell
j}(gq)^{j}.$ (37)

We observe that the ground state solution in general behaves as

$\displaystyle u_{k}(q)$
$\displaystyle=\frac{u_{k,-k}(q)}{g^{k}}+\cdots=\sum_{n=-k}^{\infty}g^{n}u_{kn}(q),$
(38) $\displaystyle\epsilon_{k}$
$\displaystyle=\frac{\epsilon_{k,-2}}{g^{2}}+\cdots=\sum_{n=-1}^{\infty}g^{2n}\epsilon_{k,2n},$

where $u_{kn}(q)$ is a polynomial of degree $3n+4k$. Under this assumption, we
can easily compute $\epsilon_{k}$ perturbatively in $g$.

##  IV Examples

In this section, we apply our formalism to various examples.

###  IV.1 A toy model: the Rosen-Morse potential

We demonstrate that the idea in Section II actually works in the QNM problem
for a simple exactly solvable toy model. What we consider is the so-called
Rosen-Morse potential, which is regarded as an integrable deformation of the
Pöschl-Teller potential. The Rosen-Morse potential was studied in the context
of the quasinormal modes in massive scalar perturbations ohashi2004 . We
revisit the same model to validate our framework. This model is given by

$\displaystyle\left(\frac{d^{2}}{dx^{2}}+\omega^{2}-V_{\text{RM}}(x)\right)\phi(x)=0,$
(39) $\displaystyle
V_{\text{RM}}(x)=\frac{1}{2\cosh^{2}x}+\mu^{2}\frac{1+\tanh x}{2}.$

where $\mu$ is a deformation parameter. If $\mu=0$, the potential reduces to
the well-known Pöschl-Teller potential. The Rosen-Morse potential (39) for
$\mu\neq 0$ is very similar to the potential for the spherically symmetric
black hole in the massive scalar perturbation ohashi2004 . We will see it in
the next subsection. We treat this system as a perturbation in the parameter
$\mu$.

We first show that this system is in fact exactly solvable. To do so, we
perform a change of variables and a transformation of the wave function by

$\displaystyle z=\frac{1}{2}(1+\tanh
x),\qquad\phi(x)=z^{-i\omega/2}(1-z)^{-i\sqrt{\omega^{2}-\mu^{2}}/2}y(z).$
(40)

Then, the new function $y(z)$ satisfies the standard hypergeometric equation:

$\displaystyle z(1-z)y^{\prime\prime}(z)+[c-(a+b+1)z]y^{\prime}(z)-aby(z)=0,$
(41)

where

$\displaystyle a$
$\displaystyle=\frac{1}{2}-\frac{i}{2}(\omega+\sqrt{\omega^{2}-\mu^{2}}+1),$
(42) $\displaystyle b$
$\displaystyle=\frac{1}{2}-\frac{i}{2}(\omega+\sqrt{\omega^{2}-\mu^{2}}-1),$
$\displaystyle c$ $\displaystyle=1-i\omega.$

For a given $\mu$, we impose the QNM-like boundary condition:

$\displaystyle\lim_{x\to-\infty}\phi(x)\sim e^{-i\omega
x},\qquad\lim_{x\to+\infty}\phi(x)\sim e^{+i\sqrt{\omega^{2}-\mu^{2}}\,x},$
(43)

where we have to choose a branch of the square root so that $\sqrt{z^{2}}=z$
for $z\in\mathbb{C}$ in order to match the boundary condition for $\mu=0$. In
terms of $y(z)$, this boundary condition is translated into the regularity
condition both at $z=0,1$ simultaneously. The regular solution at $z=0$ is
given by the Gauss hypergeometric function

$\displaystyle y(z)=F(a,b;c;z).$ (44)

Using the well-known analytic connection formula of the hypergeometric
function:

$\displaystyle
F(a,b,c;z)=\frac{\Gamma(c)\Gamma(c-a-b)}{\Gamma(c-a)\Gamma(c-b)}F(a,b,a+b-c+1;1-z)$
(45)
$\displaystyle\qquad\qquad+\frac{\Gamma(c)\Gamma(a+b-c)}{\Gamma(a)\Gamma(b)}(1-z)^{c-a-b}F(c-a,c-b,c-a-b+1;1-z),$

the regularity condition at $z=1$ requires

$\displaystyle\frac{1}{\Gamma(a)\Gamma(b)}=0.$ (46)

Therefore we obtain $a=-n$ or $b=-n$ for $n=0,1,2,\dots$. This condition leads
to the following exact spectrum:

$\displaystyle\omega^{(n,\pm)}=\pm\left(\frac{1}{2}+\mu^{2}\frac{1}{4(2n^{2}+2n+1)}\right)-i\left(n+\frac{1}{2}-\mu^{2}\frac{2n+1}{4(2n^{2}+2n+1)}\right).$
(47)

We have two symmetric branches of the spectra. The exact eigenfunction is also
given by

$\displaystyle\phi^{(n,\pm)}(x)=\left(\frac{1+\tanh
x}{2}\right)^{-i\omega^{(n,\pm)}}\left(\frac{1-\tanh
x}{2}\right)^{-i\sqrt{\omega^{(n,\pm)2}-\mu^{2}}/2}$ (48) $\displaystyle\times
F\left(-n,-n\mp i;1-i\omega^{(n,\pm)};\frac{1+\tanh x}{2}\right)$

Note that for a non-negative integer $n$, the hypergeometric function in this
equation is a polynomial of degree $n$. For simplicity, we consider the case
of $b=-n$, and abbreviate the upper index in these expressions. For the lowest
overtone number $n=0$, we have

$\displaystyle\omega$ $\displaystyle=\frac{1-i}{2}+\mu^{2}\frac{1+i}{4},$ (49)
$\displaystyle\phi(x)$ $\displaystyle=\left(\frac{1+\tanh
x}{2}\right)^{-i\omega/2}\left(\frac{1-\tanh
x}{2}\right)^{-i\sqrt{\omega^{2}-\mu^{2}}/2}.$

In the small $\mu$ limit, we have

$\displaystyle\omega^{2}$
$\displaystyle=\mathcal{E}_{0}+\mu^{2}\mathcal{E}_{1}+\mu^{4}\mathcal{E}_{2}=-\frac{i}{2}+\frac{\mu^{2}}{2}+\frac{i\mu^{4}}{8},$
(50) $\displaystyle\phi(x)$
$\displaystyle=\phi_{0}(x)+\mu^{2}\phi_{1}(x)+\mu^{4}\phi_{2}(x)+\mathcal{O}(\mu^{6}),$

where

$\displaystyle\phi_{0}(x)$ $\displaystyle=\left(\frac{1}{2\cosh
x}\right)^{-i\omega_{0}},$ (51) $\displaystyle\phi_{1}(x)$
$\displaystyle=\frac{1-i}{4}x\left(\frac{1}{2\cosh x}\right)^{-i\omega_{0}},$
$\displaystyle\phi_{2}(x)$
$\displaystyle=-\frac{i}{16}x^{2}\left(\frac{1}{2\cosh
x}\right)^{-i\omega_{0}},$

and $\omega_{0}=(1-i)/2$. These functions satisfy the same boundary condition:

$\displaystyle\lim_{x\to-\infty}\phi_{k}(x)\sim
e^{-i\omega_{0}x},\qquad\lim_{x\to+\infty}\phi_{k}(x)\sim
e^{+i\omega_{0}x},\qquad k=0,1,2,\dots.$ (52)

Note that this boundary condition is slightly different from the true QNM
boundary condition (43), but after resumming the perturbative series it is
reproduced correctly.

Now we confirm this result from perturbation theory. We consider the
perturbation in $\mu^{2}$:

$\displaystyle V_{\text{RM}}(x)$ $\displaystyle=V_{0}(x)+\mu^{2}V_{1}(x),$
(53) $\displaystyle V_{0}(x)$ $\displaystyle=\frac{1}{2\cosh^{2}x},\qquad
V_{1}(x)=\frac{1+\tanh x}{2}.$

At the lowest order, we of course obtain the Pöschl-Teller potential:

$\displaystyle\left(\frac{d^{2}}{dx^{2}}+\mathcal{E}_{0}-V_{0}(x)\right)\phi_{0}(x)=0,$
(54)

Its eigenvalue and the eigenfunction for the fundamental QNM are exactly given
by the zeroth order in (50) and (51). At the first and the second orders, we
have

$\displaystyle\left(\frac{d^{2}}{dx^{2}}+\mathcal{E}_{0}-V_{0}(x)\right)\phi_{1}(x)+(\mathcal{E}_{1}-V_{1}(x))\phi_{0}(x)=0,$
(55)
$\displaystyle\left(\frac{d^{2}}{dx^{2}}+\mathcal{E}_{0}-V_{0}(x)\right)\phi_{2}(x)+(\mathcal{E}_{1}-V_{1}(x))\phi_{1}(x)+\mathcal{E}_{2}\phi_{0}(x)=0.$

We would like to solve these inhomogeneous equations under the boundary
condition (52). Instead, it is sufficient to confirm that the functions in
(50) and (51) satisfy these differential equations. One can immediately check
it.

For higher overtone modes, since the hypergeometric function in (48) does not
change the asymptotic behavior of the solution, the same structure holds.

###  IV.2 Massive scalar perturbations

The simplest example in black hole problems is a massive scalar perturbation
of the Schwarzschild geometry. The functions in the master equation are given
by

$\displaystyle f(r)=1-\frac{2M}{r},\qquad
V(x)=f(r)\left(\frac{\ell(\ell+1)}{r^{2}}+\frac{2M}{r^{3}}+\mu^{2}\right).$
(56)

As in the Rosen-Morse potential, we regard the scalar mass square $\mu^{2}$ as
a deformation parameter: $\alpha=\mu^{2}$. Note that the function $f(r)$ does
not receive any correction. The explicit relation between $r$ and $x$ is given
by

$\displaystyle x=r+2M\log\left(\frac{r}{2M}-1\right).$ (57)

We regard $r$ as a function of $x$. The unperturbed system is just the
massless scalar case:

$\displaystyle
V_{0}(x)=f(r)\left(\frac{\ell(\ell+1)}{r^{2}}+\frac{2M}{r^{3}}\right).$ (58)

The correction in the potential is

$\displaystyle V_{1}(x)=f(r),\qquad V_{k\geq 2}(x)=0.$ (59)

The QNM frequency receives the perturbative corrections in $\mu^{2}$. To keep
the generality of $M$, we write the perturbative series as the dimensionless
form

$\displaystyle M\omega=\sum_{k=0}^{\infty}(M\mu)^{2k}w_{k},$ (60)

where the correction coefficients $w_{k}$ does not depend on $M$. Our task is
to compute $w_{k}$ order by order. We can apply the method in Section III.

Let us briefly see the boundary condition. In the case of (56), the total
boundary condition for the QNM is

$\displaystyle\lim_{x\to-\infty}\Phi(r)\sim e^{-i\omega
x},\qquad\lim_{x\to+\infty}\Phi(r)\sim e^{+i\sqrt{\omega^{2}-\mu^{2}}\,x},$
(61)

If $\mu$ is small, the boundary condition at infinity is expanded as

$\displaystyle e^{+i\sqrt{\omega^{2}-\mu^{2}}\,x}=e^{+i\omega
x}\left(1-\frac{ix}{2\omega}\mu^{2}-\frac{(i+\omega
x)x}{8\omega^{3}}\mu^{4}+\mathcal{O}(\mu^{6})\right).$ (62)

This is indeed consistent with our requirement (7).

To show an explicit result, we focus on the cases of $\ell=2,3$.888As
explained in hatsuda2020 , the Bender-Wu approach works well for larger
$\ell$. This is why we consider $\ell=2,3$ rather than $\ell=0,1$. It is
desirable to solve (12) in other approaches. It is sufficient for us to
compute the coefficients in (60) for the case of $M=1$ actually. The zeroth
order frequency for the lowest overtone number999The reader should not confuse
the subscript index here with the overtone number. is well-known:

$\displaystyle w_{0}^{\ell=2}=0.4836438722-0.0967587760i,\quad
w_{0}^{\ell=3}=0.6753662325-0.0964996277i.$ (63)

We have computed the numerical values of the perturbative coefficients $w_{k}$
up to $k=40$. The first six values are shown in Table 1. In this table, we
showed stable digits in our numerical computations. The leading and next-to-
leading corrections are consistent with the early results in cardoso2019 ;
mcmanus2019 .

Table 1: The first six perturbative corrections to the fundamental QNM frequency (60) with $\ell=2,3$ in the massive scalar perturbation. $k$ | $w_{k}^{\ell=2}$ | $4^{k}w_{k}^{\ell=3}$
---|---|---
$0$ | $0.4836438722-0.0967587760i$ | $0.6753662325-0.0964996277i$
$1$ | $0.3156326579+0.1081551348i$ | $0.9437297621+0.2278771948i$
$2$ | $0.03541170393+0.02620890155i$ | $0.2263735226+0.1075217988i$
$3$ | $0.01199156679+0.02204684913i$ | $0.2085153094+0.1986390780i$
$4$ | $0.00092115819+0.02209374509i$ | $0.2333370885+0.4509679860i$
$5$ | $-0.01001596605+0.02211024342i$ | $0.1500437709+1.0963976002i$
$6$ | $-0.02390151862+0.01898789685i$ | $-0.580414699+2.681826119i$

What do we learn about from these perturbative data? The most basic question
would be whether the perturbative series (60) is convergent or not. To see it,
we show the behavior of the ratio $w_{k-1}/w_{k}$ up to $k=40$ in figure 1.
The ratio seems to converge to a finite value, but the convergence is slow.
Using basic knowledge of complex analysis, we can estimate the radius of
convergence in a different way. The radius of convergence is determined by the
nearest singular point from the origin. In our framework, we have only the
finite number of $w_{k}$. We would like to decode the singularity structure
from these data. Probably the best tool to do so is Padé approximants.

Figure 1: To see whether the perturbative series (60) is convergent or not, we
plot the ratio $|w_{k-1}|/|w_{k}|$ for $1\leq k\leq 40$. It looks to converge
to a finite value.

Padé approximants tell us the analytic structure of a given power series. In
particular, it gives us information on sigularity structure on the original
function. See Appendix C in hatsuda2021a , for instance. Since we have the
perturbative data of (60) up to $(M\mu)^{80}$, we can construct its diagonal
Padé approximant $M\omega^{[40/40]}$. We read off the zeros and the poles of
this approximant. The results are illustrated in Figure 2. These figures imply
that the perturbative series (60) is likely a convergent series. One can
estimate its radius of convergence by computing the distance to the nearest
singular point. In this computation, one has to be care about “false” singular
points of Padé approximants. These singular points disappear if orders of Padé
approximants are changed. These are artifacts in the approximant, while the
“true” singular points are stable for Padé orders. In Figure 2, we observe
that the black dashed circle is expected to be the convergence circle. The
estimation of the radius of convergence $R$ for (60) in the complex
$M\mu$-plane is approximately given by

$\displaystyle R_{\text{fund.}}^{\ell=2}\approx 0.643,\qquad
R_{\text{fund.}}^{\ell=3}\approx 0.900.$ (64)

We do not a clear physical meaning of this radius so far. It would be
interesting to understand it.

Figure 2: The singularity structure of the $[40/40]$ Padé approximant of (60)
for $\ell=2$ (Left) and $\ell=3$ (Right) in the complex $M\mu$-plane. We show
its zeros by the blue points and poles by the orange points. The dashed curve
is a conjectural convergence circle of the perturbative series (60). Note that
the zeros and the poles inside the circle disappear when the degrees of the
Padé approximant are varied. These are artifacts for the $[40/40]$ Padé
approximant.

By using the Padé approximants, we finally extrapolate our perturbative
results to the finite parameter region, as shown in Figure 3.

Figure 3: The mass dependence of the $\ell=2$ fundamental QNM frequency for
the massive scalar perturbation. The (red) points represent the numerical
values. The (orange) dashed line and the (blue) solid line are the
perturbative series (60) up to $k=40$ and its diagonal Padé approximant,
respectively. The Padé approximant is extrapolated beyond the radius of
convergence.

###  IV.3 Slowly rotating black holes

Another simple application is the Kerr geometry. We regard the angular
momentum as a deformation parameter. Here we consider the slow rotation limit.
We briefly explain how to get the slow rotation expansion of the QNM frequency
reported in hatsuda2021 .

The perturbation of the rotating black holes are governed by the Teukolsky
equation teukolsky1972 . In hatsuda2021 , an isospectral equation to the
Teukolsky equation was proposed. This isospectral equation is much more useful
for our purpose in this paper. We start with the radial master equation

$\displaystyle\left(\frac{d^{2}}{dx^{2}}+(2M\omega)^{2}-V(x)\right)\Phi(x)=0,$
(65)

where

$\displaystyle V(x)$
$\displaystyle=f(z)\biggl{[}4c^{2}+\frac{4c(m-c)}{z}+\frac{{}_{s}A_{\ell
m}(c)+s(s+1)-c(2m-c)}{z^{2}}-\frac{s^{2}-1}{z^{3}}\biggr{]},$ (66)
$\displaystyle f(z)$ $\displaystyle=1-\frac{1}{z},\qquad x=z+\log(z-1)$

and $c=a\omega$ is related to the rotation parameter $a$. For the notational
detail, see hatsuda2021 . Of course, the slow rotation limit corresponds to
the small $c$ limit. The separation constant ${}_{s}A_{\ell m}(c)$ is
determined by the regularity condition of the angular master equation at
$\xi=\pm 1$:

$\displaystyle\biggl{[}\frac{d}{d\xi}(1-\xi^{2})\frac{d}{d\xi}+(c\xi)^{2}-2cs\xi+{}_{s}A_{\ell
m}(c)+s-\frac{(m+s\xi)^{2}}{1-\xi^{2}}\biggr{]}{}_{s}S_{\ell m}(\xi)=0.$ (67)

To compute the small $c$ expansion of the potential, we need the perturbative
series of ${}_{s}A_{\ell m}(c)$. This can be done as follows. In $c\to 0$, the
angular master equation can be solved exactly. The regular solution at
$\xi=\pm 1$ exists only for the discrete eigenvalue

$\displaystyle{}_{s}A_{\ell m}(0)$ $\displaystyle=\ell(\ell+1)-s(s+1),$ (68)

and the exact eigenfunction is given by

$\displaystyle{}_{s}S_{\ell m}^{(c=0)}(\xi)$
$\displaystyle=(1-\xi)^{-\frac{m+s}{2}}(1+\xi)^{\frac{m-s}{2}}P_{\ell+s}^{(-m-s,m-s)}(\xi),$
(69)

where $P_{n}^{(\alpha,\beta)}(z)$ is the Jacobi polynomial. We have assumed
$\ell\geq|s|$ and $|m|\leq\ell$. As in the very similar treatment in the
Bender-Wu approach, the eigenvalue ${}_{s}A_{\ell m}(c)$ and the eigenfunction
${}_{s}S_{\ell m}(\xi)$ admit the perturbative series in $c$:

$\displaystyle{}_{s}A_{\ell m}(c)=\sum_{k=0}^{\infty}c^{k}{}_{s}A_{\ell
m}^{(k)},\qquad{}_{s}S_{\ell m}(\xi)=\sum_{k=0}^{\infty}c^{k}{}_{s}S_{\ell
m}^{(k)}(\xi).$ (70)

The crucial step is to find the following general structure of the regular
function ${}_{s}S_{\ell m}^{(k)}(\xi)$:

$\displaystyle{}_{s}S_{\ell
m}^{(k)}(\xi)=(1-\xi)^{-\frac{m+s}{2}}(1+\xi)^{\frac{m-s}{2}}{}_{s}Q_{\ell
m}^{(k)}(\xi),$ (71)

where ${}_{s}Q_{\ell m}^{(k)}(\xi)$ is a polynomial of degree $\ell+s+k$ in
$\xi$. From the differential equation (67), we can fix all the coefficients in
the polynomial ${}_{s}Q_{\ell m}^{(k)}(\xi)$ and ${}_{s}A_{\ell m}^{(k)}$
order by order. This method allows us to compute the exact value of
${}_{s}A_{\ell m}^{(k)}$ up to very high orders for given $s$, $\ell$ and $m$.
We have confirmed that the first few coefficients indeed agree with the
results in berti2006 ; berti2006a .

Once we know the small $c$ expansion of ${}_{s}A_{\ell m}(c)$, we obtain the
perturbative expansion of the potential $V(x)$. Then we can apply the method
in Section II. The result is given by the following small $c$ expansion:

$\displaystyle M{}_{s}\omega_{\ell m}=\sum_{k=0}^{\infty}c^{k}{}_{s}v_{\ell
m}^{(k)}.$ (72)

However, we are interested in the perturbative expansion in terms of the
rotation parameter $a$ rather than $c=a\omega$. This expansion is easily
obtained by plugging (72) into $c=a\omega$ and by inversely expanding $c$ in
$a/M$. We finally obtain the following perturbative series

$\displaystyle M{}_{s}\omega_{\ell
m}=\sum_{k=0}^{\infty}\left(\frac{a}{M}\right)^{k}{}_{s}w_{\ell m}^{(k)},$
(73)

where the explicit values of ${}_{s}w_{\ell m}^{(k)}$ for
$(s,\ell,m)=(-2,2,0),(-2,2,1),(-2,2,2)$ up to $k=12$ are found in Table 1 in
hatsuda2021 .

###  IV.4 Almost asymptotically flat black holes

We can also apply our formalism to asymptitically non-flat geometries. We
focus on the Schwarzschild de Sitter black holes. In this case, the functions
in the minimally coupled massless scalar/vector/odd-parity gravitational
perturbations are all given by

$\displaystyle f(r)$ $\displaystyle=1-\frac{2M}{r}-\frac{\Lambda r^{2}}{3},$
(74) $\displaystyle V(x)$
$\displaystyle=f(r)\left(\frac{\ell(\ell+1)}{r^{2}}+(1-s^{2})\left(\frac{2M}{r^{3}}-\frac{4-s^{2}}{6}\Lambda\right)\right),$

where $s=0,1,2$ denotes the spin-weight of the perturbation fields, and
$\Lambda$ is the cosmological costant. We regard $\Lambda$ as a deformation
parameter. In contrast to the previous examples, the function $f(r)$ depends
on $\Lambda$. The explicit relation between $r$ and $x$ is now quite
complicated. As discussed in Section II, we have to use the relation (14) to
eliminate $r$. This can be done at least perturbatively with respect to
$\Lambda$. After this prescription, the potential in terms of $x$ receives an
infinite number of perturbative corrections. We apply the Bender-Wu approach
for such a perturbative series of the potential. In the Bender-Wu approach, we
need the Taylor series of the perturbative corrections to the potential around
the extremal point $x=\bar{x}$ of the zeroth potential. This can be done
systematically.

We expand the frequency as

$\displaystyle M\omega=\sum_{k=0}^{\infty}(9M^{2}\Lambda)^{k}w_{k}.$ (75)

The numerical values of $w_{k}$ for the fundamental mode with $\ell=2$ in the
gravitational perturbation ($s=2$) up to $k=8$ are given in Table 2.

Table 2: The first eight perturbative corrections to the fundamental QNM frequency (75) with $\ell=2$ in the odd-parity gravitational perturbation for the asymptotically dS black holes. It turns out that the same values are also obtained by the even-parity perturbation. $k$ | $w_{k}$
---|---
$0$ | $0.3736716844-0.0889623157i$
$1$ | $-0.1864855559+0.0372042528i$
$2$ | $-0.04819480629+0.01428258071i$
$3$ | $-0.02302643485+0.00713463072i$
$4$ | $-0.01415049627+0.00398414719i$
$5$ | $-0.010032759238+0.002550521089i$
$6$ | $-0.007668666891+0.001893042626i$
$7$ | $-0.006085692144+0.001548612387i$
$8$ | $-0.004939500648+0.001314426006i$

A non-trivial test of our result is to check the isospectrality between the
odd-parity and even-parity gravitational perturbations. The potential in the
even-parity gravitational perturbation is

$\displaystyle
V^{\text{even}}(x)=f(r)\frac{2}{r^{3}}\frac{9M^{3}+3\lambda^{2}Mr^{2}+\lambda^{2}(1+\lambda)r^{3}+9M^{2}\lambda
r-3M^{2}\Lambda r^{3}}{(3M+\lambda r)^{2}},$ (76)

where $\lambda=(\ell-1)(\ell+2)/2$. It is well-known that the QNM spectra in
the odd/even-parity perturbations are exactly same. The reason behind this
remarkable fact is a supersymmetric structure. See appendix A in berti2009 .
Our formalism is also applicable to this potential, and we have checked that
the isospectrality indeed holds at the perturbative level at least up to
$k=8$:

$\displaystyle w_{k}^{\text{odd}}=w_{k}^{\text{even}}.$ (77)

This is an evidence of the validity of our method.

Let us discuss the extrapolation of (75) to finite $\Lambda$. We first observe
that the perturbative series is likely convergent, but it is hard to guess the
radius of convergence from the coefficient $w_{k}$. We consider the $[4/4]$
Padé approximant by using the values in Table 2. The Padé approximant
$\omega^{[4/4]}$ for $(s,\ell)=(2,2)$ has four poles at

$\displaystyle M^{2}\Lambda=0.101-0.0134i,\quad 0.142+0.00389i,\quad
0.323+0.0678i,\quad 2.45+0.687i,$ (78)

where the first pole is relatively close to $M^{2}\Lambda=1/9$, at which the
event horizon and the de Sitter horizon coincide. It is expected that higher-
order Padé approximants capture this observation more precisely, but it is
technically difficult to check it at the moment. This observation implies that
the radius of convergence of (75) is just $|M^{2}\Lambda|=1/9$.

The extrapolation of (75) by its Padé approximant is compared to the numerical
value of the QNM frequency directly computed from (74). For
$M^{2}\Lambda=0.06$, we have

$\displaystyle M\omega^{[4/4]}_{s=2,\ell=2}(M^{2}\Lambda=0.06)\approx
0.2533-0.06304i,$ (79)

which agrees with the WKB result in zhidenko2003  and also a recent high-
precision computation in hatsuda2020d .

Figure 4: The cosmological constant dependence of the fundamental QNM
frequency for the asymptotically Schwarzschild de Sitter black holes. The
(red) points represent the numerical values, while the (blue) solid line
represents the $[4/4]$ Padé approximant.

We should note that the QNM spectral problem becomes quite different for
$\Lambda>0$ (dS) and $\Lambda<0$ (AdS). The boundary condition in the AdS case
is much more involved than the dS case cardoso2001 ; berti2009 . In this
paper, we restrict ourselves to the dS case for simplicity. It would be
interesting to clarify a physical meaning of a naïve continuation to
$\Lambda<0$ of our result. Another perturbative treatment for the (A)dS
spectral problem will be also found in Aminov:2023jve .

###  IV.5 Reissner-Nordström black holes

The spectrum for the Reissner-Nordström black holes are more involved. The
master equation in the odd-parity gravitational perturbation consists of

$\displaystyle f(r)$ $\displaystyle=1-\frac{2M}{r}+\frac{Q^{2}}{r^{2}},$ (80)
$\displaystyle V(x)$
$\displaystyle=f(r)\left(\frac{\ell(\ell+1)}{r^{2}}-\frac{q}{r^{3}}+\frac{4Q^{2}}{r^{4}}\right),$

where

$\displaystyle q=3M+\sqrt{9M^{2}+4Q^{2}(\ell-1)(\ell+2)}.$ (81)

We have two characteristic regimes: $Q=0$ and $Q=M$. We discuss perturbative
series around these two points.

####  IV.5.1 Almost chargeless limit

The first is the small charge expansion. In this case, $Q^{2}$ is a natural
deformation parameter. We write the perturbative QNM frequency as

$\displaystyle M\omega=\sum_{k=0}^{\infty}\left(\frac{Q}{M}\right)^{2k}w_{k}.$
(82)

The potential receives an infinite number of corrections. The strategy is the
same as that in the previous subsection. We show the numerical values of the
perturbative coefficients $w_{k}$ for the fundamental QNM frequency with
$\ell=2$ up to $k=4$ in Table 3. The quadratic correction $w_{1}$ matches well
with cardoso2019 .

Table 3: The low order corrections to the fundamental QNM frequency for $\ell=2$ in the Reissner-Nordström gravitational perturbation. We consider the two distinct perturbative series (82) and (83). $k$ | $w_{k}$ | $w_{k}^{\text{ext}}$
---|---|---
$0$ | $0.3736716844-0.0889623157i$ | $0.4313408007-0.0834603151i$
$1$ | $0.02581767285-0.00282403214i$ | $-0.2070138464-0.0853606869i$
$2$ | $0.02518778870+0.00020532453i$ | $0.2543444995+0.4939946909i$
$3$ | $-0.004748170246+0.002508402108i$ | $0.758606111-1.429576400i$
$4$ | $0.01557265014+0.00041287974i$ | $-6.158687644+0.575432188i$

####  IV.5.2 Almost extremal limit

We can also consider another limit $Q\to M$. In this case, $1-Q/M$ is a good
parameter. Therefore we write the frequency as

$\displaystyle
M\omega=\sum_{k=0}^{\infty}\alpha^{k}w_{k}^{\text{ext}},\qquad\alpha:=1-\frac{Q}{M}.$
(83)

Now we have

$\displaystyle
f(r)=\left(1-\frac{M}{r}\right)^{2}-\alpha\frac{2M^{2}}{r^{2}}+\alpha^{2}\frac{M^{2}}{r^{2}}.$
(84)

We also expand the potential perturbatively with respect to $\alpha$. The QNM
frequencies in the strictly extremal case ($\alpha=0$) can be computed by the
Bender-Wu approach hatsuda2020 . We do the same computation for high-order
corrections. The numerical values of $w_{k}^{\text{ext}}$ for the fundamental
QNM frequency with $\ell=2$ up to $k=4$ are shown in Table 3. The zeroth order
coefficients $w_{0}^{\text{ext}}$ agrees with the early result onozawa1996 .
We do not find any references on the perturbative corrections near the
extremal limit.

####  IV.5.3 An interpolating function

We have the two perturbative expansions of the same spectrum in the different
regimes. In each regime, we determine its Padé approximant, and can
extrapolate it to the other regime. However, to know the global behavior,
there is a better approximation, called multi-point Padé approximants
bender1978 ; bakerjr.1996 . Let us consider a rational function

$\displaystyle
M\omega^{[p/q]}=\frac{a_{0}+a_{1}Q/M+\cdots+a_{p}(Q/M)^{p}}{1+b_{1}Q/M+\cdots+b_{q}(Q/M)^{q}}.$
(85)

We fix the coefficients $a_{n}$ and $b_{n}$ so that the rational function
reproduces the both perturbative expansions around $Q/M=0$ and $Q/M=1$. For
instance to get the rational function $M\omega^{[4/4]}$ we totally need nine
data in (82) and in (83). A balanced choice is to take $w_{k}$ ($0\leq k\leq
2$) in (82) and $w_{k}^{\text{ext}}$ ($0\leq k\leq 3$) in (83). Recall the
expansion (82) has no odd-order terms. We can use this information to fix
$a_{n}$ and $b_{n}$. The explicit values of $a_{n}$ and $b_{n}$ in this case
are shown in Table 4. The interpolating function remarkably reproduces the
numerical values in the whole regime $0\leq Q/M\leq 1$, as shown in Figure 5.

Interpolating functions will be improved if one considers further perturbative
expansions around other points in the middle region. For instance, a
perturbative expansion around $Q/M\sim 0.8$ will provide us an important
information on the global structure of the imaginary part of the QNM frequency
for $\ell=2$. We do not compute it in this work, but expect that our method is
still applicable in such situations.

Table 4: The nine coefficients in the rational approximation $M\omega^{[4/4]}$ for the $\ell=2$ fundamental mode. $n$ | $a_{n}$ | $b_{n}$
---|---|---
$0$ | $0.3736716844-0.0889623157i$ |
$1$ | $-0.349769907+0.062882011i$ | $-0.92374126-0.051639318i$
$2$ | $-0.342112665-0.038824176i$ | $-0.91011322-0.313017895i$
$3$ | $0.492170748-0.023942699i$ | $1.32244473+0.24735505i$
$4$ | $-0.169504965+0.033347865i$ | $-0.454637541-0.004795289i$

Figure 5: The (red) points represent the numerical values of the QNM frequency
of the RN black holes. The (blue) solid curve is the graph of the rational
function (85) for $p=q=4$ with the coefficients in Table 4. The (orange)
dashed and (black) dotted lines represent the perturbative expansions (82) and
(83) up to $k=4$, respectively.

###  IV.6 Parameterized black hole QNMs

Recently, a simple and effective way to compute perturbative corrections was
proposed in cardoso2019 ; mcmanus2019 ; kimura2020 . We refer to it as the
parameterized QNM approach. As one can see in the previous examples, most
deformation terms in the potential take the form as linear combinations of
$1/r^{j}$ with integral $j$. At the first order in the perturbation,
corrections to the QNM frequencies are the same linear combinations of the
potential. See (17) and (18). The main idea of the parameterized QNM approach
is the following. We make a list of corrections generated by only the
$1/r^{j}$-deformations beforehand, and use it for a more complicated potential
to which corrections are linear combinations of the $1/r^{j}$-deformations.
The extension to high-order corrections is straightforward mcmanus2019 .
Physical applications of the parameterized QNM approach have been shown in
Tattersall:2019nmh ; hatsuda2020b ; deRham:2020ejn ; Volkel:2022aca ;
Volkel:2022khh ; Franchini:2022axs ; Lahoz:2023csk ; Ghosh:2023etd ;
Mukohyama:2023xyf ; Franchini:2023eda . (See also Appendix. A for
complementary discussion.)

At the technical level, it is not so easy to compute the precise values of the
quadratic corrections. In cardoso2019 ; mcmanus2019 , the authors used
numerical fittings. Since our formalism is easily applied to the setup of the
parameterized QNM approach, we re-evaluate the corrections up to the quadratic
order. This re-evaluation played an important role in the computation of
perturbative corrections for slowly rotating black holes hatsuda2020b . We
keep at least ten-digit precision for all the corrections listed in this
section. We focus on deformation of the odd-parity gravitational perturbation
of the Schwarzschild black holes. The computations for the other cases are
straightforward. The potential is

$\displaystyle
V_{0}(x)=f(r)\left(\frac{\ell(\ell+1)}{r^{2}}-\frac{3r_{H}}{r^{3}}\right),\quad
V_{1}(x)=\frac{f(r)}{r_{H}^{2}}\left(\frac{r_{H}}{r}\right)^{j},\quad V_{k\geq
2}(x)=0,$ (86)

where $r_{H}$ is the location of the event horizon and $j=0,1,2,\dots$. For
this deformation, the spectrum receives the corrections:

$\displaystyle\omega=\omega_{0}+\sum_{k=1}^{\infty}\alpha^{k}e_{j}^{(k)}.$
(87)

For $\ell=2$, we show the numerical values of $e_{j}^{(k)}$ ($0\leq j\leq 8$,
$k=1,2$) in Table 5.

Table 5: The one-parameter corrections up to the second order for $\ell=2$ in the parameterized QNM approach. $k$ | $j$ | $r_{H}e_{j}^{(k)}$
---|---|---
| $0$ | $0.2472519654+0.0926430738i$
| $1$ | $0.1598547870+0.0182084818i$
| $2$ | $0.09663224013-0.00241549645i$
| $3$ | $0.05849078501-0.00371786129i$
$1$ | $4$ | $0.03667943678-0.00043869695i$
| $5$ | $0.02403794775+0.00273079314i$
| $6$ | $0.01634281096+0.00484267168i$
| $7$ | $0.011363575081+0.006013991932i$
| $8$ | $0.007951997735+0.006536996457i$
| $0$ | $0.002868401222-0.001011345890i$
| $1$ | $-0.01439027937-0.00572350838i$
| $2$ | $-0.005756554781+0.000336740545i$
| $3$ | $-0.0006273259154-0.0004693348600i$
$2$ | $4$ | $0.0007234494450-0.0011595941966i$
| $5$ | $0.000987182421-0.001122519006i$
| $6$ | $0.0010046849768-0.0008403243677i$
| $7$ | $0.0009526541187-0.0005456646402i$
| $8$ | $0.0008715569057-0.0003017937415i$

To make a list at the quadratic order, we also have to consider two-parameter
perturbations (17) with

$\displaystyle V_{1}^{\alpha}(x)$
$\displaystyle=\frac{f(r)}{r_{H}^{2}}\left(\frac{r_{H}}{r}\right)^{i},\qquad
V_{k\geq 2}^{\alpha}(r)=0,$ (88) $\displaystyle V_{1}^{\beta}(x)$
$\displaystyle=\frac{f(r)}{r_{H}^{2}}\left(\frac{r_{H}}{r}\right)^{j},\qquad
V_{k\geq 2}^{\beta}(r)=0.$

For this perturbation, the frequency receives the corrections:

$\displaystyle\omega=\omega_{0}+\sum_{k=1}^{\infty}\sum_{\ell=0}^{k}\alpha^{\ell}\beta^{k-\ell}e_{ij}^{(\ell,k-\ell)}.$
(89)

where we have $e_{ij}^{(k,0)}=e_{i}^{(k)}$ and $e_{ij}^{(0,k)}=e_{j}^{(k)}$ by
construction. Therefore at the second order, the only unknown coefficient is
$e_{ij}^{(1,1)}$. This can be evaluated by the trick explained in (18). The
numerical values are shown in Table 6. We compare these results with
cardoso2019 ; mcmanus2019 , and found that there are significant differences.

For the error estimation of the coefficients in Tables 5 and 6, we use the
recursion relations among the coefficients in Eqs. (159) and (160). We checked
that (159) is satisfied at ${\cal O}(10^{-15})$ for linear coefficients, and
(160) is satisfied at ${\cal O}(10^{-11})$ for quadratic coefficients. This
also shows that the perturbative approach developed in the present paper works
well.

Table 6: The off-diagonal quadratic corrections for $\ell=2$ in the two-parameter perturbation. $i$ | $j$ | $r_{H}e_{ij}^{(1,1)}$
---|---|---
$0$ | $1$ | $-0.02588238896-0.02792966573i$
| $2$ | $-0.03870432587-0.02320896618i$
| $3$ | $-0.03739171923-0.01523959074i$
| $4$ | $-0.03119143980-0.01062473399i$
| $5$ | $-0.02473633363-0.00886735210i$
| $6$ | $-0.01939362275-0.00853499059i$
| $7$ | $-0.01523819641-0.00870770164i$
$1$ | $2$ | $-0.02293084111-0.00341311941i$
| $3$ | $-0.01688392216-0.00102025764i$
| $4$ | $-0.01249473743-0.0009507495878i$
| $5$ | $-0.009533459569-0.001537281111i$
| $6$ | $-0.007497937650-0.002167588509i$
| $7$ | $-0.006026376454-0.002674270477i$
$2$ | $3$ | $-0.005785247726+0.0002460429730i$
| $4$ | $-0.003236295992-0.0006934041512i$
| $5$ | $-0.002075023229-0.001296233004i$
| $6$ | $-0.001466727846-0.001581112859i$
| $7$ | $-0.001083345654-0.001682173863i$
$3$ | $4$ | $0.000315183631-0.001771852361i$
| $5$ | $0.000806605055-0.002028059473i$
| $6$ | $0.000954987015-0.001956886197i$
| $7$ | $0.001002044048-0.001756226616i$
$4$ | $5$ | $0.001737194187-0.002338806036i$
| $6$ | $0.001773835947-0.002098671851i$
| $7$ | $0.001741580709-0.001777161054i$
$5$ | $6$ | $0.001993021672-0.001958873499i$
| $7$ | $0.001947091748-0.001620453116i$
$6$ | $7$ | $0.001959730533-0.001367519282i$

###  IV.7 Series expansion method

As an application of our perturbative framework based on a method other than
the Bender-Wu approach, we study the series expansion method known as Leaver’s
method leaver1985 ; konoplya2011 . We consider the system with the
parameterized QNM potential in Eqs (111)-(112) with a single correction term

$\displaystyle\delta
V=\alpha\frac{f}{r_{H}^{2}}\left(\frac{r_{H}}{r}\right)^{j},$ (90)

where $f$ is given by $f=1-r_{H}/r$. We assume the following series expansion
of the wave function as

$\displaystyle\Phi=e^{i\omega r_{*}}\sum_{k=0}^{\infty}a_{k}f^{k+n},$ (91)

where the characteristic exponent $n$ is given by

$\displaystyle n=-2ir_{H}\omega,$ (92)

so that the QNM boundary condition at $r=r_{H}$ is satisfied. After some
calculations, we obtain recursion relations for $a_{k}$

$\displaystyle
A_{k}a_{k-1}+B_{k}a_{k}+C_{k}a_{k+1}+\alpha\sum_{m=0}^{j-2}D_{m}a_{k-m}=0,$
(93)

where coefficients $A_{k},B_{k},C_{k}$ and $D_{m}$ are given by

$\displaystyle A_{k}$ $\displaystyle=(k-2-2ir_{H}\omega)(k+2-2ir_{H}\omega),$
(94) $\displaystyle B_{k}$
$\displaystyle=3-2k(1+k)-\ell(\ell+1)+4ir_{H}\omega(1+2k)+8r_{H}^{2}\omega^{2},$
(95) $\displaystyle C_{k}$ $\displaystyle=(1+k)(1+k-2ir_{H}\omega),$ (96)
$\displaystyle D_{m}$ $\displaystyle=\frac{(-1)^{m+1}(j-2)!}{m!(j-2-m)!}.$
(97)

The coefficients $a_{k}$ with large $k$ take exponentially small value only
for the wave function with the appropriate QNM boundary condition at
$r\to\infty$. Thus, we can calculate the approximate QNM frequency by setting

$\displaystyle a_{k_{\rm max}}=0,$ (98)

with a large integer $k_{\rm max}$. However, directly solving Eq. (98)
numerically is very difficult, and then we usually use Leaver’s continued
fraction method leaver1985 ; konoplya2011  whose basic equation is
mathematically same as Eq. (98). In this section, we study this problem based
on our perturbative approach.

Expanding the coefficients $a_{k}$ and the QNM frequency $\omega$ as

$\displaystyle a_{k}$ $\displaystyle=a_{k}^{(0)}+\alpha
a_{k}^{(1)}+\alpha^{2}a_{k}^{(2)}+\cdots,$ (99) $\displaystyle\omega$
$\displaystyle=\omega_{0}+\alpha\omega_{1}+\alpha^{2}\omega_{2}+\cdots,$ (100)

the coefficients $A_{k},B_{k},C_{k}$ become

$\displaystyle A_{k}$
$\displaystyle=A_{k}^{(0)}+\alpha\omega_{1}A_{k}^{(1)}+\alpha^{2}\omega_{1}^{2}A_{k}^{(2,0)}+\alpha^{2}\omega_{2}A_{k}^{(0,1)}+\cdots,$
(101) $\displaystyle B_{k}$
$\displaystyle=B_{k}^{(0)}+\alpha\omega_{1}B_{k}^{(1)}+\alpha^{2}\omega_{1}^{2}B_{k}^{(2,0)}+\alpha^{2}\omega_{2}B_{k}^{(0,1)}+\cdots,$
(102) $\displaystyle C_{k}$
$\displaystyle=A_{k}^{(0)}+\alpha\omega_{1}C_{k}^{(1)}+\alpha^{2}\omega_{1}^{2}C_{k}^{(2,0)}+\alpha^{2}\omega_{2}C_{k}^{(0,1)}+\cdots,$
(103)

where the coefficients in RHS depend only on $\omega_{0}$. The recursion
relations in Eq. (93) at each order become

$\displaystyle{\cal O}(\alpha^{0}):\quad$ $\displaystyle
A_{k}^{(0)}a_{k-1}^{(0)}+B_{k}^{(0)}a_{k}^{(0)}+C_{k}^{(0)}a_{k+1}^{(0)}=0,$
(104) $\displaystyle{\cal O}(\alpha^{1}):\quad$ $\displaystyle
A_{k}^{(0)}a_{k-1}^{(1)}+B_{k}^{(0)}a_{k}^{(1)}+C_{k}^{(0)}a_{k+1}^{(1)}$
$\displaystyle+\omega_{1}\Big{[}A_{k}^{(1)}a_{k-1}^{(0)}+B_{k}^{(1)}a_{k}^{(0)}+C_{k}^{(1)}a_{k+1}^{(0)}\Big{]}+\sum_{m=0}^{j-2}D_{m}a_{k-m}^{(0)}=0,$
(105) $\displaystyle{\cal O}(\alpha^{2}):\quad$ $\displaystyle
A_{k}^{(0)}a_{k-1}^{(2)}+B_{k}^{(0)}a_{k}^{(2)}+C_{k}^{(0)}a_{k+1}^{(2)}$
$\displaystyle+\omega_{1}\Big{[}A_{k}^{(1)}a_{k-1}^{(1)}+B_{k}^{(1)}a_{k}^{(1)}+C_{k}^{(1)}a_{k+1}^{(1)}\Big{]}$
$\displaystyle+\omega_{1}^{2}\Big{[}A_{k}^{(2,0)}a_{k-1}^{(0)}+B_{k}^{(2,0)}a_{k}^{(0)}+C_{k}^{(2,0)}a_{k+1}^{(0)}\Big{]}$
$\displaystyle+\omega_{2}\Big{[}A_{k}^{(0,2)}a_{k-1}^{(0)}+B_{k}^{(0,2)}a_{k}^{(0)}+C_{k}^{(0,2)}a_{k+1}^{(0)}\Big{]}+\sum_{m=0}^{j-2}D_{m}a_{k-m}^{(1)}=0.$
(106)

We note that these equations correspond to the perturbative equations in Eqs
(10) and (12).

First, at ${\cal O}(\alpha^{0})$, we obtain $\omega_{0}$ using Leaver’s
continued fraction method by setting a large integer $k_{\rm max}$. Next, at
${\cal O}(\alpha^{1})$, we solve the equation

$\displaystyle a_{k_{\rm max}}^{(1)}=0,$ (107)

directly with respect to $\omega_{1}$. For this purpose, we rewrite $a_{k_{\rm
max}}^{(1)}$ as a function of $\omega_{0},\omega_{1},a_{0}^{(0)},a_{0}^{(1)}$
by using Eqs. (104)-(105) recursively, then $a_{k_{\rm max}}^{(1)}$ depends on
$\omega_{1}$ linearly. This implies that we obtain a unique $\omega_{1}$ if we
fix the value of $\omega_{0}$. In a similar way, we can solve the equation

$\displaystyle a_{k_{\rm max}}^{(2)}=0,$ (108)

directly with respect to $\omega_{2}$. In the calculation, we can set
$a_{0}^{(0)}=a_{0}^{(1)}=a_{0}^{(2)}=1$ without loss of generality. We have
confirmed that this method can reproduce a consistent result with Table 5. We
finally note that we do not need to perform the Gaussian elimination to obtain
the three term recursion relations at ${\cal O}(\alpha^{1})$ and higher order
analysis unlike usual Leaver’s continued fraction method leaver1985 ;
konoplya2011 , and this is also one of the advantage of our perturbative
approach.

##  V Outlook

In this paper, we proposed a systematic way to compute high-order perturbative
corrections to black hole quasinormal mode frequencies with continuous
deformation parameters. Our method is widely applicable to many situations,
and allows to compute the high-order corrections very accurately. We showed
various explicit examples. In particular, for the Reissner-Nordström black
holes, we can expand the quasinormal mode frequency not only around the
chargeless limit but also around the extremal limit.

There are several future directions. It is interesting to consider the near
extremal expansion of the Kerr black holes. It was argued that the QNM
frequencies in the extremal Kerr geometry have an interesting behavior in
Sasaki:1989ca . It is also interesting to develop the perturbative expansion
of rotating black holes in modified gravity theories Srivastava:2021imr ;
Pierini:2021jxd ; Pierini:2022eim ; Wagle:2021tam ; Li:2022pcy ; Cano:2023tmv
; Cano:2020cao ; Cano:2023jbk ; Cano:2021myl . In this case, the full analytic
solution with the general rotating parameter has not yet known. We inevitably
have to restrict ourselves to the perturbative treatment in terms of the
rotating parameter. We would like to extend our framework to coupled master
equations. Typically, the master equations in general relativity are
decoupled, but in modified gravity theories, they are sometimes coupled
Molina:2010fb ; Sarbach:2001mc ; Cardoso:2018ptl ; McManus:2019ulj ;
Nomura:2021efi ; Cano:2021myl ; Hui:2022vov . Therefore if we consider
perturbative expansions of modified parameters, it is desirable to generalize
our formalism to such a situation.

###### Acknowledgements.

This research is supported by JSPS KAKENHI Grant Nos. JP22K03641 (YH) and
JP22K03626 (MK).

##  Appendix A Recursion relations among coefficients in parameterized QNM
approach

When the master equation is given in a series expansion of a small parameter,
there is an ambiguity of the effective potential due to the choise of the
master variable. In this appendix, we first give a general discussion of the
ambiguity of effective potential by extending the result in kimura2020 . This
ambiguity leads to recursion relations among coefficients in the parameterized
QNM approach.

###  A.1 Parameterized QNM approach

We consider the case with $f=f_{0}=1-r_{H}/r$, and the master equation is
given by

$\displaystyle
f\frac{d}{dr}\left(f\frac{d\Phi}{dr}\right)+(\omega^{2}-V)\Phi=0,$ (109)

with

$\displaystyle V$ $\displaystyle=V_{0}+\delta V,$ (110) $\displaystyle\delta
V$
$\displaystyle=\frac{f}{r_{H}^{2}}\sum_{j=0}^{\infty}\alpha_{j}\left(\frac{r_{H}}{r}\right)^{j},$
(111)

where $V_{0}$ is the effective potential for non-perturbative case and
$\alpha_{j}$ denote the small parameters which can be written as series of a
single parameter $\alpha$

$\displaystyle\alpha_{j}=\sum_{i=1}^{\infty}\alpha^{i}A^{(i)}_{j}.$ (112)

We note that many systems can be written in this form of the master equation
cardoso2019 ; mcmanus2019 ; hatsuda2020b . The QNM frequency behaves

$\displaystyle\omega=\omega_{0}+\sum_{j=0}^{\infty}\alpha_{j}e_{j}+\sum_{j,k=0}^{\infty}\alpha_{j}\alpha_{k}e_{j,k}+\cdots.$
(113)

where $e_{j},e_{j,k},\cdots$ are model independent coefficients in
parameterized QNM approach.

When $V_{0}$ is the Regge-Wheeler potential for the odd parity gravitational
perturbation, the coefficients are related to the coefficients appearing in
subsection IV.6 as

$\displaystyle e_{j}=e_{j}^{(1)},\qquad e_{j,j}=e_{j}^{(2)},$ (114)
$\displaystyle e_{j,k}=e_{k,j}=\frac{e_{j,k}^{(1,1)}}{2}\quad(j<k),$ (115)

where numerical values of $e_{j}^{(1)},e_{j}^{(2)},e_{j,k}^{(1,1)}$ can be
seen in Tables 5 and 6.

###  A.2 Ambiguity of effective potential

In this subsection, we use the coordinate $x$ defined by $dx/dr=1/f$. The
master equation Eq. (109) in this coordinate becomes

$\displaystyle\frac{d^{2}\Phi}{dx^{2}}+(\omega^{2}-V)\Phi=0.$ (116)

We introduce a new variable $\Psi$ as101010 Note that signatures of $X$ and
$Y$ are opposite from kimura2020 .

$\displaystyle\Psi=\left(1+X\right)\Phi+Y\frac{d\Phi}{dx},$ (117)

where $X$ and $Y$ are ${\cal O}(\alpha)$ functions of $x$. If $X$ and $Y$
satisfy the relation

$\displaystyle-Y^{2}\frac{dV}{dx}+Y\left(2(\omega^{2}-V)\frac{dY}{dx}-\frac{d^{2}X}{dx^{2}}\right)+(1+X)\left(2\delta\frac{dX}{dx}+\frac{d^{2}Y}{dx^{2}}\right)=0,$
(118)

$\Psi$ satisfies an equation

$\displaystyle\frac{d^{2}\Psi}{dx^{2}}+(\omega^{2}-V-\delta W)\Psi=0,$ (119)

where $\delta W$ is given by111111 $\delta W$ also can be written in the form
$\delta W=(2dX/dx+d^{2}Y/dx^{2})/Y$.

$\displaystyle\delta W$
$\displaystyle=\frac{1}{1+X}\left(Y\frac{dV}{dx}-2(\omega^{2}-V)\frac{dY}{dx}+\frac{d^{2}X}{dx^{2}}\right),$
(120)

and this denotes the ambiguity of effective potential. We can regard that the
effective potential changes

$\displaystyle V\to V+\delta W$ (121)

due to the change of the master variable, and then the small parameters
$\alpha_{j}$ in Eq (111) are also changed. Eq. (118) can be integrated as

$\displaystyle
2C+Y\left((V-\omega^{2})Y+\frac{dX}{dx}\right)-\frac{dY}{dx}-2X-X\left(X+\frac{dY}{dx}\right)=0,$
(122)

where $C$ is the constant of integration. If we expand

$\displaystyle X$ $\displaystyle=\sum_{i=1}^{\infty}\alpha^{i}X_{i},$ (123)
$\displaystyle Y$ $\displaystyle=\sum_{i=1}^{\infty}\alpha^{i}Y_{i},$ (124)
$\displaystyle V$ $\displaystyle=V_{0}+\delta
V=\sum_{i=0}^{\infty}\alpha^{i}V_{i},$ (125) $\displaystyle\omega^{2}$
$\displaystyle=\sum_{i=0}^{\infty}\alpha^{i}{\cal E}_{i},$ (126)
$\displaystyle C$ $\displaystyle=\sum_{i=1}^{\infty}\alpha^{i}{\cal C}_{i},$
(127)

Eq. (122) can be solved order by order as

$\displaystyle X_{i}$
$\displaystyle=C_{i}-\frac{1}{2}Y_{i}^{\prime}-\frac{1}{2}\sum_{k=1}^{i-1}\sum_{j=0}^{i-k-1}({\cal
E}_{j}-V_{j})Y_{k}Y_{i-k-j}+\frac{1}{2}\sum_{k=1}^{i-1}\left(Y_{i-k}X_{k}^{\prime}-X_{i-k}Y_{k}^{\prime}-X_{i-k}X_{k}\right).$
(128)

If we also expand $\omega=\sum_{i=0}^{\infty}\alpha^{i}\omega_{i}$, ${\cal
E}_{i}$ is given by

$\displaystyle{\cal E}_{i}=\sum_{j=0}^{i}\omega_{i-j}\omega_{j}.$ (129)

Substituting the result (128) into Eq. (120), we can calculate the deformation
of the effective potential $\delta W$ as the series of $\alpha$

$\displaystyle\delta W$ $\displaystyle=\sum_{i=1}^{\infty}\alpha^{i}W_{i}.$
(130)

From Eq (120), we can write $W_{i}$ as

$\displaystyle
W_{i}=\frac{d^{2}X_{i}}{dx^{2}}+\sum_{j=0}^{i-1}\left(Y_{i-j}\frac{dV_{j}}{dx}-2({\cal
E}_{j}-V_{j})\frac{dY_{i-j}}{dx}\right)-\sum_{j=1}^{i-1}W_{i-j}X_{j}.$ (131)

For lower $i$, the explicit forms are

$\displaystyle X_{1}$ $\displaystyle=C_{1}-\frac{1}{2}\frac{dY_{1}}{dx},$
(132) $\displaystyle X_{2}$ $\displaystyle=C_{2}-\frac{1}{2}({\cal
E}_{0}-V_{0})Y_{1}^{2}+\frac{1}{8}\left(\frac{dY_{1}}{dx}\right)^{2}-\frac{1}{2}\left(C_{1}^{2}+\frac{dY_{2}}{dx}\right)-\frac{1}{4}Y_{1}\frac{d^{2}Y_{1}}{dx^{2}},$
(133)

and

$\displaystyle W_{1}$ $\displaystyle=Y_{1}\frac{dV_{0}}{dx}-2({\cal
E}_{0}-V_{0})\frac{dY_{1}}{dx}-\frac{1}{2}\frac{d^{3}Y_{1}}{dx^{3}},$ (134)
$\displaystyle W_{2}$
$\displaystyle=(Y_{2}-C_{1}Y_{1})\frac{dV_{0}}{dx}-2({\cal
E}_{0}-V_{0})\frac{d(Y_{2}-C_{1}Y_{1})}{dx}-\frac{1}{2}\frac{d^{3}(Y_{2}-C_{1}Y_{1})}{dx^{3}}$
(135)
$\displaystyle+\frac{Y_{1}}{2}\left(2\frac{dV_{1}}{dx}+Y_{1}\frac{d^{2}V_{0}}{dx^{2}}\right)-({\cal
E}_{0}-V_{0})\left[2\left(\frac{dY_{1}}{dx}\right)^{2}+Y_{1}\frac{d^{2}Y_{1}}{dx^{2}}\right]$
(136) $\displaystyle+\frac{dY_{1}}{dx}\left(-2({\cal
E}_{1}-V_{1})+\frac{5Y_{1}}{2}\frac{dV_{0}}{dx}-\frac{1}{2}\frac{d^{3}Y_{1}}{dx^{3}}\right)-\frac{Y_{1}}{4}\frac{d^{4}Y_{1}}{dx^{2}}.$
(137)

We note that $W_{i}$ contains arbitrary functions $Y_{1},Y_{2},\cdots$. If we
set $V_{i}=0$ for $i\geq 1$, the system is just a non-perturbative case whose
effective potential is $V_{0}$. Nevertheless, there is an ambiguity of
effective potential due to the change of the master variable. In this case,
the ambiguity of effective potential does not change the QNM spectrum, and we
can obtain recursion relations among coefficients in parameterized QNM
approach by setting the functions $Y_{i}$ appropriately as shown in the next
subsection.

###  A.3 Recursion relations for odd parity case

####  A.3.1 Recursion relations from the Regge-Wheeler potential

As an example, we consider the odd parity case

$\displaystyle
V=V_{0}=f_{0}\left(\frac{\ell(\ell+1)}{r^{2}}-\frac{3r_{H}}{r^{3}}\right).$
(138)

In this case, ${\cal E}_{1}=0$ because there is no correction term in the
effective potential $V$, i.e., $V_{i}=0$ for $i\geq 1$. Setting121212 From the
degrees of freedom of $Y_{2}$, we can obtain the same relation as the first
order relation among $e_{j}$. Also, $C_{1}$ does not affect the result. Thus,
we can set $Y_{2}=0$ and $C_{1}=0$.

$\displaystyle Y_{1}$
$\displaystyle=y_{j}\left(\frac{r_{H}}{r}\right)^{j}+y_{k}\left(\frac{r_{H}}{r}\right)^{k},$
(139) $\displaystyle Y_{2}$ $\displaystyle=0,$ (140) $\displaystyle C_{1}$
$\displaystyle=0,$ (141)

where $j,k\geq-1$ are integers and $y_{j},y_{k}$ are constants, Eqs
(130)-(137) lead to

$\displaystyle\delta V+\delta W=\alpha
y_{j}f_{0}\left(\frac{r_{H}}{r}\right)^{j}\bigg{[}\frac{2j{\cal
E}_{0}}{r}+\frac{(j+1)(j-2\ell)(j+2\ell+2)}{2r^{3}}$
$\displaystyle-\frac{(2j+3)r_{H}\left(j(j+3)-2\left(\ell^{2}+\ell+3\right)\right)}{2r^{4}}+\frac{(j-2)(j+2)(j+6)r_{H}^{2}}{2r^{5}}\bigg{]}+(j\leftrightarrow
k)$
$\displaystyle+\alpha^{2}y_{j}^{2}f_{0}\left(\frac{r_{H}}{r}\right)^{2j}\bigg{[}-\frac{j(3j+1){\cal
E}_{0}}{r^{2}}+\frac{j(3j+2)r_{H}{\cal
E}_{0}}{r^{3}}-\frac{3(j+1)^{2}(j-2\ell)(j+2\ell+2)}{4r^{4}}$
$\displaystyle+\frac{(3j+4)r_{H}\left(3j^{3}+12j^{2}-j(8\ell(\ell+1)+1)-2(5\ell(\ell+1)+9)\right)}{4r^{5}}$
$\displaystyle-\frac{(3j+5)r_{H}^{2}\left(3j^{3}+15j^{2}-j(4\ell(\ell+1)+7)-6\left(\ell^{2}+\ell+7\right)\right)}{4r^{6}}$
$\displaystyle+\frac{3(j-2)(j+2)^{2}(j+6)r_{H}^{3}}{4r^{7}}\bigg{]}+(j\leftrightarrow
k)$
$\displaystyle+\alpha^{2}y_{j}y_{k}f_{0}\left(\frac{r_{H}}{r}\right)^{j+k}\bigg{[}-\frac{{\cal
E}_{0}(j^{2}+4jk+j+k^{2}+k)}{r^{2}}+\frac{{\cal
E}_{0}r_{H}(j^{2}+j(4k+2)+k(k+2))}{r^{3}}$
$\displaystyle+\frac{1}{4r^{4}}\Big{(}j^{2}(-6k+4\ell(\ell+1)-11)-2j^{3}(k+3)-j^{4}+4(k(k+6)+6)\ell-k(k+1)(k+2)(k+3)$
$\displaystyle+2j(4(2k+3)\ell^{2}+4(2k+3)\ell-k(k(k+3)+4)-3)+4(k(k+6)+6)\ell^{2}\Big{)}$
$\displaystyle+\frac{r_{H}}{4r^{5}}\Big{(}j^{2}(24k-8\ell(\ell+1)+47)+6j^{3}(k+4)+3j^{4}+k^{2}(47-8\ell(\ell+1))+3k^{4}+24k^{3}$
$\displaystyle-2k(31\ell(\ell+1)+29)-16(5\ell(\ell+1)+9)+j(6k^{3}+24k^{2}-4k(8\ell(\ell+1)+1)-62\ell(\ell+1)-58)\Big{)}$
$\displaystyle+\frac{r_{H}^{2}}{4r^{6}}\Big{(}j^{2}(4(\ell^{2}+\ell-17)-30k)-6j^{3}(k+5)-3j^{4}+4k^{2}(\ell^{2}+\ell-17)-3k^{4}-30k^{3}$
$\displaystyle+k(38\ell(\ell+1)+161)+60(\ell^{2}+\ell+7)+j(-6k^{3}-30k^{2}+4k(4\ell(\ell+1)+7)+38\ell(\ell+1)+161)\Big{)}$
$\displaystyle+\frac{r_{H}^{3}(2j^{3}(k+6)+4j^{2}(3k+8)+j^{4}+2j(k+6)(k^{2}-8)+k^{2}(k+4)(k+8)-96(k+3))}{4r^{7}}\bigg{]}$
$\displaystyle+{\cal O}(\alpha^{3}),$ (142)

where we used the relation $d/dx=fd/dr$. From this result, we can read
$\alpha_{i}$ for $\delta V+\delta W$. We decompose the coefficients
$\alpha_{i}=A^{(1)}_{i}\alpha+A^{(2)}_{i}\alpha^{2}+{\cal O}(\alpha^{3})$ in
Eq. (112) as

$\displaystyle A^{(1)}_{i}$
$\displaystyle=y_{j}\partial_{y_{j}}A^{(1)}_{i}+y_{k}\partial_{y_{k}}A^{(1)}_{i}$
(143) $\displaystyle A^{(2)}_{i}$
$\displaystyle=\frac{y_{j}^{2}}{2}\partial^{2}_{y_{j}}A^{(2)}_{i}+\frac{y_{k}^{2}}{2}\partial^{2}_{y_{k}}A^{(2)}_{i}+y_{j}y_{k}\partial_{y_{j}}\partial_{y_{k}}A^{(2)}_{i}.$
(144)

Introducing $\partial_{y_{j}}A^{(1)}_{i}=r_{H}^{-1}B^{(1)}_{i}$,
$\partial_{y_{j}}\partial_{y_{k}}A^{(2)}_{i}=r_{H}^{-2}B^{(2)}_{i}$, then one
can see that the relations

$\displaystyle\partial_{y_{k}}A^{(1)}_{i}$
$\displaystyle=r_{H}^{-1}B^{(1)}_{i}|_{j\to k},$ (145)
$\displaystyle\partial^{2}_{y_{j}}A^{(2)}_{i}$
$\displaystyle=\frac{r_{H}^{-2}}{2}B^{(2)}_{i}|_{k\to j},$ (146)
$\displaystyle\partial^{2}_{y_{k}}A^{(2)}_{i}$
$\displaystyle=\frac{r_{H}^{-2}}{2}B^{(2)}_{i}|_{j\to k}$ (147)

hold from the expression of Eq. (142). The explicit forms of $B^{(1)}_{i}$ and
$B^{(2)}_{i}$ become

$\displaystyle B^{(1)}_{j+1}$ $\displaystyle=2jr_{H}^{2}{\cal E}_{0},$ (148)
$\displaystyle B^{(1)}_{j+3}$
$\displaystyle=\frac{1}{2}(j+1)(j-2\ell)(j+2\ell+2),$ (149) $\displaystyle
B^{(1)}_{j+4}$
$\displaystyle=-\frac{1}{2}(2j+3)\left(j(j+3)-2\left(\ell^{2}+\ell+3\right)\right),$
(150) $\displaystyle B^{(1)}_{j+5}$
$\displaystyle=\frac{1}{2}(j-2)(j+2)(j+6),$ (151)

and

$\displaystyle B^{(2)}_{j+k+2}$
$\displaystyle=-(j^{2}+4jk+j+k^{2}+k)r_{H}^{2}{\cal E}_{0},$ (152)
$\displaystyle B^{(2)}_{j+k+3}$
$\displaystyle=(j^{2}+j(4k+2)+k(k+2))r_{H}^{2}{\cal E}_{0},$ (153)
$\displaystyle B^{(2)}_{j+k+4}$
$\displaystyle=\frac{1}{4}\Big{(}j^{2}(-6k+4\ell(\ell+1)-11)-2j^{3}(k+3)-j^{4}+4(k(k+6)+6)\ell$
$\displaystyle-k(k+1)(k+2)(k+3)+2j(4(2k+3)\ell^{2}+4(2k+3)\ell$
$\displaystyle-k(k(k+3)+4)-3)+4(k(k+6)+6)\ell^{2}\Big{)},$ (154)
$\displaystyle B^{(2)}_{j+k+5}$
$\displaystyle=\frac{1}{4}\Big{(}j^{2}(24k-8\ell(\ell+1)+47)+6j^{3}(k+4)+3j^{4}$
$\displaystyle+k^{2}(47-8\ell(\ell+1))+3k^{4}+24k^{3}-2k(31\ell(\ell+1)+29)-16(5\ell(\ell+1)+9)$
$\displaystyle+j(6k^{3}+24k^{2}-4k(8\ell(\ell+1)+1)-62\ell(\ell+1)-58)\Big{)},$
(155) $\displaystyle B^{(2)}_{j+k+6}$
$\displaystyle=\frac{1}{4}\Big{(}j^{2}(4(\ell^{2}+\ell-17)-30k)-6j^{3}(k+5)-3j^{4}+4k^{2}(\ell^{2}+\ell-17)$
$\displaystyle-3k^{4}-30k^{3}+k(38\ell(\ell+1)+161)+60(\ell^{2}+\ell+7)$
$\displaystyle+j(-6k^{3}-30k^{2}+4k(4\ell(\ell+1)+7)+38\ell(\ell+1)+161)\Big{)},$
(156) $\displaystyle B^{(2)}_{j+k+7}$
$\displaystyle=\frac{1}{4}\Big{(}2j^{3}(k+6)+4j^{2}(3k+8)+j^{4}+2j(k+6)(k^{2}-8)$
$\displaystyle+k^{2}(k+4)(k+8)-96(k+3)\Big{)}.$ (157)

Because ${\cal E}_{1}=0$ and then $\omega=\omega_{0}$, from Eq. (113), we
obtain a relation

$\displaystyle\sum_{j=0}^{\infty}\alpha_{j}e_{j}+\sum_{j,k=0}^{\infty}\alpha_{j}\alpha_{k}e_{j,k}=0.$
(158)

From ${\cal O}(\alpha)$ and ${\cal O}(\alpha^{2})$ terms in Eq. (158), we
obtain independent recursion relations among $e_{j}$ and $e_{j,k}$

$\displaystyle 0$ $\displaystyle=\sum_{a=1}^{5}B_{j+a}^{(1)}e_{j+a}$
$\displaystyle=B_{j+1}^{(1)}e_{j+1}+B_{j+3}^{(1)}e_{j+3}+B_{j+4}^{(1)}e_{j+4}+B_{j+5}^{(1)}e_{j+5},$
(159)

and

$\displaystyle 0$
$\displaystyle=\sum_{a,b=1}^{5}B_{j+a}^{(1)}B_{k+b}^{(1)}e_{j+a,k+b}+\frac{1}{2}\sum_{a=2}^{7}B_{j+k+a}^{(2)}e_{j+k+a}$
$\displaystyle=B_{j+1}^{(1)}B_{k+1}^{(1)}e_{j+1,k+1}+B_{j+1}^{(1)}B_{k+3}^{(1)}e_{j+1,k+3}+B_{j+1}^{(1)}B_{k+4}^{(1)}e_{j+1,k+4}+B_{j+1}^{(1)}B_{k+5}^{(1)}e_{j+1,k+5}$
$\displaystyle+B_{j+3}^{(1)}B_{k+1}^{(1)}e_{j+3,k+1}+B_{j+3}^{(1)}B_{k+3}^{(1)}e_{j+3,k+3}+B_{j+3}^{(1)}B_{k+4}^{(1)}e_{j+3,k+4}+B_{j+3}^{(1)}B_{k+5}^{(1)}e_{j+3,k+5}$
$\displaystyle+B_{j+4}^{(1)}B_{k+1}^{(1)}e_{j+4,k+1}+B_{j+4}^{(1)}B_{k+3}^{(1)}e_{j+4,k+3}+B_{j+4}^{(1)}B_{k+4}^{(1)}e_{j+4,k+4}+B_{j+4}^{(1)}B_{k+5}^{(1)}e_{j+4,k+5}$
$\displaystyle+B_{j+5}^{(1)}B_{k+1}^{(1)}e_{j+5,k+1}+B_{j+5}^{(1)}B_{k+3}^{(1)}e_{j+5,k+3}+B_{j+5}^{(1)}B_{k+4}^{(1)}e_{j+5,k+4}+B_{j+5}^{(1)}B_{k+5}^{(1)}e_{j+5,k+5}$
$\displaystyle+\frac{1}{2}\Big{[}B_{j+k+2}^{(2)}e_{j+k+2}+B_{j+k+3}^{(2)}e_{j+k+3}+B_{j+k+4}^{(2)}e_{j+k+4}$
$\displaystyle+B_{j+k+5}^{(2)}e_{j+k+5}+B_{j+k+6}^{(2)}e_{j+k+6}+B_{j+k+7}^{(2)}e_{j+k+7}\Big{]}.$
(160)

We note again that $e_{j}=e_{j}^{(1)},e_{j,j}=e_{j}^{(2)}$ and
$e_{j,k}=e_{j,k}^{(1,1)}/2$ for $j\neq k$, where numerical values of
$e_{j}^{(1)},e_{j}^{(2)},e_{j,k}^{(1,1)}$ can be seen in Tables 5 and 6. Using
the first order recursion relation in Eq. (159), $e_{j}$ with higher $j$ can
be written only from those with a few lower $j$, i.e., $e_{0},e_{2}$ and
$e_{7}$ kimura2020 . However, this is not the case for the second order
recursion relation in Eq. (160). In fact, to calculate $e_{j,k}$ with higher
$j,k$ using Eq. (160), we need the values of
$e_{j,0},e_{j,2},e_{j,7},e_{k,0},e_{k,2},e_{k,7}$. To improve this point, we
study the case with the potential which contains first order correction terms
in the next subsection.

####  A.3.2 Improved recursion relation for $e_{j,k}$

We consider the Regge-Wheeler potential with first order correction terms

$\displaystyle V$ $\displaystyle=V_{0}+\delta V$
$\displaystyle=f_{0}\left(\frac{\ell(\ell+1)}{r^{2}}-\frac{3r_{H}}{r^{3}}\right)+\frac{\alpha
f_{0}}{r_{H}^{2}}\left[v_{j}\left(\frac{r_{H}}{r}\right)^{j+5}+v_{k}\left(\frac{r_{H}}{r}\right)^{k+5}\right],$
(161)

where $j,k\geq-1$ are integers and $v_{j},v_{k}$ are constants. We also assume
that $j\neq 2$ and $k\neq 2$. In this case, the QNM frequency behaves

$\displaystyle\omega=\omega_{0}+\alpha\omega_{1}+\alpha^{2}\omega_{2},$ (162)

with

$\displaystyle\omega_{1}$ $\displaystyle=v_{j}e_{j+5}+v_{k}e_{k+5},$ (163)
$\displaystyle\omega_{2}$
$\displaystyle=v_{j}^{2}e_{j+5,j+5}+2v_{j}v_{k}e_{j+5,k+5}+v_{k}^{2}e_{k+5,k+5}.$
(164)

${\cal E}_{1}=2\omega_{0}\omega_{1}$ becomes

$\displaystyle{\cal E}_{1}=2v_{j}e_{j+5}\omega_{0}+2v_{k}e_{k+5}\omega_{0}.$
(165)

For this potential $V=V_{0}+\delta V$, we set

$\displaystyle Y_{1}$
$\displaystyle=y_{j}\left(\frac{r_{H}}{r}\right)^{j}+y_{k}\left(\frac{r_{H}}{r}\right)^{k},$
(166) $\displaystyle Y_{2}$ $\displaystyle=0,$ (167) $\displaystyle C_{1}$
$\displaystyle=0,$ (168)

with

$\displaystyle y_{j}$ $\displaystyle=-\frac{2v_{j}r_{H}}{(j-2)(j+2)(j+6)},$
(169) $\displaystyle y_{k}$
$\displaystyle=-\frac{2v_{k}r_{H}}{(k-2)(k+2)(k+6)}.$ (170)

Then, Eqs (130)-(137) lead to

$\displaystyle\delta V+\delta W=\alpha
y_{j}f_{0}\left(\frac{r_{H}}{r}\right)^{j}\bigg{[}\frac{2j{\cal
E}_{0}}{r}+\frac{(j+1)(j-2\ell)(j+2\ell+2)}{2r^{3}}$
$\displaystyle-\frac{(2j+3)r_{H}\left(j(j+3)-2\left(\ell^{2}+\ell+3\right)\right)}{2r^{4}}\bigg{]}+(j\leftrightarrow
k)+\alpha^{2}{\cal
E}_{1}\left[y_{j}f_{0}\left(\frac{r_{H}}{r}\right)^{j}\frac{2j}{r}+y_{k}f_{0}\left(\frac{r_{H}}{r}\right)^{k}\frac{2k}{r}\right]$
$\displaystyle+\alpha^{2}y_{j}^{2}f_{0}\left(\frac{r_{H}}{r}\right)^{2j}\bigg{[}-\frac{j(3j+1){\cal
E}_{0}}{r^{2}}+\frac{j(3j+2)r_{H}{\cal
E}_{0}}{r^{3}}-\frac{3(j+1)^{2}(j-2\ell)(j+2\ell+2)}{4r^{4}}$
$\displaystyle+\frac{(3j+4)r_{H}\left(3j^{3}+12j^{2}-j(8\ell(\ell+1)+1)-2(5\ell(\ell+1)+9)\right)}{4r^{5}}$
$\displaystyle-\frac{(3j+5)r_{H}^{2}\left(j^{3}+3j^{2}+j(1-4\ell(\ell+1))-6(\ell^{2}+\ell-1)\right)}{4r^{6}}$
$\displaystyle-\frac{3(j-2)(j+2)^{2}(j+6)r_{H}^{3}}{4r^{7}}\bigg{]}+(j\leftrightarrow
k)$
$\displaystyle+\alpha^{2}y_{j}y_{k}f_{0}\left(\frac{r_{H}}{r}\right)^{j+k}\bigg{[}-\frac{{\cal
E}_{0}(j^{2}+4jk+j+k^{2}+k)}{r^{2}}+\frac{{\cal
E}_{0}r_{H}(j^{2}+j(4k+2)+k(k+2))}{r^{3}}$
$\displaystyle+\frac{1}{4r^{4}}\Big{(}j^{2}(-6k+4\ell(\ell+1)-11)-2j^{3}(k+3)-j^{4}+4(k(k+6)+6)\ell-k(k+1)(k+2)(k+3)$
$\displaystyle+2j(4(2k+3)\ell^{2}+4(2k+3)\ell-k(k(k+3)+4)-3)+4(k(k+6)+6)\ell^{2}\Big{)}$
$\displaystyle+\frac{r_{H}}{4r^{5}}\Big{(}j^{2}(24k-8\ell(\ell+1)+47)+6j^{3}(k+4)+3j^{4}+k^{2}(47-8\ell(\ell+1))+3k^{4}+24k^{3}$
$\displaystyle-2k(31\ell(\ell+1)+29)-16(5\ell(\ell+1)+9)+j(6k^{3}+24k^{2}-4k(8\ell(\ell+1)+1)-62\ell(\ell+1)-58)\Big{)}$
$\displaystyle-\frac{r_{H}^{2}}{4r^{6}}\Big{(}2j^{2}(3k-2(\ell^{2}+\ell-4))+2j^{3}(k+4)+j^{4}+j(2k^{3}+6k^{2}-16k\ell(\ell+1)+4k-38\ell(\ell+1)+23)$
$\displaystyle-4k^{2}(\ell^{2}+\ell-4)+k^{4}+8k^{3}+k(23-38\ell(\ell+1))-60(\ell^{2}+\ell-1)\Big{)}$
$\displaystyle-\frac{r_{H}^{3}(2j^{3}(k+6)+4j^{2}(3k+8)+j^{4}+2j(k+6)(k^{2}-8)+k^{2}(k+4)(k+8)-96(k+3))}{4r^{7}}\bigg{]}$
$\displaystyle+{\cal O}(\alpha^{3}).$ (171)

We note that the above potential at ${\cal O}(\alpha)$ does not have terms
with $(r_{H}/r)^{j+5}$ and $(r_{H}/r)^{k+5}$ unlike Eq. (142). Similar to the
discussion in the previous subsection, we can read the coefficients
$B_{i}^{(1)}$ and $B_{i}^{(2)}$ as

$\displaystyle B^{(1)}_{j+1}$ $\displaystyle=2jr_{H}^{2}{\cal E}_{0},$ (172)
$\displaystyle B^{(1)}_{j+3}$
$\displaystyle=\frac{1}{2}(j+1)(j-2\ell)(j+2\ell+2),$ (173) $\displaystyle
B^{(1)}_{j+4}$
$\displaystyle=-\frac{1}{2}(2j+3)\left(j(j+3)-2\left(\ell^{2}+\ell+3\right)\right),$
(174)

and

$\displaystyle B^{(2)}_{j+1}$
$\displaystyle=-2j(k-2)(k+2)(k+6)r_{H}^{2}\omega_{0}e_{k+5},$ (175)
$\displaystyle B^{(2)}_{k+1}$
$\displaystyle=-2k(j-2)(j+2)(j+6)r_{H}^{2}\omega_{0}e_{j+5},$ (176)
$\displaystyle B^{(2)}_{j+k+2}$
$\displaystyle=-(j^{2}+4jk+j+k^{2}+k)r_{H}^{2}{\cal E}_{0},$ (177)
$\displaystyle B^{(2)}_{j+k+3}$
$\displaystyle=(j^{2}+j(4k+2)+k(k+2))r_{H}^{2}{\cal E}_{0},$ (178)
$\displaystyle B^{(2)}_{j+k+4}$
$\displaystyle=\frac{1}{4}\Big{(}j^{2}(-6k+4\ell(\ell+1)-11)-2j^{3}(k+3)-j^{4}+4(k(k+6)+6)\ell$
$\displaystyle-k(k+1)(k+2)(k+3)+2j(4(2k+3)\ell^{2}+4(2k+3)\ell$
$\displaystyle-k(k(k+3)+4)-3)+4(k(k+6)+6)\ell^{2}\Big{)},$ (179)
$\displaystyle B^{(2)}_{j+k+5}$
$\displaystyle=\frac{1}{4}\Big{(}j^{2}(24k-8\ell(\ell+1)+47)+6j^{3}(k+4)+3j^{4}$
$\displaystyle+k^{2}(47-8\ell(\ell+1))+3k^{4}+24k^{3}-2k(31\ell(\ell+1)+29)-16(5\ell(\ell+1)+9)$
$\displaystyle+j(6k^{3}+24k^{2}-4k(8\ell(\ell+1)+1)-62\ell(\ell+1)-58)\Big{)},$
(180) $\displaystyle B^{(2)}_{j+k+6}$
$\displaystyle=-\frac{1}{4}\Big{(}2j^{2}(3k-2(\ell^{2}+\ell-4))+2j^{3}(k+4)+j^{4}$
$\displaystyle+j(2k^{3}+6k^{2}-16k\ell(\ell+1)+4k-38\ell(\ell+1)+23)$
$\displaystyle-4k^{2}(\ell^{2}+\ell-4)+k^{4}+8k^{3}+k(23-38\ell(\ell+1))-60(\ell^{2}+\ell-1)\Big{)},$
(181) $\displaystyle B^{(2)}_{j+k+7}$
$\displaystyle=-\frac{1}{4}\Big{(}(2j^{3}(k+6)+4j^{2}(3k+8)+j^{4}+2j(k+6)(k^{2}-8)$
$\displaystyle+k^{2}(k+4)(k+8)-96(k+3))\Big{)}.$ (182)

Then, the QNM frequency can be calculated from Eq. (112), and it should be
same as Eq. (162) with Eqs (163) and (164). From this condition, we obtain
independent recursion relations at ${\cal O}(\alpha^{2})$ as

$\displaystyle\frac{1}{2}(j-2)(j+2)(j+6)(k-2)(k+2)(k+6)e_{j+5,k+5}$
$\displaystyle=2\sum_{a,b=1}^{4}B_{j+a}^{(1)}B_{k+b}^{(1)}e_{j+a,k+b}+\sum_{a=2}^{7}B_{j+k+a}^{(2)}e_{j+k+a}+B_{j+1}^{(2)}e_{j+1}+B_{k+1}^{(2)}e_{k+1}.$
(183)

We note again that $j,k\geq-1$ and $j\neq 2,k\neq 2$ in the above equation.

In fact, we can obtain further independent recursion relations for $e_{j,k}$.
We consider the potential in Eq. (161) with $j\geq-1,j\neq 2$ and $k\geq-5$.
Setting

$\displaystyle Y_{1}$ $\displaystyle=y_{j}\left(\frac{r_{H}}{r}\right)^{j},$
(184) $\displaystyle Y_{2}$ $\displaystyle=0,$ (185) $\displaystyle C_{1}$
$\displaystyle=0,$ (186)

with

$\displaystyle y_{j}$ $\displaystyle=-\frac{2v_{j}r_{H}}{(j-2)(j+2)(j+6)},$
(187)

we can calculate $\delta V+\delta W$ from Eqs (130)-(137), and derive the
recursion relations similar to the above discussion. Here, we only show the
result:

$\displaystyle 0$
$\displaystyle=(j-2)(j+2)(j+6)e_{j+5,k+5}-(2j+3)(j(j+3)-2(\ell^{2}+\ell+3))e_{j+4,k+5}$
$\displaystyle+(j+1)(j-2\ell)(j+2\ell+2)e_{j+3,k+5}+4jr_{H}^{2}\mathcal{E}_{0}e_{j+1,k+5}$
$\displaystyle+4jr_{H}^{2}\omega_{0}e_{j+1}e_{k+5}-(2j+k+5)e_{j+k+6}+(2j+k+6)e_{j+k+7}.$
(188)

Using Eqs. (183) and (188), the second order coefficients $e_{j,k}$ with
higher $j,k$ can be written by those with $j,k\leq 7$ and the first order
coefficients $e_{j}$.131313 Some of coefficients $e_{j,k}$ with $j,k\leq 7$
are not independent. For example, we can choose
$e_{0,0},e_{1,0},e_{1,1},e_{2,0},e_{2,1},e_{2,2},e_{3,0},e_{3,1},e_{3,2},e_{3,3},e_{7,0},e_{7,1},e_{7,2},e_{7,3}$
and $e_{7,7}$ as independent $e_{j,k}$, then, the other $e_{j,k}$ can be
written by these.  We note that we can derive recursion relations for higher
order $\alpha$ from a straightforward extension of the above discussion.

###  A.4 Reduction of the effective potential

Using the ambiguity of effective potential, we can reduce the the effective
potential so that $\delta V$ only has lower order coefficients $\alpha_{j}$.
In kimura2020 , the first order case is discussed, but in fact, the discussion
holds even for the higher order case. For the linear order case, we can reduce
effective potential by using $O(\alpha)$ ambiguity according to kimura2020 .
For the quadratic order case, setting $Y_{1}=0$ and $Y_{2}=y_{j}(r_{H}/r)^{j}$
for the odd parity perturbation, the form of the ambiguity of the effective
potential at $O(\alpha^{2})$ becomes same as the linear order case. Then, from
the same discussion as linear case in kimura2020 , we can reduce the effective
potential at $O(\alpha^{2})$ so that $\delta V$ only has
$\alpha_{0},\alpha_{1},\alpha_{2}$ and $\alpha_{7}$ terms. Repeating this
process to higher order, we can reduce $O(\alpha^{n})$ effective potential.

## References

  * (1) T. Nakamura, K. Oohara and Y. Kojima, General Relativistic Collapse to Black Holes and Gravitational Waves from Black Holes, Prog. Theor. Phys. Suppl. 90 (1987) 1–218.
  * (2) K. D. Kokkotas and B. G. Schmidt, Quasinormal modes of stars and black holes, Living Rev. Rel. 2 (1999) 2, arXiv:gr-qc/9909058.
  * (3) H.-P. Nollert, TOPICAL REVIEW: Quasinormal modes: the characteristic ‘sound’ of black holes and neutron stars, Class. Quant. Grav. 16 (1999) R159–R216.
  * (4) E. Berti, V. Cardoso and A. O. Starinets, Quasinormal modes of black holes and black branes, Class. Quantum Grav. 26 (2009) 163001.
  * (5) R. A. Konoplya and A. Zhidenko, Quasinormal modes of black holes: From astrophysics to string theory, Rev. Mod. Phys. 83 (2011) 793–836.
  * (6) V. Ferrari and L. Gualtieri, Quasi-Normal Modes and Gravitational Wave Astronomy, Gen. Rel. Grav. 40 (2008) 945–970, arXiv:0709.0657 [gr-qc].
  * (7) E. W. Leaver, Spectral decomposition of the perturbation response of the Schwarzschild geometry, Phys. Rev. D 34 (1986) 384–408.
  * (8) N. Andersson, Evolving test fields in a black hole geometry, Phys. Rev. D 55 (1997) 468–479, arXiv:gr-qc/9607064.
  * (9) N. Andersson, Excitation of Schwarzschild black hole quasinormal modes, Phys. Rev. D 51 (1995) 353–363.
  * (10) H.-P. Nollert and B. G. Schmidt, Quasinormal modes of Schwarzschild black holes: Defined and calculated via Laplace transformation, Phys. Rev. D 45 (1992) 2617.
  * (11) E. Berti and V. Cardoso, Quasinormal ringing of Kerr black holes. I. The Excitation factors, Phys. Rev. D 74 (2006) 104020, arXiv:gr-qc/0605118.
  * (12) P. T. Leung, Y. T. Liu, W.-M. Suen, C. Y. Tam and K. Young, Quasinormal Modes of Dirty Black Holes, Phys. Rev. Lett. 78 (1997) 2894–2897.
  * (13) P. T. Leung, Y. T. Liu, W. M. Suen, C. Y. Tam and K. Young, Logarithmic perturbation theory for quasinormal modes, J. Phys. A: Math. Gen. 31 (1998) 3271–3286.
  * (14) P. T. Leung, Y. T. Liu, W. M. Suen, C. Y. Tam and K. Young, Perturbative approach to the quasinormal modes of dirty black holes, Phys. Rev. D 59 (1999) 044034.
  * (15) 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 (2019) 104077.
  * (16) 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 (2019) 044061.
  * (17) Y. Hatsuda, Quasinormal modes of black holes and Borel summation, Phys. Rev. D 101 (2020) 024008.
  * (18) E. W. Leaver, An analytic representation for the quasi-normal modes of Kerr black holes, Proc. R. Soc. Lond. A. 402 (1985) 285–298.
  * (19) S. Chandrasekhar and S. L. Detweiler, The quasi-normal modes of the Schwarzschild black hole, Proc. Roy. Soc. Lond. A 344 (1975) 441–452.
  * (20) A. Jansen, Overdamped modes in Schwarzschild-de Sitter and a Mathematica package for the numerical computation of quasinormal modes, Eur. Phys. J. Plus 132 (2017) 546, arXiv:1709.09178 [gr-qc].
  * (21) C. M. Bender and T. T. Wu, Anharmonic Oscillator, Phys. Rev. 184 (1969) 1231–1260.
  * (22) D. S. Eniceicu and M. Reece, Quasinormal modes of charged fields in Reissner-Nordstrom backgrounds by Borel-Pade summation of Bender-Wu series, Phys. Rev. D 102 (2020) 044015.
  * (23) H.-J. Blome and B. Mashhoon, Quasi-normal oscillations of a schwarzschild black hole, Phys. Lett. A 100 (1984) 231–234.
  * (24) V. Ferrari and B. Mashhoon, Oscillations of a Black Hole, Phys. Rev. Lett. 52 (1984) 1361–1364.
  * (25) V. Ferrari and B. Mashhoon, New approach to the quasinormal modes of a black hole, Phys. Rev. D 30 (1984) 295–304.
  * (26) B. Mashhoon, Quasi-Normal Modes of a Black Hole, Proceeding of the Third Marcel Grossmann Meeting on General Relativity 18 (1983) 599–608.
  * (27) B. F. Schutz and C. M. Will, Black hole normal modes - A semianalytic approach, Astrophys. J. 291 (1985) L33–L36.
  * (28) T. Sulejmanpasic and M. Ünsal, Aspects of perturbation theory in quantum mechanics: The BenderWu Mathematica package, Comput. Phys. Commun. 228 (2018) 273–289.
  * (29) J. Matyjasek and M. Opala, Quasinormal modes of black holes: The improved semianalytic approach, Phys. Rev. D 96 (2017) 024011.
  * (30) R. A. Konoplya, A. Zhidenko and A. F. Zinhailo, Higher order WKB formula for quasinormal modes and grey-body factors: Recipes for quick and accurate calculations, Class. Quantum Grav. 36 (2019) 155002.
  * (31) J. Matyjasek and M. Telecka, Quasinormal modes of black holes. II. Padé summation of the higher-order WKB terms, Phys. Rev. D 100 (2019) 124006.
  * (32) A. Ohashi and M.-a. Sakagami, Massive quasi-normal mode, Class. Quantum Grav. 21 (2004) 3973–3984.
  * (33) Y. Hatsuda and M. Kimura, Spectral Problems for Quasinormal Modes of Black Holes, Universe 7 (2021) 476.
  * (34) Y. Hatsuda, An alternative to the Teukolsky equation, Gen Relativ Gravit 53 (2021) 93.
  * (35) S. A. Teukolsky, Rotating Black Holes: Separable Wave Equations for Gravitational and Electromagnetic Perturbations, Phys. Rev. Lett. 29 (1972) 1114–1118.
  * (36) E. Berti, V. Cardoso and M. Casals, Eigenvalues and eigenfunctions of spin-weighted spheroidal harmonics in four and higher dimensions, Phys. Rev. D 73 (2006) 024013.
  * (37) E. Berti, V. Cardoso and M. Casals, Erratum: Eigenvalues and eigenfunctions of spin-weighted spheroidal harmonics in four and higher dimensions [Phys. Rev. D 73, 024013 (2006)], Phys. Rev. D 73 (2006) 109902.
  * (38) A. Zhidenko, Quasi-normal modes of Schwarzschild–de Sitter black holes, Class. Quantum Grav. 21 (2003) 273–280.
  * (39) Y. Hatsuda, Quasinormal modes of Kerr-de Sitter black holes via the Heun function, arXiv:2006.08957.
  * (40) V. Cardoso and J. P. S. Lemos, Quasinormal modes of Schwarzschild–anti-de Sitter black holes: Electromagnetic and gravitational perturbations, Phys. Rev. D 64 (2001) 084017.
  * (41) G. Aminov, P. Arnaudo, G. Bonelli, A. Grassi and A. Tanzini, Black hole perturbation theory and multiple polylogarithms, arXiv:2307.10141 [hep-th].
  * (42) H. Onozawa, T. Mishima, T. Okamura and H. Ishihara, Quasinormal modes of maximally charged black holes, Phys. Rev. D 53 (1996) 7033–7040.
  * (43) C. M. Bender and S. A. Orszag, Advanced mathematical methods for scientists and engineers, McGraw-Hill, 1978.
  * (44) G. A. Baker, Jr. and P. Graves-Morris, Pade Approximants: Encyclopedia of Mathematics and It’s Applications, Vol. 59, Cambridge University Press, January 1996.
  * (45) M. Kimura, Note on the parametrized black hole quasinormal ringdown formalism, Phys. Rev. D 101 (2020) 064031.
  * (46) O. J. Tattersall, Quasi-Normal Modes of Hairy Scalar Tensor Black Holes: Odd Parity, Class. Quant. Grav. 37 (2020) 115007, arXiv:1911.07593 [gr-qc].
  * (47) Y. Hatsuda and M. Kimura, Semi-analytic expressions for quasinormal modes of slowly rotating Kerr black holes, Phys. Rev. D 102 (2020) 044032\.
  * (48) C. de Rham, J. Francfort and J. Zhang, Black Hole Gravitational Waves in the Effective Field Theory of Gravity, Phys. Rev. D 102 (2020) 024079, arXiv:2005.13923 [hep-th].
  * (49) S. H. Völkel, N. Franchini and E. Barausse, Theory-agnostic reconstruction of potential and couplings from quasinormal modes, Phys. Rev. D 105 (2022) 084046, arXiv:2202.08655 [gr-qc].
  * (50) 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 (2022) 124036, arXiv:2209.10564 [gr-qc].
  * (51) N. Franchini and S. H. Völkel, A parametrized quasi-normal mode framework for non-Schwarzschild metrics, arXiv:2210.14020 [gr-qc].
  * (52) S. S. Lahoz and J. Noller, Testing the Speed of Gravity with Black Hole Ringdown, arXiv:2301.10272 [gr-qc].
  * (53) R. Ghosh, N. Franchini, S. H. Völkel and E. Barausse, Quasi-normal modes of non-separable perturbation equations: the scalar non-Kerr case, arXiv:2303.00088 [gr-qc].
  * (54) S. Mukohyama, K. Takahashi, K. Tomikawa and V. Yingcharoenrat, Quasinormal Modes from EFT of Black Hole Perturbations with Timelike Scalar Profile, arXiv:2304.14304 [gr-qc].
  * (55) N. Franchini and S. H. Völkel, Testing General Relativity with Black Hole Quasi-Normal Modes, arXiv:2305.01696 [gr-qc].
  * (56) M. Sasaki and T. Nakamura, Gravitational Radiation from Extreme Kerr Black Hole, Gen. Rel. Grav. 22 (1990) 1351–1366.
  * (57) 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 (2021) 064034, arXiv:2106.06209 [gr-qc].
  * (58) 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 (2021) 124017, arXiv:2103.09870 [gr-qc].
  * (59) 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 (2022) 104009, arXiv:2207.11267 [gr-qc].
  * (60) P. Wagle, N. Yunes and H. O. Silva, Quasinormal modes of slowly-rotating black holes in dynamical Chern-Simons gravity, Phys. Rev. D 105 (2022) 124003, arXiv:2103.09913 [gr-qc].
  * (61) D. Li, P. Wagle, Y. Chen and N. Yunes, Perturbations of Spinning Black Holes beyond General Relativity: Modified Teukolsky Equation, Phys. Rev. X 13 (2023) 021029, arXiv:2206.10652 [gr-qc].
  * (62) 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 (2023) 024040, arXiv:2304.02663 [gr-qc].
  * (63) P. A. Cano, K. Fransen and T. Hertog, Ringing of rotating black holes in higher-derivative gravity, Phys. Rev. D 102 (2020) 044047, arXiv:2005.03671 [gr-qc].
  * (64) P. A. Cano, K. Fransen, T. Hertog and S. Maenaut, Quasinormal modes of rotating black holes in higher-derivative gravity, arXiv:2307.07431 [gr-qc].
  * (65) P. A. Cano, K. Fransen, T. Hertog and S. Maenaut, Gravitational ringing of rotating black holes in higher-derivative gravity, Phys. Rev. D 105 (2022) 024064, arXiv:2110.11378 [gr-qc].
  * (66) C. Molina, P. Pani, V. Cardoso and L. Gualtieri, Gravitational signature of Schwarzschild black holes in dynamical Chern-Simons gravity, Phys. Rev. D 81 (2010) 124021, arXiv:1004.4007 [gr-qc].
  * (67) O. Sarbach and E. Winstanley, On the linear stability of solitons and hairy black holes with a negative cosmological constant: The Odd parity sector, Class. Quant. Grav. 18 (2001) 2125–2146, arXiv:gr-qc/0102033.
  * (68) V. Cardoso, M. Kimura, A. Maselli and L. Senatore, Black Holes in an Effective Field Theory Extension of General Relativity, Phys. Rev. Lett. 121 (2018) 251105, arXiv:1808.08962 [gr-qc].
  * (69) 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 (2019) 044061, arXiv:1906.05155 [gr-qc].
  * (70) K. Nomura and D. Yoshida, Quasinormal modes of charged black holes with corrections from nonlinear electrodynamics, Phys. Rev. D 105 (2022) 044006, arXiv:2111.06273 [gr-qc].
  * (71) L. Hui, A. Podo, L. Santoni and E. Trincherini, An analytic approach to quasinormal modes for coupled linear systems, JHEP 03 (2023) 060, arXiv:2210.10788 [gr-qc].