Spectral Quasinormal-Mode Solvers

Spectral collocation can compute many quasinormal frequencies at once, including cases where continued fractions are unavailable. It also produces many discretization eigenvalues at once, some of which have no continuum counterpart. The scientific problem is thus classification rather than diagonalization: identify modes that survive resolution, precision, basis, and boundary changes, then test their differential-equation residual independently. This page develops that numerical evidence without redefining the continuum resonance spectrum developed in Black-Hole Perturbation Theory.

Required background. Ordinary Differential Equations supplies boundary-value eigenproblems, conditioning, and residual analysis.

Helpful background. Formulations of General Relativity fixes the evolution variables, Black-Hole Perturbation Theory develops quasinormal boundary conditions, and Gravitational Radiation and Waveforms supplies the observable mode conventions.

Compactify the resonance problem

Factor the ingoing and outgoing asymptotics from the master field and map the exterior to $x\in[-1,1]$. For a Chebyshev–Lobatto grid $x_j=\cos(j\pi/N)$, replace derivatives by differentiation matrices and enforce the regular endpoint equations on the factored field. If the frequency enters polynomially, the result is

$$\left(A_0+\omega A_1+\omega^2A_2+\cdots\right)\mathbf{u}=0.$$

A quadratic problem may be linearized into a larger generalized pencil $\mathcal A\mathbf v=\omega\mathcal B\mathbf v$. Linearization preserves exact eigenvalues but changes numerical conditioning and introduces infinite eigenvalues when $\mathcal B$ is singular. The chosen factorization, map, grid, boundary rows, and pencil construction are therefore part of the method.

For coupled systems, scale variables so that one component does not dominate the matrix norm. Boundary conditions must replace, not merely supplement, the correct number of equations. Jansen 2017, §§2.2–2.4 and §3 gives an explicit collocation construction and discusses polynomial eigenvalue formulations suitable for black-hole spectra.

A concrete factored Schwarzschild operator

The continued-fraction construction defines $x=1-1/\widehat r$, $\widehat r=r/(2GM_{\rm phys})$, $\rho=-2iGM_{\rm phys}\omega$, and the factored remainder $y(x)$. In this subsection $x\in[0,1]$, with the horizon at zero. Substituting that factor into the axial gravitational equation gives

$$\begin{aligned} 0={}&x(1-x)^2y'' +\bigl[(1+2\rho)-(4+8\rho)x+(3+4\rho)x^2\bigr]y'\\ &+\bigl[3-\ell(\ell+1)-8\rho^2-4\rho +(4\rho^2+4\rho-3)x\bigr]y. \end{aligned}$$

Derivatives here are with respect to $x$. For a candidate remainder with finite derivatives, the limiting equations are

$$\begin{aligned} x=0:\quad &(1+2\rho)y'(0) +[3-\ell(\ell+1)-8\rho^2-4\rho]y(0)=0,\\ x=1:\quad &-2\rho y'(1)-[\ell(\ell+1)+4\rho^2]y(1)=0. \end{aligned}$$

These limits provide endpoint rows without evaluating divergent exponential factors. They also give a direct algebra check: inserting $y=\sum a_kx^k$ recovers the recurrence in Leaver 1985, § 2, Eqs. 3–8, p. 287. Separate the constant, linear, and quadratic powers of $\rho$ before assembling the pencil, and translate its eigenvalue back to $GM_{\rm phys}\omega=i\rho/2$.

Endpoint regularity is necessary but must not silently replace the outgoing branch. At asymptotically flat infinity, the unwanted incoming solution after this factorization contains $e^{-2i\omega r_*}$; for damped frequencies it can be exponentially small and smooth after compactification. A finite endpoint value or a small endpoint-row residual does not by itself exclude it. The retarded branch is fixed by analytic continuation from the frequency half-plane where the outgoing solution is defined. For this benchmark, the independently derived minimal-solution continued fraction supplies a particularly useful branch check. The zero-frequency branch point also requires separate treatment.

Physical modes versus spectral pollution

At a single $N$, plotting all eigenvalues is not evidence. Match candidates between $N$ and $2N$ by complex distance and eigenfunction overlap, then repeat at higher precision. A physical candidate should form a stable sequence while the cloud of spurious roots drifts. If $L(\omega)=\sum_{q=0}^{d}\omega^qL_q$, evaluate every $L_q u$ on a denser check grid that was not used in the solve. With positive quadrature weights $w_j$, define $\lVert z\rVert_w^2=\sum_jw_j\lvert z_j\rvert^2$ and report the dimensionless backward residual

$$\rho_w= \frac{\lVert L(\omega)u\rVert_w} {\displaystyle\sum_{q=0}^{d}\lvert\omega\rvert^q \lVert L_q u\rVert_w},$$

and give separately normalized residuals for every endpoint equation. The denominator must be nonzero; a zero computed denominator signals a degenerate normalization rather than a successful solve. This weighting and operator scaling make comparisons across check grids and field normalizations meaningful. They detect an eigenpair that solves the matrix pencil but poorly represents the differential equation.

Nonselfadjoint operators can have large eigenvalue condition numbers: a tiny coefficient perturbation may move a mode substantially even when $\rho_w$ is small. Report left and right eigenvector overlap or a condition estimate, and probe the smallest singular value of $L(\omega)$ around important roots. Trefethen and Embree 2005, §§1–2, pp. 1–19 and §§32–34, pp. 287–328 explain why residuals alone do not control eigenvalue error for nonnormal operators. Jansen 2017, §§4.1–4.3 and Appendix A, Figs. 11–12 demonstrates resolution tracking and precision sensitivity; it does not license accepting every stable-looking point.

A discriminating benchmark

Let $M_{\rm phys}$ denote the physical black-hole mass and $M_{\rm geom}=G M_{\rm phys}$ its geometric length. Use the Schwarzschild axial gravitational $(\ell,n)=(2,0)$ mode,

$$G M_{\rm phys}\omega=M_{\rm geom}\omega \simeq0.3736715-0.0889625i,$$

with $e^{-i\omega t}$. The displayed precision is the conversion of Leaver 1985, Table 1, p. 293, which uses $2M_{\rm geom}=1$. Run $N$, $2N$, and a third resolution, each in ordinary and higher precision. Record the matched eigenvalue, off-grid residual, endpoint residuals, condition estimate, and error against the reference. Independently compute the same root with a continued fraction or two-sided shooting Wronskian. The independent method should share only the continuum equation and convention package, not the spectral matrices or root-polishing code.

Vary the compactification and asymptotic factorization. True eigenvalues remain within the uncertainty; spurious roots often follow the changed grid distribution. For a smooth factored eigenfunction, exponential convergence may precede a roundoff plateau. If the background coefficients or eigenfunction have limited differentiability, algebraic convergence is expected and should not be disguised by fitting too short a range.

The benchmark establishes the method near one low overtone. Highly damped modes, near-extremal branches, and coupled fields can be much more ill-conditioned. A solver must restrict its claim to the tested region rather than infer uniform spectral accuracy.

Environment and rerun interface

Pin the symbolic coefficient source, angular convention, map, basis, matrix assembly, eigensolver and linear-algebra library, balancing options, precision backend, matching algorithm, and code revision. A suitable interface is qnm-spectral --config schwarzschild-l2.toml --output run/. Machine-readable outputs should contain the full unfiltered spectrum, eigenvectors or checksums, candidate matches, residual maps, precision and resolution sweeps, and environment metadata. This is a specification for reproduction, not a statement that a supplied program has been executed.

Reject candidates that move under basis or precision changes, have large off-grid residuals, appear only after an arbitrary matching tolerance, or fail the independent continued-fraction benchmark. Also reject a claimed digit when the condition estimate predicts greater sensitivity than that digit. Matrix balancing can improve the computation, but results must be transformed back and checked against the original unbalanced operator.

Diagnostic

A candidate changes by only $10^{-9}$ between $N=80$ and $N=100$, but its off-grid differential residual is $10^{-3}$ and it moves by $10^{-2}$ when precision rises from 53 to 100 bits. Is resolution stability sufficient?

Solution

No. The large independent residual and precision movement reject the claimed accuracy. They are compatible with roundoff amplification, inaccurate matrix assembly, or an unresolved mode; they do not uniquely diagnose a pseudospectral artifact. Reassemble at higher precision, track the candidate across at least three resolutions and an alternative map, and require agreement with an independent continued-fraction or shooting calculation.

References

• Jansen, A. (2017). “Overdamped modes in Schwarzschild–de Sitter and a Mathematica package for the numerical computation of quasinormal modes.” European Physical Journal Plus 132, 546. arXiv.
• Leaver, E. W. (1985). “An analytic representation for the quasi-normal modes of Kerr black holes.” Proceedings of the Royal Society A 402, 285–298. DOI.
• Trefethen, L. N., and Embree, M. (2005). Spectra and Pseudospectra: The Behavior of Nonnormal Matrices and Operators. Princeton University Press. Publisher.