Bernstein spectral method for quasinormal modes and other eigenvalue problems

Spectral methods are now common in the solution of ordinary differential eigenvalue problems in a wide variety of fields, such as in the computation of black hole quasinormal modes. Most of these spectral codes are based on standard Chebyshev, Fourier, or some other orthogonal basis functions. In this work we highlight the usefulness of a relatively unknown set of non-orthogonal basis functions, known as Bernstein polynomials, and their advantages for handling boundary conditions in ordinary differential eigenvalue problems. We also report on a new user-friendly package, called \texttt{SpectralBP}, that implements Berstein-polynomial-based pseudospectral routines for eigenvalue problems. We demonstrate the functionalities of the package by applying it to a number of model problems in quantum mechanics and to the problem of computing scalar and gravitational quasinormal modes in a Schwarzschild background. We validate our code against some known results and achieve excellent agreement. Compared to continued-fraction or series methods, global approximation methods are particularly well-suited for computing purely imaginary modes such as the algebraically special modes for Schwarzschild gravitational perturbations.


Introduction
Black holes in general relativity are simple spacetime objects, fully specified by only a handful of constants.When the spacetime around black holes is disturbed by surrounding complex distributions of matter and fields, as they are found in nature, these spacetime disturbances generically evolve in the form of damped oscillations known as quasinormal modes (QNMs).
Quasinormal modes are the characteristic ringing of spacetime around black holes.They are independent of the initial excitation that generated them, dependent only on parameters of the black hole.A wealth of information can be extracted from the quasinormal mode spectrum of a black hole, so they serve as probes for the validity of general relativity and its extensions in the strong gravity regime.Two excellent reviews on the topic with an emphasis on astrophysics can be found in [1] and [2].A review on higher dimensional black holes and their connections to strongly coupled quantum fields can be found in [3].
In general, the quasinormal mode spectrum of a black hole comes from solving an ordinary differential equation (ODE) eigenvalue problem.These usually take the form of a Schrödinger-like equation, where r * is called a tortoise coordinate.
Send offprint requests to: Various numerical methods have been developed to solve (1), such as the WKB approach, shooting methods, continued-fraction methods, and the use of Pöschl-Teller potentials.A review article with an emphasis on this topic can be found in [4].In this paper, we shall be solving (1) using a pseudospectral method.
The use of spectral and pseudospectral methods in gravitational problems is well-established [5,6], and have been applied to numerous numerical experiments such as [7,8,9] to name a few.Here we extend this library of methods to include the Bernstein polynomial basis, which has particular properties that lend its use to mixed-type boundary-value problems.
Likewise, Bernstein polynomials have been used as a function basis in the numerical solution of various differential [10,11,12,13,14], fractional differential [15], integral [16,17,18,19], integro-differential [20,21] and fractional integro-differential [22] equations.Multiple methods have been deployed in this context, such as the Bernstein-Petrov-Galerkin (BPG) method, the collocation method, operational matrices and direct integration.Our work extends the range of the Bernstein basis by exploring its use in ODE eigenvalue problems.
The aim of this paper is two-fold.First, it is a primer on how Bernstein polynomials (BPs) may be used for boundary-value problems in a general relativity setting.Second, it is an introduction to a Mathematica package we call SpectralBP that implements the pseudospectral method based on Bernstein polynomials.For examples and benchmarks, we have applied SpectralBP to a selection of eigenvalue problems in quantum mechanics and general relativity: the infinite square well, harmonic and anharmonic oscillators, and quasinormal modes of various fields in a Schwarzschild black hole.Particularly noteworthy is that SpectralBP is able to find eigenvalues with modest resources where other numerical methods find with difficulty, such as the algebraically special modes for gravitational perturbations of the Schwarzschild geometry [23,24].As will be explained below, this should be expected of any spectral method for eigenvalue problems.
The method introduced in this paper, and the accompanying Mathematica package, has seen use in general relativity [25,26,27,28,29,30,31,32,33] and quantum mechanics [34].It has been particularly useful in finding purely imaginary quasinormal modes [32,28], finding new branches of solutions hitherto unknown [29,26] and show novel and critical behaviors like spectrum bifurcation [27] and instability [30].In quantum mechanics applications, it has been shown to generate exceedingly accurate solutions where other methods require vast resources in memory and compute time [34].
Consider an n × n matrix of linear differential operators L(u, ω) dependent on a single independent variable u and polynomial in the eigenvalue ω of some maximal integer order m, Li,j (u, ω) = fi,j,0 + ω fi,j,1 and let Φ(u) be a vector of n functions dependent on u Φ(u) = (ϕ 1 (u), ϕ 2 (u), . . ., ϕ n (u)) T . ( We wish to solve the following eigenvalue problem for ω, provided the problem satisfies the following criteria: 1.The domain of the solution is compact and analytic over its whole domain.(u ∈ [a, b]) 2. The boundary conditions for all eigenfunctions ψ i (u) specifies that lim u→a ψ i (u) ∼ (u − a) q and lim u→b ψ i (u) ∼ (b − u) r for some q, r ≥ 0. 3. The eigenvalues of ω form a discrete spectrum.
The calculation of the bound state energies of quantum mechanical particles and the quasinormal modes of black hole spacetimes are examples of such a problem.
To solve (4) we use a pseudospectral method, in which the solution of the differential equation is approximated as a weighted sum of a set of basis functions, say {ϕ i (r)}, as in, This renders the initial differential problem into a system of algebraic equations the set of expansion coefficients {C i } must satisfy.Since (4) is linear, these algebraic equations can be cast as a matrix equation generically of the form of a generalized eigenvalue problem (GEP), We have developed a Mathematica package we call SpectralBP, written to streamline the numerical solution of ODE eigenvalue problems.The package utilizes the Bernstein polynomials, and the properties which make them particularly powerful in the context of boundary value problems.
A similar Mathematica package can be found in [35].It is a pseudospectral method which uses a Chebyschev polynomial basis, called QNMSpectral.This open-source package served as the initial inspiration for our work, and so the two codes unavoidably overlap in some of their functionality.We developed SpectralBP to be a superset of QNMSpectral, with the intent of developing a spectral solver not just specifically tailored to quasinormal mode calculations.It also serves to introduce the Bernstein method to the general relativity community.Aside from methods specifically tied to the Bernstein basis, SpectralBP also implements a novel algorithm for efficiently tackling transcendental and polynomial eigenvalue problems that we shall discuss in detail in a future paper [36].
This paper is organized in two parts.We first establish how the Bernstein polynomial basis may be used in ODE eigenvalue problems with boundary conditions.In Section 2, we fix our notation and enumerate the properties of the Bernstein basis relevant to the method.In Section 3, we explain how the Bernstein basis is appropriate in handling boundary conditions.In Section 4, we review standard methods for translating (4) into a generalized eigenvalue problem using a collocation method.We then enumerate various positives and negatives the Bernstein polynomial basis has compared to other bases like Fourier or Chebyschev in Section 5.
The rest of the paper involves the implementation and application of SpectralBP.Section 6 introduces the SpectralBP package and its general features.We then show in detail how SpectralBP can be used in Section 7 and Section 8, introducing functionalities of the package by working out some model problems in quantum mechanics and calculating quasinormal modes respectively.In Section 9, we look at the algebraically special modes of the Regge-Wheeler equation.In the final section, we show miscellaneous details implemented in SpectralBP: closed-form expressions of the spectral matrices, matrix inversion, and eigenfunction calculation and manipulation.

Bernstein polynomials
We review some of the key properties of Bernstein polynomials.We shall not be exhaustive and select only those properties useful to the development of SpectralBP.This section shall also fix our notation for the rest of the paper.A useful reference can be found in [12], which describes all of the properties listed here using a Bernstein basis over the interval The Bernstein basis of degree N defined over the interval For convenience, we also set B N k (u) = 0 and The Bernstein basis of degree 10 is shown in Figure 1.It is clear that at the boundaries u = a and u = b, Bernstein polynomials satisfy The derivative of a Bernstein polynomial of degree N can be expressed in terms of Bernstein polynomials of degree N − 1, satisfying the following recurrence relation, Repeated differentiation also gives A Bernstein polynomial of degree N can be expressed as a sum of Bernstein polynomials of a higher degree [37], The integral of each basis polynomial in a Bernstein basis of degree Finally, the product between two Bernstein polynomials can be expressed as single Bernstein polynomial of higher degree,

Boundary conditions
When using the Bernstein basis in mixed-type boundary-value problems, we shall see that the boundary conditions act only on a subset of the Bernstein basis.This lets us independently solve the boundary conditions separately, making the Bernstein basis particularly useful in mixed-type boundary-value problems.For the particular boundary-value problem described in Section 1, the Bernstein method reduces to a form in which each basis function satisfies the boundary conditions.We begin by approximating the solution ϕ(u) as a weighted sum of Bernstein polynomials, Let there be q boundary conditions on u = a and r boundary conditions on u = b of the following form, These constants may be interrelated.A common example would be a two-point boundary value problem of a second-order differential equation subject to mixed linear boundary conditions, which fixes a 0 , a 1 , b 0 , and b 1 .Combining ( 10) and ( 14), the m th derivative of ϕ(u) is given by We use (8) to simplify evaluating ϕ(u) at the boundaries.At u = a and u = b, we get and Thus, the boundary conditions act only first q and last r of the Bernstein basis, whose expansion coefficients are fixed via the matrix equations where and When the differential operator is linear, the modified ODE eigenvalue problem determines the rest of the expansion coefficients, where the residual function g(u, ω) is given by We consider the case where g(u, ω) vanishes, or equivalently We arrive at an ODE eigenvalue problem identical to the one we started with, but over a smaller set of basis functions It should be noted that for more standard basis functions, imposing the boundary conditions considered in (15) would involve the entire basis set.To determine the expansion coefficients, the differential equations and the boundary conditions must be solved simultaneously.In the Bernstein basis, the boundary conditions act only on the first q and last r basis polynomials, and we get their corresponding expansion coefficients for free even before considering the ODE.Though we do not prove that this advantage is unique to the Bernstein basis, we believe that any other basis must behave like Bernstein polynomials to enjoy it.That is, the nth basis function of a basis of size N must asymptote to (u − a) n towards the lower boundary and to (b − u) N −n towards the upper boundary.We express a similar sentiment for other basis functions where the condition (25) would make the residual function vanish.In the Bernstein basis, the problem is simplified since each basis polynomial satisfies the boundary conditions exactly.
Finally, we note that when the differential operator is not dependent on ω, equation ( 23) serves as a general recipe for solving boundary value problems using Bernstein polynomials.One may modify the many methods found in Section 1 to solve for the remaining undetermined coefficients.

Pseudospectral method
In this section, we review how one starts with the ODE eigenvalue problem in (4) and end up with the generalized eigenvalue problem in (6).We derive a general recipe for mapping a differential operator and function pair to a matrix and vector pair ( M(ω), C) via a collocation method in the Bernstein basis, whose closed form can be found in the last section.In the context of Chebyschev basis polynomials and Fourier basis functions, the standard reference is [38].
We start with a linear eigenvalue ODE, then show how it can be extended to polynomial eigenvalue ODEs.We extend this to include problems involving a set of dependent functions.We elaborate on special cases in A, used to convert the polynomial generalized eigenvalue problem to an eigenvalue problem.

Linear eigenvalue problem
Consider the ODE eigenvalue problem in (26), specifically of the form To arrive at a spectral matrix of size N +1, we expand the basis degree to A straightforward implementation of the collocation method would be to define a grid of Since the first q or last r Bernstein basis functions dominate the behaviour of the solution near the boundaries, we propose instead to select collocating points in the region dominated by the basis functions whose weights are still unknown.
As an illustrative example, consider the case when N = 10, q = 30 and r = 30.One can imagine rescaling a solution ϕ(u) finite at both boundaries via the transformation, The basis of ϕ(u) is in Figure 1 while the basis of φ(u) is in Figure 2. Its derivatives are similarly localized.We construct and then select grid points q through N + q + 1.
Let us now endeavour to convert the differential operator and function pair ( f (u), ψ(u)) into a matrix and vector pair A generic term in the f (u)ψ(u) is of the form f n (u) d n ψ du n .Combining (10) and (28), we get We may assign a vector to each term in f (u)ψ(u) with the condition that the differential operator is satisfied at each collocation point, where C k = C k+q and the matrix components of T (n) are given by for j, k ∈ {0, 1, . . ., N }.
To use (34), each Bernstein basis polynomial of degree N max − n through N max must be evaluated at each collocation point.Since in many applications, n max ≪ N , it would be numerically cost efficient to use (11) and rewrite (34) in terms of a single Bernstein basis degree, as in By choosing this degree to be N max , only a subset of the Bernstein basis needs to be evaluated at each collocation point−specifically those indexed in the range [q − min(n max , q), N + q + min(n max , r)].Thus, The ODE linear eigenvalue problem in (27) may be written as a generalized eigenvalue problem,

Polynomial eigenvalue problem
Consider a polynomial eigenvalue problem of order m., Using the recipe discussed in the previous section, this corresponds to an eigenvalue problem of a matrix pencil of order m, We linearize the matrix pencil by defining the following matrices, and the vector, This transforms the matrix pencil (39) to another GEP, For clarity, we typeset matrices and vectors generated from linearizing a matrix pencil by a calligraphic typeface.Generalized eigenvalue problems are more difficult to solve than regular eigenvalue problems.We describe a method to transform the above GEP to an EP in A.2 contingent on the invertibility of either M 0 or M m , leading to a modest improvement in speed.

Polynomial eigenvalue problem over several dependent functions
Consider the full problem in Section 1.In matrix form, this becomes the set of simultaneous equations, where each matrix M j,k (ω) is constructed by linearizing the matrix pencil of the kth dependent function of the jth equation, as in The set of simultaneous equations can be written as a single matrix equation by defining the following matrices, and vector, We arrive at the GEP of the full problem introduced in Section 1, The GEP of the full problem is much more complicated.Unlike the GEP of the previous subsection, it can be shown that M′ is always singular, as we show in A.2.

Advantages and disadvantages of the Bernstein basis
Having elaborated on how the Bernstein basis fits into solving a partial differential problem like (4), we discuss in this section what these properties cost and afford us, and how they compare to more standard basis functions.We also discuss some results that may be found in A.4.One may read through that Section first, and then return here.
1. Bernstein polynomials are not orthogonal.This follows from (12) and (13).This complicates an extension of the current method to partial differential equations, where the weights may be made to vary in time.2. The Bernstein basis polynomials depends on the basis degree.We cannot naively apply derivatives without costing us additional numerical resources.We need to fold in an application of (11) so that we remain in a single common basis degree.There is no operation similar to (11) for classical orthogonal polynomials, because those basis functions do not depend on the size of the basis.8), (10) and ( 12) are well-known for classical orthogonal polynomials and the Fourier basis.A simple product formula like (13) exists for Chebyschev and Fourier basis.
5. The specific form of these properties gives the Bernstein basis an advantage over other basis functions when dealing with mixed boundary value problems outlined in Section 3.
In the Bernstein basis, the boundary constraints only act on a subset of the basis set, whose weights can be fully determined independently of the differential equation.Such a luxury is not enjoyed by more standard basis functions.
For classical orthogonal polynomials and the Fourier basis, imposing the boundary conditions (15) would involve the entire basis set.Solving the differential equations and the boundary conditions must be done simultaneously.6.The lack of a residual term in (26) and the lack of additional constraints on the expansion coefficients lets us write down the algebraic equations these expansion coefficients must satisfy as a generalized eigenvalue problem in Section 4. 7.There are manipulations which can only be easily done in the Bernstein basis, discussed in A.4.For example, a tau method using Chebyschev polynomials can impose the boundary condition lim u→a ψ(u) ∼ (u − a) exactly.However, one cannot naively divide out a (u − a) term-by-term, since each Chebyschev polynomial is finite at the lower boundary.Such a rescaling can be exactly carried out in the Bernstein basis, as shown in A.4.This lets us calculate the weighted L 2 -norm of a function in the Bernstein basis in closed-form, even in cases where the weight has a pole of integer degree at the boundaries.This is useful, for example, when normalizing wavefunctions in a compactified coordinate system, as in Section 7.2.8.The numerical convergence of the Bernstein basis has been established in the context of other differential problems.Interestingly, in some cases, the Bernstein method would outperform other basis functions (including Chebyshev and Fourier) in terms of numerical cost or the accuracy of the solutions [15,17,18,19,21].We do not perform a similarly comprehensive analysis here, concentrating instead on general ideas on how the Bernstein basis may be adapted to ODE eigenvalue problems and introducing the package SpectralBP.Though we do demonstrate numerical convergence for some of the cases we tackle below.

The SpectralBP package
The SpectralBP package uses the properties of the Bernstein basis, written to streamline the calculation of the eigenvalues and eigenfunctions of (4).It is primarily distributed as a Mathematica paclet and is publicly available [39].SpectralBP commands are documented, and the package is bundled with two tutorial notebooks.After installation, the details and options of each command may be explored by prefixing a command with a question mark, as in ?GetModes, similar to built-in commands in Mathematica.
There are three types of commands in SpectralBP: Get commands, Compare commands and Print commands.The basic work flow is as follows.
1. Begin with some ODE eigenvalue problem which may not satisfy the 3 properties required in Section 1. 2. If the domain of the eigenfunctions ψ ′ i (x) is not compact, define an invertible change of variables f (x) = u so that the domain in u is compact.3.If the resulting eigenfunctions are non-analytic, one may rescale as in so that the resulting eigenfunctions ψ i (u) are analytic.
One also defines f i (u) so that all eigenfunctions ψ i (u) satisfies the same boundary conditions.The result should be an eigenvalue problem described in Section 1. 4. Use Get commands to calculate eigenvalues and eigenfunctions at different BP orders.5. Use Compare commands to filter out spurious eigenvalues and eigenfunctions.6. Use Print commands to quickly glean off information from the prior calculations.
We will discuss each command type in the following subsection before going into applications.Example notebooks can be found in the next two Sections.

Get commands
The first input of a Get command is a list of differential equations.The command automatically identifies the dependent functions, the independent variable and the eigenvariable.The command halts whenever it identifies more than one independent variable or eigenvariable, or whenever the number of dependent functions underdetermine or overdetermine the problem.
There are three Get commands, 1. GetModes[eqn,N]: Calculates the eigenvalues of the ODE eigenvalue problem stored in eqn using a basis degree of N.

GetEigenfunctions[eqn,modes,N]:
Calculates the eigenvectors corresponding to each eigenvalue in the list modes, using a basis degree of N. As discussed in the Appendix, we advise that N be identical to be basis degree the eigenvalues in modes were calculated.3. GetAccurateModes[eqn,N1,N2]: Calculates the eigenvalues using basis degrees of N1 and N2, then applies a CompareModes command to filter the eigenvalues.
By replacing the basis degree inputs with a pair of numbers, which we call a basis tuple of the form {N,prec}, eigenvalues are calculated using a basis degree of N using prec-precision numbers.That is, an alternative input scheme for the above commands is given by, In calculating the eigenvalues and eigenvectors, Get commands must be supplied with the correct domain and boundary conditions.These are controlled by 4 options, 1. LowerBound and UpperBound: defines the domain LBPower and UBPower: defines the leading polynomial power q, r at each boundary, which defaults to q = 0 and r = 0.
The option Normalization lets one choose how eigenfunctions are normalized.The option may have 4 values, 1. "UB": the coefficient of the leading polynomial expansion of the eigenfunctions at b to 1.

"LB": the coefficient of the leading polynomial expansion
of the eigenfunctions at a to 1. 3. "L2Norm": the L 2 -norm of the eigenfunctions to 1. 4. {"L2Norm",{A,B,C}}: the L 2 -norm of the eigenfunctions to 1, with a weight function underneath the integral of the form The option FinalAsymptotics lets one change the outputted eigenfunctions' asymptotics, according to manipulations detailed in A.4.

Compare commands
The spectrum calculated from a finite basis degree will be filled with either eigenvalues that have not converged or spurious eigenvalues.We have provided two ways to filter these out.These are the two Compare commands, 1. CompareModes[modes1,modes2]: Checks whether eigenvalues in the two spectra inputted share common digits, then keeps only eigenvalues that share at least 3 digits.2. CompareEigenfunctions[eqn,{modes1, modes2},{N1,N2}]: Calculates the eigenfunctions of the eigenvalues approximately common to modes1 and modes2 using a basis degree of N1 and N2 respectively.If the L 2 -norm of their difference is less than 10 −3 , the eigenvalues are kept.
There are two relevant options, 1. Cutoff: controls the minimum number of common digits for eigenvalues to be kept, which defaults to 3. 2. L2Cutoff: controls the maximum difference between two eigenfunctions, of the form 10 −n , for their corresponding eigenvalues to be kept, which defaults to n = 3.
We call eigenvalues of different spectra that share a Cutoff-number of common digits approximately common.

Print commands
There are four Print commands, 1. PrintFrequencies[modes]: plots the eigenvalues in modes on the complex plane.2. PrintEigenfunctions[eqn,modes,N]: plots the real and imaginary parts of the corresponding eigenfunctions.

PrintTable[convergedmodes]: generates a table
of eigenvalues, categorizing them into purely real, purely imaginary, and complex eigenvalues.Groups together eigenvalues satisfying ω * = ω and ω * = −ω.The input must be a pair of lists of approximately common eigenvalues, usually coming from the output of a CompareModes command.4. PrintAll[eqn,convergedmodes,N]: a shortcut to do the previous three commands in a single command.
There are three relevant options, 1. FreqName: specifies the symbol for the eigenvariable, which defaults to ω. 2. NSpectrum: specifies how many eigenvalues would be plotted, which defaults to plotting everything.3. NEigenFunc: specifies how many eigenfunctions would be plotted, which defaults to plotting everything.
The PrintTable command automatically only prints out significant digits, defined to be the digits common to both the spectra inputted.When the inputted spectra comes from two adjacent basis degrees, say N and N + 1, the right-most digits of the output may be incorrect.This occurs because the absolute error of the two spectra overlap.
We recommend using basis degrees that are far apart in the sense that the absolute error of the higher basis degree spectrum is much smaller than the absolute error of the lower basis degree spectrum.Although the practice would be numerically more costly, in this way we increase our chances that the right-most significant digit outputted is correct.

Summary of implementations
In Table 1, we summarize the different inputs needed to solve the ODE eigenvalue problems that we shall look at in the succeeding Sections.Hopefully, in the examples considered in the proceeding Sections, one is left with an impression of the general-purpose applicability and ease-of-use of SpectralBP.As shall be demonstrated, three lines of code can yield a wealth of information about the considered ODE eigenvalue problem.The difference between the examples given is just swapping in and out of differential equations, applying certain change of variables in cases where the domain is infinite, and specifying the necessary boundary conditions.

Applications in Quantum Mechanics
We first illustrate how SpectralBP is used by working through standard problems in quantum mechanics.We solve Table 1.Input scheme for the various eigenvalue ODE problems discussed in Section 7, Section 8 and Section 9.The potential V * was chosen to be (53) for the base infinite square well problem, and (56) for the perturbed infinite square well problem.The potential V † was chosen to be (58) for the quantum harmonic oscillator problem.The potential V ‡ was chosen to be (71) as the PT -symmetric anharmonic potential for specific values of λ, and (72) as the quartic anharmonic potential for specific values of β.The different variables used mark certain coordinate transformations effected to compactify an infinite domain.

Problem eqn LBPower UBPower
Infinite square well Harmonic oscillator Notebook 1 : A simple Mathematica notebook implementation of SpectralBP for the infinite square well problem.for the eigenenergies and eigenfunctions of the infinite square well and quantum harmonic potentials numerically in the first two subsection.Calculations are compared with well-known analytic results, as can be found in standard quantum mechanics textbooks like [40].
For the last two subsections, we compute the eigenenergies of the anharmonic potentials considered in [41] and [42].We compare ground state eigenenergies calculated with SpectralBP to the results of the aforementioned papers, which were both calculated perturbatively using a combination of Padé approximation and Stieltjes series.In [42], Milne's method [43] was used as an independent test.

Infinite square well
Consider the time-independent Schrödinger equation 1 2 For the infinite square well, the potential is chosen to be Its eigenenergies are The domain of solutions is the interval [0, 1] with boundary conditions, 7.1.1SpectralBP−basic implementation.
A simple implementation to solve the infinite square well problem is schematically found in Notebook 1. Lines 2 and 3 solves the ODE eigenvalue problem (52) with potential (53) using basis degrees 50 and 80 respectively.
The boundary conditions (55) are set by the option values which must be specified whenever eigenvalues and eigenvectors are calculated.Line 4 selects eigenvalues that are approximately common to modes50 and modes80.As described in the Section 6.3, this may serve as input for the PrintTable command in line 7.We have chosen to rescale the eigenenergies in lines 5 and 7 so that the output would be the first 10 perfect squares.
Line 6 plots the eigenfunctions of the inputted spectrum of the lowest three eigenvalues of modes50 using a basis degree of 50.The Print commands found in the last 3 lines output Figure 3.
As described in Section 6.3, the command PrintTable only prints out significant digits.As an illustrative example, consider the lowest rescaled eigenenergies.The absolute error for modes50 is 3.27 × 10 −22 and the absolute error for modes80 is 4.97 × 10 −31 .The PrintTable compares the two eigenvalues and detects a difference of ∼ 10 −22 , and prints out the eigenvalue up to the 21 st decimal place.

A note on machine precision.
As described in Section 6.1, one may use arbitrary precision numbers by inputting a basis tuple of the form {N,prec}.This would calculate eigenvalues using a basis degree of N with prec-precision numbers, as in Notebook 3.
Notebook 3 : An implementation demonstrating the use of arbitrary precision numbers in SpectralBP.
The PrintTable command in line 3 outputs Figure 4.The number of common modes remain at 28 (not shown), but there are more significant digits for the lowest eigenenergies.This is because the error due to floating point arithmetic at machine precision is generally small enough to resolve approximately common eigenenergies between basis degrees.When higher precision numbers are used, this error is pushed down further and may reveal more significant digits.The absolute error from approximating the solution space in a finite polynomial basis eventually dominates, and may only be corrected by using higher and higher basis degrees.
Briefly, increasing machine precision increases the significant digits (up to a point) while increasing the Bernstein The first 10 eigenenergies calculated using a basis degree of 50, plotted on the complex plane.Middle: (PrintEigenfunctions) The eigenfunctions of the first 3 eigenenergies, calculated using a basis degree of 50, normalized according to their L 2 -norm.Bottom: (PrintTable) Rescaled eigenvalues common to basis degrees of 50 and 80.There are 28 eigenenergies that share a minimum of 3 significant digits (not shown).We tabulate only the lowest 10.The spectrum calculated is in excellent agreement with (54).For completeness, let us explore the case when the exact solution is non-analytic.Suppose we perturb the potential by lifting half of the infinite square well, Fig. 4. Calculated eigenvalues common to basis tuples {50, 50} and {80, 80} (described in Section 6.1).There are 28 eigenenergies that share a minimum of 3 significant digits (not shown)−similar to Figure 3−while the number of significant digits for the lower eigenvalues have increased.
The exact solution can be derived by starting with a pair of free particle solutions at 0 ≤ x ≤ 1/2 and 1/2 ≤ x ≤ 1, then imposing the correct boundary conditions at the walls of the infinite square well and continuity relations x = 1/2.One then finds that for the boundary conditions and the continuity relations to be satisfied, the eigenenergies must be solutions to the transcendental equation, Exact solutions are non-analytic since they are not twice differentiable at x = 1/2.On the other hand, one may simply swap in the potential (56) and use a GetAccurateModes to numerically solve for these eigenenergies.We benchmark SpectralBP against the Mathematica in-built function NSolve in Table 2.
NSolve is a zero-finding algorithm, which we use to find solutions to (57).There is great agreement between the two methods.
The non-analyticity of the solutions has adversely affected how quickly the eigenenergies converge to the correct values, which is expected from a spectral method.SpectralBP was able to find all eigenenergies below 1000.On the other hand, NSolve will not find the eigenenergies indicated by *'s by default.These roots are sensitive since one must start close to the them so that NSolve can find them.The eigenenergies indicated by *'s were found by sampling the range [0,1000] with a resolution of 0.01.
We note that we have chosen odd basis tuples in the calculation so that the corresponding collocation grids avoids the point x = 1/2.Choosing even basis tuples degrades the accuracy of odd-numbered eigenenergies, and one would need to reach a basis degree of around 400 to determine the ground state energy accurate to 3 digits.

Quantum harmonic oscillator
Consider the harmonic oscillator potential, Its eigenenergies are, The domain of the solutions is the entire real line (−∞, ∞) with boundary conditions

Compactification and boundary conditions.
One may swap in the harmonic oscillator potential in the example notebooks we have presented to calculate eigenenergies and eigenfunctions, except one must include an additional step of compactifying the domain.Let us compare the spectrum calculated using two different ways of compactifying the interval (−∞, ∞).The first, has a domain of [0, 1].Some comments are in order.First, note that the exact solution in both compactified coordinates is flat at both Table 3.Comparison between compactifying using (61) and (62), using Bernstein tuples {50, 50} and {100, 100} (described in Section 6.1).For conciseness we indicate eigenergies found using (61) with a dagger † , and mark in square brackets the additional significant digits calculated using (62).Compactifying using (62) performs better, which finds more eigenvalues with more significant digits.
n En where a k , b k are the corresponding boundary locations for k = 1, 2. Second, note that the potential is singular at the boundaries in both compactified coordinates, with (64) A consequence of using the collocation grid we proposed in Section 4.1 is that we have avoided evaluating at these singularities by expanding the Bernstein basis order and choosing collocation points in the interior of the relevant domain.
Finally, we observe a dependence on the rate of convergence of the method with respect to different coordinate transformations, as can be seen in Table 3.We may attribute this discrepancy on how features such as maxima and nodes of higher energy eigenfunctions are distributed on the compactified coordinates in relation to how the collocation points are distributed on the compactified coordinates.
Consider the distance of the right-most maxima or node relative to the upper bound of a high energy eigenfunction for either transformation, (65) That is, in proportion to the length of the interval, these features are closer to the upper bound of the interval with (61) than in The same is true for the lower bound.Thus, a collocation grid defined on v 1 is unable to resolve higher energy eigenfunctions compared to v 2 since the collocation points are less densely located on where the maxima and nodes are expected to appear -ie., closer to the edge in proportion to the length of the interval for (61) than in (62) We note that a transformation such as 'spreads' these features further away from the upper bound and lower bound.An identical calculation on v 3 yields more accurate eigenenergies than on v 2 as well as finding more eigenenergies (upto 26).

Eigenfunctions−normalization and manipulation.
Consider the eigenfunctions calculated from (61) and (62).To properly normalize the eigenfunctions in the original coordinates x, one must introduce a weight function underneath the integral of their L 2 -norms in the compactified coordinates, respectively of the form As described in Section 6.1, the option value for Normalization should be {"L2Norm",{1,-1,-1}} for both v 1 and v 2 .The eigenfunctions of the three lowest eigenenergies in Table 3 may be calculated using the GetEigenfunctions command.The output is a Bernstein polynomial in the compactified variable v 2 , which may reverted to the original uncompactified coordinates by a change of variables.The eigenfunctions in x are plotted in Figure 5 together with their absolute error compared with the exact eigenfunctions.The absolute error is bounded from above, with a maximum difference between 10 −9 − 10 −11 .

Anharmonic potentials
We now benchmark SpectralBP against other numerical methods, here in the context of anharmonic potentials.We perform test calculations also done in [41] and [42], in which the time-independent Schrödinger equation has been rescaled such that, and the anharmonic potentials, were considered.In the papers cited, Padé approximation and Milne's method [43] were used to calculate the ground state energies.
The potential (71) is interesting.Although the corresponding Hamiltonian, isn't hermitian, its eigenenergies remain real and positive.This is because of its underlying PT symmetry [44], in which combining parity, P : p → −p and x → −x, and time reversal, T : p → −p, x → x, and i → −i, transformations leaves H invariant.
For both potentials, we compactify our domain via the transformation in (62).To recreate Table II of [41], we set λ = 1/7 and β = 40/49 and use basis tuples {250, 250} and {300, 300} (described in 6.1).The spectra of both potentials are found in Table 4.For a more direct comparison to Table II of [41], we use Equations ( 8) and ( 9) of [41] to calculate P (λ 2 ) and P (β) for the ground state energy.Comparing the two values coming from both basis tuples for significant digits, and we arrive at the expressions P (λ 2 ) = 5.524167213060 [22] where the last expression goes on for 21 more digits.These values are in excellent agreement with the values calculated in [41].The digits enclosed in square brackets are additional significant digits calculated by SpectralBP.
The anharmonic potential (72) was used in [42], but for different values of β.We calculated spectra using basis tuples of {150,150} and {200,200}, keeping only eigenvalues with at least 5 significant digits.In Table 5, we show only the ground state energies for a direct comparison of Table II and Table IV of [42].
The results are in great agreement with the "Exact" values calculated in [42], which were calculated using Milne's method [43].At the digits where they differ, which we have indicated in square brackets, the difference is within the error bars in both tables.The calculations took an average of 68 seconds each, running in a single 2.50 GHz Intel i5 Core with 8.00GB RAM.
With modest resources, we are able to calculate the ground state energies to high precision.This is simultaneous with an abundance of excited state energies; the calculation at β = 1/10 yielded 47 eigenenergies with at least 5 significant digits, while the calculation at β = 100 yielded 69 eigenenergies with at least 5 significant digits.

Applications in Quasinormal Modes
In general relativity, spacetime itself is treated as a dynamical entity, interacting with the matter that is placed within it.Thus, black holes found in nature are always interacting with complex distributions of matter and fields around them.In active galactic nuclei, accretion disks transport matter inward and transport angular momentum outward, heating the accretion disk into a hot plasma and immersing the black hole in a complex gravitational and electromagnetic system.Even in the absence of matter and fields, the black hole interacts with the vacuum around it, slowly evaporating due to Hawking radiation.
The standard treatment is to decompose the spacetime as in where the metric g 0 µν is that of an unperturbed black hole, such as the Schwarzschild or Kerr solution.In the linear approximation δg µν ≪ g 0 µν (so called because the perturbing metric δg µν does not back react with the background metric), these small perturbations generically take the form of damped oscillations known as quasinormal modes.When g 0 µν is spherically-symmetric, the equations for δg µν reduce to one-dimensional wave equations in certain potentials.These are the famous Regge-Wheeler and Zerilli equations for oddand even-parity perturbations, respectively.
Quasinormal modes arise as the characteristic ringing of spacetime as it is perturbed by some external field.For a given external field, these oscillations are independent of the initial excitation, their vibrations and damping specified solely by the mass, spin and charge of the black hole.As such, quasinormal modes are used as probes for the validity of general relativity in the strong gravity regime.
From a more theoretical perspective, quasinormal modes provides a test for the linear stability of more exotic Table 4. Spectra for anharmonic potentials found in (71) and (72), with λ = 1/7 and β = 40/49, calculated using basis tuples {250, 250} and {300, 300} (described in Section 6.1).Only common eigenvalues with at least 5 significant digits were kept.For (72), there are 79 such eigenvalues.We have chosen to show only the lowest 10 eigenvalues up to 40  when all quasinormal modes are damped (Im(ω) ≤ 0), the spacetime is linearly stable.In the context of AdS/CFT duality, the onset of instability of the AdS spacetime corresponds to a thermodynamic phase transition in CFT.
Review articles on quasinormal modes in an astrophysical setting -black holes, stars, and other such compact objectswe cite [1] and [2].An emphasis on higher dimensional black holes and their connection to strongly coupled quantum fields is in [3], while [4] emphasizes on the various numerical and analytical techniques that have been developed to calculate quasinormal modes.The papers [5,6] focus on the application of spectral and pseudospectral methods in gravity, of which SpectralBP is an example of.

Regge-Wheeler equation
In Section 6, we described a general work flow starting from an ODE eigenvalue problem.In this subsection we go through the first 3 steps of this work flow, starting from a standard ODE eigenvalue problem for quasinormal modes.We focus on the Regge-Wheeler equation as an illustrative example; a treatment of the Zerilli equation would proceed in a similar manner.The Regge-Wheeler equation describes axial or odd-parity perturbations of the Schwarzschild metric of mass M linearly coupled to a perturbing field of spin s and angular momentum l, where r * = r + r s ln(r/r s − 1).We are interested in solutions of the form Φ(t, r * ) = R(r) exp(−iϵt).This then turns (76) into the ODE eigenvalue problem of the form, with ϵ = 2M ω.The domain of the solutions relevant to us is non-compact, stretching from the black hole horizon at r = 1 to spatial infinity at r = ∞.Note also that the solutions are non-analytic.The coordinate singularity at r = 0 and the black hole horizon at r = 1 are both regular singular points of the ODE, while spatial infinity r = ∞ is an irregular singular point of the ODE.We may peel away the non-analytic parts by rescaling out the asymptotic behaviour of R(r) at the black hole horizon and at spatial infinity.The asymptotic behaviour of R(r) at r = ∞ can be easily determined to be R out (r) ∼ r iω exp(iωr) R in (r) ∼ r −iω exp(−iωr), (78) where we have indicated in superscript which solution is outgoing or ingoing at spatial infinity when the time dependence is restored.
Since the singularity at r = 1 is regular, we may write an indicial equation f (x) = 0 at r = 1.This can be shown to be simply which defines two solutions around r = 1, where we have indicated in subscript which solution is outgoing or ingoing at the black hole horizon when the time dependence is restored.We expect a perturbation to come from a finite location outside the black hole.As this perturbation propagates, we expect it to either fall into the black hole or out into spatial infinity.This defines the behaviour of the causal solution, and corresponds to the quasinormal mode boundary conditions An acausal solution would contain parts that are either propagating out of the black hole, or propagating in from spatial infinity.We rescale out the non-analytic parts of the desired solution, Table 5. Ground state energies calculated using the anharmonic potential (72) for different values of β, using basis tuples {150, 150} and {200, 200} (described in Section 6.1).For conciseness, we have enclosed in square brackets additional significant digits calculated by SpectralBP compared to an application of Milne's method in [42].leaving us with a differential equation in ϕ(r).We note that the additional factor of r iω is there to cancel out the asymptotic behavior of (r − 1) −iω around spatial infinity.
Explicitly, the rescaled solution at the boundaries have the following behaviours: For generic values of ω, these four solutions have very distinct behaviours.Consider the acausal solutions near their corresponding limits, When Im(ω) = 0, both solutions are highly oscillatory.Thus, the boundary conditions, filters out both undesired acausal solutions, since these solutions cannot be approximated in the Bernstein basis of finite degree.Thus, with the boundary conditions in (87), we may identify our solutions to correspond to quasinormal mode eigenfunctions, ϕ(r) = ϕ out in (r).
(88) Finally, we compactify the region [1, ∞) to [0, 1] via the change of variables r → 1/u, leaving us with This change of variables moves the regular singularity at r = 0 to u = ∞ and the irregular singularity at r = ∞ to u = 0.
We may use equation (89) as the ODE eigenvalue problem we feed into SpectralBP.However, we may improve our calculations with the transformation ω → iλ, which yields an ODE eigenvalue problem whose coefficient functions are all real, The spectral matrices constructed from (90) are strictly real.This has two consequences.First, the calculation of the spectra is quicker, which is demonstrated in Figure 6.Solving a generalized eigenvalue problem with matrices that are strictly real is computationally cheaper compared when the matrices involved are complex.Second, the calculated eigenvalues come in only two flavours: real eigenvalues, or complex conjugate pairs.Their eigenvectors are similarly real, or come in complex conjugate pairs.
When we return the imaginary number i, the eigenvalues ω are expected to be strictly imaginary or come in pairs satisfying ω = −ω * .In the proceeding subsections, we calculate all eigenvalues and eigenfunctions using (90), and then multiplying the resulting spectra with i to retrieve the spectrum of (89).

Scalar perturbations
We now calculate the quasinormal modes of a scalar perturbation (s = 0) for l = 3.A simple Mathematica implementation is in Notebook 4.
The spectrum derived from using a basis tuple of {50, 50} (described in Section 6.1) is plotted on the complex plane in Re ω Im ω Fig. 7. Calculated spectrum of a scalar field in a Schwarzschild spacetime for l = 3 using the basis tuple {50,50} (described in Section 6.1), many of which are spurious.There are eigenvalues distributed along the negative-imaginary axis because of the existence of a continuum of eigenvalues that is present there.
Figure 7. Since the ODE eigenvalue problem is quadratic in ω, there are 102 eigenvalues as follows from the discussion in Section 4.2.

Filtering spurious modes.
In Section 6.2, we described two ways to filter out spurious eigenvalues: the CompareModes command and the CompareEigenfunctions command.In Section 7, the CompareModes command on a pair of spectra was sufficient to filter out all the spurious modes.
In the current case the CompareModes command at line 6 is not sufficient.Its output in Table 6 (a) includes purely imaginary modes, which are well-known not to exist for scalar perturbations given the boundary conditions we have chosen [45].
Recall that equation (77) comes from choosing a stationary ansatz for (76).It has been shown that the retarded Green function of this wave equation possesses a branch cut on the negative-imaginary axis [46,47].It is the 'shadow' of this continuum of eigenvalues which SpectralBP feels, as can be observed in Figure 7.
To filter these modes out, we demonstrate two solutions in the Notebook 4. These can be found in lines 8 and 11.
The first method is straightforward: calculate the spectrum of a third basis tuple and select eigenvalues common to all three spectra.We have chosen {100,100} as our third basis tuple, and the corresponding output is in Table 6 (b).The purely imaginary modes are successfully filtered out.
The second method would be to compare eigenfunctions between two basis tuples.This is the purpose of the CompareEigenfunctions command, whose output on line 8 is an empty set.This confirms that these modes are indeed spurious; their eigenfunctions are not approximately equal.One is then justified to filter out the purely imaginary modes in Table 6 (a).
The calculation of a third spectrum may be numerically prohibitive, especially when only a small subset of eigenvalues are suspected to be spurious.This consideration would favour one method over the other.In this case testing only the eigenfunctions of the suspected spurious eigenvalues, as filtered in line 10, should be favoured over the former method.
This second filter works because the rescaling in equation (82) keeps other valid solutions of our ODE eigenvalue problem non-analytic.In the case of the branch cut eigenvalues, their corresponding eigenfunctions remains singular at the cosmologcal horizon after rescaling [48].Thus, the approximation of these eigenfunctions in a Bernstein basis would fail to converge near the cosmological horizon.This idea is explored further in Section 8.2.2.
This failure to converge is shown explicitly in Figure 8, where we compare the eigenfunctions of the spurious eigenvalue −18.67i and the non-spurious eigenvalue ±1.3507 . . .−0.1930 . . .i.
Using a GetEigenfunctions command, we plotted the absolute difference between the eigenfunctions of approximately common eigenvalues for two spectral basis orders.The maximum error for the spurious eigenvalue is indicative of the presence of a singularity in the eigenfunction, Table 6.Result of a CompareModes command on 2 and 3 basis tuples (discussed in Section 6.1).(a) The filtered spectrum for the duo basis tuples include purely imaginary modes, which we know to be spurious.These modes may be filtered out using a CompareEigenfunctions command.(b) The filtered spectrum for the trio of basis tuples do not include purely imaginary modes.We have printed here significant digits shared by basis tuples {80,80} and {100,100} We echo an idea from [35].One must be careful in rescaling so that boundary conditions are still capable of the undesired solutions.For example, there are instances when peeling off an extra (r − 1) −1 term so that ϕ(r) ∼ (r − 1) is desirable.This boundary condition would fail to filter out the acausal solution at the black hole horizon, since both the acausal and causal solutions vanish at r = 1.The spectral method would then try to solve for solutions of the form, which generally is a mixture of causal and acausal parts at the black hole horizon.The ultimate consequence is that the boundary-value problem no longer has a discrete spectrum of eigenvalues.Continuing to calculate the spectrum using {50, 50} and {80, 80} would result in Figure 9.As expected, SpectralBP is unable to find the desired discrete spectrum.

Algebraically special modes
It is well-known that the standing wave equation for odd-and even-parity gravitational perturbations (s = 2) has an exact solution at what is called by Chandresekhar as the algebraically special mode.It is a purely imaginary frequency which appears to separate two different branches of the quasinormal mode spectrum: a lower branch that spirals towards the imaginary axis and an upper branch corresponding to an asymptotic high-damping regime.
It is a curious mode, whose frequencies can be shown analytically [49,50,51] to be and whose corresponding eigenfunctions, with singularities properly scaled out, can be expressed analytically as a truncated polynomial.For example, for l = 2, Various numerical investigations [23,24] are hard-pressed to converge towards this exact result.It has been argued [50] that the discrepancy can be traced to two explanations: (1) the algebraically special mode is sensitive to the exact form of the gravitational potential (affecting WKB and Pöschl-Teller potential fitting) and ( 2) the sensitivity of a method to a properly defined mode number (affecting the continued fraction methods by Leaver).
In fact, numerical methods that are able to find eigenvalues on the complex plane do not generally work when those  eigenvalues are located exactly on the imaginary axis.For example, the continued fraction method is not convergent for modes on the imaginary axis [24,52,53].This disputes previous analytic and numerical results concerning Kerr QNMs on the negative-imaginary axis.One can, however, deduce the existence of these modes by finding 'mode sequences' that arbitrarily get close to the negative-imaginary axis, including the special algebraic mode [54,52].How these modes move around the negative-imaginary axis is not accessible to Leaver's method.
With respect to this, spectral methods enjoy a significant advantage over Leaver's method: an algorithm such as SpectralBP is capable of finding eigenvalues on the imaginary axis.Unlike Leaver's method, which is based on a local power series expansion at one of the horizons, spectral methods find solutions globally.This has been reported before in [35], where the spectral algorithm QNMspectral finds a novel infinite set of purely imaginary modes for massless scalar perturbations in a Schwarzschild-de Sitter background.Because the spectral method is able to find these overdamped modes, one is able to observe complex bifurcation events in which quasinormal modes sink into, move along and emerge out of the negative imaginary axis where two QNMs collide.We have also used SpectralBP to uncover an interesting scenario that occurs in a Schwarzschild AdS background [27].

Algebraically special eigenvalues
We now solve (90) for s = 2 and for l = 2, 3, 4, 5, and reverse the transformation ω → iλ to retrieve the eigenvalues of (89).We have used basis tuples of {350,350} and {400,400} (described is Section 6.1) for all calculations, and we have filtered out spurious eigenvalues on the negative-imaginary axis using CompareEigenfunctions.The resulting spectra can be seen in Table 7 and Table 8.We show only the 10 lowest damping eigenvalues, using Mathematica's notation for significant digits for the purely imaginary eigenvalues.
Table 7. Gravitational perturbations with l = 2 and l = 3, calculated using basis tuples {350,350} and {400,400}.The special algebraic modes have 295 and 227 significant digits respectively.In units where the horizon is at r = 1, we have M = 1/2, so that ω2 = −4i and ω3 = −20i according to (93).Our numerical results show agreement up to 295 and 227 significant digits, respectively.The coincidence of the calculated numerically purely imaginary mode ω ′ l with the algebraically special mode ω l is very strong.The coincidence when calculating ω ′
As a additional check, we have verified that the eigenfunction solved by SpectralBP using a basis tuple of {400,400} and l = 2 is found consistent with (94) to within and error of 10 −250 .The eigenfunctions for l = 3, 4 are also truncated polynomials, of expected degrees 41 and 121 respectively.One might need the use of higher precision Table 8.Gravitational perturbations with l = 4 and l = 5, calculated using basis tuples {350,350} and {400,400}.The special algebraic modes have 137 and 115 significant digits respectively.In units where the horizon is at r = 1, we have M = 1/2, so that ω4 = −60i and ω5 = −140i according to (93).Our numerical results show agreement up to 137 and 115 significant digits, respectively.numbers to confirm that the degree of the l = 5 eigenfunction is of degree 281.

Schwarzschild
As we have described in Section 8.1, the eigenvalues of (90) are either purely real or come in complex conjugate pairs.As a consequence of this, when we transform back to ω from λ the calculated purely imaginary eigenvalues have exactly no real part.This avoids a criticism on numerical calculations which finds a single mode near the ASM with a finite real part whose symmetric pair ω = −ω * is unexpectedly not found.
The main lesson here is that SpectralBP manages exceptionally well to find eigenvalues on the negative-imaginary axis while filtering out spurious overdamped modes, as would other spectral or pseudospectral methods.This is in contrast with continued fraction methods, which cannot converge when the real part of the eigenvalue vanishes.
As a final note, and to illustrate the resources required to calculate one of the tables in this section, a single spectrum calculation for a basis tuple of {400, 400} takes around 1 hour each, running in a single 2.50 GHz Intel i5 Core with 8.00GB RAM.

Boundary behavior of the eigenfunctions
For completeness, we give warning when labelling solutions found by spectral methods as bonafide quasinormal modes whenever the eigenvalues calculated imply that the indicial equation (79) at one or more of the singularities are non-generic.This may affect whether or not the solution found satisfies the quasinormal mode boundary conditions.
For example, the finiteness of the eigenfunctions of the special algebraic modes at the boundaries can be folded back into (82), seemingly then implying that the quasinormal mode boundary conditions are satisfied and that these imaginary frequencies correspond to quasinormal modes.
The indicial equation ( 79) is said to be generic when its two solutions, ±iω, do not differ by an integer.This is manifestly true for general complex values of ω.In this case, the power series expansion at u = 1 of the rescaled function ϕ(r) in (82) converges, whether dominant or subdominant.At the algebraically special mode, however, the indicial equation is non-generic.From (93) and M = 1/2, the solution of the indicial equation are both integers, In this case, only one power series expansion of ϕ(r) is assured to converge, corresponding to the dominant solution.
For the subdominant, say φ(r), two things may happen.First, the subdominant solution may diverge logarithmically, of the form φ(u) ∼ c 0 ϕ(u) ln(u−1)+a 0 (u−1) 8 +a 1 (u−1) 9 +. . .(96) However, a miraculous cancellation may occur [55], in which case the logarithmic term vanishes.Thus, both solutions may be expressed as a power series expansion at r = 1.It is this latter case that occurs at the algebraically special mode for the Regge-Wheeler equation.This means that the dominant and subdominant solutions, corresponding to ingoing and outgoing modes at the black hole horizon respectively, may be rescaled to have the form, For the specific case of the ASM, the following two statements are then not mutually exclusive: (1) the ASM eigenfunction, properly rescaled, has a regular, well-behaved Frobenius expansion in powers of (u − 1) and (2) it is an inextricable mixture of the two linearly independent solutions at the black hole horizon, corresponding to a causal ingoing mode and an acausal outgoing mode.The reconciliation between the analytic and numerical results is thus simple but subtle; there is no contradiction.While SpectralBP has indeed found an eigenvalue-eigenfunction pair of the Regge-Wheeler equation, this solution is an inseperable mixture of both ingoing and outgoing solutions at the black hole horizon, and therefore is not a quasinormal mode.
In summary, SpectralBP picks up the special algebraic frequency to an incredible degree of accuracy, but because of the peculiar nature of the special algebraic mode, the corresponding eigenfunction is one that does not satisfy quasinormal mode boundary conditions, as would be expected from [55].

Conclusion
This work makes a case for the use of Bernstein polynomials as a basis for spectral and pseudospectral methods applied to ordinary differential eigenvalue problems.A prime example of these problems is the calculation of quasinormal modes in black hole spacetimes.The Bernstein polynomials constitute a non-orthogonal spectral basis, which may explain why they are much less utilized compared to Chebyshev or Fourier basis functions.In contrast to its more popular counterparts though, a Bernstein basis allows one to decouple some of the spectral weights relevant to boundary conditions of ordinary differential eigenvalue problems.More specifically, the weights for the first q and last the r basis polynomials for free without recourse to the differential equations.For some applications, this proves to be a significant advantage.
We developed a user-friendly Mathematica package, SpectralBP, as a general spectral solver for eigenvalue problems.This package fully utilizes the properties of Bernstein polynomials and several other algorithmic enhancements (such as a novel inverse iteration method) that we shall describe in a later paper.As far as we know, SpectralBP is unique among existing spectral codes in its use of a Bernstein basis.We described its key functionalities and showcased several examples for its use.In particular, to serve both as tutorial and benchmarks, we featured applications of SpectralBP to a number of model eigenvalue problems in quantum mechanics.Most importantly, we have also applied SpectralBP to quasinormal mode problems in the Schwarzshild geometry.In all of our example cases, SpectralBP succeeded in providing very accurate results.Remarkably, with only modest resources, we are able calculate the algebraically special modes of Schwarzschild gravitational perturbations.Purely imaginary modes are notoriously difficult to calculate with more conventional numerical methods even when the solution is straightforward to calculate analytically, as in the case for the Schwarzschild ASM.To the best of our knowledge, ours is the most accurate numerical calculation of these algebraically special modes in the extant literature, agreeing with the analytical prediction to a staggering (294!) number of significant digits.We have supplemented our calculations with a discussion on the subtleties of the boundary conditions of the algebraically special mode.Moving forward, spectral methods should be a very useful tool in finding quasinormal modes on the negative imaginary axis.
Encouraged by these successes, we believe that SpectralBP may serve as a useful tool for the black-hole physics community or just about anyone seeking to solve a differential eigenvalue problem.Future work will look into applications of SpectralBP to the Kerr spacetime, as well as several algorithmic enhancements (such as a novel inverse iteration method) that we shall describe in a later paper.We have also used SpectralBP to discover new interesting properties of the quasinormal modes of Schwarzschild-anti-de Sitter spacetime, which will also be discussed in a later paper.

These manipulations give us
We may now plug-in the following equally spaced and Chebyschev grids, The corresponding matrices simplify to A.2 From GEP to EP Compared to GEPs, the methods for solving eigenvalue problems of the standard form (EPs) are more diverse and more studied.Iterative algorithms to solve either the entire set of eigenvalues and eigenvectors or its subsets are widely available for a general class of complex-valued matrices.Critically, EPs are numerically cheaper to solve than GEPs.Consider the polynomial eigenvalue found in Section 4.2.If one of the matrices in the GEP is non-singular, then the GEP can be converted into an EP.This is apparently dependent on whether the lowest or highest matrix, M 0 and M m , in the matrix pencil (39) are invertible.
The corresponding eigenvalue problems follows, where (107) As for the full GEP that arises in Section 4.3, a similar analysis leads to complications.First, it can be shown that M′ is always singular.To show this, let us assume that there exists some M ′ j,k that is invertible.This is to say that, with respect to the matrix pencil from which M ′ j,k was constructed (M j,k,0 + ωM j,k,1 + ω 2 M j,k,2 + • • • + ω m M j,k,m )C k = 0 (108) the matrix M j,k,0 is invertible.To illustrate that M′ is always singular, we rearrange our simultaneous set of ODE's such that M ′ j,k is now indexed by M ′ 1,1 , and then we decompose where and and that the product of any two matrices of this form is also such a matrix.For M′ to be invertible, the matrix D ′ − C ′ A ′ −1 B ′ must not be singular.However, as we have shown, D ′ and C ′ A ′ −1 B ′ are both matrices whose sub-blocks are of the form given in (113).Thus, the matrix formed by their difference would be singular, as all of the identity matrices cancel out leaving all except n − 1 rows to vanish.
On the other hand, the inversion of M′′ is a rather involved calculation best left for computers.

A.3 Eigenfunction calculation -inverse iteration
In this section, we describe briefly the inverse iteration method implemented in SpectralBP to calculate the eigenvectors of a matrix pencil.It has the advantage of working on the matrix pencil directly without the need of linearizing the polynomial eigenvalue problem.For a problem involving n dependent functions, a polynomial degree of m, and N collocation points, the size of the matrices involved reduce from (nmN ) 2 to (nN ) 2 .
Suppose µ l is some eigenvalue numerically calculated from the GEP in (49).That is, for some eigenvalue ω l that exactly satisfies (49), The error ϵ is sourced from finite precision arithmetic, and should be very small.By definition, ω l and its corresponding eigenvector v l should also satisfy the polynomial eigenvalue problem without linearization in Section 4.2, where The inverse iteration algorithm is described in Notebook 5.
Notebook 5 Inverse iteration 1: Calculate A(µ) −1 2: v (0) = a vector with ||v (0) ||2= 1 ▷ initialize v (0) 3: for k = 1, 2, 3, . . ., kmax do 4: w = A(µ) −1 v (k−1) 5: Exit for loop 8: end if 9: end for 10: Return v (final) Its output can be shown to be of the form, The eigenvector v can then be split apart into the n eigenfunctions in the Bernstein basis.The algorithm here is part of a more general inverse iteration algorithm that is useful in the calculation of eigenvalue-eigenvector pairs in polynomial and transcendental eigenvalue problems, which will be the subject of a future work.
It is quite sufficient to calculate eigenfunctions at the same BP order the input eigenvalues were derived from.The error of the eigenfunctions is dominated by the use of a finite polynomial basis and not by finite precision arithmetic, as should be apparent in the examples discussed in Section 7 and Section 8.

A.4 Eigenfunction manipulations
Suppose we start with an eigenfunction of the form given in (28).In the interest of brevity, we denote the expanded Bernstein basis order as N max = N + q + r.From the linearity of the problem, eigenfunctions are determined up to a normalization constant.We may choose a normalization constant A so that function ψ(u), given by satisfies some desirable property.The simplest choice is to either set the coefficient of the leading polynomial expansion at either boundary to 1.
A −1 = C q , → lim u→a ψa (u) ≈ (u − a) q (120) (122) with the condition that n ≥ −2q and m ≥ −2r so that the integral remains finite.Using properties (12) and ( 13) of the Bernstein basis, the integral (122) can be evaluated to where A third way to normalize would then be, When w(u) = 1, the resulting function is normalized such that its L 2 -norm in the interval [a, b] is unity.The weight function may be utilized to calculate the L 2 -norm in another set of coordinates.This typically arises when the eigenfunctions are calculated in a compactified set of coordinates.
As an example, consider the the coordinate transformation in (61) and (62) in solving the harmonic oscillator.To normalize the eigenfunctions in the uncompactified coordinate system, their respective weights are of the form (126) One may calculate the square difference between two eigenfunctions in this way.Suppose two eigenfunctions ψ 1 (u) and ψ 2 (u) calculated from a spectral basis of order N 1 and N 2 respectively.
(127) Let us say that N 2 ≥ N 1 .We may expand the BP basis order of ψ 1 (u) using (11), (128) Thus, we may write the difference between the two eigenfunctions as a new sum of BPs of order N 2 , where One may then calculate the L 2 -norm of (129) using (123).With the Bernstein basis, it is also quite easy to rescale our function as in so that the resulting eigenfunction satisfies different asymptotics at the boundaries of the form with the condition that n ′ ≥ −q, m ′ ≥ −r.This is so that Ψ(u) may still expressed in the Bernstein basis.The resulting expression follows from the definition of Bernstein polynomials (7), that is where

Fig. 2 .
Fig. 2. The set of 11 Bernstein basis polynomials appropriate when q = 30 and r = 30, and their derivatives.The basis functions are localized around the center of [a, b], as are their derivatives.

Fig. 3 .
Fig. 3. Output of Notebooks 1. Top: (PrintFrequencies)The first 10 eigenenergies calculated using a basis degree of 50, plotted on the complex plane.Middle: (PrintEigenfunctions) The eigenfunctions of the first 3 eigenenergies, calculated using a basis degree of 50, normalized according to their L 2 -norm.Bottom: (PrintTable) Rescaled eigenvalues common to basis degrees of 50 and 80.There are 28 eigenenergies that share a minimum of 3 significant digits (not shown).We tabulate only the lowest 10.The spectrum calculated is in excellent agreement with(54).
increases the number of converged modes (up to a point).
has a domain of [−1, 1].As described in Section 6.1, one may change the default domain of [0, 1] to [−1, 1] by setting the option value of LowerBound to -1.The second,

2 Fig. 5 .
Fig.5.The calculated eigenfunctions ϕ BP n (x) in the uncompactified coordinate system are plotted above, while the absolute difference between ϕ BP n and the exact eigenfunctions ϕ exact n are plotted below.The eigenfunctions were calculated with a basis degree of 100.

Fig. 6 .
Fig.6.Benchmarking for performance using basis tuples {N, N }.The blue line comes from (89), in which the coefficient functions are complex.The orange line effects the replacement ω → iλ, solving (90) in which the coefficient functions are real.Both are power laws of the form T (N ) ∼ N 3.2 , with the latter performing faster.Calculations were done in a single 2.50 GHz Intel i5 Core with 8.00GB RAM

Fig. 8 .Fig. 9 .
Fig.8.The absolute difference between eigenfunctions of approximately equal eigenvalues using Bernstein basis orders 50 and 80. ϕ1(u) calculates the absolute difference for the eigenvalue ω = −18.67i,while ϕ2(u) calculates the absolute difference for the eigenvalue ω = ±1.3507 . . .−0.1930 . . .i.The former indicates that the eigenfunctions does not converge to some non-singular function, while the latter indicates convergence.
each sub-block in A ′ , B ′ , C ′ , D ′ is of the form A ′ , B ′ , C ′ , D ′ ∼

Table 2 .
Comparison between NSolve (a zero-finding algorithm in Mathematica 11.3) and SpectralBP, for eigenenergies in the range 0 ≤ E ≤ 1000.Eigenenergies with *'s were found by NSolve by sampling the range [0,1000] with a resolution of 0.01.Unmarked eigenenergies can be found by default.Eigenenergies found by SpectralBP used basis tuples of {61, 61} and {101, 101} (described in Section 6.1).There is excellent agreement between the eigenenergies found by both methods.