Spectral Collocation and Direct Eigenvalue Methods

Spectral collocation replaces a differential equation by equations for a function’s values at finitely many points. For a black-hole resonance, first remove the known ingoing and outgoing factors; a polynomial can then approximate the remaining radial profile on a compact interval. A generalized eigensolver produces candidate frequencies, but a small matrix residual does not establish a quasinormal mode. One must also check the original equation, the radiation prescription, and convergence across genuinely different discretizations.

This page develops a finite, rerunnable benchmark: an exactly solvable Pöschl–Teller barrier, the fundamental axial Schwarzschild mode, and a Kerr radial calculation with a supplied angular eigenvalue. The last is not a coupled angular–radial solve. The tests support these selected calculations, not completeness of the computed spectrum or a general theorem excluding spectral pollution.

Required background. PDEs and Green Functions for Gravitation supplies differential operators, singular endpoints, and boundary-value problems.

Helpful background. Black-Hole Geometry and Mechanics supplies horizon geometry; Formulations of General Relativity supplies perturbative variables; Computational Relativity supplies numerical conditioning; Black-Hole Perturbation Theory supplies the master-equation context.

The compactified Schwarzschild equation

Let $\mu=GM$, where $M$ is physical mass. Use dimensionless radius $\rho=r/\mu$, tortoise coordinate $\rho_*=r_*/\mu$, and frequency $w=\mu\omega$. Thus every quoted black-hole frequency below is $GM\omega$; its damping time is $\mu/(-\operatorname{Im}w)$. For source-free axial gravitational perturbations of four-dimensional Schwarzschild, $\ell\geq2$,

$$\frac{d^2X}{d\rho_*^2} +[w^2-\widehat V_\ell(\rho)]X=0, \qquad \widehat V_\ell =f\left[\frac{\ell(\ell+1)}{\rho^2}-\frac6{\rho^3}\right], \quad f=1-\frac2\rho .$$

Here $\widehat V_\ell=\mu^2V_\ell$ and $d\rho_*/d\rho=f^{-1}$. The Regge–Wheeler equation derives this master problem. Its retarded resonance conditions are the analytically continued branches $X\sim e^{-iw\rho_*}$ at the future horizon and $X\sim e^{+iw\rho_*}$ at infinity, as explained in Quasinormal Boundary Conditions and Resonances.

Compactify and factor:

$$x=1-\frac2\rho,\qquad X=F(\rho,w)y(x),\qquad F=e^{iw\rho}(\rho-2)^{-2iw}\rho^{4iw}.$$

The physical exterior is $0\lt x\lt1$: the horizon is $x=0$, infinity $x=1$. All powers use the real logarithm of their positive, dimensionless bases on this interval. At the horizon $F$ has the ingoing power; at infinity it has the outgoing exponential and power. Overall constant factors have no effect on the homogeneous eigenproblem.

To reproduce the substitution, use

$$\frac{dx}{d\rho}=\frac2{\rho^2},\qquad \frac{d^2x}{d\rho^2}=-\frac4{\rho^3},\qquad \frac{F_\rho}{F} =iw-\frac{2iw}{\rho-2}+\frac{4iw}{\rho}.$$

Insert $X_\rho/F=y_x\,x_\rho+(F_\rho/F)y$ and differentiate once more before simplifying. Divide the original equation by $F x(1-x)^2$ in the open interval. With $L=\ell(\ell+1)$ and primes now denoting $x$ derivatives, the resulting equation is

$$\begin{aligned} 0={}&\frac{x(x-1)^2}{4}y'' +\left[\frac{(x-1)(3x-1)}4 +iw(-2x^2+4x-1)\right]y'\\ &+\left[-\frac{L+3x-3}{4} -2iw(x-1)-4w^2(x-2)\right]y . \end{aligned}$$

The leading coefficient vanishes at both ends. For a profile with finite first and second derivatives, its endpoint limits give

$$\begin{aligned} 0&=\left(\frac14-iw\right)y'(0) +\left(\frac{3-L}{4}+2iw+8w^2\right)y(0),\\ 0&=iw\,y'(1)+\left(4w^2-\frac L4\right)y(1). \end{aligned}$$

These are the endpoint rows of the discretization. They relate a value to a derivative; they do not set the value to zero. Appending Dirichlet conditions would change or overdetermine the problem.

Smoothness and the radiation prescription

The factor has a geometric interpretation: changing from stationary time to a horizon-penetrating, asymptotically null time coordinate rescales the frequency-domain field. The compactified equation then supplies limiting endpoint relations. See Macedo–Zenginoğlu 2025, §3.2, pp. 3–4.

However, boundedness or a leading logarithmic-derivative limit is not a unique outgoing prescription. For a damped frequency, an incoming solution can be exponentially smaller than the outgoing one at spatial infinity. After the outgoing factor is divided out, the incoming contribution behaves schematically as

$$y_{\rm in}\sim e^{-2iw\rho}\rho^{-4iw}, \qquad \rho=\frac2{1-x}.$$

When $\operatorname{Im}w\lt0$, its exponential decays faster than every power of $1-x$. Such a contribution can be flat at the endpoint: it and all its derivatives vanish there. In the elementary compact-support half-line model, this obstruction is exactly $e^{-2i\lambda/z}$ at $z=0$; it is smooth for damped $\lambda$ without being the outgoing branch. Galkowski–Zworski 2021, arXiv v1, pp. 1–2, Eq. (1.3) gives this precise example. It is not a substitute for a function-space theorem for the full Kerr problem.

Accordingly, this calculation tests candidates against the radiation branches selected by an independent continued-fraction calculation. Endpoint relations, finite-distance phase tests, and coefficient decay are useful diagnostics, not a proof that every smooth profile is a resonance. The distinction between bounded solutions and QNMs is also discussed in Macedo–Zenginoğlu 2025, §3.3, pp. 4–5.

Chebyshev approximation

Choose $N+1$ Chebyshev–Lobatto points

$$z_j=\cos\frac{j\pi}{N},\qquad u_j=\frac{z_j+1}{2},\qquad j=0,\ldots,N.$$

The indexing runs from $u_0=1$ to $u_N=0$. Let $p_N(u)$ be the degree-at-most-$N$ polynomial through the nodal values $v_j$. If $\ell_k(u)$ is the cardinal polynomial satisfying $\ell_k(u_j)=\delta_{jk}$, then

$$p_N(u)=\sum_{k=0}^N v_k\ell_k(u), \qquad (D_u)_{jk}=\ell_k'(u_j).$$

This explains why a differentiation matrix is dense: a global polynomial’s derivative at one node depends on all the nodal values. It also explains the need to test derivatives independently, for example by differentiating known low-degree polynomials.

Vary the distribution of physical nodes with a positive map parameter $c$:

$$x=\frac{u}{c+(1-c)u},\qquad x_u=\frac{c}{[c+(1-c)u]^2},\qquad x_{uu}=-\frac{2c(1-c)}{[c+(1-c)u]^3}.$$

The physical endpoints are unchanged. At $c=1$ the map is the identity; at other values the collocation space is made of functions $p_N(u(x))$, generally not polynomials in $x$. The correct nodal derivatives are

$$\begin{aligned} D_x&=\operatorname{diag}(x_u^{-1})D_u,\\ D_{xx}&=\operatorname{diag}(x_u^{-2})D_u^2 -\operatorname{diag}(x_{uu}x_u^{-3})D_u . \end{aligned}$$

The second term is essential. Dropping it discretizes a different differential operator.

For every off-grid or endpoint diagnostic, evaluate the same polynomial in its original coordinate:

$$\begin{aligned} u(x)&=\frac{cx}{1-(1-c)x},\\ u_x&=\frac{c}{[1-(1-c)x]^2},\qquad u_{xx}=\frac{2c(1-c)}{[1-(1-c)x]^3},\\ y(x)&=p_N(u(x)),\\ y_x&=p_N'(u(x))u_x,\\ y_{xx}&=p_N''(u(x))u_x^2+p_N'(u(x))u_{xx}. \end{aligned}$$

Interpolating the same nodal values as a new polynomial in $x$ would change the function being tested. A matrix could pass while this new function fails the differential equation—or the wrong function could conceal an error. The implementation uses barycentric evaluation with explicit Lobatto weights, $b_j=(-1)^j$ except for half-weight endpoints. Supplying those weights also avoids recomputing them with a randomized product order; the versioned API is documented in SciPy 1.13.1, BarycentricInterpolator, parameters and methods.

The squared comparison norm used below is the quadrature approximation to $\int_0^1 |y|^2dx$ for the factored profiles, on a common physical-$x$ grid. It is an auxiliary numerical diagnostic, not a black-hole energy norm or the square-integrability norm of the exponentially growing stationary QNM field.

The quadratic pencil and its finite eigenvalues

Each coefficient in the Schwarzschild equation is polynomial in $w$. Collocating it, including its limiting endpoint rows, gives

$$P_N(w)v=(A_0+wA_1+w^2A_2)v=0.$$

Introduce an auxiliary vector $q=wv$. Then

$$\begin{pmatrix}0&I\\-A_0&-A_1\end{pmatrix} \begin{pmatrix}v\\q\end{pmatrix} =w \begin{pmatrix}I&0\\0&A_2\end{pmatrix} \begin{pmatrix}v\\q\end{pmatrix}.$$

The first block row imposes $q=wv$; substituting it in the second recovers the original quadratic equation. This doubles the matrix dimension but allows a generalized eigensolver to find many candidates in one call. One must not invert $A_2$ without proving it invertible. A generalized pencil can contain infinite eigenvalues; the finite-frequency window does not include them. SciPy 1.13.1, linalg.eig, parameters and returns specifies the generalized problem and its homogeneous-eigenvalue representation.

A small normalized residual

$$r_{\rm pencil} =\frac{\lVert P_N(w)v\rVert_2} {(\lVert A_0\rVert_{\rm F}+|w|\lVert A_1\rVert_{\rm F} +|w|^2\lVert A_2\rVert_{\rm F})\lVert v\rVert_2}$$

checks the finite matrix calculation. The subscript $\mathrm F$ denotes the Frobenius matrix norm. It does not bound the distance to a continuum resonance, especially for an ill-conditioned or nonnormal problem. That issue leads to Pseudospectra and Nonnormal Sensitivity.

A nonconstant exactly solvable test

Before testing a black hole, use a barrier whose resonance frequencies and some factored profiles are known exactly:

$$\frac{d^2X}{d\xi^2} +[\Omega^2-V_0\operatorname{sech}^2(\alpha\xi)]X=0, \qquad x=\frac{1+\tanh(\alpha\xi)}2.$$

Here $\xi$ and $\Omega$ are measured in an arbitrary chosen length unit and its inverse; the benchmark takes $V_0=1$, $\alpha=0.2$. For $V_0\gt\alpha^2/4$, the positive-real-frequency outgoing family is

$$\Omega_n=\sqrt{V_0-\frac{\alpha^2}{4}} -i\alpha\left(n+\frac12\right), \qquad n=0,1,\ldots .$$

The two ends are radiative: $X\sim e^{-i\Omega\xi}$ on the left and $X\sim e^{+i\Omega\xi}$ on the right. The exactly soluble barrier is also used in WKB and Geometric-Optics Spectra to test approximation coefficients; here it tests the full differentiation operator.

Since $d/d\xi=2\alpha x(1-x)d/dx$, put

$$X=[x(1-x)]^{-i\Omega/(2\alpha)}y.$$

Divide the original equation by $4x(1-x)$ times this factor. Direct differentiation gives

$$\begin{aligned} 0={}&\alpha^2x(1-x)y''\\ &+\left[\alpha^2(1-2x) +i\alpha\Omega(2x-1)\right]y'\\ &+(\Omega^2+i\alpha\Omega-V_0)y. \end{aligned}$$

For a polynomial profile of degree $n$, its highest-power coefficient vanishes only if

$$\Omega^2+i\alpha(2n+1)\Omega-V_0-\alpha^2n(n+1)=0.$$

Completing the square recovers the displayed frequency family. The lower powers determine the polynomial coefficients; the first two profiles can be checked directly without a special-function library.

For $n=0$, $y_0=1$ solves the equation. This is a useful frequency test but a poor derivative test: both derivatives vanish. Reversing the real first-derivative coefficient would leave this fundamental test unchanged.

The first overtone supplies the missing check:

$$y_1=2x-1,\qquad \Omega_1=\sqrt{V_0-\alpha^2/4}-\frac{3i\alpha}{2}.$$

Substituting it in the correct equation gives zero. Replacing $\alpha^2(1-2x)$ by $\alpha^2(2x-1)$ instead produces $4\alpha^2(2x-1)$, including nonzero endpoint residuals. Thus the exact constant and nonconstant profiles test different parts of the operator. The executable checks both, including the reconstructed unfactored equation; it also checks mapped values and derivatives of a known polynomial in $u$.

Kerr with a fixed angular eigenvalue

For Kerr, let $a_*=a/\mu=J/(GM^2)$, $|a_*|\lt1$, and use the separated azimuthal factor $e^{im\phi}$. Define

$$\begin{aligned} \rho_\pm&=1\pm\sqrt{1-a_*^2},& d&=\rho_+-\rho_-,\\ \Delta_\rho&=(\rho-\rho_+)(\rho-\rho_-),& K&=(\rho^2+a_*^2)w-ma_* . \end{aligned}$$

The source-free radial Teukolsky equation, in the angular convention used on Angular Teukolsky Equation and Separation Constants, is

$$\Delta_\rho\mathcal R'' +2(s+1)(\rho-1)\mathcal R' +\left[ \frac{K^2-2is(\rho-1)K}{\Delta_\rho} +4isw\rho-\lambda \right]\mathcal R=0,$$

where the primes in this equation are $\rho$ derivatives and

$$\lambda=A+a_*^2w^2-2ma_*w.$$

Freeze $A$ only for this radial benchmark. In the coupled physical problem it is $A_{s\ell m}(a_*w)$, so the radial frequency dependence is no longer a quadratic polynomial with fixed coefficients. An angular–radial iteration must also track the angular branch.

Use

$$\begin{aligned} x&=\frac{\rho-\rho_+}{\rho-\rho_-},& \sigma_+&=\frac{2\rho_+w-ma_*}{d},\\ \mathcal R&=e^{iw\rho}(\rho-\rho_+)^{-s-i\sigma_+} (\rho-\rho_-)^{-1-s+2iw+i\sigma_+}y(x). \end{aligned}$$

Here the dimensionless Kerr tortoise coordinate satisfies

$$\frac{d\rho_*}{d\rho}=\frac{\rho^2+a_*^2}{\Delta_\rho}.$$

The factor removes the leading radiation behaviors

$$\begin{aligned} \mathcal R&\sim\Delta_\rho^{-s}e^{-i(w-m\mu\Omega_H)\rho_*} &&\text{at the horizon},\\ \mathcal R&\sim\rho^{-1-2s}e^{iw\rho_*} &&\text{at infinity}, \end{aligned}$$

with $\mu\Omega_H=a_*/(2\rho_+)$. The raw equation and factor are the starting point for the independent substitution, not coefficients inferred from a desired eigenvalue. See Leaver 1985, Eqs. (16)–(17) and (22)–(23), pp. 289–291; his geometrized parameters are converted here to $\rho,a_*,w$.

For the implemented $s=-2,m=1$ case, the transformed equation is

$$P y''+(B_0+wB_1)y'+(C_0+wC_1+w^2C_2)y=0,$$

where primes again mean $x$ derivatives and

$$\begin{aligned} P&=x(1-x)^2,\\ B_0&=(1-x)\left[3-x+\frac{2ia_*(1-x)}d\right],\\ B_1&=2i\left[dx-3x^2+4x-1-\frac{2(1-x)^2}{d}\right],\\ C_0&=(1-x)\left(1-\frac{2ia_*}{d}\right)-A,\\ C_1&=-2a_*+\frac{8a_*(x-1)}d +\frac{i}{d}\left[(d+2)^2+2(d-2)x\right],\\ C_2&=\frac{d+2}{4d}\left[d^2+14d+32(1-x)\right]. \end{aligned}$$

At each endpoint retain $(B_0+wB_1)y'+(C_0+wC_1+w^2C_2)y=0$. The degeneration as $d\to0$ signals that this fixed-subextremal representation is not an extremal-limit algorithm.

Convergence, matching, and independent diagnostics

The finite sweep uses all three cases, eight orders $N=12,20,24,28,32,36,40,44$, and all three maps $c=0.7,1,1.4$: 72 configurations. Every finite candidate in $0\lt\operatorname{Re}w\lt1$, $-1\lt\operatorname{Im}w\lt0$ is retained, using $\Omega$ in place of $w$ for the barrier.

Follow branches before assigning reference labels

A nearby frequency from any other calculation is not enough. At each fixed map, compare consecutive members of the declared $N$ sequence. At each fixed $N$, compare both alternate maps. A matched pair must be nearest in both directions, have frequency difference below $5\times10^{-6}$, and have a second-nearest distance exceeding the nearest distance by more than $5\times10^{-6}$ in both directions. An absent second candidate has no competing neighbor; an empty comparison configuration fails. This explicit separation rule rejects ambiguity; it is deliberately conservative for nearby or crossing branches.

Compare the factored eigenfunctions on a common 64-point Gauss–Legendre grid in physical $x$. With positive quadrature weights defining $\langle\cdot,\cdot\rangle_Q$, use

$$\eta(v,z) =1-\frac{|\langle y_v,y_z\rangle_Q|} {\lVert y_v\rVert_Q\lVert y_z\rVert_Q}.$$

This removes arbitrary complex normalization and phase. The connected set of accepted pairwise matches must contain at most one candidate from each configuration; otherwise branch identity is unresolved. Require $\eta\lt10^{-3}$ for each pair. An interior resolution needs both adjacent-resolution matches; the first and last resolutions have only one available neighbor. Each candidate still needs both other maps.

Only after constructing these branches does the benchmark compare them with independent reference frequencies. An unambiguous reference match at the highest order, $N=44$, and identity map, $c=1$, labels its whole connected branch. Other configurations are not selected independently by nearness to the reference. A rejected row means that this finite test did not validate it. It does not mean that an unlisted overtone, unresolved branch, or resonance outside the window has been proved nonexistent.

Test the function, not only the matrix

For each candidate, reconstruct $p_N(u(x))$ and its derivatives consistently. The saved diagnostics include:

Test	What it checks	Limitation
Quadratic-pencil residual	Accuracy of the finite generalized eigensolve	Does not test the continuum operator
Transformed off-grid residual	Differential equation at 37 points from $x=0.05$ to $0.95$	Can share an incorrect coefficient with the matrix
Unfactored off-grid residual	Original physical ODE after reconstructing the removed factor	Finite sampling is not a uniform residual bound
Limiting endpoint rows	The two analytically derived value–derivative relations	Do not uniquely impose analytic radiation selection
Horizon and infinity phase tests	Correct signs and powers near $x=10^{-7}$ and $1-10^{-7}$	Subdominant unwanted solutions can evade a leading phase limit
Last four Chebyshev coefficients	Resolution of the chosen factored profile in $u$	A finite small tail is not a regularity theorem
Resolution, map, and reference comparisons	Persistence of a selected frequency and profile	Do not establish completeness or accuracy of all digits

If $T_k(x_j)$ are the sampled contributions to an equation, the off-grid residual is

$$r_{\rm ODE} =\frac{\max_j\left|\sum_k T_k(x_j)\right|} {\max_j\sum_k|T_k(x_j)|}.$$

For the transformed equation, the three contributions are grouped by powers of frequency; for the raw equation they are the original ODE terms, after cancelling the common factor analytically. This avoids overflow from the damped mode’s large stationary profile. A zero global denominator is undefined, not an automatic pass. Normalizing over the sampled grid also avoids a pointwise $0/0$ at an eigenfunction node, such as the center of the first Pöschl–Teller overtone.

Endpoint residuals divide by $(1+|w|+|w|^2)(|y|+|y'|)$. The coefficient-tail ratio is the largest magnitude among the last four coefficients divided by the largest coefficient magnitude. Phase diagnostics compare the reconstructed logarithmic derivatives to the appropriate ingoing and outgoing phases; Kerr’s known spin-dependent powers are removed first.

The local thresholds are $10^{-10}$ for the pencil residual and $10^{-4}$ for each transformed, unfactored, endpoint, phase, and coefficient-tail diagnostic. Cross-resolution and cross-map thresholds are those above. Reference-labeled benchmark rows additionally require frequency disagreement below $5\times10^{-6}$. These thresholds are reproducibility criteria for this specified experiment; a measured disagreement below a threshold is not an error bar of that size.

Increasing $N$ indefinitely at fixed floating-point precision is not a cure. Derivative matrices amplify roundoff, and a resonance can be sensitive to perturbations of its discretization. Look for a jointly stable frequency, profile, residual, and map comparison, then repeat with a genuinely independent representation. Frobenius Series and Continued Fractions supplies the independent radial method, while the present polynomial representation supplies a different approximation.

Benchmark inputs and rerun

The Schwarzschild case is axial $\ell=2,n=0$. The Kerr radial case uses $a_*=0.4$, $(s,\ell,m)=(-2,2,1)$ and the fixed input

$$A=3.767574385986149+0.053240584492497i.$$

The reference targets are

Case	Dimensionless target frequency	Role
Pöschl–Teller, $n=0$	$0.994987437106620-0.100000000000000i$	Exact constant-profile control
Pöschl–Teller, $n=1$	$0.994987437106620-0.300000000000000i$	Exact nonconstant operator control
Axial Schwarzschild, $\ell=2,n=0$	$0.37367168441804184-0.08896231568893570i$	Independent continued-fraction comparison
Kerr radial, fixed $A$ above	$0.40797910871532318-0.08725714473585831i$	Fixed-input radial comparison only

The finite precision of the supplied $A$ and reference constants limits the significance of smaller disagreements. Independent continued fractions were evaluated at this same fixed $A$, not by recomputing an angular eigenvalue at every trial frequency. Depth and working precision were varied; the displayed roots are rounded to double-precision scale. The code’s parameter name omega denotes these dimensionless frequencies; it does not restore physical units automatically.

What the complete sweep shows

The saved run retains 536 candidates and 64 accepted fundamental-reference rows: 24 for Pöschl–Teller, 20 for Schwarzschild, and 20 for fixed-$A$ Kerr. The nonconstant Pöschl–Teller overtone is an additional exact operator control, not a fourth reference branch in the accepted table.

The following entries are the largest frequency disagreements across all three maps, rounded to three significant figures. A dash means that not all three reference rows passed at that order; it is not a zero error. The last column counts accepted rows out of the nine case–map combinations.

$N$	Pöschl–Teller	Schwarzschild	Fixed-$A$ Kerr	Accepted
12	$1.87\times10^{-15}$	—	—	3/9
20	$1.36\times10^{-14}$	—	—	7/9
24	$3.96\times10^{-14}$	$4.52\times10^{-9}$	$3.46\times10^{-9}$	9/9
28	$3.05\times10^{-14}$	$5.39\times10^{-10}$	$8.61\times10^{-10}$	9/9
32	$7.68\times10^{-14}$	$7.05\times10^{-11}$	$1.74\times10^{-9}$	9/9
36	$9.86\times10^{-14}$	$1.03\times10^{-11}$	$2.19\times10^{-9}$	9/9
40	$8.54\times10^{-14}$	$1.83\times10^{-12}$	$1.96\times10^{-8}$	9/9
44	$1.85\times10^{-13}$	$4.30\times10^{-13}$	$4.35\times10^{-8}$	9/9

For each black-hole case, the poorly resolved $N=12,c=1.4$ candidate prevents acceptance of the other two maps at that order. The $N=20,c=1.4$ candidate also fails its required comparison with $N=12$. Keeping those failures is part of the experiment: selecting only the best-looking map would hide them.

Schwarzschild improves throughout the displayed range, whereas Kerr’s frequency disagreement worsens after about $N=28$ even though its raw sampled residual at $N=44$ is at most $6.30\times10^{-13}$. That contrast illustrates why a small differential residual is not a frequency-error bound. The Pöschl–Teller fundamental is already an exact constant profile, so its small fluctuations mainly test finite arithmetic rather than approximation convergence.

No ambiguous nearest pairs or duplicate-configuration branch components occurred in this sweep. Synthetic frequency inputs test ambiguity and reciprocity; deliberate code and recorded-data mutations protect the branch-uniqueness check. The ordinary output does not demonstrate a pathology merely because the software can detect it.

Download and reproduce

The standalone downloads are the canonical Python driver, pinned requirements, candidate table, accepted-residual table, numerical controls and matching details, and rerun guide. Download them into a fresh directory. No private checkout or unlisted input is required. The reference environment is CPython 3.9.6, NumPy 2.0.2, SciPy 1.13.1, with complex double precision.

With that interpreter available, run the following in the download directory. The output names keep a rerun separate from the downloaded references.

python3 --version
python3 -m venv .venv-spectral
. .venv-spectral/bin/activate
python -m pip install --requirement requirements.txt
python spectral_qnm.py \
  --cases poschl-teller,schwarzschild,kerr-radial \
  --orders 12,20,24,28,32,36,40,44 \
  --maps 0.7,1.0,1.4 \
  --write-candidates rerun-candidates.csv \
  --write-residuals rerun-residuals.csv \
  --write-checks rerun-checks.json

Exercises

Recover the horizon row. Evaluate the compactified Schwarzschild equation at $x=0$ for $\ell=2$. Why does prescribing $y(0)=0$ not implement the ingoing condition?

Solution

The equation becomes

$$\left(\frac14-iw\right)y'(0) +\left(-\frac34+2iw+8w^2\right)y(0)=0.$$

For generic $w\ne-i/4$, it fixes $y'(0)/y(0)$ when the amplitude is nonzero. The ingoing factor has already been removed. Setting $y(0)=0$ usually also forces $y'(0)=0$, eliminating the desired normalized solution instead of selecting it. The exceptional frequency at which the derivative coefficient vanishes requires separate local analysis; dividing by that coefficient there is invalid.

Expose a mapped-interpolation error. Take $N=1$, $c=1.4=7/5$, and $p_1(u)=u$. At $x=1/2$, compare the correct $y(x)=p_1(u(x))$ and $y_x$ with the degree-one polynomial interpolated directly through the physical endpoint values.

Solution

The endpoint data are $y(0)=0,y(1)=1$ in either description. The correct mapped function is $u(x)=7x/(5+2x)$, giving $y(1/2)=7/12$ and $y_x(1/2)=35/36$. Interpolating directly in $x$ instead gives $y=x$, hence $y(1/2)=1/2$ and $y_x=1$. Agreement at every collocation node has not made them the same function. Off-grid tests must use the coordinate in which the approximation was defined.

Design a test that the fundamental cannot pass accidentally. Substitute $y_1=2x-1$ into the Pöschl–Teller equation and recover the $n=1$ frequency. Then reverse the real first-derivative coefficient. What residual remains?

Solution

Since $y_1'=2$ and $y_1''=0$, the correct equation reduces to

$$\left(\Omega^2+3i\alpha\Omega-V_0-2\alpha^2\right)(2x-1)=0.$$

Completing the square gives $(\Omega+3i\alpha/2)^2=V_0-\alpha^2/4$. The positive-real branch is the stated $\Omega_1$. Reversing the real derivative coefficient adds $4\alpha^2(2x-1)$, nonzero as a function and at both endpoints. The constant profile could not detect this change because its first derivative vanishes.

Interpret a stable matrix root. A frequency is unchanged across three orders and has pencil residual $10^{-13}$, but its unfactored off-grid residual is $10^{-2}$. Is it an accepted resonance?

Solution

No. The finite matrices have been solved accurately, but the reconstructed function fails the original equation at the tested points and exceeds the declared $10^{-4}$ tolerance. Check the raw-to-factored transformation, mapped derivatives, and endpoint implementation before increasing resolution. The observation neither validates that candidate nor proves that no nearby physical resonance exists.

References

• Galkowski, Jeffrey, and Maciej Zworski. “Outgoing Solutions Via Gevrey-2 Properties.” Annals of PDE 7 (2021), article 5. DOI. Open PDF, arXiv v1 with cited author pagination.
• Leaver, Edward W. “An Analytic Representation for the Quasi-Normal Modes of Kerr Black Holes.” Proceedings of the Royal Society of London A 402 (1985), 285–298. DOI.
• Macedo, Rodrigo Panosso, and Anıl Zenginoğlu. “Hyperboloidal Approach to Quasinormal Modes.” Frontiers in Physics 12, 1497601, published 17 January 2025. DOI.
• SciPy community. SciPy 1.13.1 Manual (2024). BarycentricInterpolator and linalg.eig.