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

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 secondorder 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 secondorder 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,  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 pseudospectral 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] where the horizon function of the Schwarzschild black hole is given by where M is the mass of the black hole, and r is the radial coordinate which, in principle, belongs to the interval r ∈ [0, ∞). 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 → ∞, 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 < ∞.
A. 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ödingerlike 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 where the potential is given by 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 ψ 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: 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 → ∞, where f (r) → 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. ψ s (r * ) = c 3 e −iωr * + c 4 e iωr * , r → ∞.
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 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 ) Thus, in terms of the radial coordinate the boundary condition at the horizon becomes (f ( In turn, the tortoise coordinate at the spatial infinite becomes while the asymptotic solution at the spatial infinite becomes 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 ≤ r < ∞. Hence, in terms of the new coordinate this region is bounded to the interval u ∈ [0, 1], and the potential becomes 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 ψ s → Φ s given by Thus, Eq. (3) becomes 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 Φ s (u) = (1 − u) α . By substituting this ansatz into (13) we get two solutions, The solution for α = 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 Φ s = u β to get the asymptotic solution close to the spatial infinite. Plugging this ansatz in (13) we get the following solution, 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, where φ s (u) is a regular function in the interval u ∈ [0, 1] by definition. Finally, the equation of motion describing spin 0, 1 and 2 perturbations is given by where we have used λ = ωM = ω/2. The final differential equation is then a quadratic eigenvalue problem in λ. It is worth also mentioning that in the limit of zero spin s → 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 A, V C, and V C, 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.  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 It is worth mentioning that we have found a typo in the definition of ∆ in [40], which must be ∆ = 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 ∈ [0, 1]. Then by setting 2M = 1 the effective potential (18) becomes 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 (20) in which the coefficients R(u), Q(u), and P (u) are given by respectively, and with λ standing for λ = M ω = ω/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 χ 2 = 1 − u. Nevertheless, the new coordinate also belongs to the interval χ ∈ [0, 1]. The differential equation (20) becomes (24) in which the coefficients R(χ), Q(χ), and P (χ) are given by In Sec. V D we solve Eq. (24) by using the pseudo-spectral and AIM methods and compare our results against the results of Refs. [40,41].

C. 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) As before, we change the radial coordinate to u = 1/r and consider 2M = 1. Thus, the potential becomes 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 where the coefficients R(χ), Q(χ), and P (χ) are given by where we have used the new coordinate χ 2 = 1 − u to avoid square roots. Again, we get a quadratic eigenvalue problem, and the function φ 3/2 (χ) is regular in the interval χ ∈ [0, 1]. In Sec. V E we solve Eq. (30) by using the pseudo-spectral and AIM methods and compare our results against the results of Refs. [41,48].

D. 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 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 F.

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 pseudospectral method, which recently has being applied to solve differential equations numerically in many problems. The fact that the coordinate domain is not periodic, u ∈ [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]).
The idea behind the pseudo-spectral method is to rewrite the regular function φ(u) in a base composed by cardinal functions C j (u), in the form 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 Note that (36) Evaluating on the grid, the polynomials of (34) become elements of a matrix c j (u i , λ, λ 2 ) = c j,0 (u i ) + λ c j,1 (u i ) + λ 2 c j,2 (u i ). Then, the matrix representation of the quadratic eigenvalue problem (34) can be written as where here D ji , D ji , and D (2) ji represent the cardinal function and its derivatives. Definingg = λg, the last equation may be written in the form 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 where we have defined the new matrices Notice that M 0 and M 1 are (N + 1) × (N + 1) matrices and 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 We call this particular choice as pseudo-spectral I. The second model considers two Chebyshev polynomials [45] 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 φ 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 ∈ [−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 O 1/N N 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 where primes denote derivatives with respect to to the variable x (that is defined over some interval that is not necessarily bounded), λ 0 (x) = 0 and s 0 (x) are C ∞ . 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 λ 0 and s 0 be functions of the variable x ∈ (a, b) that are C ∞ on the same interval. The differential equation (44) has a general solution of the form if for some n > 0 the condition α ≡ s n λ n = s n−1 λ n−1 (46) or equivalently 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 ω 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 λ 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 λ and s around a point ξ where the AIM is to be perform (we remind the reader that the results are independent of the choice of ξ), that is, where c i n and d i n are the Taylor coefficients of the expansions of λ n and s n around ξ, respectively. By plugging Eqs. (50) and (51) into Eqs. (48) and Eq. (49) one gets Finally, using Eqs. (52) and (53) the quantization condition, Eq. (47), becomes 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., 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 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. 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. 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. 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. Initialize ccda with data from icda and cdda with data from idda.

5.
Perform n AIM steps using the evolution Eqs. (55) and (56). That is, repeat the following n times: (a) Copy the content from ccda into pcda (b) Copy the content from cdda into pdda (c) Rewrite each element of ccda and cdda using Eqs. (55) and (56), respectively.
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 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.

A. Spin 0 QN frequencies
Our numerical results for the quasinormal frequencies of the spin 0 perturbation field are displayed in Table I. The first two column shows the data from the pseudospectral 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.

B. Spin 1 QN frequencies
Our numerical results for the quasinormal frequencies for spin 1 fields are displayed in Table II 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.

C. Spin 2 QN frequencies
Our numerical results for the quasinormal frequencies for spin 2 fields are displayed in Table III 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.

D. Spin 1/2 QN frequencies
Our numerical results for the quasinormal frequencies for spin 1/2 fields are displayed in Table IV 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 = 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 . Our numerical results are displayed in Table V, 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 ω = −in/4, n → ∞, obtained in Refs. [56,57], see also references therein.

E. Spin 3/2 QN frequencies
Our numerical results for the quasinormal frequencies for spin 3/2 fields are displayed in Table VI 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 V for the three routines we are working with. Such frequencies arise when investigating the quasinormal modes in the limit of large . 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, ∼ 5/4, · · · . l n Pseudo-spectral I (60 Polynomials)

F. Spin 5/2 QN frequencies
Our numerical results for the quasinormal frequencies for spin 5/2 fields are displayed in Table VII. 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 VIII. It turns out that these results satisfy the sequence 1/8, ∼ 3/8, ∼ 5/8, ∼ 7/8, · · · , in general, As in spin 1/2 and 3/2 perturbations, these frequencies arise when we were investigating the frequencies in the limit of large . 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 ∼ e iω t = e i(ω Re −i ω Im ) t = e ω Im t cos (ω Re t).

VI. QNMS FOR
1 AND n 1 It is also interesting to calculate the quasinormal frequencies in the limit of large , 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 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 and n, respectively. As it can be seen, the real part as a function of 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.
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 ω AIM Re − M ω PS Re | and |M ω AIM Im − M ω PS 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.
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. Ref. [41] Ref. [48] 0

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 spacetime. 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 I  -III. 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, † † † † † † † † † † † † † † † † † † † † † † † † † † † † † † * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * † AIM * Pseudo-Spectral  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 V. We also calculated the quasinormal frequencies for a spin 5/2 field, which are displayed in Table VII. 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 VIII. Hence, the pseudospectral 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 anal-ysis 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 = 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.
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 = 0 in Eq. (17), we do not get any solution. Then, we set = 1 and we get the solutions M ω 0 = ±0.110455 − 0.104896i and M ω 1 = ±0.086158 − 0.348079i, corresponding to n = 0 and n = 1, respectively. It is also interesting to point out that the solution for = 1 and n = 0 arises in both methods employed in this work, while the solution for = 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 = 1 represent an addition of angular momentum to the metric given by Eq. (1). Hence, such modes do not generate gravitational waves.