Overdamped modes in Schwarzschild-de Sitter and a Mathematica package for the numerical computation of quasinormal modes

We present a package for Mathematica that facilitates the numerical computation of the quasinormal mode (QNM) spectrum of a black hole/black brane. Requiring as input only the QNM equation(s), the application of a single Mathematica function will compute the spectrum efficiently, by discretizing the equation(s) and solving the resulting generalized eigenvalue equation. It is applicable to a large variety of black holes, independent of their asymptotics. The package comes fully documented and with several tutorials. Here we present a self-contained review of the method and consider several applications. We illustrate the method in the simplest case of scalar QNMs of a Schwarzschild black brane in anti-de Sitter. Then we go on to look at the scalar QNMs of the Schwarzschild black hole in de Sitter, in anti-de Sitter and in asymptotically flat spacetimes, finding a novel infinite set of purely imaginary modes in the first case. We also derive the QNM equations for a generic Einstein-Maxwell-scalar background and use these to compute the QNMs of the asymptotically anti-de Sitter Reissner-Nordstr\"{o}m black brane, as a further illustration and check of the method.


Introduction
Small perturbations of a black hole generically take the form of damped oscillations called quasinormal modes (QNMs). These perturbations occur only at a discrete set of frequencies, depending on the black hole itself. The spectrum of QNM frequencies contains a wealth of information about the black hole.
They describe the late time evolution of any dynamical process that results in a black hole at equilibrium. In particular in the black hole mergers recently observed by LIGO [2][3][4] the last part of the gravitational wave signal is described by quasinormal modes. With more data of black hole mergers expected in the future, matching the signal to the quasinormal modes will be a good test of General Relativity [5]. In particular as noted in [6] a more accurate determination of the mass and angular momentum of the black hole is necessary to rule out or constrain many alternative theories of gravity.
More precisely, quasinormal modes are linear perturbations δφ which behave in time as for some complex frequency ω, its (negative) imaginary part giving exponential damping and its real part giving an oscillation. By causality we require the perturbation to be ingoing at the black hole horizon, and either outgoing (asymptotically flat or de Sitter spacetimes) or finite (anti-de Sitter) at spatial infinity. This results in a discrete spectrum of frequencies, ω n . When there is no black hole there is nothing for the perturbation to dissipate energy into, there will be only normal modes, undamped oscillations, as we will see in 4.1.
A perturbation can also become unstable, when Im(ω) > 0, showing exponential increase instead of decay. This can be used to look for new solutions as the end point of such instabilities, see e.g. [24]. There are also cases known where the instability has no end point and the black hole continues to grow until it hits the anti-de Sitter boundary [32][33][34].
Through the fluid-gravity correspondence [13] the quasinormal modes of black branes are related to hydrodynamics. In particular, the quasinormal modes whose frequencies vanish as momentum is taken to zero encode in principle all hydrodynamic transport coefficients in their dispersion relation [9]. Higher modes decay more quickly and contain more information than hydrodynamics. Recently [15] found strong evidence indicating that these modes set a limit on the range of applicability of hydrodynamics.
Holography relates the physics of black holes to the physics of strongly coupled quantum field theories [14]. To each field on the gravity side corresponds a gauge-invariant operator on the QFT side. The quasinormal frequencies of this field are equal to the poles in the retarded Green's function of the operator in the dual QFT [9,12]. By summing quasinormal modes one can compute one loop determinants [20]. Out of equilibrium gauge theory plasmas are dual to out of equilibrium asymptotically anti-de Sitter black holes. At late times these are again described by quasinormal modes, see [16,17] for an explicit comparison.
Being so full of information, it would be desirable to be able to compute these for any black hole. While there are some cases where it can be done analytically (a recent quite involved example is [21], and see [22] for a review of analytic computations and approximations), of course the generic case can only be done numerically.
In this paper we present a Mathematica package which numerically computes quasinormal modes [1]. It is applicable to a broad class of cases: it does not matter what the asymptotics are, if the equations are coupled, if the frequency occurs to a higher power in the equation, or if the background is numerical (assuming this numerical background has been computed).
The method used was first used for general relativity in [27][28][29]. Essentially it discretizes the quasinormal mode equation(s) using spectral methods, and then directly solves the resulting generalized eigenvalue equation. Some great reviews on numerics in gravity are [23][24][25].
We have attempted to make the package easy to use, efficient in its computation, fully documented and with code that is itself easy to read and, if necessary, debug. Earlier versions were used successfully for the quasinormal mode computations in [32,35,36,56]. We encourage anyone to use it, and possibly to contribute by making improvements or adding new features.
The paper is structured as follows. In section 2 we describe the method used in the package. Then in section 3 we describe how to use it at a more practical level, illustrated with one of the simplest cases: the Schwarzschild-anti-de Sitter black brane, in 5 dimensions. The next two sections contain more examples. In section 4 we study the 4-dimensional Schwarzschild black hole for various asymptotics: again anti-de Sitter, flat space and de Sitter. In this last case we find an infinite set of purely imaginary modes that has been overlooked in the literature. In section 5 we derive the quasinormal mode equations for a generic Einstein-Maxwell scalar action (with a homogeneous and isotropic background). We use this to compute a more involved example: the 5-dimensional asymptotically anti-de Sitter Reissner-Nordström black brane. This serves to illustrate the method in the more complicated case of coupled equations while also giving a check on the numerics, since we also have some analytic control on these quasinormal modes. In the appendices we give a short report on the performance of the package and tables with numerical values for some quasinormal modes.

Method
To find the quasinormal mode spectrum, we have to solve the linear ordinary differential equation (ODE) (or coupled system of ODE's) describing a linearized perturbation on top of a black hole/black brane. At the horizon, causality requires us to choose only ingoing waves. Similarly at spatial infinity, or the cosmological horizon in the case of de Sitter, we must require only outgoing waves. In asymptotically anti-de Sitter spacetime we require the perturbation to be normalizable at the boundary.
These two boundary conditions can only be satisfied for a discrete set of frequencies ω ∈ {ω n |n = 1, 2, 3, ...} 1 , so at the same time as solving the ODE we have to solve for these ω n .
In this section we will discuss the method that is used in the package to do this.

Analytics
We will illustrate the method by the example of a massless scalar field in a 5-dimensional asymptotically anti-de Sitter Schwarzschild black brane. It turns out that Eddington-Finkelstein coordinates, are perfectly suited for this problem. For the time-independent backgrounds we will consider, these coordinates take the form Although not strictly necessary, they usually simplify the problem, we will see why below, and we will use these coordinates throughout this paper. In these coordinates the asymptotically AdS Schwarzschild black brane is where u = 0 corresponds to the boundary and u = 1 corresponds to the horizon. As long as the background is homogeneous, which in our case it is, we can write the fluctuation as a plane wave, where ω is the frequency and k is the momentum. The same relation would hold for any other fluctuation in this background. Inserting this into the Klein-Gordon equation we obtain, We have rescaled the quasinormal mode frequency ω and the momentum k to the dimensionless ratios λ ≡ ω/2πT and q ≡ k/2πT , where T = 1/(4π) is the temperature of the black brane.
Here we can already see one advantage of using Eddington-Finkelstein coordinates: the ansatz (2.1) implies that g tt = 0, and as a consequence the quasinormal mode equation 2.4 is linear in the frequency, whereas it would be quadratic in Fefferman-Graham coordinates. The advantage of this will become clear in subsection 2.3. Now we must analyze the behavior near the horizon and near the boundary to see how to deal with the boundary conditions. Starting with the horizon, by plugging in an ansatz φ(u) = (1 − u) p , we find that there are two solutions: δφ in (u) = const + O(1 − u) and Including time dependence, this last solution behaves as so as t increases, 1 − u has to increase as well to keep a constant phase, meaning u has to decrease, which means that this solution is outgoing. Hence we must make sure that we only get the other, ingoing solution. Notice that the ingoing solution is perfectly smooth near the horizon, while the outgoing one oscillates more and more rapidly as we approach the horizon. These properties will help us pick the correct solution, as we will see later.
Near the boundary there are two solutions, a non-normalizable mode δφ(u) ∝ 1 which we must discard, and a normalizable one δφ(u) ∝ u 4 . If we redefine δφ(u) ≡ u 3 δφ(u), then the normalizable solution hasφ approaching zero linearly, while the non-normalizable solution diverges. 2 Doing this rescaling, and also rescaling the equation itself so that it is finite but nontrivial at the boundary, the equation becomes: (2.6) Now the normalizable solution behaves perfectly smoothly both at the boundary and at the horizon, while the non-normalizable solution behaves pathologically, diverging and rapidly oscillating respectively. We will see in the next section how this can be used to automatically deal with the boundary conditions.

Discretization: Pseudospectral Methods
Having derived the equation (2.6) in a form without divergences, we will now discretize it in order to solve it numerically. For this we use pseudospectral methods, the standard reference on this is [26]. The pseudospectral method solves a (differential) equation by replacing a continuous variable, the radial one in our application, by a discrete set of points, also called collocation points. The collection of these points is usually called the grid.
A function can then be represented as the values the function takes when evaluated on the gridpoints. An equivalent and useful way of looking at this set of numbers representing a particular function is as coefficients of the so-called cardinal functions. The cardinal functions corresponding to the grid {x i |i = 0, ..., N } are polynomials C j (x) of degree N , with j = 0, ..., N , satisfying C j (x i ) = δ ij . The choice of a grid uniquely specifies the cardinal functions as A function f is then approximated as The expansion in terms of cardinal functions allows one to construct the matrix D (1) ij that represents the first derivative, D ij = C i (x j ), and similarly for higher derivatives. Solving the resulting linear equation will give a function which solves the original equation exactly at the collocation points. The hope is that as the number of gridpoints is increased, it will also solve the equation at other points.
For this to work, the choice of collocation points is crucial. The choice which tends to work best and is therefore most often used is the Chebyshev grid: (2.9) For this grid it can be proven that any analytic function can be approximated with exponential convergence in N . This can be understood by realizing that when we increase N , not only does the largest distance between gridpoints decrease proportional to 1/N , but at the same time the order of the cardinal functions increases, leading to a numerical error scaling as (1/N ) N . Note that while this means that we can approximate the equation with exponential convergence, this does not mean that the number of quasinormal modes we find will grow exponentially with N , for more details see Appendix A.
As given, these points lie in the interval [-1,1], but this can be rescaled and shifted to any other interval. In particular, we will want to rescale it to [0,1] if the horizon is at 1.
For this grid, the cardinal functions are linear combinations of Chebyshev polynomials T n (x): In Fig. (1) we show what these functions look like for N = 6. Of course these functions are all perfectly smooth, and at the endpoints they are either 0 or 1 (since the endpoints themselves are collocation points). With a linear combination of these smooth and finite functions we can never approximate either a function that is diverging or that is rapidly oscillating. As a consequence, we have already implicitly solved the boundary conditions by choosing these functions as a basis.

Generalized Eigenvalue Equation
Now we have turned the problem of solving a linear ODE, subject to specific boundary conditions, into solving a matrix equation, where the boundary conditions are already implicitly solved. It is still not purely numerical though, as the equation depends on the frequency, which we have to solve for as well.
This can be done by recognizing that this is a generalized eigenvalue problem. The simplest type of quasinormal mode equation is of the form where each of the c i are linear in ω: To be completely explicit, in our example Eq. (2.6), setting the momentum q to zero, this will give c 0,0 = −3 − 9u 4 , c 0,1 = 6iu, c 1,0 = u(3 − 7u 4 ), c 1,1 = 4iu 2 , c 2,0 = u(u − u 5 ) and c 2,1 = 0. Each of these coefficients c i,j (u) is turned into a vector by evaluating it on the gridpoints. These vectors are multiplied with the corresponding derivative matrices D (n) ij and the resulting matrices added, to bring the equation into the form, where the M i are now purely numerical matrices. Explicitly, ij , and similarly for M 1 . Equation (2.12) is precisely a generalized eigenvalue equation. This can be solved directly using Mathematica's built-in function Eigenvalues (or Eigensystem to get the eigenfunctions as well).

Extensions
The method presented works not just for the simple case used to illustrate it, in fact it works quite generally with a few simple modifications.
Firstly, suppose that instead of one ODE, we have a coupled system of ODE's. Say there are two of them (but the following applies generally), of the form We can discretize this in a similar way, by joining the two functions φ and ψ into a single vector (φ(u), ψ(u)). The matrix M 0 of Eq. (2.12) becomes where the 0 everywhere indicates that this is the piece of order 0 in the frequency, the second index indicates the equation number and the argument indicates the function, so thatM 0,1 (φ) is the matrix coefficient of φ in the first equation. Of course there is a similar equation for the linear term in the frequency.
Further, we can generalize to an equation depending on the frequency as an arbitrary polynomial. Whether coupled or not, the procedure above will bring such a (system of) equation(s) to the form Note that in the case of a coupled equation, the φ here would be a vector composed of several functions, as above.
Notice that for a coupled system of N eq equations with the maximal power of the frequency being p, the resulting matrix will be N eq p(N + 1) (where N is the N in Eq. (2.9)). So both extensions come at a price in computational time.
An alternative here is to solve det(M(ω)) = 0. This is prohibitively slow to do symbolically, but it can be done numerically by sweeping the complex ω-plane. This has the advantage that higher powers of ω don't require larger matrices, but the extra complexity of selecting and possibly refining a grid of points in the complex plane. 3

Using the package
In this section, we will show at a very practical level how to work with the package. Our starting point will be the properly rescaled equation for a massless scalar in an asymptotically anti-de Sitter Schwarzschild black brane, at zero momentum, Eq. (2.6).
First some practical remarks. The package can be found here [1], where installation instructions can also be found. Once installed there are several ways to get started and familiarize oneself with it. Apart from the present paper, the package also has its own documentation. This can be accessed through Mathematica by going to Help, Documentation and typing in "QNMspectral", and describes all the functions and their options. It includes some tutorials, which also go into more detail on how to obtain Eq. (2.6). We also provide a separate notebook with several examples, also found at [1].
Here we continue with Eq. (2.6). The function which implements everything mentioned in the previous section is called GetModes. Assuming Eq. (2.6) is stored in Mathematica under eq, one can compute the quasinormal modes simply as This does the computation with N = 40 (meaning with 41 gridpoints), at machine precision (machine precision has about 16 digits, anything less than that in the last argument will use machine precision). Note that it is not necessary to specify what are the frequency, the radial variable and the fluctuation, this can be determined unambiguously from the equation and is done automatically. Now the result, a list of quasinormal mode frequencies, is stored in modes and can be displayed by evaluating PlotFrequencies[modes,Name→''ω/2πT ''], producing Fig. (2).
The computation used a grid of 41 points, so there are 41 eigenvalues. Typically the quasinormal modes lie in an approximately straight line, they certainly should in this case. This means that most of the eigenvalues found are numerical artifacts, only a few are accurate. This is unavoidable and will continue to be the case if we increase N or the precision. We will get more accurate results, but also more results in general, so still only a small fraction of the total will be accurate.
To test whether a computed eigenvalue really is a quasinormal mode and not just a numerical artifact, one has to repeat the computation at different grid sizes and precisions, and look for convergence. A simple and efficient implementation of this is given in the function GetAccurateModes. Used for example as this does the computation twice, once with N = 40 at machine precision and once with N = 80 at precision 40. It returns only those modes which occur in both computations, and of the modes it returns, it takes only those digits which agree. Note that for this to actually give only converged modes, one has to be sure that the equation itself is either fully analytic, or if it is numeric that the error in the equation is not the dominant one. One also has to be sure that there is no common error in both computations, as this will not get filtered out. This can happen when taking the grid sizes too close.
We can display these results in a plot of the complex plane and a table by evaluating ShowModes[modesaccurate,Name→''ω/2πT '',Precision→Infinity], giving Fig. (3). The option Precision→Infinity is used here to show all computed digits, by default only the first 6 are shown. 4 These values can be compared with [11], where we see that indeed all displayed digits are accurate.
We see that the lowest mode has been computed quite accurately, and as the mode number increases the precision of our result decreases. Apart from the caveats mentioned above, in practice this often gives only correct results. We advise however to do more elaborate checks to ensure accuracy.
One additional check that one can do is to look at the eigenfunctions as well, and check that they are smooth and finite. By default, these are not computed. To compute them, simply add the option Eigenfunctions→True to either GetModes or GetAccurateModes. Repeating the computation (3.2) with this option, we can now plot the eigenfunctions by PlotEigenfunctions[modesaccurate], producing Fig. (4).
In this case, one can see that all eigenfunctions are perfectly smooth, and go to 0 at the boundary, as was expected from the rescaling we did. Notice that they are normalized . Eigenfunctions come in conjugate pairs, only those with positive imaginary part are shown. All are smooth, normalized to go to 1 at the horizon, and go to 0 at the boundary as was expected from the rescaling that we did.
to be 1 at the horizon.

Schwarzschild black hole
As another simple example and to illustrate the general applicability we now consider the fluctuations of a massless scalar on top of a Schwarzschild black hole in various 4dimensional spacetimes: anti-de Sitter, flat space and de Sitter. These black holes are all maximally symmetric solutions of the equations of motion coming from the action The background metric we take is where we set Newton's constant G = 1, M is the mass of the black hole, and depending on the asymptotics we have = 1, 0, −1 and Λ = −3/L 2 , 0, 3/L 2 , for anti-de Sitter, flat space and de Sitter respectively.
Apart from the different asymptotics to the previous example we take a different horizon topology: a sphere instead of a plane. This has the consequence that plane waves are no longer solutions, and we must instead use spherical harmonics, labelled by a multipole number l. The perturbation we can do is then: These perturbations have to satisfy the Klein-Gordon equation, which takes the form: In each of the cases considered, the behavior of the fluctuation near the horizon is the same as in the black brane discussed before. There is an outgoing mode δφ out (u) = (1−u) iλ , and an ingoing mode δφ in (u) = const., due to our use of Eddington-Finkelstein coordinates. This is perfect, as the numerical approximation is only able to resolve the ingoing mode. We will explore each of the cases in turn, starting with the one most similar to the previous example: anti-de Sitter.

Anti-de Sitter
Just as in the previous example, anti-de Sitter has a conformal boundary at u = 0. The scalar field contains a normalizable mode as well as a non-normalizable mode, behaving near the boundary respectively as φ ∝ u 3 and φ ∝ const. We do not want our fluctuation to change the boundary, therefore we demand the non-normalizable mode to vanish. Similar to before, we do this by redefining so that the normalizable mode approaches 0 as u → 0 and the non-normalizable mode diverges.
The blackening factor, in Eq. (4.2) with = 1, has a horizon at some u = u h , in terms of which the mass is M = (1 + L 2 u 2 h )/(2L 2 u 3 h ) and the black hole temperature is T = (3 + L 2 u 2 h )/(4πL 2 u h ). We will set u h = 1 without loss of generality. Replacing M by L using the equation above, replacing ω by λ ≡ ω/(2πT ) and doing the rescaling of δφ we obtain the quasinormal mode equation: Note that we now have the dimensionless ratio M/L as a parameter describing the relative size of the black hole and the whole anti-de Sitter universe.
In Fig. (5) we show the quasinormal modes for l = 0, 1, 2 as M/L is varied between 0 and infinity. Some of these values were computed before in [10], with which we of course agree 5 .
In the limit M/L → ∞ the Schwarzschild black brane is approached. The l-dependence drops out and all modes converge at those of the brane, indicated by the blue dots. The opposite limit M/L → 0 is that of empty anti-de Sitter. In this limit we obtain the normal modes of empty anti-de Sitter, ω n L = ±(l + 1 + 2n), n = 1, 2, ... [10,44,45]. These are pure oscillations without decay, since there is no horizon to decay into. 5 The case where the scalar also has a conformal coupling was studied in [60,61].

Flat Space
Now we discuss the case of flat asymptotics. The background is obtained by setting = 0 in Eq. (4.2), which when inserted into Eq. (4.4) gives the quasinormal mode equation for the asymtotically flat Schwarzschild black hole: The behavior of φ is singular at spatial infinity, u = 0, but we can scale this out. There are again two solutions, φ out (u) = exp (2iω/u) u 1−2iω (1 + ...) representing waves going out to infinity, and φ in (u) = u(1 + ...) representing waves coming in from infinity. We must demand the latter to vanish, so we redefine so thatδφ(u) ∝ u for the outgoing modes, while the modes coming in from spatial infinity diverge. We set the horizon u h = 1/(2M ) = 1 and define λ = ωM , to get the final equation, In Fig. 6 we show the scalar quasinormal modes for l = 0, 1, 2, 3, along with a table for the first four modes for each l. These points reproduce table 1 of [40] exactly.

De Sitter
The Schwarzschild-de Sitter black hole is the most physically rich case of all. There are two horizons: the cosmological de Sitter horizon and the black hole horizon, whose coordinates we will denote as u c and u b respectively. The blackening factor, Eq. (4.2) with = −1, can be conveniently rewritten in terms of these quantities as (4.9) Figure 6.
Scalar quasinormal modes of the Schwarzschild black hole in asymptotically flat spacetime, for multipole number l = 0 (blue circles), l = 1 (red squares), l = 2 (green diamonds) and l = 3 (orange triangles). On the right is a table with the lowest 4 modes for each l.
Note that this is symmetric in u c , u b , but for our naming to make sense we require that u c < u < u b , so in particular we have to restrict 0 < u c < u b .
These horizons have surface gravities respectively. The quasinormal mode equation in these parameters becomes, The background solution has two independent dimensionful parameters, from which we can extract a dimensionless ratio which we will take as M/L, the mass of the black hole divided by the radius of de Sitter. In terms of u c and u b this is Solving Eq. (4.10) near the cosmological horizon we again find two solutions, δφ(u) = (u − u c ) −iλc , with λ c ≡ ω/κ c , or δφ(u) → const. Restoring time dependence on the first solution we get δφ(t, u) = exp −iω(t + 1 2πTc log(u − u c )) . To keep a constant phase as t increases, we need to decrease u, which means that this solution is going into the cosmological horizon. In order to keep only this solution we need it to be smooth while the outgoing one oscillates rapidly, so we have to redefine so thatφ(u) ∝ (u − u c ) for the ingoing solution, while the outgoing solution oscillates rapidly.
Making this replacement, rescaling the radial coordinate u c < u < u b to 0 < x < 1, and setting u b = 1, we obtain the final quasinormal equation, shown in appendix B.
Before going into the numerical results it is instructive to look at the limiting cases. One limit is the extremal limit M/L → 1/(3 √ 3). In this limit u c → u b = 1, but the proper distance remains finite and one obtains the Nariai spacetime. Here the quasinormal modes can be calculated analytically [37] to be: extremal: ω n /κ c = −(n − 1/2)i + l(l + 1) − 1/4 , n = 1, 2, ... , (4.13) The other limit, M/L → 0, has two qualitatively different interpretations: it can be reached by taking M → 0 at fixed L, corresponding to the limit of empty de Sitter, or L → ∞ at fixed M , corresponding to the limit of the Schwarzschild black hole in asymptotically flat spacetime. The former has analytic quasinormal modes [41]: de Sitter: ω 1 /κ c = −li; ω n = −(l + n)i , n = 2, 3, ... . (4.14) For l = 0 this gives a zero mode, which is present for any M/L, reflecting the fact that φ can be shifted by a constant.
There is no analytic result for the quasinormal modes of Schwarzschild black hole in asymptotically flat spacetime, but of course they have been computed numerically in e.g. [40] and in the previous section.
For small M/L, we expect the equilibration process of the full space time not to be influenced by the presence of a small black hole, and conversely we do not expect the equilibration of the black hole to be significantly different from the same black hole in asymptotically flat spacetime. So we expect to see two decoupled sets of quasinormal modes which are small deformations of those of empty de Sitter and of Schwarzschild in asymptotically flat spacetime.
In Fig. (7) we show the imaginary part of the quasinormal modes, for l = 0 and for the full range of M/L. In both cases, blue lines correspond to complex frequencies, while red lines correspond to purely imaginary frequencies. On the left plot, the modes are normalized with respect to κ c , and at M/L ≈ 0 we indeed see the pure de sitter modes of Eq. (4.14). As M/L grows they move down the complex plane, but they stay on the imaginary axis for the whole range of M/L. The extremal values of Eq. (4.13) are approached by the other set of modes, which are complex in general but become purely imaginary in the extremal limit M/L → 1/(3 √ 3). In Fig. (7b) the same data is presented, but now normalized as ωM . In these units, the complex modes approach those of the Schwarzschild black hole in flat space, as can be seen by comparing with Table (6b), and as indicated by the blue dashed lines. All of the imaginary modes now disappear in the limit, simply because M goes to zero. In this limit it becomes numerically difficult to compute the complex modes, as there are so many smaller purely imaginary modes.
Note that in this asymptotically flat limit, the boundary condition at the cosmo-logical horizon approaches the boundary condition of outgoing waves in asymptotically flat spacetime continuously. This is most easily seen in different coordinates, ds 2 = −f (r)dt 2 + f (r) −1 dr 2 + r 2 dθ 2 + sin 2 (θ)dφ 2 , where f is the same as in Eq. (4.9), with r = 1/u. Defining the "tortoise" coordinate r by dr /dr = 1/f (r) and φ(t, r, θ, φ) = e −iωt r −1 Y lm (θ, φ)ψ(r) we obtain the Schrödinger form of the quasinormal mode equation, The Schrödinger potential V (r) approaches zero both at the cosmological horizon of Schwarzschild de Sitter and as r → ∞ in the asymptotically flat case, so that expanding ψ near these points we find in both cases ψ(r) = exp (±iωr ), where we have to choose the plus sign to obtain the outgoing wave.
In the anti-de Sitter case, the Schrödinger potential no longer approaches zero at the boundary, and so in the flat limit the boundary condition does not approach that of the asymptotically flat Schwarzschild black hole. This is why in this case, discussed in section 4.1, there isn't a set of modes converging to the asymptotically flat Schwarzschild values.
So the limits match, but what happens in between is also interesting. Firstly, for small M/L the purely imaginary mode coming from de Sitter is dominant, so the late time behavior of a small perturbation (that is sufficiently spread out to excite both sets of modes) is purely damped at late times. At M/L ≈ 0.051 there is a crossover, where the complex mode becomes dominant and late time behavior is a damped oscillation. This crossover occurs at M/L ≈ 0.090, 0.047, 0.032 for l = 1, 2, 3, so for l = 1 the imaginary mode dominates in the largest range, which matches the fact that the dominant mode in pure de Sitter is that of l = 1 (excluding the nondynamical zero mode at l = 0).
For larger M/L we no longer expect two decoupled sets of modes, and in fact we observe an interaction between the two sets, which produces the oscillations seen in Fig.  (7). This can be seen more clearly in Fig. (8a), where we show the spectrum in the complex plane, plotted for all M/L at the same time, with arrows indicating increasing M/L. What happens is that as a purely imaginary mode moves down the axis and a complex mode reaches a similar imaginary part, the complex mode moves towards the axis, as if attracted by the imaginary mode, and moves away again once the difference in imaginary parts increases again.
The "attraction" becomes stronger for higher complex modes, and for each it is strongest when it approaches the lowest imaginary mode. This makes intuitive sense, as the larger the mode number, the higher the value of M/L when it crosses the lowest imaginary mode.
The consequence of this interaction is that black holes of certain masses, where the imaginary parts of a complex mode and a purely imaginary mode are equal, will oscillate slightly less. It would be hard to detect this in a real evolution though, as the effect is quite small for the lower modes, and it occurs at a different mass for each mode.
Starting from the 9th lowest complex mode the complex modes actually reach the imaginary axis and hit the imaginary mode before moving away from it again, in the region indicated by the box in Fig. (8a). What happens is more clearly seen in Fig. (8b), which shows only the 9th mode in the region of this box, plotting the real and imaginary part separately as a function of M/L. Since φ is real, complex modes must occur in conjugate pairs, but once they become purely imaginary they can split. The conjugate pair moves towards the real axis and splits into two purely imaginary modes when they hit the axis. One of these imaginary modes moves up the axis towards the original purely imaginary mode, which it then pairs up with to become complex again. The other continues down the axis and stays imaginary.
For higher l these interactions are weaker, with l = 1 still showing some oscillations but from l = 2 on they appear to be absent. Apart from that the picture at larger l is qualitatively similar to what has been discussed. In Table (1) we present some numerical values. These match the values of [38], which used a sixth order WKB approximation, except that we find the additional infinite set of purely imaginary modes, which have not been reported before in the literature. The WKB approximation is known to miss overdamped modes [7], which in some cases misses the mode that is physically the most relevant.
These modes do fit nicely with [46], which did an explicit time evolution of a scalar perturbation of a black hole in de Sitter with Λ 1 (i.e. very small M/L). They observe an initial exponential decay with approximately the asymptotically flat Schwarzschild black hole QNM, and a late time exponential decay with approximately the empty de Sitter QNM. Here we showed how these QNMs are modified for larger M/L, and that this picture only holds up to some crossover value where the complex mode is dominant.
Although we did not compute electromagnetic or metric perturbations it seems clear that there also one expects purely imaginary modes as continuations of pure de Sitter modes, just as for the scalar. One might then also expect similar oscillations in the complex modes at the lowest l = s. It is therefore interesting to look at [42] in this light, which studies the high overtones of these electromagnetic and metric perturbations. It is shown that at fixed ΛM 2 = 0.02 and for large overtone n, Re(ω n ) oscillates as a function of Im(ω n ). This might be similar to the oscillation seen in Fig. (8a), since taking a snapshot at constant ΛM 2 would still produce oscilaltory behavior as a function of n, although it is difficult to give a more precise relation, since we consider only scalar modes, at much lower n.

Einstein-Maxwell-scalar model
To illustrate the method in a more complex example, we will look at the quasinormal modes of the 4 + 1-dimensional Reissner-Nordström anti-de Sitter black brane. Before we specialize to this case we will first derive the equations in a much more general setting, namely any Einstein-Maxwell-scalar black brane that is homogeneous and isotropic (and ofcourse time-independent), in any number of dimensions.
We consider the following action, where F = dA is the field strength of a U (1) gauge field A, which has its kinetic term modified by a function Z(φ) of a scalar field φ, we have a cosmological constant Λ and a potential V (φ) for the scalar field, and finally ξ is an arbitrary normalization factor for the scalar field kinetic term. We work in units where κ = 1/(8πG) = 1/2. The equations of motion derived from this action are as follows, We specialize to a homogeneous, isotropic and time-independent planar black hole background, allowing us to use the following ansatz, In these coordinates the temperature of any black hole with a horizon at u h is given by In order to bring the quasinormal mode equations into a useful form, we follow the approach of [9], which we generalize from pure gravity to this rather general Einstein-Maxwell-scalar case.
We start by fluctuating all of the fields as plane waves, with momentum q in the x direction,

(5.4)
In what follows we will denote spatial indices other than x, so transverse to the momentum, with α = β.
We can classify all of these fluctuations according to their helicity, or their transformation properties under the remaining SO(2) symmetry of rotations around the x-asix. 6 This immediately groups the fluctuation equations into three mutually decoupled groups, also called channels, with helicities 0, ±1 and ±2 and containing the fields, (helicity ± 2): h α,β , h α,α − h β,β , (helicity ± 1): a α , h t,α , h r,α , h x,α , (helicity 0): δφ, a t , a r , a x , h t,t , h t,r , h t,x , h r,r , h r,x , h x,x , h , (5.5) where for helicity 0 we defined h ≡ Σ d+1 α=4 h αα /(d − 2). Note that for d = 3 there is only one transverse spatial direction, so there are no helicity 2 fields. For any higher d the decomposition is as above.
Since the quasinormal modes that we want to calculate are gauge invariant quantities, we expect to be able to find gauge invariant equations for them. In order to do this, we have to express the fluctuations in gauge invariant combinations, under infinitesimal diffeomorphisms ξ µ (u, t, x) = e −i(ωt−qx) ξ µ (u), and infinitesimal U (1) gauge transformations λ(u, t, x) = e −i(ωt−qx) λ(u), which act on the fluctuations as where all covariant derivatives are taken only with respect to the background metric (other contributions would be second order).
By simply taking an arbitrary linear combination of all the fluctuations in a given channel, performing the above gauge transformation on this combination and demanding the variation to vanish we find the following gauge invariant variables, h .

(5.7)
In the helicity 2 case h α,α − h β,β is also gauge invariant and yields the same equation. The combinations of Eq. (5.7) are gauge invariant for any dimension d ≥ 4. For d = 3 the only difference is that there is no helicity 2 component as noted above, but the other combinations remain identical. These are the unique (up to taking linear combinations) gauge invariant combinations of the fluctuations, though if one allows for their radial derivatives as well there are more possibilities [47].
Substituting the gauge invariant fluctuations into the fluctuation equations, they can be completely decoupled from the gauge-dependent ones, again simply by taking an arbitrary linear combination of the equations in a given channel and demanding the coefficients of the gauge-dependent functions and their derivatives to vanish. We stress that it is not necessary to make any gauge choice in order to do this, although it simplifies the process.
In the general case the gauge invariant variables in the same channel do remain coupled however, resulting in a single decoupled equation for the helicity 2 channel, two coupled equations for helicity 1, and three coupled equations for helicity 0.
Very schematically, not writing any of the coefficients, the equations take the following form, for the helicity 2 channel, for the helicity 1 channel, and for the helicity 0 channel, (5.10) Each term has a coefficient that in general depends on u, the frequency ω and the momentum q, and any parameters of the background solution.
Note that even though the original coupled equations are only quadratic in ω, through the decoupling higher powers of ω arise. The equations go up to ω, ω 3 and ω 5 for helicity 2, 1 and 0 respectively, though the imposing of boundary conditions may raise this further, as we saw in sections 4.2 and 4.3.
At zero momentum, most equations decouple. The equations for E γ (γ = x, y, z) become identical to eachother, as do those for Z γ (γ = 1, 2, 3), however the fluctuation Z 2 in general still occurs in the equation for Z φ . Furthermore, if the background gauge field vanishes the equations for E γ (γ = x, y, z) decouple, since the action is quadratic in the gauge field, and if the background value of the scalar vanishes, the equation for Z φ decouples only if the potentials V (φ), Z(φ) are at least quadratic in φ.

Reissner-Nordström anti-de Sitter black brane
We now apply the previous to the asymptotically AdS 5 Reissner-Nordström black brane. We will encounter most of the numerical difficulties of the general Einstein-Maxwell scalar black brane, while still having an analytic background and even some analytic control on the quasinormal modes. This allows us to demonstrate the method in a more complicated case while also showing that it gives correct results.
The background, in terms of the ansatz of Eq. (5.3), is as follows, and ofcourse φ = 0. We fix the horizon to be at u = 1, which fixes M = 1 + Q 2 . The temperature of this brane is T = (1 − Q 2 /2)/π, so the extremal solution is Q extr = √ 2, Figure 9. Quasinormal mode spectrum of the Reissner-Nordström anti-de Sitter black brane, with charge Q/Q extr = 1/2 and momentum q/πT = 1.5.
which has a vanishing temperature. It has an entropy density s = 1/(4π) and a chemical potential µ = √ 3Q. This background has a single dimensionless ratio, which we choose to take asQ ≡ Q/Q extr .
We simply plug this background into the generally derived quasinormal mode equations. As before, the boundary condition of ingoing modes at the horizon is enforced automatically by using Eddington-Finkelstein coordinates. The only thing we still need to do is fix the boundary conditions at the boundary u = 0, where we must set the nonnormalizable mode to zero. For each gauge-independent fluctuation, the solutions near the boundary, and the rescaling we do to enforce this boundary condition, are where each fluctuation is now rescaled such that the non-normalizable mode diverges, while the normalizable mode approaches 0 linearly. Having made these redefinitions we are ready to compute the quasinormal modes. In Fig. (9) we show the full quasinormal mode spectum of all the channels, for the background charge Q/Q extr = 1/2 and momentum of the perturbations q/πT = 1.5. Modes of the different channels nearly overlap. This can be understood from the fact that, as we mentioned in the previous section, at zero momentum the equations for all the E γ become identical to eachother, giving the nearly overlapping spin 0 and spin 1 modes coming from E x and E y . The equations for all the Z γ also become identical, giving the nearly overlapping spin 0, 1 and 2 modes. It is somewhat remarkable though that the modes are still so close at the reasonably high q/πT = 1.5.
In Table (2) we give the numerical values of the lowest 10 quasinormal modes in each channel, for Q/Q extr = 1/2 and q/πT = 1.5 as in Fig. (9).
What depends much more strongly on the momentum are the three modes near ω = 0, shown enlarged in the inset. These are the hydrodynamic modes, which have the property that ω(q → 0) → 0. The spin 1 channel contains the shear mode, while the spin 0 channel contains the sound mode, and a mode governing charge diffusion.
In Fig. (10) we show the momentum-dependence of these hydrodynamic modes, again for Q/Q extr = 1/2. At low momentum these modes are given by the dispersion relations(see e.g. [18]), as: where and P are the energy and pressure, η and ζ the shear and bulk viscosity, v 2 s = ∂P/∂ is the speed of sound and D is the charge diffusion, normalized to be dimensionless and 1 at Q = 0.
The shear viscosity is universally fixed, for any two derivative gravitational action, to be η/s = 1/(4π) [19]. Furthermore the theory we are looking at is conformal, so that the bulk viscosity vanishes and v 2 s = 1/(d − 1) = 1/3. This means that the sound mode is also completely fixed. While the charge diffusion constant D is not fixed by η/s and conformality, it can be calculated analytically [48,49], so that the three hydrodynamic modes become, where we used + P = 4M = 4(1 + 2Q 2 ), following from a standard renormalization procedure [50]. In Fig. (10) Eqs. (5.14) are plotted as dashed lines going through the numerically computed points. In each case the dispersion relations describe the modes well at low momenta, although we can already see the higher order corrections around q/πT ≈ 1. We also check that the dispersion relations of Eq. (5.14) are followed for any other Q.

Discussion
The example of the Schwarzschild-de Sitter black hole illustrates one of the main advantages of the method we employed. As one of the simplest black holes, its quasinormal modes have been computed already in 1990 [39], yet all this time the infinite set of purely imaginary modes, which depending on the black hole mass may even be the dominant mode, have been missed. This is either because the approximation used misses this type of modes (e.g. the WKB approximation [51]), or in other cases. e.g. with Leaver's method [53], because the method requires an initial guess, effectively giving only what's expected. We checked that even the thirteenth order WKB approximation with Padé resummation of [52], while giving very accurate results for the complex modes, still completely misses the purely imaginary modes.
The method used here is completely a-priori, it has no guesses, no assumptions other than that the eigenfunctions be analytical, and no approximations other than the grid-size, which can easily be varied.
Two additional examples, rederiving the QNM's of [9] and [35], are included in [1]. Although we did not discuss a case where the background is known only numerically, it should be clear that as long as the numerical background is known to a high enough precision this will not give any problems, as an analytical background has to be converted to a numerical one in any case. In [32] we did use the package successfully for a numerical background, and the same method as we use here has been used with numerical backgrounds for instance in [48,54,55], and in analytic but more involved situations for instance in [65][66][67][68][69].
In an upcoming work [56] we will study the QNM's of the Reissner-Nordström black hole in asymptotically de Sitter spacetime using [1], in particular in relation to strong cosmic censorship.
There are several possible extensions to the current method. A limitation of the current approach is that the equation, after the rescaling to implicitly solve the boundary conditions, has to have a polynomial dependence on the frequency. This is not always the case, see for instance [57]. There are two main ways around this issue. The most closely related one is to take the determinant of the matrix representing the full QNM equation, as a function of the frequency, and look for zeroes by evaluating it on numerical points in the complex plane [58]. Another way, developed and used in [59], is to use the Newton-Raphson method. This can be advantageous even when not strictly necessary, as it allows one to efficiently track a single mode as a parameter is varied [62][63][64].
Another extension that is possible and would be interesting is to nonhomogeneous backgrounds. The main complication there is that the matrices will get a lot larger, by a factor of the size of the grid in the extra direction. All these extensions would fit nicely in the current framework and are worth adding in the future.
We hope that the present package will be of use to the community. We have certainly found it very convenient in our own work, and believe others could benefit from it. In addition, we encourage anyone to contribute with optimizations or bug fixes or extra features, it is hosted on GitHub, which facilitates such collaboration.
to keep in mind that for N eq coupled equations with the maximal power of the frequency being p the matrix size scales as N eq p.  Figure 11. Benchmark for performance and its dependence on the grid size N , for various fixed values of the precision p. In addition the ... line has precision p = N/2. Lines are power law fits of t(N ) = t 0 N q , finding q ≈ 2.377 using machine precision, q ≈ 3 using constant high precision and q ≈ ... using precision p = N/2. On the right a table with several absolute values, running on a single 1,7 GHz Intel i7 Core, with 8GB of RAM.
In Fig. (11) we show a log-log plot of the time t taken to compute the quasinormal modes (in seconds) as a function of the grid size N , using machine precision (blue points), a precision of 50 digits (yellow points) and a precision depending on the grid size as p = N/2 (green points). In all cases we find a power law t(N ) = t 0 N q . For the computation with machine precision, the power seems to be different for smaller grid sizes (q ≈ 2.38) and for larger grid sizes (q ≈ 3.2). For the case with 50 digits of precision, we find q ≈ 3.1, while for the one with precision increasing as N/2 we find q ≈ 3.3. Summarizing, the scaling with N is roughly as N 3 , except for grids with roughly less than 400 points where it scales as N 2.38 . The dependence on precision is very mild, as we can see by comparing the yellow and green points, except in the transition from machine precision to higher precision, which slows the computation down by roughly 2 orders of magnitude.
To the right of Fig. (11) we give a few explicit values of the time taken on a laptop. Note that for larger grid sizes and/or using higher precision, the computation of the eigenvalues takes up virtually all of the time, while only for low grid sizes the other steps in the computation, namely the construction of the matrices, become important. We check that for a grid of 40 points with machine precision, which we consider roughly the smallest grid that can still give valuable results, about 65% of the time is taken up by highly optimized built-in functions of Mathematica (which are programmed in C). These are the construction of the derivative matrices and the actual computation of the eigenvalues. This means that in the worst case, where the rest of the computation can be optimized so much as to take a negligible amount of time, the speed gain would be about 35%. However, since for example the construction of the numerical matrix, which involves computations with vectors (in particular the grid), does take some time, we believe the actual room for improvement to be much smaller.
In Fig. (12) we show the number of accurately computed modes as a function of the Note that before the different colors branch off the straight line, they follow it exactly, lying under the blue dots. This means that when it comes to the digits of precision used, for a given grid size N it needs to be above some threshold to compute the most accurate modes, but increasing it further has no effect. As long as the precision is high enough, the number of accurate modes increases linearly with the grid size, up until the point where the precision is no longer high enough, when the number of modes actually goes down a little with increasing grid size. In order to see what features we can expect to be generic, independent of the equation we are solving, and which are not, we repeat this procedure on the right-hand side of Fig.  (12) for the same scalar but now with a large momentum of q/πT = 160 (these modes were also computed in [31], we ofcourse agree on the result). This may seem like a trivial change, and on the level of the equation it is, but the solutions become highly oscillatory and thus harder to capture with a few Chebyshev polynomials.
From these two cases, we expect that the two features mentioned above are generic, at least after N is large enough to resolve the first mode.
The details are of course not generic, such as the actual value of n(N ) for a specific N , or the actual value of the precision p(N ) needed, or the coefficient of the linear growth of n(N ) (curiously, at high momentum we need less precision than at zero momentum).

B Quasinormal mode equations
In this appendix we write the quasinormal mode equations in their rescaled form which can directly be used in the package.
Quasinormal mode equations Here we present the scalar quasinormal mode equations of the Schwarzschild black hole for de Sitter and flat asymptotics, in a form that can directly be used for numerical computation.
For the Schwarzschild black hole, the anti-de Sitter case was given in Eq. (4.5) and the asymptotically flat case in Eq. (4.8). In de Sitter, using the radial variable v = (u − u c )/(u b − u c ) with the cosmological horizon at v = 0 and the black hole horizon at v = 1 , and having rescaled φ according to Eq. (4.12) to enforce the correct boundary conditions, the equation is with, , , The equations for the general Einstein-Maxwell-scalar backgrounds, and those for the specific case of the Reissner-Nordström black brane can be found in [1].

C Numerical values
In this appendix we provide some quantitative results. All complex modes come with their negative complex conjugate which we do not show.
In Table (1) we show the lowest lying quasinormal mode of the Schwarzschild-de Sitter black hole discussed in section 4.3, for l = 1, 2, and the second lowest for l = 2. We take as parameter ΛM 2 = 3M 2 /L 2 , for ease of comparison with [38] which computed the same QNMs with a sixth order WKB approximation. Results agree to a high precision, except that we find an extra set of purely imaginary modes (in the top left entry this is the dominant mode).