# Asymptotics of linear differential systems and application to quasi-normal
modes of nonrotating black holes

David Langlois  Université de Paris, CNRS, Astroparticule et Cosmologie,
F-75006 Paris, France  Karim Noui  Institut Denis Poisson (UMR 7013),
Université de Tours, Université d’Orléans, Parc de Grandmont, 37200 Tours,
France  Université de Paris, CNRS, Astroparticule et Cosmologie, F-75006
Paris, France  Hugo Roussille  Université de Paris, CNRS, Astroparticule et
Cosmologie, F-75006 Paris, France  Institut Denis Poisson (UMR 7013),
Université de Tours, Université d’Orléans, Parc de Grandmont, 37200 Tours,
France

###### Abstract

The traditional approach to perturbations of nonrotating black holes in
General Relativity uses the reformulation of the equations of motion into a
radial second-order Schrödinger-like equation, whose asymptotic solutions are
elementary. Imposing specific boundary conditions at spatial infinity and near
the horizon defines, in particular, the quasi-normal modes of black holes. For
more complicated equations of motion, as encountered for instance in modified
gravity models with different background solutions and/or additional degrees
of freedom, such a convenient Schrödinger-like reformulation might be
unavailable, even in a generalised matricial form. In order to tackle such
cases, we present a new approach that analyses directly the first-order
differential system in its original form and extracts the asymptotic behaviour
of perturbations. As a pedagogical illustration, we apply this treatment to
the perturbations of Schwarzschild black holes and then show that the standard
quasi-normal modes can be obtained numerically by solving this first-order
system with a spectral method. This new approach paves the way for a generic
treatment of the asymptotic behaviour of black hole perturbations and the
identification of quasi-normal modes in theories of modified gravity.

##  I Introduction

The oscillations of black holes (BH) have been studied theoretically for
several decades. Today, with the first observations of gravitational waves
emitted by BH mergers, one can now hope to observe directly these oscillations
via their GW signatures, especially in the ringdown phase of the signal when
the post-merger black hole relaxes to a Kerr black hole, according to General
Relativity. One of the major goals of future detections will be to check
whether the observed oscillations coincide with the predictions based on
General Relativity (see e.g. Berti:2005ys ; Berti:2018vdi ). This is also an
ideal playground to test alternative theories of gravitation. Indeed, even if
the background BH solution may coincide with that of GR, the linear
perturbations in general obey different equations of motion.

During the ringdown phase, at least in the linear regime, the GW signal is
expected to mainly consist of a superposition of the so-called quasi-normal,
or resonant, modes (QNMs) which have been excited by the merger and then decay
via GW radiation: these modes correspond to the proper oscillation modes of
the black hole and are characterised by a complex frequency $\omega$, whose
imaginary part quantifies their damping rate.

In the simplest case of nonrotating black holes, i.e. Schwarzschild black
holes, the computation of QNMs is based on the classical papers by Regge &
Wheeler Regge:1957td  and later Zerilli Zerilli:1970se , who reformulated the
linearised Einstein equations in the frequency domain, which are first-order
with respect to the radial coordinate, as a second-order Schrödinger-like
equation. This familiar equation, with a specific potential for axial and
polar metric perturbations, is the standard starting point for the numerical
calculations or semi-analytical treatments of QNMs, using for instance well-
known methods in quantum mechanics.

Understanding the asymptotic behaviour of the perturbations at the horizon and
at spatial infinity is crucial for QNMs, which are defined by very specific
boundary conditions. Indeed, they correspond to purely outgoing radiation at
spatial infinity and ingoing radiation at the horizon. Imposing these specific
boundary conditions leads to a discrete set of allowed frequencies.

When the equations of motion of the perturbations are written as a second-
order Schrödinger equation, obtaining their asymptotic behaviour is immediate,
as it simply depends on the asymptotic behaviour of the effective potential.
In the context of modified gravity however, the problem can become more
involved for several reasons. First, the background metric can differ from the
standard GR solutions, i.e. be different from Schwarzschild in the nonrotating
case. Moreover, modified theories often involve additional fields, such as
scalar fields, which increases the number of degrees of freedom and therefore
the complexity of the linear equations of motion.

In several interesting cases, the equations of motion can be rewritten as a
generalised $N$-dimensional matrix Schrödinger-like system for $N$ fields
$\Psi_{i}$, of the typical form (see e.g. McManus:2019ulj )

$f\frac{d}{dr}\left(f\frac{d\Psi_{i}}{dr}\right)+\left(\omega^{2}-f\,V_{ij}\right)\Psi_{j}=0\,,$
(1)

where $f(r)=1-r_{s}/r$ and the $N\times N$ matrix $V_{ij}$ of radial
potentials usually vanishes or becomes a constant diagonal matrix at spatial
infinity. The boundary conditions are still easy to infer from such a
differential system.

However, one could also encounter more complex situations where such a nice
reformulation of the equations of motion is not available. Specific examples
will be given in a companion paper, in the context of Degenerate Higher-Order
Scalar-Tensor (DHOST) theories Langlois:2015cwa ; Crisostomi:2016czh ;
Achour:2016rkg ; BenAchour:2016fzp  which provide the most general viable set
of scalar-tensor theories to date. In those examples, it is not clear whether
one can rewrite the polar equations of motion as a second-order Schrödinger-
like system of the form (1). For this reason, it would be very useful to be
able to analyse directly the first-order system of equations in its original
form and to extract directly from it the asymptotic behaviour of
perturbations.

The purpose of this paper is to present such a systematic treatment of a
general first-order differential system. In order to reach this goal, we use
recent developments that appeared in the mathematical literature. These
results enable us to determine, via a systematic algorithm, the asymptotic
structure of the solutions of a generic first-order differential system. For
pedagogical reasons, we use here this algorithm to recover the asymptotic
solutions for the axial, or odd-parity, modes and for the polar, or even-
parity, modes of the standard Schwarzschild solution. This paper will be
completed by a companion paper that applies the same method to a few black
hole solutions in DHOST theories.

The outline of the paper is the following. In the next section, we review the
standard derivation for the Schwarzschild perturbations, distinguishing as
usual the axial and polar modes. In section III, we present our new approach
and show explicitly how this new method enables us to recover the usual
asymptotic solution, working directly with the first order system. We also
show how the quasi-normal modes can be computed in this new perspective. We
then present, in section IV, the general algorithm, carefully listing the
various steps of the algorithm depending on the structure of the system. We
give a summary and open some perspectives in the concluding section. A few
appendices contain some additional details.

##  II A short review on Regge-Wheeler and Zerilli equations

In this section, we review the standard procedure to derive the equations of
motion for the perturbations of a Schwarzschild black hole in general
relativity, originally obtained by Regge and Wheeler Regge:1957td  for the
axial, or odd-parity, modes and Zerilli Zerilli:1970se  for the polar, or
even-parity, modes. These equations can be shown to reduce to a Schrödinger-
like equation with an effective potential characterising the “dynamics” of the
linear perturbations.

###  II.1 Linear perturbations of Einstein equations about the Schwarzschild
black hole

We start with the four-dimensional Einstein-Hilbert action in vacuum (with no
cosmological constant) for the metric $\tensor{g}{{}_{\mu}{}_{\nu}}$,

$S[g_{\mu\nu}]=\frac{1}{16\pi G_{N}}\int\differential[4]{x}\sqrt{-g}\,R\,,$
(2)

where $g\equiv\det(g_{\mu\nu})$ is the determinant of the metric, $R$ the
four-dimensional Ricci scalar and $G_{N}$ denotes Newton’s constant, which
actually will not show up in the equations of motion since we are not
considering any matter field here.

####  II.1.1 Linearised general relativity

Given any background metric $\tensor{\overline{g}}{{}_{\mu}{}_{\nu}}$ solution
to the Einstein equations, one can introduce the perturbed metric

$\tensor{g}{{}_{\mu}{}_{\nu}}=\tensor{\overline{g}}{{}_{\mu}{}_{\nu}}+\tensor{h}{{}_{\mu}{}_{\nu}}$
(3)

where the $\tensor{h}{{}_{\mu}{}_{\nu}}$ denote the linear perturbations of
the metric. In order to derive the linear equations of motion that govern the
evolution of $h_{\mu\nu}$, one expands the Einstein-Hilbert action (2) up to
the second order in $\tensor{h}{{}_{\mu}{}_{\nu}}$. The Euler-Lagrange
equations associated with the quadratic part of this expansion then provide
the linearised equations of motion for $h_{\mu\nu}$.

By expanding (2), one obtains the following quadratic action for $h_{\mu\nu}$,

$\displaystyle S_{\rm quad}[h_{\mu\nu}]$ $\displaystyle=$
$\displaystyle\frac{1}{16\pi
G_{N}}\int\differential[4]{x}\sqrt{-\overline{g}}\left\\{-\frac{1}{2}\tensor{h}{{}_{\mu}{}_{\nu}}\tensor{h}{{}^{\mu}{}^{\nu}}\bar{R}+\frac{1}{4}h^{2}\bar{R}+h\tensor{h}{{}_{\mu}{}_{\nu}}\tensor{\bar{R}}{{}^{\mu}{}^{\nu}}+4\tensor{h}{{}_{\mu}^{\rho}}\tensor{h}{{}^{\mu}{}^{\nu}}\tensor{\bar{R}}{{}_{\nu}{}_{\rho}}-2\tensor{h}{{}^{\mu}{}^{\nu}}\tensor{h}{{}^{\rho}{}^{\sigma}}\tensor{\bar{R}}{{}_{\mu}{}_{\rho}{}_{\nu}{}_{\sigma}}\right.$
(4) $\displaystyle\left.\hskip
42.67912pt+\frac{1}{2}(\bar{\nabla}_{\mu}h)(\bar{\nabla}^{\mu}h)-2(\bar{\nabla}_{\mu}\tensor{h}{{}^{\mu}_{\nu}})(\bar{\nabla}_{\rho}\tensor{h}{{}_{\nu}^{\rho}})-(\bar{\nabla}_{\mu}h)(\bar{\nabla}_{\nu}\tensor{h}{{}^{\mu}{}^{\nu}})\right.$
$\displaystyle\left.\hskip
42.67912pt+3(\bar{\nabla}_{\nu}\tensor{h}{{}_{\mu}{}_{\rho}})(\bar{\nabla}^{\rho}\tensor{h}{{}^{\mu}{}^{\nu}})-\frac{1}{2}(\bar{\nabla}_{\rho}\tensor{h}{{}_{\mu}{}_{\nu}})(\bar{\nabla}^{\rho}\tensor{h}{{}^{\mu}{}^{\nu}})\right\\},$

where $\bar{R}_{\mu\nu\rho\sigma}$, $\bar{R}_{\mu\nu}$, $\bar{R}$ and
$\bar{\nabla}_{\mu}$ are respectively the Riemann tensor, the Ricci tensor,
the Ricci scalar and the covariant derivative associated with the background
metric $\overline{g}_{\mu\nu}$. The indices are lowered or raised with
$\overline{g}_{\mu\nu}$ and $h\equiv\overline{g}^{\mu\nu}h_{\mu\nu}$ denotes
the trace of the metric perturbation. The linearised Einstein equations are
then given by the Euler-Lagrange equations of (4) and can be written in the
form

$\displaystyle{\cal E}_{\mu\nu}$ $\displaystyle\equiv$
$\displaystyle\bar{\nabla}_{\sigma}\bar{\nabla}^{\sigma}h_{\mu\nu}+\bar{\nabla}_{\mu}\bar{\nabla}_{\nu}h+(\bar{\nabla}_{\alpha}\bar{\nabla}_{\beta}h^{\alpha\beta}-\bar{\nabla}_{\sigma}\bar{\nabla}^{\sigma}h)\overline{g}_{\mu\nu}+2\bar{\nabla}_{(\mu}\bar{\nabla}_{\alpha}h^{\alpha}_{\nu)}-6\bar{\nabla}_{\alpha}\bar{\nabla}_{(\mu}h_{\nu)}^{\alpha}$
(5) $\displaystyle+$
$\displaystyle\bar{R}_{\mu\nu}h-\bar{R}h_{\mu\nu}+\frac{1}{2}\bar{R}\,\overline{g}_{\mu\nu}h+\bar{R}^{\alpha\beta}\overline{g}_{\mu\nu}h_{\alpha\beta}+8\bar{R}_{\alpha(\mu}h^{\alpha}_{\nu)}=0\,,$

where use the standard notation $A_{(\mu\nu)}\equiv(A_{\mu\nu}+A_{\nu\mu})/2$
for the symmetrisation of any rank-2 tensor $A_{\mu\nu}$.

Let us now specialise these equations to the case where the background metric
is the Schwarzschild metric, expressed as

$\tensor{\bar{g}}{{}_{\mu}{}_{\nu}}\differential\tensor{x}{{}^{\mu}}\differential\tensor{x}{{}^{\nu}}=-\left(1-\frac{r_{s}}{r}\right)\differential{t}^{2}+\left(1-\frac{r_{s}}{r}\right)^{-1}\differential{r}^{2}+r^{2}\left(\differential{\theta}^{2}+\sin^{2}\theta\differential{\varphi}^{2}\right)\,,$
(6)

where $r_{s}=2M_{s}$ is the Schwarzschild radius, $M_{s}$ being the mass of
the black hole.

Given the spherical symmetry of the background solution, it is convenient to
decompose the metric perturbations $h_{\mu\nu}$ into (scalar, vectorial and
tensorial) spherical harmonics that are defined from the standard $Y_{\ell
m}(\theta,\varphi)$ functions and their derivatives with respect to $\theta$
and $\varphi$. They are labelled by the two multipole integers $\ell$ and $m$
(with $\ell\geq 0$ and $-\ell\leq m\leq\ell$).

Furthermore, one can distinguish axial and polar modes, which behave
differently under the parity transformation $\vec{r}\rightarrow-\vec{r}$: the
polar, or even-parity, modes transform as $(-1)^{\ell}$, similarly to the
scalar spherical harmonics $Y_{\ell m}(\theta,\varphi)$, whereas the axial, or
odd-parity, modes transform as $(-1)^{\ell+1}$. These modes can be treated
separately as they are decoupled at linear order. Moreover, we consider here
only the modes $\ell\geq 2$. The particular cases of the $\ell=0$ and $\ell=1$
modes are briefly discussed in Appendix A.3.

Since the background metric is static, it is also convenient to decompose the
time dependence of the perturbations into Fourier modes,

$F(t,r)=\int_{-\infty}^{+\infty}\differential\omega\,\tilde{F}(\omega,r)e^{-i\omega
t}\,.$ (7)

In the rest of this paper, we will use the same notation for the time-
dependent function $F$ and its Fourier transform, as there will be no
ambiguity. From a practical point of view, we simply replace every time
derivative by a multiplication by $-i\omega$ in the linearised equations,
which leads to a system of ordinary differential equations with respect to the
radial variable $r$.

In both axial and polar sectors, the equations of motion can be reduced to a
system of two first order ordinary differential equations, as we will show
below.

####  II.1.2 Axial perturbations

We choose the usual Regge-Wheeler gauge Regge:1957td  to describe the axial
modes. As recalled in Appendix A.1, in this gauge the perturbations for
$\ell\geq 2$ are parametrised by three families of functions $h_{0}^{\ell m}$,
$h_{1}^{\ell m}$ and $h_{2}^{\ell m}$ according to

$\displaystyle h_{t\theta}=\frac{1}{\sin\theta}\sum_{\ell,m}h_{0}^{\ell
m}\partial_{\varphi}{Y_{\ell m}}(\theta,\varphi),\qquad
h_{t\varphi}=-\sin\theta\sum_{\ell,m}h_{0}^{\ell m}\partial_{\theta}{Y_{\ell
m}}(\theta,\varphi),$ $\displaystyle
h_{r\theta}=\frac{1}{\sin\theta}\sum_{\ell,m}h_{1}^{\ell
m}\partial_{\varphi}{Y_{\ell m}}(\theta,\varphi),\qquad
h_{r\varphi}=-\sin\theta\sum_{\ell,m}h_{1}^{\ell m}\partial_{\theta}{Y_{\ell
m}}(\theta,\varphi),$ (8)

while the other components vanish.

For these perturbations, the equations of motion (5) reduce to the following
three equations

$\displaystyle{\cal
E}_{t\theta}=2\left(\frac{r_{s}}{r}-1-\lambda\right)h_{0}(t,r)+r(r-r_{s})\partialderivative[2]{h_{0}}{r}-2(r-r_{s})\partialderivative{h_{1}}{t}-r(r-r_{s})\partialderivative{h_{1}}{t}{r}=0\,,$
(9) $\displaystyle{\cal
E}_{r\theta}=-2\lambda\,h_{1}(t,r)-\frac{2r^{2}}{r-r_{s}}\partialderivative{h_{0}}{t}+\frac{r^{3}}{r-r_{s}}\partialderivative{h_{0}}{t}{r}-\frac{r^{3}}{r-r_{s}}\partialderivative[2]{h_{1}}{t}=0\,,$
$\displaystyle{\cal
E}_{\theta\theta}=2r_{s}h_{1}(t,r)+2r(r-r_{s})\partialderivative{h_{1}}{r}-\frac{2r^{3}}{r-r_{s}}\partialderivative{h_{0}}{t}=0\,,$

where we have introduced the notation

$\displaystyle 2\lambda\;\equiv\;\ell(\ell+1)-2\,,$ (10)

as the equations ${\cal E}_{t\varphi}=0$, ${\cal E}_{r\varphi}=0$, ${\cal
E}_{\varphi\varphi}=0$ and ${\cal E}_{\theta\varphi}=0$ are identical to the
above ones.

Since there are only two independent functions, $h_{0}$ and $h_{1}$, one
expects one of the above equations to be redundant. This is indeed verified by
noting the following relation between the equations (9) and their derivatives,
written now in the frequency domain,

$\derivative{\mathcal{E}_{r\theta}}{r}+\frac{ir^{2}\omega}{(r-r_{s})^{2}}\mathcal{E}_{t\theta}+\frac{r_{s}}{r(r-r_{s})}\mathcal{E}_{r\theta}+\frac{\lambda}{r(r-r_{s})}\mathcal{E}_{\theta\theta}=0\,.$
(11)

This shows that the two equations $\mathcal{E}_{r\theta}=0$ and
$\mathcal{E}_{\theta\theta}=0$ are sufficient to fully describe the dynamics
of axial perturbations. As a consequence, the initial system (9) reduces to

$\displaystyle\derivative{Y}{r}=M(r)Y\,,\quad
M(r)=\left(\begin{array}[]{cc}{2}/{r}&2i\lambda(r-r_{s})/{r^{3}}-i\omega^{2}\\\
-{ir^{2}}/{(r-r_{s})^{2}}&-{r_{s}}/{r(r-r_{s})}\end{array}\right)\,,$ (14)

where the two components of the column vector $Y\equiv{}^{T}\\!(Y_{1},Y_{2})$
are $Y_{1}(r)\equiv h_{0}(r)$ and $Y_{2}(r)\equiv h_{1}(r)/\omega$. Notice
that we divided the variable $h_{1}(r)$ by $\omega$ in the definition of
$Y_{2}$ in order to get a system which does not involve powers of $\omega$
higher than 2, or equivalently which is at most second order in time if one
inverts the Fourier transform (7).

####  II.1.3 Polar perturbations

After fixing the gauge, polar perturbations are parametrised by four families
of functions $H_{0}^{\ell m},H_{1}^{\ell m},H_{2}^{\ell m}$ and $K^{\ell m}$
as shown in Appendix A.2. The nonvanishing metric perturbations then read

$\displaystyle h_{tt}=A(r)\sum_{\ell,m}H_{0}^{\ell m}(t,r)Y_{\ell
m}(\theta,\varphi),\quad h_{tr}=\sum_{\ell,m}H_{1}^{\ell m}(t,r)Y_{\ell
m}(\theta,\varphi),$ (16) $\displaystyle
h_{rr}=\frac{1}{A(r)}\sum_{\ell,m}H_{2}^{\ell m}(t,r)Y_{\ell
m}(\theta,\varphi),\quad h_{ab}=\sum_{\ell,m}K^{\ell m}(t,r)g_{ab}Y_{\ell
m}(\theta,\varphi),$

where $A(r)\equiv 1-r_{s}/r$ is included in the definitions for later
convenience, and the indices $a$ or $b$ in the last equation are the angles
$\theta$ or $\varphi$.

The linearised Einstein’s equations yield seven distinct equations, which can
be found in (214) of Appendix B. After a few manipulation, also discussed in
Appendix B, one finds that these equations of motion can be reduced to two
first-order equations only. In the frequency domain, they read

$\displaystyle\derivative{Y}{r}=M(r)Y\,,\quad M(r)=\frac{1}{3r_{s}+2\lambda
r}\begin{pmatrix}\frac{r_{s}(3r_{s}+(\lambda-2)r)-2r^{4}\omega^{2}}{r(r-r_{s})}&\frac{2i(\lambda+1)(r_{s}+\lambda
r)+2ir^{3}\omega^{2}}{r^{2}}\\\ \frac{ir\left(9r_{s}^{2}-8\lambda
r^{2}+8(\lambda-1)r_{s}r\right)+4ir^{5}\omega^{2}}{2(r-r_{s})^{2}}&\frac{2r^{4}\omega^{2}-r_{s}(3r_{s}+3\lambda
r+r)}{r(r-r_{s})}\\\ \end{pmatrix}\ $ (17)

where now the two components of $Y$ are defined by $Y_{1}(r)\equiv K(r)$ and
$Y_{2}(r)\equiv H_{1}(r)/\omega$. Similarly to the axial sector, the
definition of $Y_{2}$ is motivated by the fact that the resulting system
involves at most $\omega^{2}$ terms.

###  II.2 Schrödinger-like equation and effective potential

In both axial and polar sectors, the equations of motion have been recast in
the form of a system consisting simply of two first-order differential
equations (with respect to the radial variable), namely (14) for axial
perturbations and (17) for polar perturbations. In both cases, we now recall
how this system can be rewritten as a Schrödinger-like equation.

####  II.2.1 From the first order system to the Schrödinger-like equation

As shown in Regge:1957td  and Zerilli:1970se , one can rewrite these systems
as a single second order (in radial derivatives) Schrödinger-like equation for
a unique dynamical variable. Reformulating a first order system of this kind
as a Schrödinger equation is, in general, not an easy task because one has to
ensure that the Schrödinger equation is second order in time and in space. It
requires, in particular, a decoupling of the dynamical variables involved in
the original first order system and a “clever” choice for the dynamical
variable that should satisfy the second order Schrödinger equation.

In this section, we will describe how this works for the two systems (14) and
(17) which take the general form

$\displaystyle\derivative{Y}{r}=M(r)Y\,,$ (18)

where the coefficients of the matrix $M$ are polynomials (of degree at most 2)
in $\omega$ and rational functions in $r$.

First, we consider the general (linear) change of vector

$\displaystyle Y(r)=P(r)\hat{Y}(r)\,,$ (19)

where $\hat{Y}$ is a new column vector and the two dimensional invertible
matrix $P$ has not been fixed at this stage. We also define a new radial
coordinate $r_{*}$ and introduce the “Jacobian” of the transformation
$n(r)\equiv{{\rm d}r}/{{\rm d}r_{*}}$. Now, the idea is to show that it is
possible to find a matrix $P$ such that the new system satisfied by $\hat{Y}$
takes the canonical form

$\derivative{\hat{Y}}{r_{*}}=\begin{pmatrix}0&1\\\
V(r)-\omega^{2}&0\end{pmatrix}\hat{Y}\,,$ (20)

where the potential $V(r)$ depends on $r$, but not on $\omega$. Somehow, the
first component $\hat{Y}_{1}$ plays the role of the “momentum” conjugate to
the second component $\hat{Y}_{2}$ which would immediately implies that
$\hat{Y}_{1}$ is the “canonical” variable satisfying the required Schrödinger-
like equation

$\displaystyle\derivative[2]{\hat{Y}_{1}}{r_{*}}+\left(\omega^{2}-V(r)\right)\hat{Y}_{1}=0\,.$
(21)

####  II.2.2 Axial modes

Applying this procedure to the system (14) for the axial perturbations is
rather simple111When one changes variables according to (19), the new variable
$\hat{Y}$ satisfies the differential equation
$\displaystyle\derivative{\hat{Y}}{r_{*}}=\hat{M}\hat{Y}\,,\qquad\hat{M}\equiv
n(r)(P^{-1}MP-P^{-1}P^{\prime})\,,$ (22) where $P^{\prime}$ is the derivative
of $P$ with respect to $r$, $M$ is the matrix introduced in (14) while
$\hat{M}$ is the matrix entering in the system (20). They take a similar form
$M=M_{[0]}+\omega^{2}M_{[2]}$ and
$\hat{M}=\hat{M}_{[0]}+\omega^{2}\hat{M}_{[2]}$ where the expressions of
$M_{[0]}$, $M_{[2]}$, $\hat{M}_{[0]}$ and $\hat{M}_{[2]}$ are trivially
obtained. As $P$ does not depend on $\omega$, the relation between $M$ and
$\hat{M}$ translates into the two matricial relations
$\hat{M}_{[2]}=n(r)P^{-1}M_{[2]}P$ and
$\hat{M}_{[0]}=n(r)(P^{-1}M_{[0]}P-P^{-1}P^{\prime})\,$ which can be viewed as
8 equations for the 6 unknowns $n(r)$, $V(r)$ together with the four
components of $P$. Interestingly, the system is not overdetermined and admits
a solution for $P$ (23), for the potential $V(r)$ (25) and for the function
$n(r)$ which can be shown to be associated with the tortoise coordinate (24).
Details can be found in the Appendix D of the companion paper.. Indeed, the
appropriate transition matrix is given by

$\displaystyle P(r)=\begin{pmatrix}1-{r_{s}}/{r}&r\\\
-{ir^{2}}/{(r-r_{s})}&0\end{pmatrix}\,,$ (23)

while $n(r)=1-{r_{s}}/{r}$, which means that $r_{*}$ coincides with the
“tortoise” coordinate,

$r_{*}\equiv\int\frac{dr}{1-{r_{s}}/{r}}=r+r_{s}\ln(r/r_{s}-1)\,.$ (24)

Finally the effective potential $V_{\text{odd}}(r)$ for the axial
perturbations takes the form

$V_{\text{odd}}(r)=\left(1-\frac{r_{s}}{r}\right)\frac{2(\lambda+1)r-3r_{s}}{r^{3}}.$
(25)

Note that this potential vanishes both at spatial infinity
($r\rightarrow+\infty$) and at the horizon ($r\rightarrow r_{s}$).

####  II.2.3 Polar modes

The case of polar perturbations is slightly more involved. Starting from the
system (17), we find that the transition matrix leading to a canonical form
(20) is given by222 We follow the same method as the one described in the
previous footnote for the axial mode.

$\displaystyle P=\begin{pmatrix}\frac{3r_{s}^{2}+3\lambda
r_{s}r+2r^{2}\lambda(\lambda+1)}{2r^{2}(3r_{s}+2\lambda r)}&1\\\
-i+\frac{ir_{s}}{2(r-r_{s})}+\frac{3ir_{s}}{3r_{s}+2\lambda
r}&-\frac{ir^{2}}{r-r_{s}}\end{pmatrix}\,,$ (26)

with, in addition, $n(r)=1-{r_{s}}/{r}$, which means that $r_{*}$ is still the
tortoise coordinate (24). Finally, the corresponding potential
$V_{\text{even}}(r)$ reads

$V_{\text{even}}(r)=\left(1-\frac{r_{s}}{r}\right)\frac{9r_{s}^{3}+18r_{s}^{2}r\lambda+12r_{s}r^{2}\lambda^{2}+8r^{3}\lambda^{2}(1+\lambda)}{r^{3}(3r_{s}+2r\lambda)^{2}}\,.$
(27)

Figure 1: Illustration of the effective potentials (for axial and polar modes)
for a Schwarzschild black hole. The parameters are such that $r_{s}=2$ (i.e.
the mass of the black is $1$ in natural units) and $\ell=2$ here.

Despite their different analytic forms, we notice in Fig.(1) that the
potentials $V_{\text{odd}}(r)$ and $V_{\text{even}}(r)$ are quite similar,
although distinct. In fact, there exists an underlying symmetry between these
two potentials, further explained in Chandrasekhar:1985kt , leading to the
isospectrality theorem which states that the spectra of axial and polar
perturbations are exactly the same.

####  II.2.4 Quasi-normal modes and boundary conditions

Finding quasi-normal modes requires to impose the appropriate boundary
conditions: the modes must be outgoing at infinity and ingoing at the horizon.

Since both $V_{\text{odd}}$ and $V_{\text{even}}$ go to zero at infinity and
at the horizon, equation (21) becomes asymptotically

$\derivative[2]{\hat{X}_{1}}{r_{*}}+\omega^{2}\hat{X}_{1}\approx
0\qquad(r_{*}\rightarrow\pm\infty)\,,$ (28)

where $\approx$ is an equality up to sub-leading corrections333Near the
horizon, $V={\cal O}(r-r_{s})$ for both potentials, hence we assume
$r_{s}^{2}\omega^{2}\gg r/r_{s}-1$. At infinity, $V={\cal O}(1/r^{2})$ for
both potentials as well, hence we assume $\omega^{2}r^{2}\gg 1$ in this
limit.. Therefore, at both boundaries, the function $\hat{X}_{1}$ behaves like

$\displaystyle\hat{X}_{1}(r)\approx{\cal A}\,e^{i\omega r_{*}}+{\cal
B}\,e^{-i\omega r_{*}}\,,$ (29)

where ${\cal A}$ and ${\cal B}$ are integration constants which take different
values at the horizon and at infinity.

The physical interpretation of these modes is more transparent if we include
their time dependence explicitly, which gives

$\left\\{\begin{aligned} \hat{X}_{1}(t,r)&\approx{\cal
A}_{\text{hor}}\,e^{-i\omega(t-r_{*})}+{\cal
B}_{\text{hor}}\,e^{-i\omega(t+r_{*})}\quad\text{when}\,\,{r\longrightarrow
r_{s}},\\\ \hat{X}_{1}(t,r)&\approx{\cal
A}_{\infty}\,e^{-i\omega(t-r_{*})}+{\cal
B}_{\infty}\,e^{-i\omega(t+r_{*})}\quad\text{when}\,\,{r\longrightarrow\infty}\,.\end{aligned}\right.$
(30)

We can interpret each term as a radially propagating wave: the terms
proportional to ${\cal A}_{\text{hor}}$ and ${\cal A}_{\infty}$ are outgoing
while the terms proportional to ${\cal B}_{\text{hor}}$ and ${\cal
B}_{\infty}$ are ingoing. Imposing a purely outgoing behaviour at infinity and
a purely ingoing behaviour at the horizon, i.e. such that ${\cal
A}_{\text{hor}}=0$ and ${\cal B}_{\infty}=0$ severely restricts the possible
values of $\omega$. These values can be found numerically by integrating the
Schrödinger-like equation (see Leaver:1985ax  and the reviews Kokkotas:1999bd
; Nollert:1999ji ; Berti:2009kk ; Konoplya:2011qq ).

Finally, one can easily deduce the asymptotic expansion of the original
gravitational perturbations using the transformations (19). For the axial
modes, the leading order terms at infinity are thus given by

$\displaystyle h_{0}(r)\approx i\omega r\left({\cal A}_{\infty}e^{i\omega
r_{*}}-{\cal B}_{\infty}e^{-i\omega r_{*}}\right),\quad
h_{1}(r)\approx-i\omega r\left({\cal A}_{\infty}e^{i\omega r_{*}}+{\cal
B}_{\infty}e^{-i\omega r_{*}}\right)\,,$ (31)

while the leading order terms at the horizon read

$\displaystyle h_{0}(r)\approx i\omega r_{s}\left({\cal
A}_{\text{hor}}e^{i\omega r_{*}}-{\cal B}_{\text{hor}}e^{-i\omega
r_{*}}\right)\,,\quad h_{1}(r)\approx-\frac{i\omega
r_{s}^{2}}{\varepsilon}\left({\cal A}_{\text{hor}}e^{i\omega r_{*}}+{\cal
B}_{\text{hor}}e^{-i\omega r_{*}}\right)\,,$ (32)

where we have introduced the variable $\varepsilon\equiv r-r_{s}$ which
satisfies $\varepsilon\ll r_{s}$ near the horizon.

For the polar modes, the leading order terms at infinity are

$\displaystyle K(r)\approx i\omega\left({\cal A}_{\infty}e^{i\omega
r_{*}}-{\cal B}_{\infty}e^{-i\omega r_{*}}\right)\,,\quad H_{1}(r)\approx
r\omega^{2}\left({\cal A}_{\infty}e^{i\omega r_{*}}-{\cal
B}_{\infty}e^{-i\omega r_{*}}\right)\,,$ (33)

while the leading terms at the horizon are a bit more involved and read

$\displaystyle K(r)$ $\displaystyle\approx$
$\displaystyle\frac{\lambda+1+2i\omega r_{s}}{r_{s}}{\cal
A}_{\text{hor}}e^{i\omega r_{*}}+\frac{\lambda+1-2i\omega r_{s}}{r_{s}}{\cal
B}_{\text{hor}}e^{-i\omega r_{*}}\,,$ (34) $\displaystyle H_{1}(r)$
$\displaystyle\approx$ $\displaystyle\frac{ir_{s}\omega(1-2i\omega
r_{s})}{2\varepsilon}{\cal A}_{\text{hor}}e^{i\omega
r_{*}}+\frac{ir_{s}\omega(1+2i\omega r_{s})}{2\varepsilon}{\cal
B}_{\text{hor}}e^{-i\omega r_{*}}\,.$ (35)

In the next section, we will recover these asymptotic behaviours in a
completely different way.

##  III First order approach to Schwarzschild perturbations

As we have seen in subsection II.2, finding a (second-order) Schrödinger-like
equation for the metric perturbations starting from the Einstein equations
requires some manipulations of the equations of motion and an appropriate
choice of the function that verifies the Schrödinger-like equation.

The rest of this paper will be devoted to obtaining the asymptotic behaviours
of the perturbations by using a different method. Although this is of course
not necessary for the perturbations of Schwarzschild in General Relativity,
our method may prove to be very useful in situations where a Schrödinger-like
system is not obvious to find or even impossible to reach. In such a case, one
would need an alternative method to determine the asymptotic limits of the
solutions of the system, and from them, to compute the quasi-normal modes.

The general method will be described in a systematic way in the next section.
As the general procedure is somewhat tedious, we have preferred to present it
first, in a pedestrian way, for the perturbations of Schwarzschild. A more
mathematically-minded reader might prefer to jump directly to the next section
and later come back to this section to find a particular application of the
general method.

###  III.1 Method

Ignoring the traditional Schrödinger reformulation, we now go back to the
original first-order system given in (14) or (17). Schematically, we thus have
a first-order system of the form

$\displaystyle\derivative{{Y}}{r}=M(r)Y\,,\qquad$ (36)

where $Y(r)$ is a column vector and $M(r)$ a square matrix. In order to study
the system at spatial infinity, say, i.e. when $r\rightarrow\infty$, one can
expand the matrix $M(r)$ in powers of $r$,

$M(r)=M_{p}\,r^{p}+\dots+M_{0}+M_{-1}\,\frac{1}{r}+{\cal O}(\frac{1}{r^{2}})$
(37)

where all the matrix coefficients $M_{i}$ are $r$-independent. We stop here
the expansion at order $1/r$, which is sufficient for the simplest cases, but
higher orders might be needed in general.

If all matrices $M_{i}$ are diagonal, it is immediate to integrate the
truncated system, which then consists of $n$ ordinary differential equations
of the form

$y^{\prime}(r)=\left(\lambda_{p}r^{p}+\dots+\lambda_{0}+\frac{\mu}{r}\right)y(r)\,,$
(38)

whose solution is

$y(r)=y_{0}\,e^{q(r)}r^{\mu}\,,\quad
q(r)=\frac{\lambda_{p}}{p+1}r^{p+1}+\dots+\lambda_{0}r\,.$ (39)

Putting together these $n$ solutions, we thus get the solution to the system
(36), assuming all matrices $M_{i}$ in (37) are diagonal, in the form

$Y(r)=e^{\mathbf{\Upsilon}(r)}r^{\mathbf{\Delta}}\mathbf{F}(r)Y_{0}$ (40)

where $Y_{0}$ is a constant vector, corresponding to the $n$ integration
constants, $\mathbf{\Upsilon}$ is a diagonal matrix whose coefficients are
polynomials of degree at most $p+1$, $\mathbf{\Delta}$ is a constant diagonal
matrix and $\mathbf{F}(r)$ is a matrix which is regular at infinity (i.e.
whose limit is finite).

Of course, in general, the matrices $M_{i}$ are not diagonal but, remarkably,
it is always possible to transform the truncated system into a fully diagonal
system, in a finite number of steps following an algorithm introduced in
wasow_asymptotic_1965 ; balser_computation_1999 ; barkatou_algorithm_1999 ;
pflugel_root-free_2019 , which we will present in full details in the next
section.

At each step in the algorithm, one introduces a new vector $\tilde{Y}$,
related to the vector $Y$ of the previous step by

$Y=P\tilde{Y}\,,$

where $P$ is an invertible matrix so that the previous system (36) is
transformed into a new, but equivalent, system of the form

$\displaystyle\derivative{\tilde{Y}}{r}=\tilde{M}(r)\tilde{Y},\qquad\tilde{M}(r)\equiv
P^{-1}MP-P^{-1}\derivative{P}{r}\,.$ (41)

The idea is then to choose an appropriate transition matrix $P$ at each step
in order to diagonalise, order by order, the matrices that appear in the
expansion of $M$. Once all the matrices are diagonalised, one can integrate
directly the diagonal system, as we have seen earlier, and obtain the general
asymptotic solution of the system.

For the asymptotic behaviour near the horizon, one proceeds in the same way by
noting that the variable $z=1/(r-r_{s})$ goes to infinity when $r\rightarrow
r_{s}$. In the rest of this section, we will illustrate the algorithm by
considering in turn the asymptotic behaviours of the axial and polar modes.

###  III.2 Axial modes

The analysis of the asymptotic behaviour of the first order system (14) is
relatively simple and instructive. We recall that the system is of the form

$\displaystyle\derivative{{Y}}{r}={M}(r){Y}\,,$ (42)

with

$\displaystyle Y(r)\equiv\left(\begin{array}[]{c}h_{0}(r)\\\
h_{1}(r)/\omega\end{array}\right)\,,\qquad
M(r)\equiv\left(\begin{array}[]{cc}{2}/{r}&2i\lambda(r-r_{s})/{r^{3}}-i\omega^{2}\\\
-{ir^{2}}/{(r-r_{s})^{2}}&-{r_{s}}/{r(r-r_{s})}\end{array}\right)\,.$ (47)

####  III.2.1 Asymptotic analysis at spatial infinity

We first study the asymptotic behaviour at spatial infinity, i.e. when
$r\rightarrow\infty$. The asymptotic expansion of the matrix $M(r)$ at large
$r$ reads

$\displaystyle M(r)=M_{0}+\frac{1}{r}M_{-1}+{\cal
O}\left(\frac{1}{r^{2}}\right)\,,\qquad
M_{0}\equiv-i\left(\begin{array}[]{cc}0&\omega^{2}\\\
1&0\end{array}\right)\,,\quad M_{-1}\equiv{2}\left(\begin{array}[]{cc}1&0\\\
-ir_{s}&0\end{array}\right)\,.$ (52)

The leading term $M_{0}$ is diagonalisable and one can go to a basis where it
is diagonal, by introducing the new vector $Y^{(1)}$ defined by

$\displaystyle Y\equiv P_{(1)}{Y}^{(1)}\,,\qquad
P_{(1)}=\left(\begin{array}[]{cc}\omega&-\omega\\\ 1&1\end{array}\right)\,.$
(55)

According to (41), this gives the new system

$\displaystyle\derivative{{Y}^{(1)}}{r}={M}^{(1)}{Y}^{(1)}\,,\quad{M}^{(1)}(r)={M}^{(1)}_{0}+\frac{1}{r}{M}^{(1)}_{-1}+{\cal
O}\left(\frac{1}{r^{2}}\right)\,,$ (56)

with

$\displaystyle{M}^{(1)}_{0}\equiv\left(\begin{array}[]{cc}-i\omega&0\\\
0&i\omega\end{array}\right)\,,\quad{M}^{(1)}_{-1}\equiv\left(\begin{array}[]{cc}-i\omega
r_{s}+1&i\omega r_{s}-1\\\ -i\omega r_{s}-1&i\omega
r_{s}+1\end{array}\right).$ (61)

We need some extra work to diagonalise the next-to-leading order matrix
${M}^{(1)}_{-1}$ while keeping the leading order matrix diagonal.

This can be achieved by introducing a new vector $Y^{(2)}$ defined by

$\displaystyle Y^{(1)}\equiv P_{(2)}Y^{(2)}\,,\qquad
P_{(2)}\;=I+\frac{1}{r}\Xi\,,$ (62)

where $I$ is the identity matrix and $\Xi$ a constant matrix. Indeed, it is
immediate to see that such a change of variable leads to the equivalent
differential system,

$\displaystyle\derivative{{Y}^{(2)}}{r}={M}^{(2)}{Y}^{(2)}\,,\quad{M}^{(2)}(r)={M}^{(2)}_{0}+\frac{1}{r}{M}^{(2)}_{-1}+{\cal
O}\left(\frac{1}{r^{2}}\right)\,,$ (63)

with

$\displaystyle{M}^{(2)}_{0}={M}^{(1)}_{0}\,,\quad{M}^{(2)}_{-1}={M}^{(1)}_{-1}+[{M}^{(1)}_{0},\Xi]\,.$
(64)

The leading matrix remains unchanged while one can easily find a matrix $\Xi$
so that ${M}^{(2)}_{-1}$ is diagonal. Notice that $\Xi$ appears in (64) only
in a commutator with the diagonal matrix $M_{0}^{(1)}$, hence the diagonal
part of $\Xi$ is irrelevant and we can already fix the diagonal terms of $\Xi$
to $0$. In this case, the solution to (64) with ${M}^{(2)}_{-1}$ diagonal is
unique and given by

$\displaystyle\Xi\;=\;\frac{1}{2i\omega}\left(\begin{array}[]{cc}0&i\omega
r_{s}-1\\\ {i\omega r_{s}+1}&0\end{array}\right)\,.$ (67)

We have thus managed to obtain a fully diagonalised system, up to order $1/r$,
with the matrix

$\displaystyle{M}^{(2)}(r)=\left(\begin{array}[]{cc}-i\omega&0\\\
0&i\omega\end{array}\right)+\frac{1}{r}\left(\begin{array}[]{cc}1-i\omega
r_{s}&0\\\ 0&1+i\omega r_{s}\end{array}\right)+{\cal
O}\left(\frac{1}{r^{2}}\right).$ (72)

This system can be immediately integrated in the form (40), and the asymptotic
solution reads

$\displaystyle{Y}^{(2)}(r)\;=\;\left(1+{\cal
O}\left({1}/{r}\right)\right)\left(\begin{array}[]{c}c_{-}\,e^{-i\omega
r}r^{1-i\omega r_{s}}\\\ c_{+}\,e^{+i\omega r}r^{1+i\omega
r_{s}}\end{array}\right)\,,$ (75)

where $c_{\pm}$ are integration constants. Taking into account the time
dependency $e^{-i\omega t}$ of the modes, the two components $Y^{(2)}_{\mp}$
of $Y^{(2)}$ are of the form

$\displaystyle e^{-i\omega t}\,Y^{(2)}_{\mp}(r)=\left(1+{\cal
O}\left({1}/{r}\right)\right)c_{\mp}re^{-i\omega(t\pm(r+r_{s}\ln
r))}=c_{\mp}\left(r+{\cal O}(1)\right)e^{-i\omega(t\pm r_{*})}\,,$ (76)

where it is convenient to use the “tortoise” coordinate $r_{*}$, introduced in
(24), noting that

$\displaystyle r_{*}=r+r_{s}\ln(r/r_{s}-1)=r+r_{s}\ln r+{\cal O}(1)\,.$ (77)

As a consequence, one can identify $Y^{(2)}_{-}$ as an ingoing mode and
$Y^{(2)}_{+}$ as an outgoing mode at spatial infinity.

Finally, we can return to the original vector $Y$ thanks to the transformation

$\displaystyle
Y=P_{(1)}P_{(2)}Y^{(2)}=\left(\begin{array}[]{cc}\omega&-\omega\\\
1&1\end{array}\right)\left(1+\frac{\Xi}{r}\right)Y^{(2)}\,,$ (80)

in order to obtain the asymptotic expansion of the two original gravitational
perturbations $h_{0}$ and $h_{1}$ at spatial infinity,

$\displaystyle h_{0}(r)$ $\displaystyle=$
$\displaystyle\omega\left(c_{-}e^{-i\omega r_{*}}-c_{+}e^{+i\omega
r_{*}}\right)\left(r+{\cal O}(1)\right)\,,$ (81) $\displaystyle h_{1}(r)$
$\displaystyle=$ $\displaystyle\omega\left(c_{-}e^{-i\omega
r_{*}}+c_{+}e^{+i\omega r_{*}}\right)\left(r+{\cal O}(1)\right)\,.$ (82)

One can immediately check that these expressions agree with the asymptotic
expansion (31) obtained from the Schrödinger-like equation (with
$c_{-}=-i{\cal B}_{\infty}$ and $c_{+}=-i{\cal A}_{\infty}$).

####  III.2.2 Asymptotic analysis near the black hole horizon

Let us now study the behaviour of the axial modes near the horizon. In this
case, it is convenient to introduce the new radial variable $\varepsilon\equiv
r-r_{s}$ and expand the matrix $M$ for the system (47) in powers of
$\varepsilon$. One finds444Note that $\varepsilon$ goes to zero here, in
contrast to the previous case where the variable $r$ was going to infinity.
One could work in a fully analogous system by using the variable
$z=1/\varepsilon$, with the system
$\frac{dY}{dz}=\tilde{M}(z)Y\,,\quad\tilde{M}=-\frac{1}{z^{2}}M(z^{-1})=-M_{2}-M_{1}\frac{1}{z}-M_{0}\frac{1}{z^{2}}\,.$
(83) In the present case, one must push the expansion up to order $1/z^{2}$
because the leading matrix $M_{2}$ is nilpotent.

$\displaystyle
M(\varepsilon)=\frac{1}{\varepsilon^{2}}M_{2}+\frac{1}{\varepsilon}M_{1}+M_{0}+{\cal
O}(\varepsilon)\,,$ (84)

with the matrix coefficients

$\displaystyle M_{2}\equiv\left(\begin{array}[]{cc}0&0\\\
-ir_{s}^{2}&0\end{array}\right)\,,\quad
M_{1}\equiv\left(\begin{array}[]{cc}0&0\\\
-2ir_{s}&-1\end{array}\right)\,,\quad
M_{0}\equiv\left(\begin{array}[]{cc}2/r_{s}&-i\omega^{2}\\\
-i&1/r_{s}\end{array}\right)\,.$ (91)

An important difference with the previous situation is that the leading term
$M_{2}$ is no longer diagonalisable but nilpotent instead. We thus need to
first perform a transformation that yields a diagonalisable leading matrix,
taking advantage of the derivative term in (41). This can be done with the
transformation

$\displaystyle Y\equiv P_{(1)}Y^{(1)}\,,\quad
P_{(1)}(\varepsilon)\equiv\left(\begin{array}[]{cc}1&0\\\
0&1/\varepsilon\end{array}\right)\,,$ (94)

leading to the new system

$\displaystyle\derivative{{Y}^{(1)}}{\varepsilon}={M}^{(1)}{Y}^{(1)}\,,\quad{M}^{(1)}(\varepsilon)=-\frac{1}{\varepsilon}\left(\begin{array}[]{cc}0&i\omega^{2}\\\
ir_{s}^{2}&0\end{array}\right)+{\cal O}(1)\,.$ (97)

The transformation (94) has eliminated the term in $1/\varepsilon^{2}$ in the
expansion and the leading term $M^{(1)}_{1}$ is now diagonalisable, so that
only the expansion of $M^{(1)}$ up to order $1/\varepsilon$ is required (see
discussion in the footnote). It is worth noticing that $M^{(1)}_{1}$ receives
contributions from $M_{2}$, $M_{1}$ and $M_{0}$. In particular, some of its
coefficients involve the frequency $\omega$ which is originally present only
in $M_{0}$.

The final step of the analysis consists in diagonalising the system (97), via
the transformation

$\displaystyle Y^{(1)}=P_{(2)}Y^{(2)}\,,\quad
P_{(2)}\equiv\left(\begin{array}[]{cc}\omega&-\omega\\\
r_{s}&r_{s}\end{array}\right)\,,$ (100)

leading to

$\displaystyle\derivative{{Y}^{(2)}}{\varepsilon}=M^{(2)}Y^{(2)}\,,\qquad
M^{(2)}(\varepsilon)\equiv\frac{1}{\varepsilon}\left(\begin{array}[]{cc}-i\omega
r_{s}&0\\\ 0&i\omega r_{s}\end{array}\right)+{\cal O}(1)\,.$ (103)

Integrating this equation yields

$\displaystyle{Y}^{(2)}(\varepsilon)=(1+{\cal
O}(\varepsilon))\left(\begin{array}[]{c}c_{-}\varepsilon^{-i\omega r_{s}}\\\
c_{+}\varepsilon^{+i\omega r_{s}}\end{array}\right)=(1+{\cal
O}(\varepsilon))\left(\begin{array}[]{c}c_{-}e^{-i\omega r_{*}}\\\
c_{+}e^{+i\omega r_{*}}\end{array}\right)\,,$ (108)

where we have again expressed the result in terms of the tortoise coordinate
$r_{*}$, which behaves as $r_{*}=r_{s}\ln\varepsilon+{\cal O}(1)$ near the
horizon. One can immediately recognize the ingoing and outgoing modes at the
horizon.

Finally, one can return to the original functions, via
$Y=P_{(1)}P_{(2)}Y^{(2)}$, and derive the expressions

$\displaystyle h_{0}(r)$ $\displaystyle=$
$\displaystyle\omega\left(c_{-}e^{-i\omega r_{*}}-c_{+}e^{+i\omega
r_{*}}\right)\left(1+{\cal O}(\varepsilon)\right)\,,$ (109) $\displaystyle
h_{1}(r)$ $\displaystyle=$ $\displaystyle\frac{\omega
r_{s}}{\varepsilon}\left(c_{-}e^{-i\omega r_{*}}+c_{+}e^{+i\omega
r_{*}}\right)\left(1+{\cal O}(\varepsilon)\right)\,,$ (110)

which coincide with the asymptotic expansions (32) obtained from the
Schrödinger-like equation (with $c_{-}=-ir_{s}{\cal B}_{\text{hor}}$,
$c_{+}=-ir_{s}{\cal A}_{\text{hor}}$).

###  III.3 Polar modes

The dynamics of the polar perturbations is described by the first-order system
(17), of the form

$\displaystyle\derivative{{Y}}{r}={M}(r){Y}\,,\qquad\text{with}\qquad
Y(r)\equiv\left(\begin{array}[]{c}K(r)\\\
H_{1}(r)/\omega\end{array}\right)\,,$ (113)

and the matrix

$\displaystyle M(r)=\frac{1}{3r_{s}+2\lambda
r}\begin{pmatrix}\frac{r_{s}(3r_{s}+(\lambda-2)r)-2r^{4}\omega^{2}}{r(r-r_{s})}&\frac{2i(\lambda+1)(r_{s}+\lambda
r)+2ir^{3}\omega^{2}}{r^{2}}\\\ \frac{ir\left(9r_{s}^{2}-8\lambda
r^{2}+8(\lambda-1)r_{s}r\right)+4ir^{5}\omega^{2}}{2(r-r_{s})^{2}}&\frac{2r^{4}\omega^{2}-r_{s}(3r_{s}+3\lambda
r+r)}{r(r-r_{s})}\\\ \end{pmatrix}\,.$ (114)

####  III.3.1 Asymptotic analysis at spatial infinity

Expanding (114) in powers of $r$, one gets

$\displaystyle M(r)$ $\displaystyle=$ $\displaystyle\begin{pmatrix}0&0\\\
\frac{i\omega^{2}}{\lambda}&0\end{pmatrix}r^{2}+\begin{pmatrix}-\frac{\omega^{2}}{\lambda}&0\\\
\frac{ir_{s}\omega^{2}(4\lambda-3)}{2\lambda^{2}}&\frac{\omega^{2}}{\lambda}\end{pmatrix}r+\begin{pmatrix}-\frac{(2\lambda-3)r_{s}\omega^{2}}{2\lambda^{2}}&\frac{i\omega^{2}}{\lambda}\\\
-2i+\frac{3i\left(4\lambda^{2}-4\lambda+3\right)r_{s}^{2}\omega^{2}}{4\lambda^{3}}&\frac{(2\lambda-3)r_{s}\omega^{2}}{2\lambda^{2}}\\\
\end{pmatrix}$ (115) $\displaystyle+$
$\displaystyle\frac{1}{r}\begin{pmatrix}-\frac{\left(4\lambda^{2}-6\lambda+9\right)r_{s}^{2}\omega^{2}}{4\lambda^{3}}&-\frac{3ir_{s}\omega^{2}}{2\lambda^{2}}\\\
\frac{i\left(8(1-2\lambda)\lambda^{3}r_{s}-(27-4\lambda(\lambda(8\lambda-9)+9))r_{s}^{3}\omega^{2}\right)}{8\lambda^{4}}&\frac{\left(4\lambda^{2}-6\lambda+9\right)r_{s}^{2}\omega^{2}}{4\lambda^{3}}\end{pmatrix}+\mathcal{O}\left(\frac{1}{r^{2}}\right)\,.$

In contrast with the axial modes at spatial infinity, the leading matrix is of
order $r^{2}$ and is nilpotent. So, in principle, one needs to apply a
procedure similar to the near-horizon analysis of axial modes, which will be
presented in full generality in the next section, and then diagonalise in turn
all subsequent orders. All this involves many steps which are straightforward
but rather tedious to describe.

To shorten our discussion, we provide directly the transformation that
combines all these intermediate steps, given by

$Y=P\tilde{Y}\,,\qquad P=\begin{pmatrix}{\cal S}+{\cal T}&{\cal S}-{\cal T}\\\
{\cal U}-{\cal V}&{\cal U}+{\cal V}\end{pmatrix},$ (116)

with the functions

$\displaystyle{\cal S}(r)\equiv\frac{i(r-r_{s})((2\lambda-3)r_{s}+4\lambda
r)}{4\lambda r}+\frac{i\lambda}{2r\omega^{2}}\,,\quad{\cal
T}(r)\equiv\frac{(1-2\lambda)r_{s}+2(1+2\lambda)r}{4r\omega}\,,$
$\displaystyle{\cal U}(r)\equiv
r^{2}+\frac{2\lambda-3}{4\lambda}r_{s}r\,,\quad{\cal
V}(r)\equiv\frac{ir}{2\omega}\,.$ (117)

This leads to the new system

$\derivative{{\tilde{Y}}}{r}=\tilde{M}(r)\tilde{Y}\,,\qquad\tilde{M}(r)=\begin{pmatrix}i\omega&0\\\
0&-i\omega\\\ \end{pmatrix}+\begin{pmatrix}-1+i\omega r_{s}&0\\\ 0&-1-i\omega
r_{s}\\\ \end{pmatrix}\frac{1}{r}+\mathcal{O}\left(\frac{1}{r^{2}}\right),$
(118)

which is diagonal and whose solution is

$\tilde{Y}(r)\;=\;\left(\begin{array}[]{c}c_{-}\,e^{-i\omega r}r^{-1-i\omega
r_{s}}\\\ c_{+}\,e^{+i\omega r}r^{-1+i\omega
r_{s}}\end{array}\right)\left(1+{\cal
O}\left({1}/{r}\right)\right)=\frac{1}{r}\left(\begin{array}[]{c}c_{-}\,e^{-i\omega
r_{*}}\\\ c_{+}\,e^{+i\omega r_{*}}\end{array}\right)\left(1+{\cal
O}\left({1}/{r}\right)\right)\,.$ (119)

This result is very similar to that obtained for axial perturbations (75),
even though the asymptotic expansion of the matrix $M$ is rather different. In
terms of the original functions, we find

$\displaystyle K(r)=\frac{i}{\omega}H_{1}(r)=i(c_{-}\,e^{-i\omega
r_{*}}+c_{+}\,e^{+i\omega r_{*}})\left(1+{\cal
O}\left({1}/{r}\right)\right)\,,$ (120)

which agree with (33) (with $c_{-}=-\omega{\cal B}_{\infty}$ and
$c_{+}=\omega{\cal A}_{\infty}$).

####  III.3.2 Asymptotic analysis at the black hole horizon

We finally turn to the near-horizon behaviour of polar modes. The expansion of
the matrix (114) in terms of the small parameter $\varepsilon\equiv r-r_{s}$
yields

$\displaystyle
M(\varepsilon)=\frac{1}{\varepsilon^{2}}M_{2}+\frac{1}{\varepsilon}M_{1}+M_{0}+{\cal
O}(\varepsilon),\quad M_{2}=\begin{pmatrix}0&0\\\
\gamma_{2}&0\end{pmatrix},\;M_{1}=\begin{pmatrix}\alpha_{1}&0\\\
\gamma_{1}&\delta_{1}\end{pmatrix},\;M_{0}=\begin{pmatrix}\alpha_{0}&\beta_{0}\\\
\gamma_{0}&\delta_{0}\end{pmatrix},$ (121)

where only a few of the coefficients $\alpha_{I}$, $\beta_{I}$ and
$\gamma_{I}$ will be needed explicitly.

Once more, the dominant $M_{2}$ is a nilpotent matrix and, as in the axial
case, we use the transformation

$\displaystyle Y=P_{(1)}Y^{(1)}\quad\text{with}\quad
P_{(1)}(\varepsilon)\equiv\left(\begin{array}[]{cc}1&0\\\
0&1/\varepsilon\end{array}\right)\,,$ (124)

which gives the new system

$\displaystyle\derivative{{Y}^{(1)}}{\varepsilon}={M}^{(1)}{Y}^{(1)}\,,\quad{M}^{(1)}(\varepsilon)=\frac{1}{\varepsilon}\left(\begin{array}[]{cc}\alpha_{1}&\beta_{0}\\\
\gamma_{2}&1+\delta_{1}\end{array}\right)+{\cal O}(1)\,,$ (127)

with the coefficients

$\displaystyle\alpha_{1}=-(1+\delta_{1})=\frac{1+\lambda-2r_{s}^{2}\omega^{2}}{3+2\lambda}\,,\quad\beta_{0}=\frac{2i}{r_{s}^{2}}\frac{(\lambda+1)^{2}+r_{s}^{2}\omega^{2}}{3+2\lambda}\,,\quad\gamma_{2}=\frac{ir_{s}^{2}}{2}\frac{1+4r_{s}^{2}\omega^{2}}{3+2\lambda}\,.$
(128)

The leading matrix can now be diagonalised via the transformation

$\displaystyle Y^{(1)}=P_{(2)}Y^{(2)}\,,\quad\text{with}\quad
P_{(2)}=\begin{pmatrix}\alpha-\beta&\alpha+\beta\\\
1&1\end{pmatrix}\quad\text{and}\quad\alpha=\frac{\alpha_{1}}{\gamma_{2}}\,,\quad\beta=\frac{i\omega
r_{s}}{\gamma_{2}}\,.$ (129)

leading to the system

$\derivative{Y^{(2)}}{\varepsilon}=M^{(2)}Y^{(2)}\,,\quad
M^{(2)}=\frac{1}{\varepsilon}\begin{pmatrix}-i\omega r_{s}&0\\\ 0&i\omega
r_{s}\end{pmatrix}+\mathcal{O}(1)\,.$ (130)

Note that this expression is extremely simple and does not involve $\lambda$,
as expected, even though it appears explicitly in $M^{(1)}$. We obtain
immediately the asymptotic behaviour of $X^{(2)}$ near the horizon

$Y^{(2)}(\varepsilon)=(1+\mathcal{O}(\varepsilon))\begin{pmatrix}c_{-}\,e^{-i\omega
r_{*}}\\\ c_{+}\,e^{+i\omega r_{*}}\end{pmatrix}\,,$ (131)

which reproduces the same result as for the axial mode (108). In terms of the
original gravitational functions $H_{1}(r)$ and $K(r)$, using the
transformation $Y=P_{(1)}P_{(2)}Y^{(2)}$, we recover the result (35), with

$c_{+}=\frac{i}{2}r_{s}(1-2i\omega r_{s}){\cal A}_{\text{hor}}\qquad
c_{-}=-\frac{i}{2}{r_{s}(1+2i\omega r_{s})}{\cal B}_{\text{hor}}\,.$ (132)

This completes our study of all asymptotic behaviours of Schwarzschild
perturbations, demonstrating that one can recover the standard results
directly from the linearised Einstein’s equations, without resorting to the
Schrödinger-like reformulation of the system.

###  III.4 Quasi-Normal modes

In this section, we show how the Schwarzschild quasinormal modes can be
recovered numerically, using directly the first-order system instead of the
Schrödinger equation (21). We restrict ourselves to the polar modes and
consider the system (113-114). The computation of the axial quasinormal modes
would be completely similar.

By definition of the quasi-normal modes, we impose that the solutions are
outgoing at spatial infinity and ingoing at the horizon, which means, using
the results of subsection III.2, that the two components of the vector $Y$
satisfy

$\displaystyle Y_{1}(r)\equiv K(r)$ $\displaystyle=$ $\displaystyle
K_{\infty}(r)e^{i\omega r_{*}}=\tilde{K}_{\infty}(r)e^{i\omega r}r^{i\omega
r_{s}}$ (133) $\displaystyle=$ $\displaystyle K_{\text{h}}(r)e^{-i\omega
r_{*}}=\tilde{K}_{\text{h}}(r)(r-r_{s})^{-i\omega r_{s}}\,,$ (134)

where $K_{\infty}$ (and $\tilde{K}_{\infty}$) is finite at infinity while
$K_{\text{h}}$ (and $\tilde{K}_{\text{h}}$) is finite at the horizon, and also

$\displaystyle Y_{2}(r)\equiv H_{1}(r)/\omega$ $\displaystyle=$ $\displaystyle
H_{\infty}(r)re^{i\omega r_{*}}=\tilde{H}_{\infty}(r)e^{i\omega r}r^{1+i\omega
r_{s}}$ (135) $\displaystyle=$ $\displaystyle
H_{\text{h}}(r)\varepsilon^{-1}e^{-i\omega
r_{*}}=\tilde{H}_{\text{h}}(r)(r-r_{s})^{-1-i\omega r_{s}}\,,$ (136)

where again $H_{\infty}$ (and $\tilde{H}_{\infty}$) is finite at infinity
while $H_{\text{h}}$ (and $\tilde{H}_{\text{h}}$) is finite at the horizon.

Therefore, we look for solutions of (113-114) using the ansatz

$\displaystyle K(r)=e^{i\omega r}r^{i\omega
r_{s}}\left(\frac{r-r_{s}}{r}\right)^{-i\omega r_{s}}f_{K}(r)\,,\quad
H(r)=e^{i\omega r}r^{1+i\omega
r_{s}}\left(\frac{r-r_{s}}{r}\right)^{-1-i\omega r_{s}}f_{H}(r)\,,$ (137)

where the functions $f_{K}$ and $f_{H}$ are supposed to be finite (hence
bounded) both at the horizon and at spatial infinity, in agreement with the
required boundary conditions. Furthermore, we introduce the new variable

$u=\frac{2r_{s}}{r}-1\,,$ (138)

so that the black hole horizon is located at $u=1$ and spatial infinity at
$u=-1$. Each function entering in the equations (137) is now treated as a
function of $u$ and the system of equations (113-114) can be expressed in the
form

$\displaystyle{\cal P}_{11}(u)f_{K}(u)+{\cal P}_{12}(u)f_{H}(u)+{\cal
Q}_{1}(u)f_{K}^{\prime}(u)=0\,,$ (139) $\displaystyle{\cal
P}_{21}(u)f_{K}(u)+{\cal P}_{22}(u)f_{H}(u)+{\cal
Q}_{2}(u)f_{H}^{\prime}(u)=0\,,$

where a prime denotes here a derivative with respect to $u$, and the functions
${\cal P}_{ij}$ and ${\cal Q}_{i}$ are polynomials in $u$. This is possible
because the matrix $M$ given in (114) contains only rational fractions of $r$.

In order to solve the system (139) numerically, we adapt the spectral method
presented in Jansen:2017oag  and we decompose $f_{K}(u)$ and $f_{H}(u)$ onto a
basis of Chebyshev polynomials. The facts that the functions ${\cal P}_{ij}$
and the ${\cal Q}_{i}$ are polynomials (hence $C^{\infty}$-functions) and the
Chebyshev polynomials are bounded at the boundaries ensure the boundedness of
$f_{K}(u)$ and $f_{H}(u)$ which is sufficient to enforce the required boundary
conditions. This is called a “behavioural” boundary condition Boyd1989 .

Then, any smooth and continuous complex-valued function $g(u)$ defined on the
interval $[-1,1]$ can be written as an infinite sum of Chebyshev polynomials
$T_{n}(u)$,

$g(u)=\sum_{n=0}^{\infty}g_{n}T_{n}(u)\,,$ (140)

where $g_{n}$ are complex coefficients. We can approximate the function $g$ by
truncating this series at a given order $N$, the approximation getting better
as $N$ is increased. Hence, we decompose the two functions $f_{K}$ and $f_{H}$
as follows,

$\displaystyle f_{K}(u)\approx\sum_{n=0}^{N}\alpha_{n}T_{n}(u)\,,\qquad
f_{H}(u)\approx\sum_{n=0}^{N}\beta_{n}T_{n}(u)\,,$ (141)

where $\alpha_{n}$ and $\beta_{n}$ are complex coefficients. Notice that the
symbol $\approx$ means that we truncated the series at an order $N$, then the
equality is not exact.

The next step is to express the differential system (139) as a linear system
for the coefficients $\alpha_{n}$ and $\beta_{n}$, which is always possible
due to fundamental relations satisfied by Chebyshev polynomials555The
Chebyshev polynomials satisfy the properties  $\displaystyle
T^{\prime}_{n}(u)=\sum_{m-n=2k+1}\frac{2m}{1+\delta_{n0}}\,T_{m}(u)\,,\quad(uT_{n})(u)=\sum_{m}\frac{1}{2}((1+\delta_{n,1})\delta_{n-1,m}+\delta_{n+1,m})\,T_{m}(u)\,,$
(142) where $\delta_{m,n}$ is the Kronecker symbol and $k\in\mathbb{N}$ in the
first sum. . As a consequence, the differential system (139) can be recast as
the following system of algebraic equations

$M_{N}(\omega)V_{N}(\alpha_{n},\beta_{n})=0\,,$ (143)

where $M_{N}$ is a $2(N+1)\times 2(N+1)$ matrix whose expansion in powers of
$\omega$ reads

$M_{N}(\omega)=M_{N[0]}+M_{N[1]}\omega+M_{N[2]}\omega^{2}\,,$ (144)

while the $2(N+1)$-dimensional vector $V_{N}(\alpha_{n},\beta_{n})$ is such
that

${}^{T}V_{N}(\alpha_{n},\beta_{n})\equiv\begin{pmatrix}\alpha_{0},&\cdots,&\alpha_{N},&\beta_{0},&\cdots,&\beta_{N}\end{pmatrix}\,.$
(145)

Following Jansen:2017oag , we can reformulate this system as

$\tilde{M}_{N}(\omega)\tilde{V}_{N}(\alpha_{n},\beta_{n})=0\,,$ (146)

where the matrix $\tilde{M}_{N}$ is now of dimension $4(N+1)$ and defined by

$\tilde{M}_{N}=\tilde{M}_{N[0]}+\tilde{M}_{N[1]}\omega\quad\text{and}\quad\tilde{M}_{N[0]}=\begin{pmatrix}M_{N[0]}&M_{N[1]}\\\
0&I\end{pmatrix}\,,\quad\tilde{M}_{N[1]}=\begin{pmatrix}0&M_{N[2]}\\\
-I&0\end{pmatrix}\,.$ (147)

Finding the values of $\omega$ such that the system (146) is nontrivial is
called a _generalised eigenvalues problem_ and can be done by a numerical
engine such as Mathematica or Scipy. In practice, we have computed the
eigenvalues for two different values of $N$ and kept only the ones that are
present in the two calculations, which allows us to get rid of spurious modes.
We plot the quasinormal modes detected in Fig.(2). They coincide with the
well-known first quasi-normal modes of Schwarzschild.

Figure 2: Quasinormal modes numerically found by Mathematica for $r_{s}=1$ and
$\ell=2$ ($\lambda=2$). The blue dots are generalised eigenvalues for $N=30$,
the orange dots generalised eigenvalues for $N=60$, and the red dots are the
modes detected (eigenvalues that change by a factor $10^{-3}$ or less). We
observe a symmetry with respect to the imaginary axis. The positions of the
first modes are $\omega_{0}=\pm 0.747-0.178i$, $\omega_{1}=\pm 0.693-0.548i$
and $\omega_{2}=\pm 0.602-0.957i$.

##  IV General analysis

As we have seen in the previous section, it is possible to compute the quasi-
normal modes of black holes in general relativity without reformulating the
linearised Einstein equations in terms of a Schrödinger-like equation. The
advantage of this method is that it can be straightforwardly generalised to
the study of black holes in theories of modified gravity where it might be
difficult or impossible to reduce the linearised equations to a Schrödinger-
like form.

In this section, we present a systematic algorithm for a generic first-order
system of the form (42), which has been developed in the mathematics
literature, first in wasow_asymptotic_1965  and more recently in
balser_computation_1999 ; pflugel_root-free_2019 ; barkatou_algorithm_1999 ;
barkatou_algorithm_1997 ; abdelaziz_hilali_solutions_1987 . The various steps
of the algorithm presented in this section are summarised in the flowchart
diagram depicted in Appendix C.

###  IV.1 Asymptotic solution: overview

We consider a general system of first-order ordinary differential equations of
the form

$\derivative{Y}{z}=M(z)Y\,,$ (148)

where $Y$ is a $n$-dimensional column vector, $M$ an $n\times n$-dimensional
matrix and $z$ a real variable defined in some interval. In the following, we
will consider only the asymptotic behaviour when $z\rightarrow+\infty$, but it
is straightforward to extend the algorithm near a finite value $z_{0}$ where
the system is singular, by a suitable change of the variable $z$.

We then assume that one can expand $M$ in powers of $z$, up to some order
(depending on the required precision of the asymptotic expansion) as follows,

$M(z)=M_{r}z^{r}+\dots+M_{0}+\dots M_{r-f}z^{r-f}+{\cal
O}(z^{r-f-1})=z^{r}\sum_{k=0}^{f}M_{r-k}z^{-k}+{\cal O}(z^{r-f-1})\,,$ (149)

where the integer $r$ is called the Poincaré rank of the system, and the
$M_{i}$ are $z$-independent matrices. In most cases666Note that, in some
cases, the variable $z$ in the expression (150) differs from the variable $z$
in the original system (148), because a change of variable is necessary, as
will be discussed around Eq. (167). Morever, the special case where
$M(z)=M_{-1}/z+{\cal O}(z^{-2})$ with $M_{-1}$ nilpotent leads to a $\ln z$
behaviour at large $z$, as discussed at the end of subsection IV.3., the
general solution to the system (148) admits an asymptotic expansion of the
form wasow_asymptotic_1965

$Y(z)=e^{\mathbf{\Upsilon}(z)}r^{\mathbf{\Delta}}\,\mathbf{F}(z)Y_{0}\,,$
(150)

where $Y_{0}$ is a constant vector, corresponding to $n$ integration
constants, $\mathbf{\Upsilon}$ is a diagonal matrix whose coefficients are
polynomials of degree at most $r+1$, $\mathbf{\Delta}$ is a constant diagonal
matrix and $\mathbf{F}(z)$ is a matrix which is regular at infinity.

The goal of the algorithm presented below is to determine explicitly the
expression (150) up to some irrelevant sub-leading terms. As we have already
seen in the previous section, the guiding principle in order to obtain this
expression is to fully diagonalise the differential system, up to the
appropriate order, by using iteratively transformations of the vector $Y$ into
a new vector $\tilde{Y}$, of the form

$Y(z)=P(z)\tilde{Y}(z)\,,$

where $P$ is an invertible matrix. The system (148) is then transformed into a
new but equivalent differential system, given by

$\displaystyle\derivative{\tilde{Y}}{z}=\tilde{M}(z)\tilde{Y},\qquad\tilde{M}(z)\equiv
P^{-1}MP-P^{-1}\derivative{P}{z}\,.$ (151)

The end point of this procedure is a system where the matrix coefficients in
the expansion of the form (149) are diagonal at each order. It is then
immediate to integrate the system and to find the solution in the form (150),
as discussed in subsection III.1.

In the following subsections, we describe the algorithm step by step. We have
also inserted two subsections that contain examples chosen to illustrate some
of the finer points of the algorithm. The algorithm contains several branches,
depending on whether the leading term $M_{r}$ in the expansion of $M(z)$ is
diagonalisable or not.

###  IV.2 Case 1: the leading term is diagonalisable

The simplest situation is when the leading matrix $M_{r}$ is diagonalisable.
In this case, one first uses the transformation $Y=P_{(1)}{Y}^{(1)}$ where
$P_{(1)}$ is a constant matrix that diagonalises $M_{r}$, which gives the new
system

$\displaystyle\derivative{{Y}^{(1)}}{z}={M}^{(1)}{Y}^{(1)},\qquad{M}^{(1)}(z)=D_{r}z^{r}+{M}^{(1)}_{r-1}z^{r-1}+\cdots+{M}^{(1)}_{0}+{M}^{(1)}_{-1}\frac{1}{z}+{\cal
O}\left(\frac{1}{z^{2}}\right)\,,$ (152)

where the matrix $D_{r}$ is diagonal.

One then seeks to transform the next-to-leading matrix ${M}^{(1)}_{r-1}$ into
a diagonal matrix (if it is not already) without affecting the diagonal form
of the leading order. This can be accomplished with a new transformation

$Y^{(1)}=P_{(2)}{Y}^{(2)}\,,\qquad P_{(2)}(z)=I+\frac{1}{z}\Xi_{(2)}\,,$ (153)

where $\Xi_{(2)}$ is a constant matrix. Indeed, this yields the new system

$\displaystyle\derivative{{Y}^{(2)}}{z}={M}^{(2)}{Y}^{(1)},\quad{M}^{(2)}(z)=D_{r}z^{r}+D_{r-1}z^{r-1}+{M}^{(2)}_{r-2}z^{r-2}+\cdots+{M}^{(2)}_{-1}\frac{1}{z}+{\cal
O}\left(\frac{1}{z^{2}}\right)\,,$ (154)

with

$\displaystyle D_{r-1}\;=\;M_{r-1}^{(1)}+[D_{r},\Xi^{(2)}]\,,$ (155)

which is imposed to be diagonal via an appropriate choice777To find $\Xi$ such
that the matrix $\tilde{D}=M+[D,\Xi]$ is diagonal, $M$ being arbitrary and $D$
diagonal, one notices that $[D,\Xi]_{ij}=(d_{i}-d_{j})\Xi_{ij}$ where $d_{i}$
are the eigenvalues of $D$. Consequently, $\tilde{D}$ is given by the diagonal
component of $M$ and the coefficients of $\Xi$ satisfy
$(d_{i}-d_{j})\Xi_{ij}+M_{ij}=0$, which always admit at least one solution for
each $\Xi_{ij}$. for $\Xi_{(2)}$. Furthermore, $D_{r-1}$ is the diagonal part
of $M_{r-1}^{(1)}$.

One can proceed similarly to “diagonalise” all the other terms, order by
order, until one gets a system of the form888Note that we could have proceeded
in a single step by introducing the new variable $\tilde{Y}$ defined by
$Y=P(z)\tilde{Y}$ with
$P(z)=P_{0}+\frac{1}{z}P_{1}+\cdots+\frac{1}{z^{r+1}}P_{r+1}$ and determining
the constant matrices $P_{j}$ so that $\tilde{M}(z)$ is equal to (156). The
calculation we have just done proves this is possible with
$\tilde{Y}=Y^{(r+2)}$.

$\displaystyle\derivative{{Y}^{(r+2)}}{z}={M}^{(r+2)}{Y}^{(r+2)}\,,\quad
M^{(r+2)}(z)=D_{r}z^{r}+\cdots+D_{0}+D_{-1}\frac{1}{z}+{\cal
O}\left(\frac{1}{z^{2}}\right)\,,$ (156)

where all matrices are diagonal up to order $1/z$. The system can then be
immediately integrated, to yield

$\displaystyle
Y^{{(r+2)}}(z)=e^{\mathbf{\Upsilon}(z)}z^{\mathbf{\Delta}}\,\mathbf{F}(z)\,Y_{0}\,,\quad\mathbf{\Delta}\equiv
D_{-1}\,,\quad\mathbf{\Upsilon}(z)\equiv
D_{r}\frac{z^{r+1}}{r+1}+\cdots+D_{0}z\,,$ (157)

where $Y_{0}$ is a constant vector.

The asymptotic expansion of the original vector $Y$ can be simply deduced from
the combined transformations, i.e.

$\displaystyle Y=P_{(1)}P_{(2)}\cdots P_{(r+2)}Y^{(r+2)}\,.$ (158)

Since the $P_{(j)}$ are polynomials of $1/z$, $Y$ has exactly the same
exponential behaviour (in its asymptotic expansion) as $Y^{(r+2)}$.

###  IV.3 Case 2: the leading term is non-diagonalisable, similar to a single-
block Jordan matrix

Solving asymptotically a system where the dominant term $M_{r}$ is not
diagonalisable is more challenging. The basic idea consists in finding a
transformation where the leading term of the new matrix becomes
diagonalisable. This can be done by reducing progressively the Poincaré rank
of the system until the leading term is diagonalisable, in which case the
procedure of the previous subsection becomes applicable. If the leading term
never gets diagonalisable down to the rank $r=-1$, then the general formula
(150) for the asymptotic expansion is not valid but the system can
nevertheless be integrated explicitly.

The reduction of the Poincaré rank together with the diagonalisation of the
leading term is done in different steps, which we now describe, first when the
leading term is similar to a Jordan matrix with a single block. The case of a
Jordan matrix with several blocks will be discussed later, in subsection IV.5.

Step 1. Transformation to a Jordan block

Starting from the asymptotic expansion (149) of the matrix $M$, we use the
transformation $X=P_{(1)}X^{(1)}$ to write
$M^{(1)}_{r}=P_{(1)}^{-1}M_{r}P_{(1)}$ in a Jordan canonical form (although
with a lower triangular matrix). We assume here that $M^{(1)}_{r}$ contains a
single (lower triangular) Jordan block with eigenvalue $\lambda$, i.e. of the
form

$\displaystyle M^{(1)}_{r}=\begin{pmatrix}\lambda&0&\cdots&&\\\
1&\lambda&0&\cdots&\\\ 0&1&\lambda&0&\cdots\\\
\vdots&&&\end{pmatrix}\equiv\lambda I+J(n)\,,$ (159)

where $J(n)$ has the property to be nilpotent (we recall that $n$ is the
dimension of the matrix).

Step 2. Transformation to a nilpotent matrix

We then apply the transformation

$\displaystyle Y^{(1)}=P_{(2)}{Y}^{(2)}\,,\qquad
P_{(2)}(z)\equiv\exp\left(\frac{\lambda}{r+1}z^{r+1}\right)I\,,$ (160)

which renders the leading term nilpotent999This follows from the relation
$\displaystyle P_{(2)}^{-1}\left(z^{r}(\lambda
I+J(n)\right)P_{(2)}-P_{(2)}^{-1}\derivative{P_{(2)}}{z}=z^{r}(M_{r}^{(1)}-\lambda
I)=z^{r}J(n)\,.$ (161)

$\displaystyle
M^{(2)}(z)=J(n)z^{r}+M^{(2)}_{r-1}z^{r-1}+\cdots+M^{(2)}_{0}+M^{(2)}_{-1}\frac{1}{z}+{\cal
O}\left(\frac{1}{z^{2}}\right)\,.$ (162)

Step 3. Normalisation and reduction of the Poincaré rank

The next step consists in reducing the Poincaré rank of the system by using
the transition matrix

$\displaystyle P(z)=D({n},z)\equiv\begin{pmatrix}1&0&0&\cdots&\cdots&0\\\
0&z&0&\cdots&\cdots&0\\\ 0&0&z^{2}&0&\cdots&0\\\ \vdots&&&\ddots&&\vdots\\\
0&\cdots&&&\cdots&z^{n-1}\end{pmatrix}\,,$ (163)

which satisfies the useful property

$\displaystyle P^{-1}J(n)P=\frac{1}{z}J(n)\,.$ (164)

A transformation with the above $P$ will thus reduce the order of the leading
term $J(n)z^{r}$, but will also affect the sub-dominant terms in the expansion
(162) of $M^{(2)}$, in particular $M^{(2)}_{r-1}$ which could generate terms
whose order is higher than $r-1$ in the new matrix.

To avoid this situation, we need first to “normalise” the system, with the
transformation

$\displaystyle{P}_{(3)}(z)=I+\frac{1}{z}\Lambda_{(3)}\,,$ (165)

where $\Lambda_{(3)}$ is a constant matrix, chosen such that such that the
next-to-leading order matrix $M^{(3)}_{r-1}$ in the new matrix expansion
contains only zeros except possibly in the first row. Let us stress that this
transformation leaves the leading term of the expansion unchanged. The new
system associated with ${M}^{(3)}$ is said to be normalised.

One can then perform the transformation generated by the transition matrix

$P_{(4)}(z)=D({n},z)\,,$ (166)

which, in most cases, gives a reduced Poincaré rank. There are however a few
exceptions where the reduction does not work. These special cases require a
more general transformation, with a transition matrix of the form

$\displaystyle P_{(4)}(z)=D({n},z^{p/q})\,,\qquad(1\leq p\leq q\leq n)$ (167)

where $p$ and $q$ are co-prime integers. For example, when $n=4$, the possible
choices are $\\{1/4,1/3,1/2,2/3,3/4,1\\}$, where the last value corresponds to
the generic case (166). To identify the appropriate value of $p/q$, one must
test successively the possible values, in decreasing order, until the
transformation (167) effectively leads to a system with a lower Poincaré rank.
The corresponding value of $p/q$ is said to be “admissible”. In practice, this
can be understood as a change of variable101010In this case, the asymptotic
expansion of the solution may have an exponential behaviour where the argument
$Q(z)$ is not a polynomial of $z$ but rather a polynomial of $z^{1/q}$., $z$
being replaced by $u=z^{p/q}$.

Step 4. Diagonalisable or not diagonalisable?

The next step depends on the nature of the system $(Y^{(4)},M^{(4)})$, which
possesses a lower Poincaré rank than the initial system. If the leading term
of $M^{(4)}$ is diagonalisable, one proceeds as in subsection IV.2.

If $M^{(4)}$ is not diagonalisable, one needs to reduce again the Poincaré
rank of the system, unless one has already reached $r=-1$, in which case one
can jump directly to the next paragraph. Otherwise, one must distinguish the
following different cases.

  * •

If the leading term is similar to a single-block Jordan matrix and we took
$p/q=1$ in the previous step, we repeat the procedure of this subsection.

  * •

If the leading term is similar to a single-block Jordan matrix but we took
$p/q<1$ in the previous step, we discard the last step, and start again with
the normalised system $M^{(3)}$. However, this time, we normalise the system
up to second order: after having normalised $M_{-1}$, we repeat the procedure
with $z^{2}$ instead of $z$ in ${P}_{(3)}$ (165) and require that $M_{-2}$ has
a specific form. Details can be found in balser_computation_1999 . If
necessary, one can pursue the normalisation to higher orders.

  * •

If the Jordan canonical form of the leading term contains several blocks, we
go to subsection IV.5.

Eventually we obtain either a system with a diagonalisable leading term, which
can be solved following subsection IV.2, or a system of Poincaré rank $r=-1$
with a nilpotent leading term. In the latter case, the solution is equivalent
to a polynomial of $\ln z$ at large $z$. Indeed, a system of the form

$\displaystyle\derivative{Y}{z}=\frac{\mu_{0}}{z}\begin{pmatrix}0&0&\cdots&&\\\
1&0&0&\cdots&\\\ 0&1&0&0&\cdots\\\ \vdots&&&\end{pmatrix}Y\,,$ (168)

where $\mu_{0}$ is an arbitrary constant, is easily integrated. The components
$Y_{i}$ (for $1\leq i\leq n$) are obtained iteratively and are given by
$Y_{1}(z)=\xi_{1}$, $Y_{2}(z)=\xi_{1}\ln z+\xi_{2}$ and more generally,

$\displaystyle Y_{i}(z)=\sum_{j=1}^{i}\frac{\xi_{j}}{(i-j)!}(\mu_{0}\ln
z)^{i-j}\,,$ (169)

where the $\xi_{i}$ are $n$ constants of integration. All the components of
$Y$ are thus polynomials of $\ln z$ at large $z$.

###  IV.4 An example with a nilpotent leading term

Let us give a concrete example of the procedure used for systems with a
nilpotent leading term. We consider the two-dimensional system defined by

$\displaystyle\derivative{Y}{z}=M(z)Y\,,\qquad M(z)=\begin{pmatrix}0&1\\\
0&0\end{pmatrix}z^{2}+\begin{pmatrix}1&0\\\ 0&-1\end{pmatrix}\,,$ (170)

and let us determine its asymptotic solution at large $z$, following the
algorithm described above.

We first put the leading term in its lower triangular Jordan form:

$P_{(1)}=\begin{pmatrix}0&1\\\ 1&0\end{pmatrix}\quad\implies\quad
M^{(1)}(z)=\begin{pmatrix}0&0\\\ 1&0\end{pmatrix}z^{2}+\begin{pmatrix}-1&0\\\
0&1\end{pmatrix}\,.$ (171)

Since the leading term is already nilpotent, step 2 is irrelevant. Moreover,
the system is already normalised since the next-to-leading order term
vanishes.

We can thus move directly to the reduction of the order of the system and
consider the transformation of the form (163):

$P_{(2)}(z)=\begin{pmatrix}1&0\\\ 0&z\end{pmatrix}\quad\implies\quad
M^{(2)}(z)=\begin{pmatrix}0&0\\\ 1&0\end{pmatrix}z+\begin{pmatrix}-1&0\\\
0&1\end{pmatrix}\,.$ (172)

The order has been reduced but the leading term is still nilpotent. Since the
reduction was obtained via a transformation with $p/q=1$, we continue the
process by doing a new iteration of the algorithm. We first normalise the
system with a transformation of the form (165),

$P_{(3)}(z)=I+\frac{1}{z}\begin{pmatrix}0&-1\\\
0&0\end{pmatrix}\quad\implies\quad M^{(3)}(z)=\begin{pmatrix}0&0\\\
1&0\end{pmatrix}z+\begin{pmatrix}0&1\\\
0&-1\end{pmatrix}\frac{1}{z}+\begin{pmatrix}0&-2\\\
0&0\end{pmatrix}\frac{1}{z^{2}}\,,$ (173)

and again reduce the order of the system with the transformation

$P_{(4)}(z)=\begin{pmatrix}1&0\\\ 0&z\end{pmatrix}\quad\implies\quad
M^{(4)}(z)=\begin{pmatrix}0&1\\\ 1&0\end{pmatrix}+\begin{pmatrix}0&-2\\\
0&-2\end{pmatrix}\frac{1}{z}\,.$ (174)

The leading term is now diagonalisable. We diagonalise it explicitly, via

$P_{(5)}=\begin{pmatrix}-1&1\\\ 1&1\end{pmatrix}\quad\implies\quad
M^{(5)}(z)=\begin{pmatrix}-1&0\\\ 0&1\end{pmatrix}+\begin{pmatrix}0&0\\\
-2&-2\end{pmatrix}\frac{1}{z}\,,$ (175)

then we diagonalise the next-to-leading term, with a transformation of the
form (153),

$P_{(6)}(z)=\begin{pmatrix}1&0\\\ 1/z&1\end{pmatrix}\quad\implies\quad
M^{(6)}(z)=\begin{pmatrix}-1&0\\\ 0&1\end{pmatrix}+\begin{pmatrix}0&0\\\
0&-2\end{pmatrix}\frac{1}{z}+\mathcal{O}\left(\frac{1}{z^{2}}\right)\,.$ (176)

We have thus managed to fully diagonalise the system, which immediately gives
us the asymptotic solution

$Y^{(6)}(z)=\left(1+{\cal
O}\left({1}/{z}\right)\right)\begin{pmatrix}\exp(-z)&0\\\
0&\frac{1}{z^{2}}\exp(z)\end{pmatrix}Y_{0},\qquad
Y_{0}\equiv\begin{pmatrix}\xi_{1}\\\ \xi_{2}\end{pmatrix}\,,$ (177)

where $Y_{0}$ is a constant column vector. As a consequence, to obtain the
behaviour of $Y$ in the original system, we use the combined transformations

$Y=\left(\prod_{j=1}^{6}P_{(j)}\right)Y^{(6)},$ (178)

which implies

$Y(z)=\left(1+{\cal
O}\left({1}/{z}\right)\right)\begin{pmatrix}\xi_{1}\exp(-z)z^{2}+\xi_{2}\exp(z)\\\
-2\xi_{1}\exp(-z)\end{pmatrix}.$ (179)

For this particular example, it turns out that the original system (170) can
be solved exactly, with the solution

$Y(z)=\begin{pmatrix}\frac{1}{2}\xi_{1}\exp(-z)\left(1+2z+2z^{2}\right)+\xi_{2}\exp(z)\\\
-2\xi_{1}\exp(-z)\end{pmatrix}.$ (180)

One can thus check that the asymptotic solution (179) agrees with the
asymptotic behaviour of the exact solution.

###  IV.5 Case 3: $M_{r}$ is similar to a Jordan matrix with several blocks

We now briefly discuss (without entering into too many details, which can be
found in balser_computation_1999 ) the more general case where $M_{r}$ is
block diagonalisable and its canonical Jordan form admits several Jordan
blocks. The first two steps of subsection IV.3 still apply to this case and
one can find a transformation (with a constant matrix $P$) such that the new
system associated with $M^{(2)}$ (we use the same notation as in subsection
IV.3) has a block diagonal leading term $M^{(2)}_{r}$ with Jordan lower
triangular blocks, each block being either nilpotent or 1-dimensional:

$\displaystyle M^{(2)}_{r}=\begin{pmatrix}J(n_{1})&0&\cdots&&\\\
0&J(n_{2})&0&\cdots&\\\ \vdots&0&\ddots&0&\cdots\\\
&\vdots&0&\lambda_{1}&0&\cdots\\\ &&\vdots&0&\lambda_{2}&0\\\
&&&\vdots&0&\ddots\\\ \end{pmatrix}\,,\qquad
J(n)\equiv\begin{pmatrix}0&0&\cdots&&\\\ 1&0&0&\cdots&\\\ 0&1&0&0&\cdots\\\
\vdots&&&\end{pmatrix}.$ (181)

The Jordan form is chosen so that the blocks $J(n)$ are ordered by decreasing
size ($n_{1}\geq n_{2}\geq\cdots$). We will use this block structure as a
layout for the block structure of the other matrices that appear in the
expansion of $M^{(2)}$. And each block will be denoted by two indices, $(KL)$,
corresponding to a submatrix of dimensions $n_{K}\times n_{L}$.

The principle of the diagonalisation procedure is similar to what was done in
sections IV.2 and IV.3. However, it is now possible to have both
diagonalisable blocks and nilpotent blocks. Those must be dealt with
separately to get the full asymptotic behaviour of the system. In order to do
this, one can generalise the order-by-order procedure of subsection IV.2: this
is called the “Splitting Lemma” in balser_computation_1999 . It is not
detailed here, but can be understood by considering blocks instead of scalars
in the computations of subsection IV.2111111In the case where $M^{(2)}_{r}$
consists of a 2-block Jordan matrix, one would use a transformation of the
form  $P=\begin{pmatrix}I&\sum_{j=1}^{p}\Xi_{j}\,z^{-j}\\\
\sum_{j=1}^{p}\Lambda_{j}\,z^{-j}&I\end{pmatrix}\,,$ (182) where the $\Xi_{j}$
and $\Lambda_{j}$ are constant matrices. Such a transformation, which
generalises (153), enables us to transform each $M^{(2)}_{r-j}$ in the same
block diagonal form as $M^{(2)}_{r}$ with a convenient choice of $\Xi_{i}$ and
$\Lambda_{i}$. Therefore, the initial system gives two decoupled sub-systems
and, for each one, we proceed along the same lines as in the previous
section..

One can use this lemma to block diagonalise $M^{(2)}$, order by order : the
two global blocks considered will be the nilpotent part of $M^{(2)}_{r}$ and
its diagonalisable part. The latter can be dealt with using the procedure
given in subsection IV.2, while the former must be addressed using a
generalised version of the procedure given in subsection IV.3. We give here
more details about the last part and, in the rest of this section, assume
without loss of generality that $M^{(2)}_{r}$ contains only nilpotent blocks,
such that

$M^{(2)}_{r}=\begin{pmatrix}J(n_{1})&0&\cdots&&\\\ 0&J(n_{2})&0&\cdots&\\\
\vdots&0&J(n_{3})&0&\cdots\\\ &\vdots&&&\ddots\end{pmatrix}$ (183)

.

As a preliminary step, one must “normalise” this new system, via a
transformation

${P}_{(3)}=I+\frac{1}{z}\Lambda_{(3)}\,,$ (184)

where the constant matrix $\Lambda_{(3)}$ is chosen such that the $(KL)$
blocks of the resulting sub-leading order matrix ${M}^{(3)}_{r-1}$ has the
following structure:

  * -

either all rows are equal to zero except possibly the first one if $K\leq L$,

  * -

or all columns are equal to zero except possibly the last one if $K>L$.

Notice that the leading order term ${M}^{(3)}_{r}$ is unchanged by this
transformation while lower order terms ${M}^{(3)}_{r-j}$ are modified but
there are no conditions when $j\geq 2$. Furthermore, the transformation
$\Lambda_{(3)}$ such that the previous conditions are met can be shown to be
unique, if we require in addition that all its blocks $\Lambda_{(3)}^{KL}$
have zero last row if $K\leq L$ and zero first column if $K>L$.

We are now in the position to reduce the Poincaré rank of the system using a
transformation of the form of (163) or (167) for each block, i.e.

$P^{(4)}(z)=\begin{pmatrix}D(n_{1},z^{p/q})&0&\cdots&&\\\
0&D(n_{2},z^{p/q})&0&\cdots&\\\ \vdots&0&D(n_{3},z^{p/q})&0&\cdots\\\
&\vdots&&&\ddots\end{pmatrix}\,,$ (185)

where the matrices $D({n},z)$ have been defined in (163); $p$ and $q$ are co-
prime integers (with $1\leq p\leq q\leq n_{1}$), or equal in the case $p/q=1$.

###  IV.6 A higher dimensional example with $p/q\neq 1$

We now present a higher dimensional ($n=5$) example, adapted from
pflugel_root-free_2019 , where the dominant term in the asymptotic expansion
of the matrix $M$ has a non trivial canonical Jordan form with two Jordan
blocks. The matrix $M(z)$ is given by

$M(z)=\begin{pmatrix}0&z^{3}&-z&1&2z\\\ -z^{2}&z&0&-z&0\\\ z&1&0&z^{3}&1\\\
1&-z&1&z&z^{3}\\\ z&0&-3z&0&-1\\\ \end{pmatrix}\equiv
M_{3}z^{3}+M_{2}z^{2}+M_{1}z+M_{0}\,,$ (186)

where the leading term $M_{3}$ is nilpotent and has a 2-block Jordan
structure.

We perform a first transformation $Y=P_{(1)}Y^{(1)}$ so that the leading term
has now the following Jordan (lower triangular) canonical form (the matrix
$P_{(1)}$ can easily been deduced):

$M^{(1)}(z)=\begin{pmatrix}-1&0&-3z&0&z\\\ z^{3}&z&1&-z&1\\\ 1&z^{3}&0&1&z\\\
0&-z&0&z&-z^{2}\\\ 2z&1&-z&z^{3}&0\\\
\end{pmatrix}\,\qquad\Longrightarrow\quad
M^{(1)}_{3}=\begin{pmatrix}0&0&0&0&0\\\ 1&0&0&0&0\\\ 0&1&0&0&0\\\ 0&0&0&0&0\\\
0&0&0&1&0\\\ \end{pmatrix}\,.$ (187)

The block structure of $M^{(1)}_{3}$ defines the layout that we will be using
to compute the asymptotic expansion of the solution.

We notice that the next-to-leading term $M^{(1)}_{2}$ in the expansion of
$M^{(1)}$ is already normalised. Therefore, we can immediately try to reduce
the order of the system thanks to a new transformation
$Y^{(1)}=P_{(2)}Y^{(2)}$,

$P_{(2)}=\begin{pmatrix}1&0&0&0&0\\\ 0&z&0&0&0\\\ 0&0&z^{2}&0&0\\\
0&0&0&1&0\\\ 0&0&0&0&z\\\ \end{pmatrix}\quad\implies\quad
M^{(2)}=\begin{pmatrix}-1&0&-3z^{3}&0&z^{2}\\\ z^{2}&z-\frac{1}{z}&z&-1&1\\\
\frac{1}{z^{2}}&z^{2}&-\frac{2}{z}&\frac{1}{z^{2}}&1\\\ 0&-z^{2}&0&z&-z^{3}\\\
2&1&-z^{2}&z^{2}&-\frac{1}{z}\\\ \end{pmatrix}\,.$ (188)

However, we immediately see that the order of the system has not diminished.
This example falls in the cases where we need to change the variable $z$ or,
equivalently, i.e. to make a transformation of the form (167) for each Jordan
block, We must therefore cancel the previous transformation (188) and instead
consider $Y^{(1)}=\tilde{P}_{(2)}\tilde{Y}^{(2)}$, with

$\displaystyle\tilde{P}_{(2)}(z)=\begin{pmatrix}1&0&0&0&0\\\
0&z^{p/q}&0&0&0\\\ 0&0&z^{2p/q}&0&0\\\ 0&0&0&1&0\\\
0&0&0&0&z^{p/q}\end{pmatrix}\,.$ (189)

Following the method described below Eq. (167), we note that the largest
Jordan block is of dimension 3, therefore we should take 2 co-prime integers
between $1$ and $3$ for $p$ and $q$ with $p\leq q$. The possible choices for
the ratio $p/q$ belong to the set $\\{1/3,1/2,2/3\\}$, since $p/q=1$ does not
work. The largest value is $p/q=2/3$, which gives for the matrix $M^{(3)}$ the
expression

$\begin{pmatrix}-1&0&-3z^{7/3}&0&z^{5/3}\\\
z^{7/3}&z-\frac{2}{3z}&z^{2/3}&-{z}^{1/3}&1\\\
\frac{1}{z^{4/3}}&z^{7/3}&-\frac{4}{3z}&\frac{1}{z^{4/3}}&{z}^{1/3}\\\
0&-z^{5/3}&0&z&-z^{8/3}\\\ 2{z}^{1/3}&1&-z^{5/3}&z^{7/3}&-\frac{2}{3z}\\\
\end{pmatrix}.$ (190)

We observe that the subdiagonal terms have order $7/3$. To keep this value of
$p/q$, we must make sure that no other term behaves like $z^{\alpha}$ with
$\alpha>7/3$. However in this case there is a $z^{8/3}$ term. Therefore, the
value $2/3$ is not admissible and we have to consider the next possible choice
which is $p/q=1/2$. Such a change of variable leads to the matrix

$\tilde{M}^{(2)}=\begin{pmatrix}-1&0&-3z^{2}&0&z^{3/2}\\\
z^{5/2}&z-\frac{1}{2z}&\sqrt{z}&-\sqrt{z}&1\\\
\frac{1}{z}&z^{5/2}&-\frac{1}{z}&\frac{1}{z}&\sqrt{z}\\\
0&-z^{3/2}&0&z&-z^{5/2}\\\ 2\sqrt{z}&1&-z^{3/2}&z^{5/2}&-\frac{1}{2z}\\\
\end{pmatrix}.$ (191)

Now, it verifies the requirements and we thus keep the value $p/q=1/2$ and
continue the process.

The previous change of variable leads to a differential system where the
coefficients of $M^{(3)}$ are non-integer powers functions of $z$. To apply
the algorithm, we have to make a change of coordinate so that the system
involves only integer powers of $z$. This can easily be done by introducing
the new coordinate $u$ defined by $z=u^{2}$. As a consequence, the new
differential system is now given by

$\derivative{Y^{(3)}}{u}=M^{(3)}(u)Y^{(3)}\,,\qquad
M^{(3)}(u)=\begin{pmatrix}-2u&0&-6u^{5}&0&2u^{4}\\\
2u^{6}&\frac{2u^{4}-1}{u}&2u^{2}&-2u^{2}&2u\\\
\frac{2}{u}&2u^{6}&-\frac{2}{u}&\frac{2}{u}&2u^{2}\\\
0&-2u^{4}&0&2u^{3}&-2u^{6}\\\ 4u^{2}&2u&-2u^{4}&2u^{6}&-\frac{1}{u}\\\
\end{pmatrix},$ (192)

where $Y^{(3)}(u)\equiv\tilde{Y}^{(2)}(z)$ and $M^{(3)}(u)\equiv
2u\tilde{M}^{(2)}(z)$ with $z=u^{2}$. As the leading term is not nilpotent, we
keep the value of $p/q$. If it had been nilpotent, we would have had to go
back one step and normalise up to the next order.

We can continue the algorithm with this new system: we will to do a new change
of variables, reduce the order, and decouple the system… We will not present
more steps as the rest of the computations is similar to what was done here
and in previous sections. Nonetheless, for the sake of completeness, we give
the final result. We show that, after enough steps of the algorithm, the
initial system can be equivalently reformulated as

$\displaystyle\derivative{Y^{(4)}}{w}=M^{(4)}(w)Y^{(4)}\,,$ (193)

where $w=z^{1/6}$ and $M^{(4)}(w)$ is the following diagonal matrix

$\displaystyle M^{(4)}(w)=\text{Diag}[$ $\displaystyle
3^{4/3}(1-i\sqrt{3})w^{19}+2w^{11},-2\times 3^{4/3}w^{19}+2w^{11},$ (194)
$\displaystyle
3^{4/3}(1+i\sqrt{3})w^{19}+2w^{11},6iw^{20}+3w^{11},-6iw^{20}+3w^{11}]+\mathcal{O}(w^{9})\,,$

up to order $\mathcal{O}(w^{9})$. Integrating such a system is immediate and
yields the leading orders of the asymptotic expansion of $Y^{(4)}$ from which
we can extract the asymptotic expansion of the original variable $Y$.

##  V Conclusion

In this work, we have studied the asymptotic behaviours, both at spatial
infinity and near the horizon, of the linear perturbations about Schwarzschild
black holes. Instead of following the traditional approach that consists in
rewriting the equations of motion in the form of a stationary Schrödinger-like
equation, which is second-order with respect to the radial coordinate, we have
worked directly with the first-order equations of motion (in the frequency
domain). For this direct approach to the asymptotic behaviour, we have used an
algorithm that has been developed in several recent articles published in
mathematical journals.

The principle of this algorithm is to transform the differential system, via
successive changes of functions, until it can be written in an explicitly
diagonal form, up to the required order (in the small parameter characterising
the asymptotic regime). This procedure automatically provides the combination
of the metric perturbations that encapsulates the physical degree of freedom
in this asymptotic region and enables one to separate the ingoing and outgoing
physical modes. Although we have worked in the standard Regge-Wheeler gauge,
the same approach would work similarly for any other gauge choice.

Beyond its application to the perturbations of black holes, this systematic
approach to the asymptotic behaviour could be very useful for many problems in
physics. This is why we have devoted the last part of this paper to a
pedagogical presentation of the algorithm, with a few illustrative examples.

For black holes, the knowledge of the asymptotic behaviour of the
perturbations is an indispensable first step in the determination of the
quasi-normal modes. Indeed, these modes are characterised by the following
boundary conditions: a purely outgoing behaviour at spatial infinity and
purely ingoing behaviour at the horizon. Imposing these boundary conditions,
we have shown that the quasi-normal modes can be found numerically, without
resorting to the Schrödinger-like formulation, providing an alternative
approach to the standard method.

This novel approach could be especially useful in the context of generalised
black hole solutions, for instance in modified gravity theories, where the
equations of motion for the perturbations are different and extra fields can
be present. In a companion paper, we have applied the same algorithm for a few
black holes solutions within scalar-tensor theories that belong to the most
general known family: DHOST (Degenerate Higher-Order Scalar-Tensor) theories.
For the polar modes, we were not able to rewrite the equations of motion in
the form of a Schrödinger-like system and this new approach thus becomes
invaluable in order to identify the corresponding quasi-normal modes. The same
method could be applied to the study of other types of black holes, or even
completely different physical systems.

###### Acknowledgements.

We would like to thank Oleg Lisovyi for instructive discussions and guidance
on the mathematics literature, as well as Emanuele Berti and Vitor Cardoso for
very useful correspondence on quasi-normal modes. We have used Mathematica for
many of the calculations involved in this work. KN acknowledges the support
from the CNRS grant 80PRIME and thanks the Laboratory of Physics at the ENS in
Paris for its hospitality.

##  Appendix A Gauge transformations

For completeness, we summarise in this Appendix the gauge fixing procedure for
polar and axial perturbations about a Schwarzschild black hole in General
Relativity, as originally discussed in Regge:1957td  and Zerilli:1970se .

Due to the invariance of the theory under space-time diffeomorphisms, the
metric perturbations are not completely determined $h_{\mu\nu}$. Indeed, any
infinitesimal change of coordinates $x^{\mu}\rightarrow x^{\mu}+\xi^{\mu}$
induces the transformation

$\displaystyle h_{\mu\nu}\rightarrow
h_{\mu\nu}+\nabla_{\mu}\xi_{\nu}+\nabla_{\nu}\xi_{\mu}\,$ (195)

at the linear level. These transformations can be “projected” in the axial or
polar sectors, which we examine in turn.

###  A.1 Axial perturbations: Regge-Wheeler gauge

Before gauge fixing, axial perturbations are parametrised by three families of
functions $h_{0}^{\ell m}$, $h_{1}^{\ell m}$ and $h_{2}^{\ell m}$ of the
variables $(r,t)$, according to

$\displaystyle h_{t\theta}=\frac{1}{\sin\theta}\sum_{\ell,m}h_{0}^{\ell
m}(t,r)\partial_{\varphi}{Y_{\ell m}}(\theta,\varphi),\qquad
h_{t\varphi}=-\sin\theta\sum_{\ell,m}h_{0}^{\ell
m}(t,r)\partial_{\theta}{Y_{\ell m}}(\theta,\varphi),$ $\displaystyle
h_{r\theta}=\frac{1}{\sin\theta}\sum_{\ell,m}h_{1}^{\ell
m}(t,r)\partial_{\varphi}{Y_{\ell m}}(\theta,\varphi),\qquad
h_{r\varphi}=-\sin\theta\sum_{\ell,m}h_{1}^{\ell
m}(t,r)\partial_{\theta}{Y_{\ell m}}(\theta,\varphi),$ (196) $\displaystyle
h_{ab}=\sin\theta\sum_{\ell,m}h_{2}^{\ell
m}(t,r)\epsilon_{c(a}D^{c}\partial_{b)}Y_{\ell m}(\theta,\varphi)\,,$

where, in the last equation, the indices $a$ and $b$ belong to the set
$\\{\theta,\varphi\\}$, $\epsilon_{ab}$ is the totally antisymmetric symbol
such that $\epsilon_{\theta\varphi}=+1$ and $D_{a}$ is the 2-dimensional
covariant derivative associated with the metric of the 2-sphere
$d\theta^{2}+\sin^{2}\theta\,d\varphi^{2}$. More explicitely, the angular
components of the metric can be written

$\displaystyle h_{\theta\theta}$
$\displaystyle=\sum_{\ell,m}\frac{1}{\sin\theta}h_{2}^{\ell
m}(t,r)\left(\partial_{\theta}\partial_{\varphi}-\operatorname{cotan}\theta\,\partial_{\varphi}\right)Y_{\ell
m}(\theta,\varphi)\,,$ (197) $\displaystyle h_{\theta\varphi}$
$\displaystyle=h_{\varphi\theta}=-\sum_{\ell,m}\sin\theta\,h_{2}^{\ell
m}(t,r)\left(\frac{\ell(\ell+1)}{2}+\partial_{\theta}^{2}\right)Y_{\ell
m}(\theta,\varphi)\,,$ (198) $\displaystyle h_{\varphi\varphi}$
$\displaystyle=-\sum_{\ell,m}h_{2}^{\ell
m}(t,r)\sin\theta\left(\partial_{\theta}\partial_{\varphi}-\operatorname{cotan}\theta\,\partial_{\varphi}\right)Y_{\ell
m}(\theta,\varphi)\,.$ (199)

All the other components of the axial perturbations vanish.

In the axial sector, the nonzero components of the generator $\xi^{\mu}$ that
preserves the odd parity of the perturbations can be decomposed into spherical
harmonics as follows,

$\displaystyle\xi_{\theta}=\sum_{\ell,m}\xi^{\ell
m}(t,r)\partial_{\theta}Y_{\ell,m}(\theta,\varphi)\,,\qquad\xi_{\varphi}=\sum_{\ell,m}\xi^{\ell
m}(t,r)\partial_{\varphi}Y_{\ell,m}(\theta,\varphi)\,,$ (200)

and the induced gauge transformations on the functions $h_{0}$, $h_{1}$ and
$h_{2}$ are given, according to (195), by

$\displaystyle h_{0}\rightarrow h_{0}-\dot{\xi}\,,\qquad h_{1}\rightarrow
h_{1}-\xi^{\prime}+\frac{2}{r}\xi\,,\quad h_{2}\;\rightarrow\;h_{2}-2\xi\,,$
(201)

where we have dropped the indices $(\ell m)$ for simplicity. A dot and a prime
denote a derivative with respect to $t$ and $r$, respectively.

As a consequence, one can always choose a gauge in which $h_{2}^{\ell m}=0$
which is the well-known Regge-Wheeler gauge for the axial perturbations
Regge:1957td . Notice that this gauge choice is possible for $\ell\geq 2$ only
(the cases $\ell=0$ and $\ell=1$ will be discussed later below).

###  A.2 Even-parity or polar perturbations: Zerilli gauge

Before gauge fixing, polar perturbations of the metric are parametrised by
seven families of functions $H_{0}^{\ell m},H_{1}^{\ell m},H_{2}^{\ell m}$,
$\alpha^{\ell m}$, $\beta^{\ell m}$, $K^{\ell m}$ and $G^{\ell m}$ of the
variables $(r,t)$ which appear in the components of the metric perturbations
as follows,

$\displaystyle h_{tt}=A(r)\sum_{\ell,m}H_{0}^{\ell m}(t,r)Y_{\ell
m}(\theta,\varphi),\quad h_{tr}=\sum_{\ell,m}H_{1}^{\ell m}(t,r)Y_{\ell
m}(\theta,\varphi),$ (202) $\displaystyle
h_{rr}=\frac{1}{A(r)}\sum_{\ell,m}H_{2}^{\ell m}(t,r)Y_{\ell
m}(\theta,\varphi),$ (203) $\displaystyle h_{ta}=\sum_{\ell,m}\beta^{\ell
m}(t,r)\partial_{a}Y_{\ell m}(\theta,\varphi),\quad
h_{ra}=\sum_{\ell,m}\alpha^{\ell m}(t,r)\partial_{a}Y_{\ell
m}(\theta,\varphi),$ (204) $\displaystyle h_{ab}=\sum_{\ell,m}K^{\ell
m}(t,r)g_{ab}Y_{\ell m}(\theta,\varphi)+\sum_{\ell,m}G^{\ell
m}(t,r)D_{a}D_{b}Y_{\ell m}(\theta,\varphi)\,.$ (205)

More precisely, the angular part of the metric can be written as

$\displaystyle h_{\theta\theta}$ $\displaystyle=\sum_{\ell,m}K^{\ell
m}(t,r)Y_{\ell m}(\theta,\varphi)+\sum_{\ell,m}G^{\ell
m}(t,r)\partial_{\theta}^{2}Y_{\ell m}(\theta,\varphi)\,,$ (206)
$\displaystyle h_{\theta\varphi}$
$\displaystyle=h_{\varphi\theta}=-\sum_{\ell,m}G^{\ell
m}(t,r)\operatorname{cotan}\theta\,\partial_{\varphi}Y_{\ell
m}(\theta,\varphi)\,,$ (207) $\displaystyle h_{\varphi\varphi}$
$\displaystyle=\sum_{\ell,m}\sin^{2}\theta\,K^{\ell m}(t,r)Y_{\ell
m}(\theta,\varphi)+\sum_{\ell,m}G^{\ell
m}(t,r)\left(\partial_{\varphi}^{2}+\sin\theta\cos\theta\,\partial_{\theta}\right)Y_{\ell
m}(\theta,\varphi)\,.$ (208)

Similarly to the axial sector, this parametrisation is redundant and can be
simplified by gauge fixing. Now, linear diffeomorphisms which preserve even-
parity of the metric components are generated by vector fields $\xi$ whose
components decompose into spherical harmonics as follows,

$\displaystyle\xi_{t}$ $\displaystyle=\sum_{\ell,m}T^{\ell m}(t,r)Y_{\ell
m}(\theta,\varphi)\,,\quad\xi_{r}=\sum_{\ell,m}R^{\ell m}(t,r)Y_{\ell
m}(\theta,\varphi)\,,$ (209) $\displaystyle\xi_{\theta}$
$\displaystyle=\sum_{\ell,m}\Theta^{\ell m}(t,r)\partial_{\theta}Y_{\ell
m}(\theta,\varphi)\,,\quad\xi_{\varphi}=\sum_{\ell,m}\Theta^{\ell
m}(t,r)\partial_{\varphi}Y_{\ell m}(\theta,\varphi)\,.$

Here $T^{\ell m}$, $R^{\ell m}$ and $\Theta^{\ell m}$ are arbitrary functions
of $(t,r)$. These linear diffeomorphisms induce gauge transformations on the
functions that parametrise metric perturbations according to

$\displaystyle\begin{aligned} H_{0}^{\ell m}(t,r)&\longrightarrow H_{0}^{\ell
m}(t,r)+\frac{2}{A(r)}\dot{T}^{\ell m}(t,r)+A^{\prime}(r)R^{\ell m}(t,r),\\\
H_{1}^{\ell m}(t,r)&\longrightarrow H_{1}^{\ell m}(t,r)+\dot{R}^{\ell
m}(t,r)+T^{\prime\ell m}(t,r)+\frac{A^{\prime}(r)}{A(r)}T^{\ell m}(t,r),\\\
H_{2}^{\ell m}(t,r)&\longrightarrow H_{2}^{\ell m}(t,r)+2A(r)R^{\prime\ell
m}(t,r)-A^{\prime}(r)R_{\ell m}(t,r),\\\ \beta^{\ell
m}(t,r)&\longrightarrow\beta^{\ell m}(t,r)+T^{\ell m}(t,r)+\dot{\Theta}^{\ell
m}(t,r),\\\ \alpha^{\ell m}(t,r)&\longrightarrow\alpha^{\ell m}(t,r)+R^{\ell
m}(t,r)+\Theta^{\prime\ell m}(t,r)-\frac{2}{r}\Theta^{\ell m}(t,r),\\\ K^{\ell
m}(t,r)&\longrightarrow K^{\ell m}(t,r)+\frac{2A(r)}{r}R^{\ell m}(t,r),\\\
G^{\ell m}(t,r)&\longrightarrow G^{\ell m}(t,r)+2\Theta^{\ell
m}(t,r)\,.\end{aligned}$ (210)

An immediate consequence of the gauge transformations is that one can choose
the gauge parameter $\xi$ such that $G^{\ell m}=0$ by fixing $\Theta^{\ell
m}$, then $\alpha^{\ell m}=0$ and $\beta^{\ell m}=0$ by fixing $R^{\ell m}$
and $T^{\ell m}$ respectively, in the case where $\ell\geq 2$. This gauge is
known as the Zerilli gauge Zerilli:1970se  (see Kobayashi:2014wsa  for a
recent presentation in the context of modified gravity).

###  A.3 Monopole and dipole perturbations

We consider here the special cases $\ell=0$ and $\ell=1$.

####  A.3.1 Axial modes

For the axial modes, the components $h_{ab}$ vanish identically for $\ell=1$
(axial perturbations do not have $\ell=0$ components) which means that $h_{2}$
does not show up in the components of the metric. Hence, when $\ell=1$, it is
necessary to make a different gauge choice. In general, one chooses $h_{1}=0$
which fixes the gauge parameter $\xi$ up to a function of the form
$C(t)r^{2}$. Therefore, $h_{0}$ inherits a residual gauge invariance given by
$h_{0}\rightarrow h_{0}+F(t)r^{2}$ where $F(t)$ is an arbitrary function. Then
$h_{0}$ can be shown to satisfy the equation of motion,

$\displaystyle 2h_{0}(r)-rh_{0}^{\prime}(r)=0\,.$ (211)

Therefore, the mode $h_{0}$ is not propagating.

####  A.3.2 Polar modes

Let us now turn to polar perturbations. In the case $\ell=0$, $H_{0}$,
$H_{1}$, $H_{2}$ and $K$ are the only non-vanishing components of the metric
perturbations whereas $T$ and $R$ are the only non-vanishing components of the
gauge parameter (so that the gauge transformation preserves the monopole). As
in the general case, one can choose $R$ to fix $K=0$. Then, one can in
principle make use of $T$ to get rid of $H_{1}$ (we could have also set
$H_{0}=0$). Finally, we are left with only two non-vanishing functions which
are either $H_{2}$ or $H_{0}$ and we will compute the corresponding equations
of motion in the next section.

The main difference, concerning the gauge fixing, between the general case and
the case $\ell=1$ lies in the fact that, in the latter, $h_{ab}$ can be shown
to depend on the difference $G-K$ only, so that one can fix $K=0$ without loss
of generality. Furthermore, one can make the gauge fixing $G=0$ by an
appropriate choice of $\Theta$. Then, one makes use of $T$ to fix $\beta=0$.
Finally, one uses the remaining free gauge function $R$ to fix $\alpha=0$. At
the end, we are left with the three non-vanishing functions $H_{0}$, $H_{1}$
and $H_{2}$. The dynamics of these three free parameters will be studied in
the next section as well.

Concerning the monopole $(\ell=0)$, we showed in subsubsection II.1.3 that its
dynamics is fully described in terms of the functions $H_{0}$ and $H_{2}$
only, as all the others can be sent to $0$ by gauge fixing. Thus the equations
of motion simplifies drastically and, after some calculations, give

$\displaystyle H_{0}(r)-H_{2}(r)=0\,,\qquad
H_{2}(r)+(r-r_{s})H_{2}^{\prime}(r)=0\,.$ (212)

The solution reads $H_{2}(r)=C/(r-r_{s})$ and the mode is not propagating.

Finally, the dynamics of the polar dipole $(\ell=1)$ is described by the three
non-vanishing functions $H_{0}$, $H_{1}$ and $H_{2}$ which satisfy the three
independent equations,

$\displaystyle 2H_{2}(r)+(r-r_{s})H_{2}^{\prime}(r)=0\,,\quad H_{1}(r)+i\omega
H_{2}(r)=0\,,$ $\displaystyle H_{0}(r)+(r_{s}-r)H_{0}^{\prime}(r)-2ir\omega
H_{1}(r)+H_{2}(r)=0\,.$ (213)

Indeed, the full set of the original Einstein equations is equivalent to this
one which can easily be solved explicitly but its solution is not relevant for
our purpose. Nonetheless, we see immediately from the equations that, like the
monopole, the polar dipole does not propagate. This is why we do not consider
it in the rest of the paper.

##  Appendix B Equations of motion for the polar perturbations

In this appendix, we present the equations of motion satisfied by the polar
perturbations and show how the system (17) is obtained. The Euler-Lagrange
equations equations of motion (5) yield, in the polar sector,

$\displaystyle{\cal E}_{tt}=$
$\displaystyle-2(\lambda+2)\left(1-\frac{r_{s}}{r}\right)H_{2}(t,r)-2\lambda\left(1-\frac{r_{s}}{r}\right)K(t,r)-\frac{2}{r}(r-r_{s})^{2}\partialderivative{H_{2}}{r}$
(214)
$\displaystyle+\left(6r-11r_{s}+\frac{5r_{s}^{2}}{r}\right)\partialderivative{K}{r}+2(r-r_{s})^{2}\partialderivative[2]{K}{r}=0\,,$
$\displaystyle{\cal E}_{tr}=$
$\displaystyle-2(\lambda+1)H_{1}(t,r)-2r\partialderivative{H_{2}}{t}+r\frac{2r-3r_{s}}{r-r_{s}}\partialderivative{K}{t}+2r^{2}\partialderivative{K}{t}{r}(t,r)=0\,,$
$\displaystyle{\cal E}_{rr}=$
$\displaystyle-2\frac{\lambda+1}{1-r_{s}/r}H_{0}(t,r)+\frac{2}{1-r_{s}/r}H_{2}(t,r)+\frac{2\lambda}{1-r_{s}/r}K(t,r)+2r\partialderivative{H_{0}}{r}-r\frac{2r-r_{s}}{2(r-r_{s})}\partialderivative{K}{r}$
$\displaystyle-\frac{4r^{2}}{r-r_{s}}\partialderivative{H_{1}}{t}+\frac{2r^{4}}{(r-r_{s})^{2}}\partialderivative[2]{K}{t}=0\,,$
$\displaystyle{\cal E}_{t\theta}=$
$\displaystyle-\frac{r_{s}}{r}H_{1}(t,r)-(r-r_{s})\partialderivative{H_{1}}{r}+r\partialderivative{H_{2}}{t}+r\partialderivative{K}{t}=0\,,$
$\displaystyle{\cal E}_{r\theta}=$
$\displaystyle\frac{2r-3r_{s}}{2(r-r_{s})}H_{0}(t,r)-\frac{2r-r_{s}}{2(r-r_{s})}H_{2}(t,r)-r\partialderivative{H_{0}}{r}+r\partialderivative{K}{r}+\frac{r^{2}}{r-r_{s}}\partialderivative{H_{1}}{t}=0\,,$
$\displaystyle{\cal E}_{\theta\theta}=$
$\displaystyle\frac{2r+r_{s}}{2}\partialderivative{H_{0}}{r}+\frac{2r-r_{s}}{2}\partialderivative{H_{2}}{r}-(2r-r_{s})\partialderivative{K}{r}+r(r-r_{s})\partialderivative[2]{H_{0}}{r}-r(r-r_{s})\partialderivative[2]{K}{r}$
$\displaystyle-r\frac{2r-r_{s}}{r-r_{s}}\partialderivative{H_{1}}{t}-2r^{2}\partialderivative{H_{1}}{t}{r}+\frac{r^{3}}{r-r_{s}}\partialderivative[2]{H_{2}}{t}+\frac{r^{3}}{r-r_{s}}\partialderivative[2]{K}{t}=0\,,$
$\displaystyle{\cal E}_{\theta\varphi}=$ $\displaystyle
H_{0}(t,r)-H_{2}(t,r)=0\,.$

The equations of motion ${\cal E}_{t\varphi}=0$, ${\cal E}_{r\varphi}=0$ and
${\cal E}_{\varphi\varphi}=0$ are identical to ${\cal E}_{t\theta}=0$, ${\cal
E}_{r\theta}=0$ and ${\cal E}_{\theta\theta}=0$, respectively.

We can immediately solve the last equation of the system (214) and replace
$H_{2}$ by $H_{0}$ in all the other equations. We thus get six equations for
only three independent functions $K$, $H_{0}$ and $H_{1}$, and we want to
extract three “simple” independent equations out of them. One can then note
that the combination

$\mathcal{E}\equiv\frac{ir_{s}}{4\omega
r(r-r_{s})}\mathcal{E}_{tr}+\frac{1}{2}\mathcal{E}_{rr}+\mathcal{E}_{r\theta}$
(215)

is purely algebraic, i.e. it does not involve any derivatives of the
functions. Moreover, we find that the system $\mathcal{E}_{tr}$,
$\mathcal{E}_{t\theta}$, $\mathcal{E}_{r\theta},\mathcal{E}$ enables us to
recover $\mathcal{E}_{tt}$ and $\mathcal{E}_{\theta\theta}$ so that we can
restrict immediately to the system formed by these four equations which, after
some simple calculations, are given by the system of differential equations

$\displaystyle K^{\prime}(r)-\frac{1}{r}H_{0}(r)-\frac{i(\lambda+1)}{\omega
r^{2}}H_{1}(r)+\frac{1}{r}\frac{2r-3r_{s}}{2(r-r_{s})}K(r)=0\,,$ (216)
$\displaystyle H_{1}^{\prime}(r)+\frac{i\omega
r}{r-r_{s}}H_{0}(r)+\frac{r_{s}}{r(r-r_{s})}H_{1}(r)+\frac{i\omega
r}{r-r_{s}}K(r)=0\,,$ $\displaystyle
H_{0}^{\prime}(r)-K^{\prime}(r)+\frac{r_{s}}{r(r-r_{s})}H_{0}(r)+\frac{i\omega
r}{r-r_{s}}H_{1}(r)=0\,,$

together with the algebraic equation

$\displaystyle\left(\frac{3r_{s}}{r}+2\lambda\right)H_{0}(r)+\left(\frac{ir_{s}(\lambda+1)}{\omega
r^{2}}-2i\omega r\right)H_{1}(r)+\frac{3r_{s}^{2}+2r_{s}(2\lambda-1)r-4\lambda
r^{2}+4\omega^{2}r^{4}}{2r(r-r_{s})}K(r)=0\,.$

One equation is still redundant. However, we can solve the algebraic equation
for $H_{0}$ and substitute its expression into the first three equations. This
shows that the third is not independent from the first two. Finally, we obtain

$\displaystyle K^{\prime}(r)=$
$\displaystyle\frac{3r_{s}^{2}+r_{s}(\lambda-2)r-2\omega^{2}r^{4}}{r(r-r_{s})(3r_{s}+2\lambda
r)}K(r)+\frac{i}{\omega
r^{2}}\left(\lambda+1+\frac{-r_{s}(\lambda+1)+2\omega^{2}r^{3}}{3r_{s}+2\lambda
r}\right)H_{1}(r)\,,$ (217) $\displaystyle H_{1}^{\prime}(r)=$
$\displaystyle\frac{ir(9r_{s}^{2}+8r_{s}(\lambda-1)r-8\lambda
r^{2}+4\omega^{2}r^{4})}{2(r-r_{s})^{2}(3r_{s}+2\lambda r)}\omega
K(r)-\frac{3r_{s}^{2}+r_{s}(1+3\lambda)r-2\omega^{2}r^{4}}{r(r-r_{s})(3r_{s}+2\lambda
r)}H_{1}(r)\,,$

and we obtain the required form (17) with the definitions $X_{1}(r)\equiv
K(r)$ and $X_{2}(r)\equiv H_{1}(r)/\omega$.

##  Appendix C Flowchart for the algorithm

In this appendix, we draw a flowchart to illustrate the algorithm that we are
using to compute the asymptotic behaviour of a solution of a first order
system.

It should be noted that, in principle, one can skip the first question “Is the
leading term diagonalisable?” and put directly the leading order term in its
Jordan form. Indeed, when the leading term is diagonalisable, putting it into
its Jordan form is equivalent to diagonalising it and the resulting Jordan
matrix is made of $d$ one-dimensional blocks where $d$ is the dimension of the
system, thus of the matrix. Therefore, the procedure for splitting the system
into several subsystems described in subsection IV.5 is in this case
equivalent to the procedure described in subsection IV.2 where we are treating
several blocks.

Start: Consider the leading term in the expansion of $M(z)$ Is the leading
term diagonalisable? Diagonalise step by step with transformations in
subsection IV.2 and Stop Put it under the form (181) Does it have
diagonalisable blocks? Split the system into its nilpotent and diagonalisable
part (Splitting Lemma of IV.5) and proceed with each part separately Is the
Poincaré rank $r$ equal to -1? Stop: compute the logarithmic asymptotic
behaviour (see 169) Normalise the next-to-leading order with transformations
(184) Use a transformation of the form (163) for each block ((185) with
$p/q=1$) to try to reduce the order Normalise the next order Is the order
reduced? Try again with a change of variables $z\rightarrow z^{p/q}$ (see Eq.
185) Is the leading order still nilpotent? noyesyesnopt(all blocksare
nilpotent)noyesnoyesno

## References

  * (1) E. Berti, V. Cardoso, and C. M. Will, “On gravitational-wave spectroscopy of massive black holes with the space interferometer LISA,” Phys. Rev. D 73 (2006) 064030, gr-qc/0512160.
  * (2) E. Berti, K. Yagi, H. Yang, and N. Yunes, “Extreme Gravity Tests with Gravitational Waves from Compact Binary Coalescences: (II) Ringdown,” Gen. Rel. Grav. 50 (2018), no. 5 49, 1801.03587.
  * (3) T. Regge and J. A. Wheeler, “Stability of a Schwarzschild singularity,” Phys. Rev. 108 (1957) 1063–1069.
  * (4) F. J. Zerilli, “Effective potential for even parity Regge-Wheeler gravitational perturbation equations,” Phys. Rev. Lett. 24 (1970) 737–738.
  * (5) 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), no. 4 044061, 1906.05155.
  * (6) D. Langlois and K. Noui, “Degenerate higher derivative theories beyond Horndeski: evading the Ostrogradski instability,” JCAP 1602 (2016), no. 02 034, 1510.06930.
  * (7) M. Crisostomi, K. Koyama, and G. Tasinato, “Extended Scalar-Tensor Theories of Gravity,” JCAP 1604 (2016), no. 04 044, 1602.03119.
  * (8) J. Ben Achour, D. Langlois, and K. Noui, “Degenerate higher order scalar-tensor theories beyond Horndeski and disformal transformations,” Phys. Rev. D93 (2016), no. 12 124005, 1602.08398.
  * (9) J. Ben Achour, M. Crisostomi, K. Koyama, D. Langlois, K. Noui, and G. Tasinato, “Degenerate higher order scalar-tensor theories beyond Horndeski up to cubic order,” JHEP 12 (2016) 100, 1608.08135.
  * (10) S. Chandrasekhar, “The mathematical theory of black holes,” in Oxford, UK: Clarendon (1992) 646 p., Oxford, UK: Clarendon (1985) 646 P., 1985\.
  * (11) E. W. Leaver, “An Analytic representation for the quasi normal modes of Kerr black holes,” Proc. Roy. Soc. Lond. A402 (1985) 285–298.
  * (12) K. D. Kokkotas and B. G. Schmidt, “Quasinormal modes of stars and black holes,” Living Rev. Rel. 2 (1999) 2, gr-qc/9909058.
  * (13) H.-P. Nollert, “Topical Review: Quasinormal modes: the characteristic ‘sound’ of black holes and neutron stars,” Class. Quant. Grav. 16 (1999) R159–R216.
  * (14) E. Berti, V. Cardoso, and A. O. Starinets, “Quasinormal modes of black holes and black branes,” Class. Quant. Grav. 26 (2009) 163001, 0905.2975.
  * (15) R. A. Konoplya and A. Zhidenko, “Quasinormal modes of black holes: From astrophysics to string theory,” Rev. Mod. Phys. 83 (2011) 793–836, 1102.4014.
  * (16) W. Wasow, Asymptotic Expansions for Ordinary Differential Equations.  Courier Dover publications, 1965.
  * (17) W. Balser, “Computation of formal fundamental solutions,” Linear Algebra and its Applications 288 (Feb., 1999) 293–312.
  * (18) M. Barkatou and E. Pflügel, “An Algorithm Computing the Regular Formal Solutions of a System of Linear Differential Equations,” Journal of Symbolic Computation 28 (Oct., 1999) 569–587.
  * (19) E. Pflügel, “A Root-Free Splitting-Lemma for Systems of Linear Differential Equations,” arXiv:1911.05837 [cs] (Nov., 2019). arXiv: 1911.05837.
  * (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), no. 12 546, 1709.09178.
  * (21) J. P. Boyd, Chebyshev & Fourier Spectral Methods.  Springer Berlin Heidelberg, 1989.
  * (22) M. A. Barkatou, “An algorithm to compute the exponential part of a formal fundamental matrix solution of a linear differential system,” Applicable Algebra in Engineering, Communication and Computing 8 (Jan., 1997) 1–23.
  * (23) Abdelaziz Hilali, Solutions formelles de systèmes différentiels linéaires au voisinage d’un point singulier.  PhD thesis, Université Joseph-Fourier - Grenoble I, 1987.
  * (24) T. Kobayashi, H. Motohashi, and T. Suyama, “Black hole perturbation in the most general scalar-tensor theory with second-order field equations II: the even-parity sector,” Phys. Rev. D89 (2014), no. 8 084042, 1402.6740.