# Revisiting the quasinormal modes of the Schwarzschild black hole: Numerical
analysis

Luis A. H. Mamani<EMAIL_ADDRESS>Centro de Ciências Exatas
Naturais e Tecnológicas,
Universidade Estadual da Região Tocantina do Maranhão,
Rua Godofredo Viana 1300, 65901- 480, Imperatriz, Maranhão, Brazil  Angel D.
D. Masa<EMAIL_ADDRESS>Centro de Ciências Naturais e Humanas,
Universidade Federal do ABC, Avenida dos Estados 5001, 09210-580 Santo André,
São Paulo, Brazil  Lucas Timotheo Sanches<EMAIL_ADDRESS>Centro de
Ciências Naturais e Humanas, Universidade Federal do ABC, Avenida dos Estados
5001, 09210-580 Santo André, São Paulo, Brazil  Vilson T. Zanchin
<EMAIL_ADDRESS>Centro de Ciências Naturais e Humanas, Universidade
Federal do ABC, Avenida dos Estados 5001, 09210-580 Santo André, São Paulo,
Brazil

###### Abstract

We revisit the problem of calculating the quasinormal modes of spin $0$,
$1/2$, $1$, $3/2$, $2$, and spin $5/2$ fields in the asymptotically flat
Schwarzschild black hole spacetime. Our aim is to investigate the problem from
the numerical point of view, by comparing some numerical methods available in
the literature and still not applied for solving the eigenvalue problems
arising from the perturbation equations in the Schwarzschild black hole
spacetime. We focus on the pseudo-spectral and the asymptotic iteration
methods. These numerical methods are tested against the available results in
the literature, and confronting the precision between each other. Besides
testing the different numerical methods, we calculate higher overtones
quasinormal frequencies for all the investigated perturbation fields in
comparison with the known results. In particular, we obtain purely imaginary
frequencies for spin $1/2$ and $3/2$ fields that are in agreement with
analytic results reported previously in the literature. The purely imaginary
frequencies for the spin $1/2$ field are exactly the same as the frequencies
obtained for the spin $3/2$ field. In turn, the quasinormal frequencies for
the spin $5/2$ perturbation field are calculated for the very first time, and
purely imaginary frequencies are found also in this case. We conclude that
both methods provide accurate results and they complement each other.

##  I Introduction

Perturbation theory is a very useful theoretical toolkit for the investigation
of properties of a physical system. For example, the stability under small
perturbations. In the harmonic oscillator problem, it drives into a second-
order differential equation with Dirichlet boundary conditions whose solutions
are characterized by a set of discrete real frequencies, i.e., normal modes
(NMs). However, there are physical systems whose boundary conditions drive
solutions with complex frequencies, i.e., quasinormal modes (QNMs), for
example, a harmonic oscillator into a dissipative medium, see for instance
[1]. Thus, the investigation of the quasinormal (QN) frequencies and its
mathematical properties become fascinating subject which may shed light into
the understanding of universal properties of the physical system under
investigation.

In the context of gravity theories, perturbation theory is important for
several reasons. One of the motivations is the investigation of gravitational
waves spectroscopy [2]. One also may use perturbation theory to investigate
the stability under small perturbations of determined background solutions. It
has been shown that the perturbation equations may be written as a second-
order differential equations allowing us to use numerical techniques
implemented in differential equations to solve them. One of the boundary
conditions considers that classically noting comes out from the black hole
interior, so that the boundary condition at the horizon are ingoing waves. In
turn, at the spatial infinite the condition are outgoing waves, because
nothing can come from outside of the spacetime. Solving the perturbation
equations under these particular boundary conditions drive to solutions with
discrete (complex) eigenvalues. Such eigenvalues are frequencies representing
the characteristic oscillations of the black hole that relaxes after being
perturbed. One very interesting property of these frequencies is that they do
not depend on the initial perturbation and are fixed completely by the
properties of the black hole under consideration.

The quasinormal modes are of particular interest in black hole astrophysics.
Direct observations showed in the coalescence of a binary system emits
gravitational waves in the form of QNMs, that is, the final system obeys the
predictions of black hole perturbation theory [3]. This fact alone is enough
to motivate the development, testing, and comparison of different methods for
finding QNMs, but teh relevance o QNMs go far further this fact. For reviews
and additional discussions about QNMs in different contexts see, for instance,
[4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24,
25, 26, 27, 28] and references therein.

In turn, the first observation of the shadow of the supermassive black hole
M87* by the Event Horizon Telescope Collaboration (EHT) [29, 30] open a new
window for the investigation of strong gravitational field phenomena. There is
also the studies proposing a connection between the real part of the
quasinormal (QN) frequency with the radius of the shadow, see the Refs. [31,
32, 33, 34, 35]. Among the new possibilities, there are researched papers
attempting to find constraints in parameters arising in different models like
in general uncertainty principle (GUP), see for instance Ref. [36, 37].

In this paper we are going to revisit the calculation of quasinormal modes for
integer spin $0$, $1$ and $2$ fields, as well as for semi-integer spin $1/2$,
$3/2$ and $5/2$ fields. Since Chandrasekhar calculated the quasinormal modes
for $s=2$ in Ref. [38], the problem of calculating QNMs for other spin fields
were previously investigated in the literature using different techniques
(numeric and analytic), see for instance Refs. [39, 40, 41, 42, 43]. However,
our approach here is from the numerical point of view, for doing so we are
going to use two numerical methods well established in the literature. The
first one is the pseudo-spectral method used to solve differential equations
expanding the solution in a base composed by special functions [44]. This
method was used to calculate the quasinomal modes of Schwarzschild black hole
for spin zero field in Ref. [45]. However, we extend the method for
calculating the QNMs for spin $1/2$, $1$, $3/2$, $2$ and $5/2$ fields. We also
use the asymptotic iteration method (AIM) proposed originally in Ref. [46].
This method was extended to solve quasinormal modes in Ref. [43]. In this
paper we review the relatively unexplored asymptotic iteration method and
apply it to the QNM problem. We also introduce a new software package that
implements the latter method for usage in general second order ODEs.

This paper is organized as follows. In Section II we write the equations of
motion describing the spin $0$, $1/2$, $1$, $3/2$, $2$ and $5/2$ fields in a
suitable form to apply the numerical methods. In Section III we review and
discuss the pseudo-spectral method, we focus in the way how this method can be
applied, by expanding the solution using one or two special functions. Section
IV is devoted to discuss about the AIM and it extension to calculate QN
frequencies. We also present a open source code that can be used freely. In
turn, in Section V we present our numerical results obtained by both methods,
we also compare against numerical results available in the literature. We
leave the discussion of the QNMs in the limit of large angular for Section VI,
were we also compare against analytic results. Finally, our main conclusions
are presented in Section VII. Additional details are presented in Appendices A
and B.

##  II Equations of motion

In this section we write the equations describing the field perturbations on
the gravitational background solution of the Einstein equations. We focus on
the metric for an spherical symmetric black hole, which is given by [47]

$ds^{2}=-f(r)\,dt^{2}+\frac{1}{f(r)}dr^{2}+r^{2}d\theta^{2}+r^{2}\sin^{2}{\theta}\,d\varphi^{2},$
(1)

where the horizon function of the Schwarzschild black hole is given by

$f(r)=1-\frac{2M}{r}.$ (2)

where $M$ is the mass of the black hole, and $r$ is the radial coordinate
which, in principle, belongs to the interval $r\in[0,\,\infty)$. The
coordinates in the metric (1) are known as Schwarzschild coordinates. As it is
well known, this metric presents an event horizon at $r=2M$ and a curvature
(physical) singularity at $r=0$. In the asymptotic region, i.e., $r\to\infty$,
the metric reduces to flat metric. As long as the quasinormal modes in this
black hole spacetime are concerned, the interesting region is the spacetime
region spanned by the radial coordinate $r$ in the interval $2M<r<\infty$.

###  II.1 Spin $0$, $1$ and $2$ perturbations

Here we revisit the study of perturbations of integer spin, such that scalar,
vector, and gravitational perturbations in the Schwarzschild black hole
peacetime. This is a long standing problem, and there are a considerable
amount of results published in the literature, what is certainly interesting
for the purpose of the present work. In fact, it was proven that the equations
of motion can be written in a compact form, the so called, Schrödinger-like
differential equations, see for instance [4]. Thus, for massless scalar
($s=0$), electromagnetic ($s=1$) and vector type gravitational perturbations
($s=2$), the Schrödinger-like equations are given by

$\frac{d^{2}\psi_{{s}}(r)}{dr_{*}^{2}}+\left[\omega^{2}-V_{{s}}(r)\right]\psi_{{s}}(r)=0.$
(3)

where the potential is given by

$V_{s}(r)=f(r)\left(\frac{\ell(\ell+1)}{r^{2}}+(1-s^{2})\frac{2M}{r^{3}}\right),$
(4)

where the tortoise coordinate is defined in terms of the areal coordinate $r$
by $dr_{*}=dr/f(r)$. So far, the problem of calculating quasinormal modes was
reduced to solve an eigenvalue problem. We will see that it is possible to
solve this problem following two approaches, one of them expanding the
function $\psi$ in a base composed by especial functions, while the other
solving directly the second-order differential equation.

On the other hand, note that the potential (4) is zero at the horizon,
$f(r_{h})=0$. Thus, the Schrödinger-like equation reduces to a single harmonic
oscillator problem, whose solutions are:

$\psi_{s}(r)=c_{1}\,e^{-i\omega r_{*}}+c_{2}\,e^{i\omega r_{*}},\quad r\to
r_{h}.$ (5)

The first of these solutions is interpreted as an ingoing wave, i.e., a wave
that travels inward and eventually falls into the black hole event horizon.
The second solution is interpreted an outgoing wave, i.e., a wave that travels
outward with respect to the black hole and can escape to space infinity. Waves
travelling as this second solution would represent waves coming from the
interior of the black hole. Since the perturbation theory is implemented using
classical assumptions, nothing is expected to come out from the black hole
interior, thus, in the following analysis we impose the first solution as
boundary condition at the horizon, we set $c_{2}=0$.

We also need to investigate the spatial infinity, at $r\to\infty$, where
$f(r)\to 1$ and the effective potential (4) also vanishes. Thus, in such a
limit the general solutions to the wave equation (3) have the same form as the
function given in Eq. 5, i.e.

$\psi_{s}(r_{*})=c_{3}\,e^{-i\omega r_{*}}+c_{4}\,e^{i\omega r_{*}},\quad
r\to\infty.$ (6)

The first solution is interpreted physically as waves coming in from outside
the universe and must be avoided setting $c_{3}=0$. In turn, the second
solution represent waves going out the universe, this is the boundary
condition at the spatial infinite. Finally, note that the boundary conditions
do not depend explicitly on the angular momentum $\ell$ neither the spin.

It is interesting to see the behavior of the tortoise coordinate, which close
to the horizon is given by ($r_{h}=2M$)

$r_{*}=\int\frac{dr}{f^{\prime}(r_{h})(r-r_{h})}\approx\frac{\ln{(r-r_{h})}}{f^{\prime}(r_{h})},\quad
r\to r_{h}$ (7)

Thus, in terms of the radial coordinate the boundary condition at the horizon
becomes ($f^{\prime}(r_{h})=1/r_{h}$)

$\psi_{s}(r)\sim
e^{-i\omega\frac{\ln{(r-r_{h})}}{f^{\prime}(r_{h})}}\sim\left(r-r_{h}\right)^{-i\frac{\omega}{f^{\prime}(r_{h})}}.$
(8)

In turn, the tortoise coordinate at the spatial infinite becomes

$r_{*}=\int\frac{dr}{f(r)}\approx r+r_{h}\,\ln{r},\quad r\to\infty$ (9)

while the asymptotic solution at the spatial infinite becomes

$\psi_{s}(r)\sim e^{i\omega(r+r_{h}\ln{r})}\sim r^{i\,r_{h}\omega}e^{i\omega
r}.$ (10)

In the following, we are going to change our strategy and use a new coordinate
defined by $u=2M/r$. This is equivalent to the choice $u=1/r$ and then
normalizing the mass $M$ to $2M=1$. The relation between the tortoise
coordinate $r_{*}$ and the new coordinate $u$ becomes $du/dr_{*}=-u^{2}f(u)$.

We also constraint our analysis to the outer region of the black hole, such
that $r_{h}\leq r<\infty$. Hence, in terms of the new coordinate this region
is bounded to the interval $u\in[0,1]$, and the potential becomes

$V_{s}(u)=f(u)\,u^{2}\left[\ell\left(1+\ell\right)+\left(1-s^{2}\right)\,u\right],$
(11)

where $f(u)=1-\,u$.

To implement the pseudo-spectral method, one must write the equations on the
background metric written terms of the Eddington-Filkenstein coordinates, see
for instance the discussion in Ref. [45]. However, we found a short way to
write the equations going directly from the Schrodinger-like equation by
implementing some transformations. At the end, this transformations lead to a
differential equation which is the same as the one obtained from the metric in
Eddington-Filkenstein coordinates. In order to write the perturbation
equations in terms of the Eddington-Filkenstein coordinate we must implement
the transformation $\psi_{s}\to\Phi_{s}$ given by

$\psi_{s}=\frac{\Phi_{s}(u)}{u}e^{-i\omega\,r_{*}(u)}.$ (12)

Thus, Eq. (3) becomes

$\begin{split}&\left[s^{2}u^{2}-\ell(\ell+1)u-2\,i\,\omega\right]\Phi_{s}(u)\\\
&-u\left(u^{2}-2\,i\,\omega\right)\Phi_{s}^{\prime}(u)+\left(1-u\right)\,u^{3}\,\Phi_{s}^{\prime\prime}(u)=0,\end{split}$
(13)

where we have set $M=1/2$ so that $r_{h}=1$. The asymptotic solutions close to
the horizon may be calculated using the ansatz $\Phi_{s}(u)=(1-u)^{\alpha}$.
By substituting this ansatz into (13) we get two solutions,

$\alpha=0,\qquad\qquad\alpha=2\,i\,\omega.$ (14)

The solution for $\alpha=0$ is interpreted physically as the ingoing waves at
the horizon, while the other one is interpreted as a wave coming out from the
black hole interior. Therefore, the second solution is neglected in the
following analysis.

In the same way, we consider the ansatz $\Phi_{s}=u^{\beta}$ to get the
asymptotic solution close to the spatial infinite. Plugging this ansatz in
(13) we get the following solution,

$\Phi_{s}(u)=c_{5}\,e^{2\,i\,\omega/u}u^{-2\,i\,\omega}+c_{6}\,u.$ (15)

We want the divergent solution, such that we set $c_{6}=0$. Then, we implement
the final transformation which takes into consideration the boundary
conditions,

$\Phi_{s}(u)=e^{2\,i\,\omega/u}u^{-2\,i\,\omega}\phi_{s}(u),$ (16)

where $\phi_{s}(u)$ is a regular function in the interval $u\in[0,\,1]$ by
definition. Finally, the equation of motion describing spin $0$, $1$ and $2$
perturbations is given by

$\begin{split}&\Big{[}\ell(\ell+1)\,u-s^{2}u^{2}-4\,i\lambda-16\,u(1+u)\lambda^{2}\Big{]}\phi_{s}(u)+\\\
&\Big{[}u^{3}+4\,i\,u\left(1-2u^{2}\right)\lambda\Big{]}\phi_{s}^{\prime}(u)-(1-u)u^{3}\phi_{s}^{\prime\prime}(u)=0,\end{split}$
(17)

where we have used $\lambda=\omega M=\omega/2$. The final differential
equation is then a quadratic eigenvalue problem in $\lambda$. It is worth also
mentioning that in the limit of zero spin $s\to 0$, Eq. (17) reduces to Eq.
(4.8) of Ref. [45]. These results just proves that the alternative way for
getting the equations for the integer spin perturbations presented here is
consistent with other approaches found in the literature.

The differential equation (17) was solved numerically by means of the pseudo-
spectral and AIM methods. The results for perturbations of spin $0$, $1$, and
$2$ are presented in Secs. V.1, V.3, and V.3, respectively, where a comparison
with the corresponding data in the literature is also performed.

In the following we extend the analysis of the present section to other kind
of perturbations.

###  II.2 Spin $1/2$ perturbations

For half-integer spin perturbations the history is different, the differential
equations are quite distinct from (17). The equation for the spin 1/2 Dirac
field as a perturbation on the Schwarzschild background was derived in Ref.
[40] by using the Newman-Penrose formalism. The analysis was generalized for
arbitrary half-integer spin in Ref. [41]. The resulting equation of motion for
the perturbations may be written in the Schrödinger-like form Eq. (3), where
the potential for the massless spin 1/2 field is given by

$V_{\scriptscriptstyle{1/2}}=\frac{\left(1+\ell\right)\sqrt{f(r)}}{r^{2}}\left[\left(1+\ell\right)\sqrt{f(r)}+\frac{3M}{r}-1\right].$
(18)

It is worth mentioning that we have found a typo in the definition of $\Delta$
in [40], which must be $\Delta=r(r-2M)$.

We then implement the same transformations done in the integer spin cases.
First, we change the radial coordinate to $u=2M/r$, which is defined in the
interval $u\in[0,1]$. Then by setting $2M=1$ the effective potential (18)
becomes

$V_{\scriptscriptstyle{1/2}}=(1+\ell)u^{2}\sqrt{f(u)}\left[\left(1+\ell\right)\sqrt{f(u)}+\frac{3}{2}u-1\right].$
(19)

It is interesting pointing out that the asymptotic solutions do not depend on
the spin of the field, for that reason the asymptotic solutions for this
problem are the same as those obtained in Eqs. (14) and (15). Then, similar
transformations as those giving in Eqs. (12) and (16) can be applied also in
the present spin 1/2 case. Thus, the differential equation to be solved is
given by

$\begin{split}&R(u)\phi_{\scriptscriptstyle{1/2}}(u)+Q(u)\phi_{\scriptscriptstyle{1/2}}^{\prime}(u)+P(u)\,\phi_{\scriptscriptstyle{1/2}}^{\prime\prime}(u)=0,\end{split}$
(20)

in which the coefficients $R(u)$, $Q(u)$, and $P(u)$ are given by

$\displaystyle R(u)=\,$ $\displaystyle
u^{3}+u(1+\ell)\left(1+\ell-\sqrt{1-u}\,\right)$ (21)
$\displaystyle+\frac{u^{2}}{2}\left[(1+\ell)\left(3\sqrt{1-u}-4\right)-2\ell^{2}\right]$
$\displaystyle-4\,i\,(1-u)\lambda-16u(1-u^{2})\lambda^{2},$ $\displaystyle
Q(u)=$ $\displaystyle
u^{3}(1-u)+4\,i\,u\,\lambda\left(1-u-2u^{2}+2u^{3}\right),$ (22)
$\displaystyle P(u)=$ $\displaystyle-u^{3}(1-u)^{2},$ (23)

respectively, and with $\lambda$ standing for $\lambda=M\omega=\omega/2$.

As it can be seen, the differential equation (20) is impregnated by square
roots that may difficult the convergence of the numerical methods. To avoid
the square roots we implement an additional change of variable $\chi^{2}=1-u$.
Nevertheless, the new coordinate also belongs to the interval
$\chi\in[0,\,1]$. The differential equation (20) becomes

$\begin{split}&R(\chi)\phi_{\scriptscriptstyle{1/2}}(\chi)+Q(\chi)\phi_{\scriptscriptstyle{1/2}}^{\prime}(\chi)+P(\chi)\,\phi_{\scriptscriptstyle{1/2}}^{\prime\prime}(\chi)=0,\end{split}$
(24)

in which the coefficients $R(\chi)$, $Q(\chi)$, and $P(\chi)$ are given by

$\displaystyle R(\chi)=\,$ $\displaystyle
2(1-\chi^{2})\left[\left(\ell+1\right)\left(1+2\ell\,\chi-3\chi^{2}\right)+2\ell\,\chi+2\chi^{3}\right]$
(25)
$\displaystyle-8\,i\,\chi\,\lambda-32\,\chi\left(2-3\chi^{2}+\chi^{4}\right)\lambda^{2},$
$\displaystyle Q(\chi)=$
$\displaystyle(\chi^{2}-1\big{)}\left[\big{(}1-\chi^{2}\big{)}^{2}-8\,i\big{(}1-4\chi^{2}+2\chi^{4}\big{)}\lambda\right],$
(26) $\displaystyle P(\chi)=$
$\displaystyle-\chi\big{(}1-\chi^{2}\big{)}^{3}.$ (27)

In Sec. V.4 we solve Eq. (24) by using the pseudo-spectral and AIM methods and
compare our results against the results of Refs. [40, 41].

###  II.3 Spin $3/2$ perturbations

As well as for spin $1/2$ perturbation, the perturbation equation for spin
$3/2$ field is different. To get the equation we are going to use the result
obtained in Ref. [41], specifically Eq. (37) of this reference, setting
$s=3/2$ on this equation we get the potential of the Schrödinger-like equation
(3)

$\begin{split}V_{\scriptscriptstyle{3/2}}&=\frac{(1+\ell)(2+\ell)(3+\ell)\sqrt{f(r)}}{\left[2M+r(1+\ell)(3+\ell)\right]^{2}}\bigg{(}\frac{2M^{2}}{r^{2}}\\\
&+(1+\ell)(3+\ell)\left[(2+\ell)\sqrt{f(r)}+\frac{3M}{r}-1\right]\bigg{)}.\end{split}$
(28)

As before, we change the radial coordinate to $u=1/r$ and consider $2M=1$.
Thus, the potential becomes

$\begin{split}V_{\scriptscriptstyle{3/2}}&=\frac{u^{2}(1+\ell)(2+\ell)(3+\ell)\sqrt{1-u}}{2\big{[}u+(1+\ell)(3+\ell)\big{]}^{2}}\bigg{(}\frac{u^{2}}{2}\\\
&+(1+\ell)(3+\ell)\left[(2+\ell)\sqrt{1-u}+\frac{3u}{2}-1\right]\bigg{)}.\end{split}$
(29)

Following the procedure implemented from Eq. (12) to Eq. (15), where we go
from the Schrödinger-like equation to a differential equation suited to the
numerical methods we are working with, the differential equation (20) becomes

$\begin{split}&R(\chi)\phi_{\scriptscriptstyle{3/2}}(\chi)+Q(\chi)\phi_{\scriptscriptstyle{3/2}}^{\prime}(\chi)+P(\chi)\,\phi_{\scriptscriptstyle{3/2}}^{\prime\prime}(\chi)=0,\end{split}$
(30)

where the coefficients $R(\chi)$, $Q(\chi)$, and $P(\chi)$ are given by

$\displaystyle R(\chi)=\,$ $\displaystyle
2\big{(}1-\chi^{2}\big{)}\Big{[}6+11\ell+6\ell^{2}+\ell^{3}+2\chi^{5}$ (31)
$\displaystyle+4\chi^{4}(2+\ell)+2\chi^{3}\left(3+4\ell+\ell^{2}\right)$
$\displaystyle\\!\\!-\chi^{2}\left(2-7\ell-6\ell^{2}-\ell^{3}\right)+2\chi(2+\ell)^{2}\left(2+4\ell+\ell^{2}\right)\Big{]}$
$\displaystyle\\!-16\chi(2+\chi+\ell)^{2}\lambda\Big{[}i+4\left(2-\chi^{2}\right)\left(1-\chi^{2}\right)\lambda\Big{]},$
$\displaystyle Q(\chi)=$
$\displaystyle-\left(1-\chi^{2}\right)\left(2+\chi+\ell\right)^{2}$ (32)
$\displaystyle\times\left[\left(1-\chi^{2}\right)^{2}-8\,i\,(1-4\chi^{2}+2\chi^{4})\lambda\right],$
$\displaystyle P(\chi)=$
$\displaystyle-\chi\left(1-\chi^{2}\right)^{3}\left(2+\chi+\ell\right)^{2}.$
(33)

where we have used the new coordinate $\chi^{2}=1-u$ to avoid square roots.
Again, we get a quadratic eigenvalue problem, and the function
$\phi_{\scriptscriptstyle{3/2}}(\chi)$ is regular in the interval
$\chi\in[0,1]$. In Sec. V.5 we solve Eq. (30) by using the pseudo-spectral and
AIM methods and compare our results against the results of Refs. [41, 48].

###  II.4 Spin $5/2$ perturbations

It is believed that the investigation of higher spin fields may shed some
light on the understanding of fundamental physics, like on new unifying
theories for the fundamental interactions, or on new phenomenology beyond the
standard model. The main motivation for investigating the spin $5/2$ field
perturbation is the Rarita-Schwinger theory. Inspired by such a theory, the
authors of Ref. [49] computed some physical observable for the spin
$\frac{5}{2}$-field. In this section we use the generic equation obtained in
Ref. [41], specifically Eq. (37), to determine the quasinormal frequencies of
this perturbation field on the Schwarzschild black hole. The differential
equation for the perturbations becomes so huge and for that reason we write it
in Appendix A. The resulting equation is solved numerically by using the
pseudo-spectral and AIM methods. The numerical results are displayed in Sec.
V.6.

##  III The pseudo-spectral method

It is well known that Fourier method is appropriate to solve periodic
problems, nevertheless it cannot be applied for nonperiodic problems due to
the Gibbs phenomenon arising at the boundaries [50]. An alternative method to
solve nonperiodic problems is the pseudo-spectral method, which recently has
being applied to solve differential equations numerically in many problems.
The fact that the coordinate domain is not periodic, $u\in[0,1]$, allows us to
use this method in our problem. Thus, the quadratic eigenvalue problem can be
written in the form (using the notation of Ref. [45]).

$c_{0}(u,\lambda,\lambda^{2})\phi(u)+c_{1}(u,\lambda,\lambda^{2})\phi^{\prime}(u)+c_{2}(u,\lambda,\lambda^{2})\phi^{\prime\prime}(u)=0.$
(34)

The coefficients of this equation can be written as
$c_{j}(u,\lambda,\lambda^{2})=c_{j,0}(u)+\lambda\,c_{j,1}(u)+\lambda^{2}\,c_{j,2}(u)$,
where $c_{j,0}(u)$, $c_{j,1}(u)$, and $c_{j,2}(u)$ are polynomials of $u$
only.

The idea behind the pseudo-spectral method is to rewrite the regular function
$\phi(u)$ in a base composed by cardinal functions $C_{j}(u)$, in the form

$\phi_{s}(u)=\sum_{j=0}^{N}g(u_{j})\,C_{j}(u),$ (35)

where $g(u)$ is a function of $u$. The next step is to evaluate the
differential equation (including these functions) on a grid or collocation
points. The best choice is the Gauss-Lobato grid given by

$u_{i}=\frac{1}{2}\left(1\pm\cos{\left[\frac{i}{N}\pi\right]}\right),\quad
i=0,1,2,\cdots,N$ (36)

Note that (36) maps the interval $[-1,1]$ into $[0,1]$.

Evaluating on the grid, the polynomials of (34) become elements of a matrix
$c_{j}(u_{i},\lambda,\lambda^{2})=c_{j,0}(u_{i})+\lambda\,c_{j,1}(u_{i})+\lambda^{2}\,c_{j,2}(u_{i})$.
Then, the matrix representation of the quadratic eigenvalue problem (34) can
be written as

$\begin{split}\left(\tilde{M}_{0}+\tilde{M}_{1}\lambda+\tilde{M}_{2}\lambda^{2}\right)g=0,\end{split}$
(37)

where

$\begin{split}(\tilde{M}_{0})_{ji}=c_{0,0}(u_{i})D_{ji}+c_{1,0}(u_{i})D^{(1)}_{ji}+c_{2,0}(u_{i})D^{(2)}_{ji},\\\
(\tilde{M}_{1})_{ji}=c_{0,1}(u_{i})D_{ji}+c_{1,1}(u_{i})D^{(1)}_{ji}+c_{2,1}(u_{i})D^{(2)}_{ji},\\\
(\tilde{M}_{2})_{ji}=c_{0,2}(u_{i})D_{ji}+c_{1,2}(u_{i})D^{(1)}_{ji}+c_{2,2}(u_{i})D^{(2)}_{ji},\end{split}$
(38)

here $D_{ji}$, $D^{(1)}_{ji}$, and $D^{(2)}_{ji}$ represent the cardinal
function and its derivatives. Defining $\tilde{g}=\lambda g$, the last
equation may be written in the form

$\begin{split}\tilde{M}_{0}\,g+\left(\tilde{M}_{1}+\tilde{M}_{2}\lambda\right)\tilde{g}=0,&\end{split}$
(39)

This is the first step to linearize the quadratic eigenvalue problem. For a
generalization of this procedure see for instance Ref. [51]. Therefore, the
matrix representation of the eigenvalue problem may be written as

$\left(M_{0}+M_{1}\,\lambda\right)\cdot\vec{g}=\mathbb{0},$ (40)

where we have defined the new matrices

$M_{0}=\begin{pmatrix}\tilde{M}_{0}&\tilde{M}_{1}\\\
\mathbb{0}&\mathbb{1}\end{pmatrix},\,M_{1}=\begin{pmatrix}\mathbb{0}&\tilde{M}_{2}\\\
-\mathbb{1}&\mathbb{0}\end{pmatrix},\,\vec{g}=\begin{pmatrix}g\\\
\tilde{g}\end{pmatrix}.$ (41)

Notice that ${M}_{0}$ and ${M}_{1}$ are $(N+1)\times(N+1)$ matrices and
$\vec{g}$ is a $(N+1)-$dimensional vector with components $g_{j}=g(u_{j})$,
$j=0,\,1,...,N$. Finally, the QN frequencies are determined solving the linear
eigenvalue problem (40). The last procedure was explained for a quadratic
eigenvalue problem. This procedure can be easily extended for arbitrary order
of the eigenvalue problem whenever the power of the frequency is an integer.
However, if the value of the potential changes at the infinite spatial, as in
the case of massive scalar field, see for instance [52], the implementation of
the pseudo-spectral method is not obvious because the power of the frequency
turns out semi-integer such that it is not possible to write the eigenvalue
problem in the form of Eq. (37).

Having described how to calculate the eigenvalues, we need to specify the
cardinal functions. We realized that these functions may depend on one or more
Chebyshev polynomials of the first kind $T_{k}(u)$. In the following, we
consider two forms for the cardinal functions. The first model considers one
Chebyshev polynomial in the form

$C_{j}(u)=T_{j}(u).$ (42)

We call this particular choice as pseudo-spectral I. The second model
considers two Chebyshev polynomials [45]

$C_{j}(u)=\frac{2}{Np_{j}}\sum_{m=0}^{N}\frac{1}{p_{m}}T_{m}(u_{j})T_{m}(u),\quad\begin{cases}p_{0}=2,\\\
p_{N}=2,\\\ p_{j}=1.\end{cases}$ (43)

We call this choice pseudo-spectral II. It is worth mentioning that the
pseudo-spectral method inevitably leads to the emergence of spurious solutions
that do not have any physical meaning. To eliminate the spurious solutions we
use the fact that the relevant QN frequencies do not depend on the number of
Chebyshev polynomials being considered. An additional check of consistency is
to plot the function $\phi_{s}(u)$, which must satisfy the boundary
conditions, i.e., regular at the horizon and divergent at the spatial
infinity. Note that the problem of calculating QN frequencies does not depend
on any initial guess, as the shooting method, for example. We get directly the
frequencies using, for instance, Mathematica’s built-in function Eigenvalues,
or Eigensystem. One may consider this fact as an advantage in relation to
other methods available in the literature.

The Chebyshev polynomials of the first kind $T_{k}(x)$ are defined in the
interval $x\in[-1,1]$, and have special properties [50]. Note also that the
collocation points (36) map this interval into the interval $[0,1]$. In turn,
the error associated with the pseudo-spectral method is of the order
$\mathcal{O}\left(1/N^{N}\right)$ for sufficiently smooth regular functions
[53]. For further discussions see for instance Ref. [45].

##  IV The Asymptotic Iteration Method

The Asymptotic Iteration Method (AIM) is a numerical method recently proposed
by Ciftci et. al. [46] for solving homogeneous second order ordinary
differential equations of the form

$y^{\prime\prime}(x)-\lambda_{0}(x)y^{\prime}(x)-s_{0}(x)y(x)=0,$ (44)

where primes denote derivatives with respect to to the variable $x$ (that is
defined over some interval that is not necessarily bounded),
$\lambda_{0}(x)\neq 0$ and $s_{0}(x)$ are $C_{\infty}$. These equations can be
found in many different areas of physics, such as the time independent
Schrödinger equation in Quantum Mechanics, or in General Relativity, such as
the differential equations for black hole perturbations Eq. (3) (where we can
restore the first derivative if the standard radial coordinate is used instead
of tortoise coordinates). Here we present a brief review concerning the AIM
and discuss some implementation details of such a method. The AIM is based
upon the following theorem:

###### Theorem 1

Let $\lambda_{0}$ and $s_{0}$ be functions of the variable $x\in(a,b)$ that
are $C_{\infty}$ on the same interval. The differential equation (44) has a
general solution of the form

$y(x)=\exp\left(-\int\alpha\mathrm{d}t\right)\times\\\
\left[C_{2}+C_{1}\int^{x}\exp\left(\int^{t}(\lambda_{0}(\tau)+2\alpha(\tau))\mathrm{d}\tau\right)\mathrm{d}t\right]$
(45)

if for some $n>0$ the condition

$\alpha\equiv\frac{s_{n}}{\lambda_{n}}=\frac{s_{n-1}}{\lambda_{n-1}}$ (46)

or equivalently

$\delta\equiv s_{n}\lambda_{n-1}-\lambda_{n}s_{n-1}=0$ (47)

is satisfied, where

$\displaystyle\lambda_{k}(x)\equiv$
$\displaystyle\lambda^{\prime}_{k-1}(x)+s_{k-1}(x)+\lambda_{0}(x)\lambda_{k-1}(x)$
(48) $\displaystyle s_{k}(x)\equiv$ $\displaystyle
s^{\prime}_{k-1}(x)+s_{0}\lambda_{k-1}(x)$ (49)

with $k$ being a integer that ranges from $1$ to $n$.

From now on, we shall refer to the condition expressed by Eq. (47) as the _AIM
quantization condition_. Provided that Theo. 1 is satisfied we can find both
the eigenvalues and eigenvectors of the second order ODE using, respectively,
Eq. (47) and Eq. (45). More specifically, the quasinormal modes of a perturbed
black hole will be the complex frequency values $\omega$ that satisfy Eq. (47)
for any value of $x$.

Despite being quite general, the method presents a computational difficulty
hidden in Eq.(48) and Eq. (49): Not only the definition of the $n$-th
coefficients are coupled and recursive, they also involve the derivatives of
previous entries. In practice this means that to compute the quantization
condition, Eq. (47), using $n$ iterations we end up computing the $n$-th
derivatives of $\lambda_{0}$ and $s_{0}$ multiple times. Depending on the size
of the original functions, the size and complexity of each coefficient can
quickly spiral out of control. To address these issues, Cho et. al. have
proposed in Ref. [43] to instead of computing these coefficients directly, use
a Taylor expansion of both $\lambda$ and $s$ around a point $\xi$ where the
AIM is to be perform (we remind the reader that the results are independent of
the choice of $\xi$), that is,

$\displaystyle\lambda_{n}(\xi)=$
$\displaystyle\sum_{i=0}^{\infty}c^{i}_{n}(x-\xi)^{i},$ (50) $\displaystyle
s_{n}(\xi)=$ $\displaystyle\sum_{i=0}^{\infty}d^{i}_{n}(x-\xi)^{i},$ (51)

where $c^{i}_{n}$ and $d^{i}_{n}$ are the Taylor coefficients of the
expansions of $\lambda_{n}$ and $s_{n}$ around $\xi$, respectively. By
plugging Eqs. (50) and (51) into Eqs. (48) and Eq. (49) one gets

$\displaystyle c^{i}_{n}=$
$\displaystyle(i+1)c^{i+1}_{n-1}+d^{i}_{n-1}+\sum_{k=0}^{i}c^{k}_{0}c^{i-k}_{n-1},$
(52) $\displaystyle d^{i}_{n}=$
$\displaystyle(i+1)d^{i+1}_{n-1}+\sum_{k=0}^{i}d^{k}_{0}c^{i-k}_{n-1}.$ (53)

Finally, using Eqs. (52) and (53) the quantization condition, Eq. (47),
becomes

$\delta\equiv d^{0}_{n}c^{0}_{n-1}-d^{0}_{n-1}c^{0}_{n}=0.$ (54)

In order to better visualize and understand the improved algorithm, it is
useful to arrange the $c^{i}_{n}$ (or $d^{i}_{n}$) coefficients as elements
$c_{i,n}$ (or $d_{i,n}$) of a matrix $C$ (or $D$), where the index $i$
indicates the matrix row and the index $n$ represents the matrix column, i.e.,

$C=\begin{pmatrix}c_{0,0}&c_{0,1}&\cdots&c_{0,n-1}&c_{0,n}\\\
c_{1,0}&c_{1,1}&\cdots&c_{1,n-1}&c_{1,n}\\\
\vdots&\vdots&\ddots&\vdots&\vdots\\\
c_{i,0}&c_{i,1}&\cdots&c_{i,n-1}&c_{i,n}\end{pmatrix}.$

Notice that, when using Eq. (54), only the last and the before last top
elements of the matrix (from left to right) are actually required, i.e., only
the elements $c_{0,n-1}$ and $c_{0,n}$ are necessary. Note also that, by using
Eqs. (52) and (53), to compute an element in row $i$ and column $n$, one needs
to have previously computed the column $n-1$ up to at least row $i+1$. In
practice this means that the matrix $C$ “grows diagonally” and in order to
compute $n$ columns one needs at least $i=n$ rows. This formulation motivates
us to view Eqs. (52) and (53) not as recursion relations, but as recipes for
iteration, i.e., given the first column of the matrix $C$, one can use Eqs.
(52) and (53) rewritten as

$\displaystyle c^{i}_{n+1}=$
$\displaystyle(i+1)c^{i+1}_{n}+d^{i}_{n}+\sum_{k=0}^{i}c^{k}_{0}c^{i-k}_{n},$
(55) $\displaystyle d^{i}_{n+1}=$
$\displaystyle(i+1)d^{i+1}_{n}+\sum_{k=0}^{i}d^{k}_{0}c^{i-k}_{n},$ (56)

to compute the next column of the matrix.

With this insight, we can now devise an algorithm that performs $n$ iterations
of the AIM. Remember that if $n$ iterations are to be performed, one needs at
least $i=n$ rows of coefficients and thus we shall truncate the Taylor
expansions at $i=n$. The algorithm steps are the following:

  1. 1.

Construct two arrays of size $n$ where the $i$-th element is $c^{i}_{0}$ (or
$d^{i}_{0}$) where $i$ ranges from zero to $n$. We shall call these icda
(initial $c$ data array) and idda (initial $d$ data array).

  2. 2.

Construct two arrays of size $n$ to contain the current column of $c$ (or $d$)
indexes. We shall call these ccda (current $c$ data array) and cdda (current
$d$ data array)

  3. 3.

Construct two arrays of size $n$ to contain the previous column of $c$ (or
$d$) indexes. We shall call these pcda (previous $c$ data array) and pdda
(previous $d$ data array).

  4. 4.

Initialize ccda with data from icda and cdda with data from idda.

  5. 5.

Perform $n$ AIM steps using the evolution Eqs. (55) and (56). That is, repeat
the following $n$ times:

    1. (a)

Copy the content from ccda into pcda

    2. (b)

Copy the content from cdda into pdda

    3. (c)

Rewrite each element of ccda and cdda using Eqs. (55) and (56), respectively.

  6. 6.

After $n$ iterations, the current and previous data array contain the sought
coefficients. Apply the quantization condition, Eq. (54), using the first
indexes of each array (as they represent the $i=0$ coefficients). Explicitly,
perform cdda[1]*pcda[1] - pdda[1]*ccda[1]111Given that Julia uses 1 based
array indexes, we are also using 1 based arrays for the algorithmic
description. This means cdda[1] refers to the first element of the cdda array
and so on so forth..

  7. 7.

Finding the roots of the resulting expression from the last step yields the
eigenvalues of the ODE (in the context of this work, the quasinormal modes).

This implementation is realized in the Julia [54] package called
QuasinormalModes.jl [55]. The implementation makes use of an additional buffer
array for each coefficient family in order to allow for thread based
parallelization to take place during the main AIM loop. All AIM numerical
results from this work were obtained with the aforementioned package.

##  V Numerical results for the QN frequencies

###  V.1 Spin 0 QN frequencies

Our numerical results for the quasinormal frequencies of the spin $0$
perturbation field are displayed in Table 1. The first two column shows the
data from the pseudo-spectral method with different numbers of interpolating
polynomials, the third column shows the results form the AIM method, while the
fourth and fifth columns are reproduction of the results from Refs. [41, 42],
respectively. Notice that in Refs. [41, 42] the WKB method was used to
calculate the quasinormal frequencies. As it can be seen, the pseudo-spectral
method I (calculated with $60$ polynomials) and pseudo-spectral method II
(calculated with $40$ polynomials) provide results which are practically the
same as those provided by the AIM, within six decimal places of precision. The
numerical methods employed here provide more accurate results than those
obtained by using the WKB approximation, and allows us to calculate additional
frequencies for spin 0 fields not reported previously in the literature.

$l$ | $n$ |  | Pseudo-spectral
---
I (60 Polynomials)
| Pseudo-spectral
---
II (40 polynomials)
| AIM
---
100 Iterations
Ref. [41] | Ref. [42]
0 | $0$ | $\pm 0.110455-0.104896i$ | $\pm 0.110455-0.104896i$ | $\pm 0.110455-0.104896i$ | $0.1046-0.1152i$ | $\pm 0.1105-0.1008i$
1 | $0$ | $\pm 0.292936-0.097660i$ | $\pm 0.292936-0.097660i$ | $\pm 0.292936-0.097660i$ | $0.2911-0.0980i$ | $\pm 0.2929-0.0978i$
| $1$ | $\pm 0.264449-0.306257i$ | $\pm 0.264449-0.306257i$ | $\pm 0.264449-0.306257i$ | — | $\pm 0.2645-0.3065i$
2 | $0$ | $\pm 0.483644-0.096759i$ | $\pm 0.483644-0.096759i$ | $\pm 0.483644-0.096759i$ | $0.4832-0.0968i$ | $\pm 0.4836-0.0968i$
| $1$ | $\pm 0.463851-0.295604i$ | $\pm 0.463851-0.295604i$ | $\pm 0.463851-0.295604i$ | $0.4632-0.2958i$ | $\pm 0.4638-0.2956i$
| $2$ | $\pm 0.430544-0.508558i$ | $\pm 0.430544-0.508558i$ | $\pm 0.430544-0.508558i$ | — | $\pm 0.4304-0.5087i$
3 | $0$ | $\pm 0.675366-0.096500i$ | $\pm 0.675366-0.096500i$ | $\pm 0.675366-0.096500i$ | $0.6752-0.0965i$ | —
| $1$ | $\pm 0.660671-0.292285i$ | $\pm 0.660671-0.292285i$ | $\pm 0.660671-0.292285i$ | $0.6604-0.2923i$ | —
| $2$ | $\pm 0.633626-0.496008i$ | $\pm 0.633626-0.496008i$ | $\pm 0.633626-0.496008i$ | $0.6348-0.4941i$ | —
| $3$ | $\pm 0.598773-0.711221i$ | $\pm 0.598773-0.711221i$ | $\pm 0.598773-0.711221i$ | — | —
4 | $0$ | $\pm 0.867416-0.096392i$ | $\pm 0.867416-0.096392i$ | $\pm 0.867416-0.096392i$ | $0.8673-0.0964i$ | —
| $1$ | $\pm 0.855808-0.290876i$ | $\pm 0.855808-0.290876i$ | $\pm 0.855808-0.290876i$ | $0.8557-0.2909i$ | —
| $2$ | $\pm 0.833692-0.490325i$ | $\pm 0.833692-0.490325i$ | $\pm 0.833692-0.490325i$ | $0.8345-0.4895i$ | —
| $3$ | $\pm 0.803288-0.697482i$ | $\pm 0.803288-0.697482i$ | $\pm 0.803288-0.697482i$ | $0.8064-0.6926i$ | —
| $4$ | $\pm 0.767733-0.914019i$ | $\pm 0.767733-0.914019i$ | $\pm 0.767733-0.914019i$ | — | —
Table 1:  Quasinormal frequencies of the spin $0$ perturbations compared
against the results of Refs. [41, 42].

###  V.2 Spin 1 QN frequencies

Our numerical results for the quasinormal frequencies for spin $1$ fields are
displayed in Table 2 compared against the results of Refs. [41, 42], where the
WKB method was used. The first two columns show the data from the pseudo-
spectral method with different numbers of interpolating polynomials, the third
column shows the results form the AIM method, while the fourth and fifth
columns are reproduction of the results from Refs. [41, 42], respectively. As
it can be seen from the table, the pseudo-spectral method I (calculated with
$60$ polynomials) and pseudo-spectral method II (calculated with $40$
polynomials) provide results which are practically identical to those provided
by the AIM. The numerical methods we are working with provide more accurate
results than those obtained by using the WKB approximation, and allows us to
present additional frequencies for spin 1 fields not reported previously in
the literature.

$l$ | $n$ |  | Pseudo-spectral
---
I (60 Polynomials)
| Pseudo-spectral
---
II (40 polynomials)
| AIM
---
100 Iterations
Ref. [41] | Ref. [42]
1 | $0$ | $\pm 0.248263-0.092488i$ | $\pm 0.248263-0.092488i$ | $\pm 0.248263-0.092488i$ | $0.2459-0.0931i$ | $\pm 0.2482-0.0926i$
| $1$ | $\pm 0.214515-0.293668i$ | $\pm 0.214515-0.293667i$ | $\pm 0.214515-0.293668i$ | — | $\pm 0.2143-0.2941i$
2 | $0$ | $\pm 0.457596-0.095004i$ | $\pm 0.457595-0.095004i$ | $\pm 0.457596-0.095004i$ | $0.4571-0.0951i$ | $\pm 0.4576-0.0950i$
| $1$ | $\pm 0.436542-0.290710i$ | $\pm 0.436542-0.290710i$ | $\pm 0.436542-0.290710i$ | $0.4358-0.2910i$ | $\pm 0.4365-0.2907i$
| $2$ | $\pm 0.401187-0.501587i$ | $\pm 0.401187-0.501587i$ | $\pm 0.401187-0.501587i$ | — | $\pm 0.4009-0.5017i$
3 | $0$ | $\pm 0.656899-0.095616i$ | $\pm 0.656899-0.095616i$ | $\pm 0.656899-0.095616i$ | $0.6567-0.0956i$ | $\pm 0.6569-0.0956i$
| $1$ | $\pm 0.641737-0.289728i$ | $\pm 0.641737-0.289728i$ | $\pm 0.641737-0.289728i$ | $0.6415-0.2898i$ | $\pm 0.6417-0.2897i$
| $2$ | $\pm 0.613832-0.492066i$ | $\pm 0.613832-0.492066i$ | $\pm 0.613832-0.492066i$ | $0.6151-0.4901i$ | $\pm 0.6138-0.4921i$
| $3$ | $\pm 0.577919-0.706331i$ | $\pm 0.577919-0.706331i$ | $\pm 0.577919-0.706330i$ | — | $\pm 0.5775-0.7065i$
4 | $0$ | $\pm 0.853095-0.095860i$ | $\pm 0.853095-0.095860i$ | $\pm 0.853095-0.095810i$ | $0.8530-0.0959i$ | —
| $1$ | $\pm 0.841267-0.289315i$ | $\pm 0.841267-0.289315i$ | $\pm 0.841267-0.289315i$ | $0.8411-0.2893i$ | —
| $2$ | $\pm 0.818728-0.487838i$ | $\pm 0.818728-0.487838i$ | $\pm 0.818728-0.487838i$ | $0.8196-0.4870i$ | —
| $3$ | $\pm 0.787748-0.694242i$ | $\pm 0.787748-0.694242i$ | $\pm 0.787748-0.694242i$ | $0.7909-0.6892i$ | —
| $4$ | $\pm 0.751549-0.910242i$ | $\pm 0.751549-0.910242i$ | $\pm 0.751549-0.910242i$ | — | —
Table 2:  Quasinormal frequencies of the spin $1$ perturbations compared
against the results of Refs. [41, 42].

###  V.3 Spin 2 QN frequencies

Our numerical results for the quasinormal frequencies for spin $2$ fields are
displayed in Table 3 compared against the results of Refs. [41, 42], where the
WKB method was employed. The first two columns show the data from the pseudo-
spectral method with different numbers of interpolating polynomials, the third
column shows the results form the AIM method, while the fourth and fifth
columns are reproduction of the results from Refs. [41, 42], respectively. As
it can be seen from the table, the pseudo-spectral method I and pseudo-
spectral method II provide results which are practically equal to those
provided by the AIM. The numerical results we are working with provide more
accurate results than those obtained by using the WKB approximation. It is
worth mentioning that, in this case, we obtain additional solutions to the
eigenvalue problem which do not represent gravitational waves, see the
discussion in Appendix B for more details on this point.

$l$ | $n$ |  | Pseudo-spectral
---
I (60 Polynomials)
| Pseudo-spectral
---
II (40 polynomials)
| AIM
---
100 Iterations
Ref. [41] | Ref. [42]
2 | $0$ | $\pm 0.373672-0.088962i$ | $\pm 0.373672-0.088962i$ | $\pm 0.373672-0.088962i$ | $0.3730-0.0891i$ | $\pm 0.3736-0.0890i$
| $1$ | $\pm 0.346711-0.273915i$ | $\pm 0.346711-0.273915i$ | $\pm 0.346711-0.273915i$ | $0.3452-0.2746i$ | $\pm 0.3463-0.2735i$
| $2$ | $\pm 0.301053-0.478277i$ | $\pm 0.301053-0.478277i$ | $\pm 0.301053-0.478277i$ | — | $\pm 0.2985-0.4776i$
3 | $0$ | $\pm 0.599443-0.092703i$ | $\pm 0.599443-0.092703i$ | $\pm 0.599443-0.092703i$ | $0.5993-0.0927i$ | $\pm 0.5994-0.0927i$
| $1$ | $\pm 0.582644-0.281298i$ | $\pm 0.582644-0.281298i$ | $\pm 0.582644-0.281298i$ | $0.5824-0.2814i$ | $\pm 0.5826-0.2813i$
| $2$ | $\pm 0.551685-0.479093i$ | $\pm 0.551685-0.479093i$ | $\pm 0.551685-0.479027i$ | $0.5532-0.4767i$ | $\pm 0.5516-0.4790i$
| $3$ | $\pm 0.511962-0.690337i$ | $\pm 0.511962-0.690337i$ | $\pm 0.511962-0.690337i$ | — | $\pm 0.5111-0.6905i$
4 | $0$ | $\pm 0.809178-0.094164i$ | $\pm 0.809178-0.094164i$ | $\pm 0.809178-0.094164i$ | $0.8091-0.0942i$ | $\pm 0.8092-0.0942i$
| $1$ | $\pm 0.796632-0.284334i$ | $\pm 0.796632-0.284334i$ | $\pm 0.796632-0.284334i$ | $0.7965-0.2844i$ | $\pm 0.7966-0.2843i$
| $2$ | $\pm 0.772710-0.479908i$ | $\pm 0.772710-0.479908i$ | $\pm 0.772710-0.479908i$ | $0.7736-0.4790i$ | $\pm 0.7727-0.4799i$
| $3$ | $\pm 0.739837-0.683924i$ | $\pm 0.739837-0.683924i$ | $\pm 0.739837-0.683924i$ | $0.7433-0.6783i$ | $\pm 0.7397-0.6839i$
| $4$ | $\pm 0.701516-0.898239i$ | $\pm 0.701516-0.898239i$ | $\pm 0.701516-0.898239i$ | — | $\pm 0.7006-0.8985i$
Table 3:  Quasinormal frequencies of spin $2$ perturbations compared against
the results of Refs. [42, 41].

###  V.4 Spin $1/2$ QN frequencies

Our numerical results for the quasinormal frequencies for spin $1/2$ fields
are displayed in Table 4 compared against results available in the literature.
The first two columns show the data from the pseudo-spectral method with
different numbers of interpolating polynomials, the third column shows the
results form the AIM method, while the fourth and fifth columns are
reproduction of the results from Refs. [40, 41], respectively. As it can be
seen, the results obtained using the pseudo-spectral I and II are in perfect
agreement with the results obtained using the AIM within the decimal places
considered. Note that the numerical methods we are working with provide more
accurate results than the results reported in Refs. [40, 41], where the
authors employed the WKB approximation. Note that we also show additional
frequencies not reported previously in the literature, for example for
$\ell=1$ and $n=1$.

It is worth pointing out that we also found purely imaginary frequencies, that
arise when investigating the quasinormal modes in the limit of large $\ell$.
Our numerical results are displayed in Table 5, where we show the first five
purely imaginary frequencies. As it is seen from the table, the agreement
between the numerical methods is perfect for low overtones but it gets worse
for higher overtones. It is worth mentioning that these results are in perfect
agreement with the analytic solution, $M\omega=-i{n}/{4},\quad n\to\infty$,
obtained in Refs. [56, 57], see also references therein.

$l$ | $n$ |  | Pseudo-spectral
---
I (60 Polynomials)
| Pseudo-spectral
---
II (40 polynomials)
| AIM
---
100 Iterations
Ref. [41] | Ref. [40]
0 | $0$ | $\pm 0.182963-0.096982i$ | $\pm 0.182963-0.096982i$ | $\pm 0.182963-0.096824i$ | — | —
1 | $0$ | $\pm 0.380037-0.096405i$ | $\pm 0.380037-0.096405i$ | $\pm 0.380037-0.096405i$ | $0.3786-0.0965i$ | $0.379-0.097i$
| $1$ | $\pm 0.355833-0.297497i$ | $\pm 0.355833-0.297497i$ | $\pm 0.355833-0.297497i$ | — | —
2 | $0$ | $\pm 0.574094-0.096305i$ | $\pm 0.574094-0.096305i$ | $\pm 0.574094-0.096305i$ | $0.5737-0.0963i$ | $0.574-0.096i$
| $1$ | $\pm 0.557015-0.292715i$ | $\pm 0.557015-0.292715i$ | $\pm 0.557015-0.292715i$ | $0.5562-0.2930i$ | $0.556-0.293i$
| $2$ | $\pm 0.526607-0.499695i$ | $\pm 0.526607-0.499695i$ | $\pm 0.526607-0.499695i$ | — | —
3 | $0$ | $\pm 0.767355-0.096270i$ | $\pm 0.767355-0.096270i$ | $\pm 0.767355-0.096270i$ | $0.7672-0.0963i$ | $0.767-0.096i$
| $1$ | $\pm 0.754300-0.290968i$ | $\pm 0.754300-0.290968i$ | $\pm 0.754300-0.290968i$ | $0.7540-0.2910i$ | $0.754-0.291i$
| $2$ | $\pm 0.729770-0.491910i$ | $\pm 0.729770-0.491910i$ | $\pm 0.729770-0.491910i$ | $0.7304-0.4909i$ | $0.730-0.491i$
| $3$ | $\pm 0.696913-0.702293i$ | $\pm 0.696913-0.702293i$ | $\pm 0.696913-0.702293i$ | — | —
4 | $0$ | $\pm 0.960293-0.096254i$ | $\pm 0.960293-0.096254i$ | $\pm 0.960293-0.096254i$ | $0.9602-0.0963i$ | $0.960-0.096i$
| $1$ | $\pm 0.949759-0.290148i$ | $\pm 0.949759-0.290148i$ | $\pm 0.949759-0.290148i$ | $0.9496-0.2902i$ | $0.950-0.290i$
| $2$ | $\pm 0.929494-0.488116i$ | $\pm 0.929494-0.488116i$ | $\pm 0.929494-0.488116i$ | $0.9300-0.4876i$ | $0.930-0.488i$
| $3$ | $\pm 0.901129-0.692520i$ | $\pm 0.901129-0.692520i$ | $\pm 0.901129-0.692520i$ | $0.9036-0.6892i$ | $0.904-0.689i$
| $4$ | $\pm 0.867043-0.905047i$ | $\pm 0.867008-0.905066i$ | $\pm 0.867043-0.905047i$ | — | —
Table 4:  Quasinormal frequencies of the spin $1/2$ perturbations compared against the results of Refs. [40, 41].  | Pseudo-spectral
---
I (60 Polynomials)
| Pseudo-spectral
---
II (40 polynomials)
| AIM
---
100 Iterations
$-0.250000i$ | $-0.250000i$ | $-0.250000i$
$-0.500000i$ | $-0.500000i$ | $-0.500000i$
$-0.750000i$ | $-0.750000i$ | $-0.750000i$
$-1.000000i$ | $-1.000000i$ | $-1.000031i$
$-1.2499998i$ | $-1.250000i$ | $-1.246550i$
Table 5:  Purely imaginary frequencies for spin $1/2$ perturbations. The
numerical values of such frequencies are exactly the same as for the purely
imaginary frequencies arising in the QNM of spin 3/2 perturbations.

###  V.5 Spin $3/2$ QN frequencies

Our numerical results for the quasinormal frequencies for spin $3/2$ fields
are displayed in Table 6 compared against results available in the literature.
The first two columns show the data from the pseudo-spectral method with
different numbers of interpolating polynomials, the third column shows the
results form the AIM method, while the fourth and fifth columns are
reproduction of the results from Refs. [41, 48], respectively. As it can be
seen, the results obtained by using the pseudo-spectral I and II are in
perfect agreement with the results obtained using the AIM. We also realize
that these results are in very good agreement with the results reported in
Ref. [48], where the authors also employed the IAM, and with the results from
Ref. [41], where the authors employed the WKB approximation.

As in the spin $1/2$ field perturbations, we also find purely imaginary
frequencies for the spin $3/2$ field. The numerical results are displayed in
Table 5 for the three routines we are working with. Such frequencies arise
when investigating the quasinormal modes in the limit of large $\ell$. Notice
that these results are also in agreement with the analytic solutions obtained
in Refs. [56, 57]. It is worth pointing out that the numerical values of these
purely imaginary frequencies are exactly the same for spin $1/2$ and $3/2$
fields. We do not have an explanation for this fact, maybe it is just a
coincidence. Notice also that these frequencies can be written as fractions,
multiples of $1/4$, i.e., $1/4$, $2/4$, $3/4$, $4/4$, $\sim 5/4$, $\cdots$.

$l$ | $n$ |  | Pseudo-spectral
---
I (60 Polynomials)
| Pseudo-spectral
---
II (40 polynomials)
| AIM
---
100 Iterations
Ref. [41] | Ref. [48]
0 | $0$ | $\pm 0.311292-0.090087i$ | $\pm 0.311292-0.090087i$ | $\pm 0.311292-0.090087i$ | — | $0.3112-0.0902i$
1 | $0$ | $\pm 0.530048-0.093751i$ | $\pm 0.530048-0.093751i$ | $\pm 0.530048-0.093751i$ | — | $0.5300-0.0937i$
| $1$ | $\pm 0.511392-0.285423i$ | $\pm 0.511392-0.285423i$ | $\pm 0.511392-0.285423i$ | — | $0.5113-0.2854i$
2 | $0$ | $\pm 0.734750-0.094878i$ | $\pm 0.734750-0.094878i$ | $\pm 0.734750-0.094878i$ | $\pm 0.7346-0.0949i$ | $0.7347-0.0948i$
| $1$ | $\pm 0.721047-0.286906i$ | $\pm 0.721047-0.286906i$ | $\pm 0.721047-0.286906i$ | $\pm 0.7206-0.2870i$ | $0.7210-0.2869i$
| $2$ | $\pm 0.695287-0.485524i$ | $\pm 0.695287-0.485524i$ | $\pm 0.695287-0.485524i$ | — | $0.6952-0.4855i$
3 | $0$ | $\pm 0.934364-0.095376i$ | $\pm 0.934364-0.095376i$ | $\pm 0.934364-0.095376i$ | $\pm 0.9343-0.0954i$ | $0.9343-0.0953i$
| $1$ | $\pm 0.923502-0.287560i$ | $\pm 0.923502-0.287560i$ | $\pm 0.923502-0.287560i$ | $\pm 0.9233-0.2876i$ | $0.9235-0.2875i$
| $2$ | $\pm 0.902599-0.483957i$ | $\pm 0.902599-0.483957i$ | $\pm 0.902599-0.483957i$ | $\pm 0.9031-0.4835i$ | $0.9025-0.4839i$
| $3$ | $\pm 0.873342-0.687024i$ | $\pm 0.873343-0.687024i$ | $\pm 0.873342-0.687024i$ | — | $0.8732-0.6870i$
4 | $0$ | $\pm 1.131530-0.095640i$ | $\pm 1.131530-0.095640i$ | $\pm 1.131530-0.095640i$ | $\pm 1.1315-0.0956i$ | $1.1315-0.0956i$
| $1$ | $\pm 1.122523-0.287908i$ | $\pm 1.122523-0.287908i$ | $\pm 1.122523-0.287908i$ | $\pm 1.1224-0.2879i$ | $1.1225-0.2879i$
| $2$ | $\pm 1.104976-0.483096i$ | $\pm 1.104976-0.483096i$ | $\pm 1.104976-0.483096i$ | $\pm 1.1053-0.4828i$ | $1.1049-0.4830i$
| $3$ | $\pm 1.079852-0.683000i$ | $\pm 1.079852-0.683000i$ | $\pm 1.079852-0.683000i$ | $\pm 1.0817-0.6812i$ | $1.0798-0.6829i$
| $4$ | $\pm 1.048599-0.889113i$ | $\pm 1.048596-0.889115i$ | $\pm 1.048599-0.889113i$ | — | $1.0484-0.8890i$
Table 6:  Quasinormal frequencies of spin $3/2$ perturbations compared against
the results of Refs. [48, 41].

###  V.6 Spin $5/2$ QN frequencies

Our numerical results for the quasinormal frequencies for spin $5/2$ fields
are displayed in Table 7. As it can be seen, the results obtained by using the
pseudo-spectral I and II are in perfect agreement with the results obtained
using the AIM within the decimal places considered. It is worth mentioning
that we present here the quasinormal frequencies for spin $5/2$ perturbation
fields for the very first time

In turn, we also found purely imaginary frequencies. Our numerical results are
displayed in Table 8. It turns out that these results satisfy the sequence
$1/8$, $\sim 3/8$, $\sim 5/8$, $\sim 7/8$, $\cdots$, in general,

$M\omega=-i\frac{(2n+1)}{8},\qquad n=0,1,2,3,\cdots.$ (57)

As in spin $1/2$ and $3/2$ perturbations, these frequencies arise when we were
investigating the frequencies in the limit of large $\ell$. As it can be seen,
the discrepancy of both methods increases as the imaginary frequency gets
negative.

Finally, differently from the quasinormal frequencies with real and imaginary
parts, the purely imaginary frequencies represent damping solutions, while the
frequencies with real and imaginary parts represent oscillatory solutions
being damped by the imaginary part of the frequency. This is easily understood
because the solution goes like $\sim
e^{i\omega\,t}=e^{i(\omega_{Re}-i\,\omega_{Im})\,t}=e^{\omega_{Im}\,t}\cos{\left(\omega_{Re}\,t\right)}$.

$l$ | $n$ |  | Pseudo-spectral
---
I (60 Polynomials)
| Pseudo-spectral
---
II (40 polynomials)
| AIM
---
100 Iterations
0 | $0$ | $\pm 0.462727-0.092578i$ | $\pm 0.462727-0.092578i$ | $0.462727-0.092577i$
1 | $0$ | $\pm 0.687103-0.094566i$ | $\pm 0.687103-0.094566i$ | $0.687103-0.094566i$
| $1$ | $\pm 0.670542-0.285767i$ | $\pm 0.670542-0.285767i$ | $0.670542-0.285767i$
2 | $0$ | $\pm 0.897345-0.095309i$ | $\pm 0.897345-0.095309i$ | $0.897345-0.095309i$
| $1$ | $\pm 0.884980-0.287266i$ | $\pm 0.884980-0.287266i$ | $0.884980-0.287266i$
| $2$ | $\pm 0.861109-0.483113i$ | $\pm 0.861109-0.483113i$ | $0.861109-0.483113i$
3 | $0$ | $\pm 1.101190-0.095648i$ | $\pm 1.101190-0.095648i$ | $1.101190-0.095648i$
| $1$ | $\pm 1.091300-0.287886i$ | $\pm 1.091300-0.287886i$ | $1.091300-0.287886i$
| $2$ | $\pm 1.071999-0.482895i$ | $\pm 1.071999-0.482895i$ | $1.071999-0.482895i$
| $3$ | $\pm 1.044272-0.682307i$ | $\pm 1.044272-0.682307i$ | $1.044272-0.682307i$
4 | $0$ | $\pm 1.301587-0.095829i$ | $\pm 1.301587-0.095829i$ | $1.301587-0.095829i$
| $1$ | $\pm 1.293328-0.288184i$ | $\pm 1.293328-0.288184i$ | $1.293328-0.288184i$
| $2$ | $\pm 1.277107-0.482604i$ | $\pm 1.277107-0.482604i$ | $1.277107-0.482604i$
| $3$ | $\pm 1.253526-0.680366i$ | $\pm 1.253526-0.680366i$ | $1.253526-0.680366i$
| $4$ | $\pm 1.223513-0.882554i$ | $\pm 1.223512-0.882553i$ | $1.223513-0.882554i$
Table 7:  Quasinormal frequencies of spin $5/2$ perturbations.  | Pseudo-spectral
---
I (60 Polynomials)
| Pseudo-spectral
---
II (40 polynomials)
| AIM
---
100 Iterations
$-0.125000i$ | $-0.125000i$ | $-0.125000i$
$-0.375602i$ | $-0.375602i$ | $-0.378659i$
$-0.626877i$ | $-0.626877i$ | $-0.623931i$
$-0.878946i$ | $-0.878948i$ | $-0.907374i$
Table 8:  Purely imaginary frequencies of spin $5/2$ perturbations.

##  VI QNMs for $\ell\gg 1$ and $n\gg 1$

It is also interesting to calculate the quasinormal frequencies in the limit
of large $\ell$, where analytic solutions are available in the literature to
compare with. The analytic solutions we obtained in Ref. [58], such that the
real and imaginary parts of the frequencies are given by

$\begin{split}M\omega_{\text{Re}}=\,&\frac{1}{3\sqrt{3}}\left(\ell+\frac{1}{2}\right),\\\
M\omega_{\text{Im}}=\,&\frac{1}{3\sqrt{3}}\left(n+\frac{1}{2}\right).\end{split}$
(58)

Before comparing our numerical results against the analytic solutions, it is
worth mentioning that these analytic solutions were obtained for integer spin
perturbations. We do not expected that these results can be applied for semi-
integer field perturbations, in principle.

Let us now compare our numerical results against the analytic solutions (58).
Our numerical results for $s=0$ field are displayed in Figs. 1 and 2, as a
function of $\ell$ and $n$, respectively. As it can be seen, the real part as
a function of $\ell$ fits very well the analytic result (58). Meanwhile, the
imaginary part of the frequency as a function of $n$ also fits very well the
analytic solution.

Figure 1: The figure shows the real and imaginary parts of the frequency as a
function of $\ell$ for $n=0$ considering $s=0$ field. Dagger represents the
AIM results, while asterisk the pseudo-spectral and diamond analytic results
of Eq. (58).

Figure 2: The figure shows the real and imaginary parts of the frequency as a
function of $n$ for $\ell=30$ considering $s=0$ field. Dagger represents the
AIM results, while asterisk the pseudo-spectral and diamond analytic results
of Eq. (58).

It is worth pointing out that the numerical results obtained using the pseudo-
spectral and AIM are in agreement as seen in these figures. To observe the
numerical difference we plotted the difference of the frequencies
$|M\omega^{\text{AIM}}_{\text{Re}}-M\omega^{\text{PS}}_{\text{Re}}|$ and
$|M\omega^{\text{AIM}}_{\text{Im}}-M\omega^{\text{PS}}_{\text{Im}}|$ in
logarithmic scale and displayed the results in Fig. 3. It is seen that the
difference between the numerical results is very small and slightly increases
with the increasing of $n$.

Figure 3: The figure shows the difference of the real and imaginary parts of
the frequency in logarithmic scale as a function of $n$ obtained in both, the
AIM and Pseudo-Scalar methods for spin 0 perturbation.

In turn, our numerical results for spin $1/2$ field are displayed in Fig 4. As
it can be seen, the numerical and analytic results are in agreement. This
means that the analytic results are also valid for perturbations for spin
$1/2$.

Figure 4: The figure shows the real and imaginary parts of the frequency as a
function of $n$ for $\ell=30$ considering $s=1/2$ field. Dagger represents the
AIM results, while asterisk the pseudo-spectral and diamond analytic results
of Eq. (58).

##  VII Discussion and Conclusions

In this paper we calculated the quasinormal frequencies for spin $0$, $1/2$,
$1$, $3/2$, $2$ and $5/2$ fields on the Schwarzschild black hole in
asymptotically flat space-time. We have employed the pseudo-spectral and AIM
methods. The main difference between these methods lies in the way they solve
the eigenvalue problem. While the pseudo-spectral method expand the solution
in a base of cardinal functions, the AIM calculates the roots of a
characteristic polynomial. In doing so, the results obtained for a spin zero
field applying the pseudo-spectral method are in good agreement with the
results obtained by applying the AIM and also in good agreement with
literature values.

We displayed results with six decimal places in Tables 1 \- 3. The results
show that both methods are in agreement at least up to the sixth decimal
place. We point out that the quasinormal frequencies obtained for spin $0$,
$1$, and $2$ fields do not bring any additional information to what we known
from the literature. In turn, for spin $1/2$ and $3/2$ fields our results
obtained for both methods are in agreement, these results are also in
agreement with results available in the literature. For these fields, we have
obtained additional frequencies which are purely imaginary, see Table 5. We
also calculated the quasinormal frequencies for a spin $5/2$ field, which are
displayed in Table 7. We have observed that the pseudo-spectral method and the
AIM show good agreement. Additional frequencies, which are purely imaginary,
were also obtained and displayed in Table 8. Hence, the pseudo-spectral method
calculated a set of quasinormal frequencies, while the AIM used a guess to
start looking for solutions. Our main conclusion is that both methods
complemented each other in the task for calculating the QNM frequencies.

Additional comments concerning the numerical analysis may be of interest. The
use of $60$ polynomials for the pseudo-spectral method I and $40$ polynomials
for the pseudo-spectral method II are the minimum necessary to reach the
precision one desires in the calculation, say 30 decimal places in such a
case. An additional issue concerning the pseudo-spectral method is related to
the reason why such a method does not yield, for instance, the 30 expected
eigen-frequencies for $\ell=30$. We tried increasing the number of polynomials
but, with fixed precision, were unable to go beyond 16 solutions. We realized
that increasing the number of polynomials and the precision we are able to
reach 30 eigen-frequencies. Similar difficulties are also present in the AIM.

In the next stage of the present work we are going to address the problem of
calculating QNMs for massive test fields, extending the analysis to other
black hole solutions like Reissnerd-Nordström, Kerr and Kerr-Newman black
holes, as well as considering the contribution of the cosmological constant.

###### Acknowledgements.

We would like acknowledge Alex S. Miranda for useful discussions along the
development of this work. L. T. S. and A. D. D. Masa are partly founded by
Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), Brazil,
Finance Code 001. L. A. H. M. is partially founded by the Universidade
Estadual da Região Tocantina do Maranhão (UEMASUL, Brazil). V. T. Z. thanks
financial support from Conselho Nacional de Desenvolvimento Científico e
Tecnológico (CNPq), Brazil, Grant No. 309609/2018-6, and from CAPES, Brazil,
Grant No. 88887.310351/2018-00.

##  Appendix A Fundamental equations for the spin $5/2$ perturbations

In this appendix we write additional details of the spin $5/2$ field. First,
we write the effective potential $V_{\scriptscriptstyle{5/2}}(u)$ that appears
in the Schrödinger-like equation when written in terms of the variable
$u=1/r$. It assumes the form

$\begin{split}V_{\scriptscriptstyle{5/2}}&=\frac{u^{2}f}{4S(u)}\bigg{\\{}4\bigg{[}75L_{2}\left(23+3L_{3}\right)u+9L_{2}^{3}u^{6}+L_{2}^{2}\left[25L_{4}-9u^{7}\Big{(}9+L_{3}\\{3+u^{2}-3u\\}+u\left\\{3+6u^{2}-13u\right\\}\Big{)}\right]\\\
&-450u^{2}\left[2u-\left(3+\ell\right)^{2}\right]\bigg{]}-15u^{2}f^{1/2}\bigg{[}2u\left(35\sqrt{7+L_{3}}-221+3L_{3}^{2}\left(\sqrt{7+L_{3}}-6u-16\right)+4L_{3}^{3}(u-1)\right.\\\
&\left.+14u\left(5\sqrt{7+L_{3}}-98\right)+L_{3}\left(22\sqrt{7+L_{3}}-189-456u+14u\sqrt{7+L_{3}}\right)\right)\\\
&+L_{2}\Big{[}12-20\sqrt{7+L_{3}}+9u^{3}(33+46u)+L_{3}\left(141u^{3}-4\sqrt{7+L_{3}}\right)\Big{]}\bigg{]}\bigg{\\}},\end{split}$
(59)

where, to simplify the expressions, we introduced the notation

$\displaystyle S(u)=\,$
$\displaystyle\left[5\left(L_{2}+6u\right)+3L_{2}u^{3}f^{3/2}(u)\right]^{2},$
$\displaystyle L_{1}=\,$ $\displaystyle 11+\ell(6+\ell),$ $\displaystyle
L_{2}=\,$ $\displaystyle(1+\ell)(5+\ell),$ $\displaystyle L_{3}=\,$
$\displaystyle\ell(6+\ell),$ $\displaystyle L_{4}=\,$
$\displaystyle(2+\ell)(4+\ell).$

Following the procedure explained in Section II.2, we transform the
Schrödinger-like equation into a equation suitable to apply the pseudo-
spectral method. Moreover, in order to avoid square roots we have used the
variable $\chi^{2}=1-u$, such that the final differential equation is

$R(\chi)\phi_{\scriptscriptstyle{5/2}}(\chi)+Q(\chi)\phi_{\scriptscriptstyle{5/2}}^{\prime}(\chi)+P(\chi)\phi_{\scriptscriptstyle{5/2}}^{\prime\prime}(\chi)=0,$
(60)

where $P(\chi)$, $Q(\chi)$, and $R(\chi)$ are given by

$\begin{split}R(\chi)=&\chi\left(1-\chi^{2}\right)\bigg{[}100L_{2}L_{4}-800L_{2}L_{4}\left(\chi^{2}-1\right)+60\Big{[}10(17+L_{3})+L_{2}\chi\left(\sqrt{7+L_{3}}(5+L_{3})-3\right)\Big{]}\left(\chi^{2}-1\right)^{2}\\\
&+15L_{2}\chi\left(305+133L_{3}\right)\left(\chi^{2}-1\right)^{5}+36L_{2}^{2}\left(10+3L_{3}\right)\left(\chi^{2}-1\right)^{7}+108L_{2}^{2}L_{3}\left(\chi^{2}-1\right)^{8}\\\
&+36L_{2}^{2}(L_{3}-10)(\chi^{2}-1)^{9}-180L_{2}^{2}(\chi^{2}-1)^{10}+18L_{2}(2(5+L_{3})^{2}-305\chi)(\chi^{2}-1)^{6}\\\
&-60\chi\left[35\sqrt{7+L_{3}}-636+L_{3}\left(7\sqrt{7+L_{3}}-208+L_{3}\left(2L_{3}-7\right)\right)\right]\left(\chi^{2}-1\right)^{4}\bigg{]}\\\
&-30\chi\left(1-\chi^{2}\right)^{4}\Big{(}240+\chi\left[35\sqrt{7+L_{3}}-221-L_{3}\left(189-22\sqrt{7+L_{3}}+L_{3}\Big{\\{}48-3\sqrt{7+L_{3}}+4L_{3}\Big{\\}}\right)\right]\Big{)}\\\
&-16\lambda\chi\Big{[}i+4\lambda\left(2-3\chi^{2}+\chi^{4}\right)\Big{]}\,t(\chi),\\\
Q(\chi)=&\,\left(\chi^{2}-1\right)\Big{[}1+\chi^{2}\left(\chi^{2}-2\right)\left(1-16i\lambda\right)-8i\lambda\Big{]}\,t(\chi)\\\
P(\chi)=&\,\chi\big{(}\chi^{2}-1\big{)}^{3}\,t(\chi),\end{split}$ (61)

where $t(\chi)$ is an auxiliary coefficient given by
$t(\chi)=\left[30\left(\chi^{2}-1\right)+L_{2}\left(3\chi^{9}-9\chi^{7}+9\chi^{5}-3\chi^{3}-5\right)\right]^{2}$.

##  Appendix B Further comments on the numerical QN frequencies

In this appendix we comment on additional details found when analysing the
quasinormal frequencies of the Schwarzschild black hole by means pseudo the
AIM and pseudo-spectral methods for integer spin perturbation fields $s=2$ and
$s=1$.

For completeness, we started our search for solutions to the eigenvalue
problem for $s=2$ by setting $\ell=0$ in Eq. (17), we do not get any solution.
Then, we set $\ell=1$ and we get the solutions $M\omega_{0}=\pm
0.110455-0.104896i$ and $M\omega_{1}=\pm 0.086158-0.348079i$, corresponding to
$n=0$ and $n=1$, respectively. It is also interesting to point out that the
solution for $\ell=1$ and $n=0$ arises in both methods employed in this work,
while the solution for $\ell=1$ and $n=1$ arises in the asymptotic iteration
method alone. As a consistency check, we solved the same problem by employing
the Leaver continued fraction method [39] and obtained the same results as
from the AIM. However, as investigated by Regge-Wheeler [59], and Zerilli [59]
modes with $\ell=1$ represent an addition of angular momentum to the metric
given by Eq. (1). Hence, such modes do not generate gravitational waves.

On the other hand, investigating the electromagnetic perturbation, $s=1$, for
$\ell=0$ we obtained

$M\omega_{n}=i\frac{(n+1)}{4},\qquad n=0,1,2,\cdots$ (62)

in both methods.

## References

  * Goldstein [1980] H. Goldstein,  _Classical Mechanics_ (Addison-Wesley, 1980).
  * Abbott  _et al._ [2017] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 119, 161101 (2017), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ [2016] B. Abbott, R. Abbott, T. Abbott, M. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. Adhikari, and et al., Physical Review Letters 116, 10.1103/physrevlett.116.221101 (2016).
  * Berti  _et al._ [2009] E. Berti, V. Cardoso, and A. O. Starinets, Classical and Quantum Gravity 26, 163001 (2009).
  * Konoplya and Zhidenko [2011] R. A. Konoplya and A. Zhidenko, Rev. Mod. Phys. 83, 793 (2011), arXiv:1102.4014 [gr-qc] .
  * York [1983] J. W. York, Jr., Phys. Rev. D 28, 2929 (1983).
  * Ferrari and Mashhoon [1984] V. Ferrari and B. Mashhoon, Phys. Rev. D 30, 295 (1984).
  * Leaver [1986] E. W. Leaver, Phys. Rev. D 34, 384 (1986).
  * Nollert and Schmidt [1992] H.-P. Nollert and B. G. Schmidt, Phys. Rev. D 45, 2617 (1992).
  * Andersson and Linnæus [1992] N. Andersson and S. Linnæus, Phys. Rev. D 46, 4179 (1992).
  * Nollert [1993] H.-P. Nollert, Phys. Rev. D 47, 5253 (1993).
  * Kokkotas and Schmidt [1999] K. D. Kokkotas and B. G. Schmidt, Living Rev. Rel. 2, 2 (1999), arXiv:gr-qc/9909058 .
  * Horowitz and Hubeny [2000] G. T. Horowitz and V. E. Hubeny, Phys. Rev. D 62, 024027 (2000), arXiv:hep-th/9909056 .
  * Nollert [1999] H.-P. Nollert, Class. Quant. Grav. 16, R159 (1999).
  * Cardoso and Lemos [2001a] V. Cardoso and J. P. S. Lemos, Phys. Rev. D 63, 124015 (2001a), arXiv:gr-qc/0101052 .
  * Cardoso and Lemos [2001b] V. Cardoso and J. P. S. Lemos, Phys. Rev. D 64, 084017 (2001b), arXiv:gr-qc/0105103 .
  * Konoplya [2002] R. A. Konoplya, Phys. Rev. D 66, 044009 (2002), arXiv:hep-th/0205142 .
  * Starinets [2002] A. O. Starinets, Phys. Rev. D 66, 124013 (2002), arXiv:hep-th/0207133 .
  * Clarkson and Barrett [2003] C. A. Clarkson and R. K. Barrett, Class. Quant. Grav. 20, 3855 (2003), arXiv:gr-qc/0209051 .
  * Kovtun and Starinets [2005] P. K. Kovtun and A. O. Starinets, Phys. Rev. D 72, 086009 (2005), arXiv:hep-th/0506184 .
  * Cardoso  _et al._ [2009] V. Cardoso, A. S. Miranda, E. Berti, H. Witek, and V. T. Zanchin, Phys. Rev. D 79, 064016 (2009), arXiv:0812.1806 [hep-th] .
  * Miranda  _et al._ [2008] A. S. Miranda, J. Morgan, and V. T. Zanchin, JHEP 2008 (11), 030, arXiv:0809.0297 [hep-th] .
  * Morgan  _et al._ [2009] J. Morgan, V. Cardoso, A. S. Miranda, C. Molina, and V. T. Zanchin, JHEP 2009 (09), 117, arXiv:0907.5011 [hep-th] .
  * Miranda  _et al._ [2009] A. S. Miranda, C. A. Ballon Bayona, H. Boschi-Filho, and N. R. F. Braga, JHEP 2009 (11), 119, arXiv:0909.1790 [hep-th] .
  * Mamani  _et al._ [2014] L. A. H. Mamani, A. S. Miranda, H. Boschi-Filho, and N. R. F. Braga, JHEP 2014 (3), 1, arXiv:1312.3815 [hep-th] .
  * Mamani  _et al._ [2018] L. A. H. Mamani, J. Morgan, A. S. Miranda, and V. T. Zanchin, Phys. Rev. D 98, 026006 (2018), arXiv:1804.01544 [gr-qc] .
  * Mamani  _et al._ [2019] L. A. H. Mamani, A. S. Miranda, and V. T. Zanchin, Eur. Phys. J. C 79, 1 (2019), arXiv:1809.03508 [hep-th] .
  * Mamani  _et al._ [2022] L. A. H. Mamani, D. Hou, and N. R. F. Braga, , 1 (2022), arXiv:2204.08068 [hep-ph] .
  * Akiyama  _et al._ [2019a] K. Akiyama _et al._ (Event Horizon Telescope), Astrophys. J. Lett. 875, L1 (2019a), arXiv:1906.11238 [astro-ph.GA] .
  * Akiyama  _et al._ [2019b] K. Akiyama _et al._ (Event Horizon Telescope), Astrophys. J. Lett. 875, L6 (2019b), arXiv:1906.11243 [astro-ph.GA] .
  * Perlick  _et al._ [2015] V. Perlick, O. Y. Tsupko, and G. S. Bisnovatyi-Kogan, Phys. Rev. D 92, 104031 (2015), arXiv:1507.04217 [gr-qc] .
  * Bisnovatyi-Kogan and Tsupko [2017] G. S. Bisnovatyi-Kogan and O. Y. Tsupko, Universe 3, 57 (2017), arXiv:1905.06615 [gr-qc] .
  * Cuadros-Melgar  _et al._ [2020] B. Cuadros-Melgar, R. D. B. Fontana, and J. de Oliveira, Phys. Lett. B 811, 135966 (2020), arXiv:2005.09761 [gr-qc] .
  * Jusufi [2021] K. Jusufi, Gen. Rel. Grav. 53, 87 (2021), arXiv:2007.16019 [gr-qc] .
  * Jusufi [2020] K. Jusufi, Phys. Rev. D 101, 124063 (2020), arXiv:2004.04664 [gr-qc] .
  * Neves [2020a] J. C. S. Neves, Eur. Phys. J. C 80, 717 (2020a), arXiv:2005.00483 [gr-qc] .
  * Neves [2020b] J. C. S. Neves, Eur. Phys. J. C 80, 343 (2020b), arXiv:1906.11735 [gr-qc] .
  * Chandrasekhar and Detweiler [1975] S. Chandrasekhar and S. Detweiler, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 344, 441 (1975).
  * Leaver [1985] E. W. Leaver, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 402, 285 (1985).
  * Cho [2003] H. T. Cho, Phys. Rev. D 68, 024003 (2003), arXiv:gr-qc/0303078 .
  * Shu and Shen [2005] F.-W. Shu and Y.-G. Shen, Phys. Lett. B 619, 340 (2005), arXiv:gr-qc/0501098 .
  * Konoplya [2004] R. A. Konoplya, J. Phys. Stud. 8, 93 (2004).
  * Cho  _et al._ [2012] H. T. Cho, A. S. Cornell, J. Doukas, T. R. Huang, and W. Naylor, Adv. Math. Phys. 2012, 281705 (2012), arXiv:1111.5024 [gr-qc] .
  * Boyd [2001] J. P. Boyd,  _Chebyshev and Fourier Spectral Methods_ , 2nd ed., Dover Books on Mathematics (Dover Publications, Mineola, NY, 2001).
  * Jansen [2017] A. Jansen, Eur. Phys. J. Plus 132, 546 (2017), arXiv:1709.09178 [gr-qc] .
  * Ciftci  _et al._ [2003] H. Ciftci, R. L. Hall, and N. Saad, Journal of Physics A: Mathematical and General 36, 11807–11816 (2003).
  * Schwarzschild [1916] K. Schwarzschild, Sitzungsber. Preuss. Akad. Wiss. Berlin (Math. Phys.) 1916, 189 (1916), arXiv:physics/9905030 .
  * Chen  _et al._ [2016] C. H. Chen, H. T. Cho, A. S. Cornell, and G. Harmsen, Phys. Rev. D 94, 044052 (2016), arXiv:1605.05263 [gr-qc] .
  * Shklyar  _et al._ [2010] V. Shklyar, H. Lenske, and U. Mosel, Phys. Rev. C 82, 015203 (2010), arXiv:0912.3751 [hep-ph] .
  * Arfken and Weber [1995] G. B. Arfken and H. J. Weber,  _Mathematical methods for physicists; 4th ed._ (Academic Press, San Diego, CA, 1995).
  * Tisseur and Göttel [2017] F. Tisseur and S. Göttel, MIMS Preprint , 95 (2017).
  * Konoplya and Zhidenko [2005] R. A. Konoplya and A. V. Zhidenko, Phys. Lett. B 609, 377 (2005), arXiv:gr-qc/0411059 .
  * Grandclement and Novak [2009] P. Grandclement and J. Novak, Living Rev. Rel. 12, 1 (2009), arXiv:0706.2286 [gr-qc] .
  * Bezanson  _et al._ [2017] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, SIAM review 59, 65 (2017).
  * Sanches [2022] L. T. Sanches, Journal of Open Source Software 7, 4077 (2022).
  * Cho [2006] H. T. Cho, Phys. Rev. D 73, 024019 (2006), arXiv:gr-qc/0512052 .
  * Khriplovich and Ruban [2006] I. B. Khriplovich and G. Y. Ruban, Int. J. Mod. Phys. D 15, 879 (2006), arXiv:gr-qc/0511056 .
  * Ferrari [1984] B. Ferrari, Valeria; Mashhoon, Physical Review Letters vol. 52 iss. 16 52, 10.1103/PhysRevLett.52.1361 (1984).
  * Zerilli [1970] F. J. Zerilli, Physical Review D vol. 2 iss. 10 2, 10.1103/PhysRevD.2.2141 (1970).