Perturbative Stability and Strong-Field Viability

A non-GR compact-object solution becomes physically credible only after every propagating sector has been tested around the same branch and with the same boundary conditions used to define the background. The first question is spectral—does any allowed mode grow exponentially?—but viability also requires a well-posed principal part, positive kinetic energy, controlled characteristic speeds, and frequencies below the effective-theory cutoff. A radial calculation can disprove stability; by itself it cannot establish complete stability.

Helpful background. Einstein Equations fixes the background equations; Gravitational Actions and Boundary Terms determines the quadratic action and surface terms; Formulations of General Relativity supplies constraints and hyperbolicity tests; Black-Hole Geometry and Mechanics supplies horizon generators and regular coordinates; Black-Hole Perturbation Theory supplies mode decompositions and radiative conditions; and Relativistic Matter and Compact Objects supplies stellar surfaces and matter modes.

What a stability claim must cover

Write each field as $\Phi=\Phi_0+\varepsilon\,\delta\Phi$ and expand the equations consistently through $O(\varepsilon)$. On a spherical background, parity and angular momentum often split perturbations into radial, axial, and polar systems. An extra scalar or vector can couple to the metric in only some sectors, and a fluid star adds acoustic and composition modes. Rotation mixes multipoles and opens superradiant channels. Consequently, “stable against radial perturbations” is a sector-qualified result, not shorthand for mode stability.

For a reduced master variable one frequently obtains

$$-\frac{d^2\Psi}{dx^2}+V(x)\Psi=\omega^2\Psi,$$

where $x$ is a tortoise coordinate. A self-adjoint stellar or boxed problem is unstable if it has $\omega^2\lt0$. A black-hole quasinormal problem instead imposes ingoing behavior at the future horizon and outgoing or normalizable behavior at the outer boundary; its operator is generally non-self-adjoint, and growth is determined by the sign of $\operatorname{Im}\omega$ under the declared $e^{-i\omega t}$ convention. These are different spectra.

A reproducible Schwarzschild benchmark

Before attacking coupled non-GR equations, test the discretization on the Schwarzschild axial equation. In this benchmark use geometrized units $G=c=1$ and let $M=GM_{\rm phys}$ denote the mass length. Radii and times are measured in $M$, and frequencies in $M^{-1}$. With $M=1$, $\ell=2$, and

$$x=r+2M\ln\!\left(\frac{r}{2M}-1\right), \qquad V_{\mathrm{RW}}=\left(1-\frac{2M}{r}\right) \left[\frac{\ell(\ell+1)}{r^2}-\frac{6M}{r^3}\right],$$

put Dirichlet walls at $x=-20$ and $x=40$. This deliberately produces a real, self-adjoint verification spectrum rather than physical Schwarzschild quasinormal frequencies. The original odd-parity reduction and potential are given by Regge–Wheeler 1957, Eqs. (23)–(25), pp. 1066–1067.

Save the Python fence as stability_spectrum.py. The script must produce three refinement rows, compare the finite-difference eigenvalue with independent DOP853 shooting, recover the exact free-box eigenvalue, and reject a negative sampled potential.

import sys
import numpy as np
from scipy.integrate import solve_ivp
from scipy.linalg import eigh_tridiagonal
from scipy.optimize import brentq
from scipy.special import lambertw

M, ELL = 1.0, 2
LEFT, RIGHT = -20.0, 40.0

def radius(x):
    return 2.0*M*(1.0 + lambertw(np.exp(x/(2.0*M) - 1.0)).real)

def potential(x):
    r = radius(x)
    f = 1.0 - 2.0*M/r
    return f*(ELL*(ELL + 1.0)/r**2 - 6.0*M/r**3)

def fd_mode(n, use_potential=True):
    h = (RIGHT - LEFT)/(n + 1)
    x = LEFT + h*np.arange(1, n + 1)
    diagonal = 2.0/h**2 + (potential(x) if use_potential
                           else np.zeros_like(x))
    off = -np.ones(n - 1)/h**2
    values, vectors = eigh_tridiagonal(
        diagonal, off, select="i", select_range=(0, 0))
    psi = vectors[:, 0]
    residual = np.max(np.abs(
        diagonal*psi + np.r_[off*psi[1:], 0.0]
        + np.r_[0.0, off*psi[:-1]] - values[0]*psi))
    return values[0], residual

def shoot(lam, rtol):
    def rhs(x, y):
        return (y[1], (potential(x) - lam)*y[0])
    sol = solve_ivp(rhs, (LEFT, RIGHT), (0.0, 1.0),
                    rtol=rtol, atol=rtol*1e-3, method="DOP853")
    return sol.y[0, -1]

def shooting_mode(guess, rtol):
    root = brentq(lambda lam: shoot(lam, rtol),
                  guess - 0.002, guess + 0.002,
                  xtol=rtol, rtol=max(rtol, 1e-14))
    return root, abs(shoot(root, rtol))

def main():
    free_box = (np.pi/(RIGHT - LEFT))**2
    print("N eigenvalue fd_shoot_difference free_box_error "
          "matrix_residual boundary_mismatch")
    previous = None
    for n, tol in ((400, 1e-7), (800, 1e-9), (1600, 1e-11)):
        eigenvalue, residual = fd_mode(n)
        free_value, _ = fd_mode(n, use_potential=False)
        shot, mismatch = shooting_mode(eigenvalue, tol)
        difference = abs(eigenvalue - shot)
        free_error = abs(free_value - free_box)
        print(f"{n:d} {eigenvalue:.12f} {difference:.3e} "
              f"{free_error:.3e} {residual:.3e} {mismatch:.3e}")
        if previous is not None and difference >= previous:
            raise SystemExit("finite-difference/shooting agreement did not converge")
        previous = difference
    sample = potential(np.linspace(LEFT, RIGHT, 10000))
    if eigenvalue <= 0.0 or free_error > 1e-9 or np.min(sample) < -1e-14:
        raise SystemExit("GR positivity benchmark failed")
    print(f"FREE_BOX_EXACT={free_box:.12f} "
          f"MIN_RW_POTENTIAL={np.min(sample):.3e}")

if __name__ == "__main__":
    main()

The reproducible environment is CPython version 3.9.6 with NumPy 2.0.2 and SciPy 1.13.1:

python3 -m pip install numpy==2.0.2 scipy==1.13.1
python3 stability_spectrum.py
python3 stability_spectrum.py > stability_spectrum.txt

The rerun gives:

grid $N$	eigenvalue $\omega_0^2$	finite-difference–shooting difference	free-box error	matrix residual	boundary mismatch
400	0.023373720225	$3.223\times10^{-7}$	$1.402\times10^{-8}$	$3.958\times10^{-15}$	$1.199\times10^{-4}$
800	0.023373960081	$8.022\times10^{-8}$	$3.514\times10^{-9}$	$1.451\times10^{-14}$	$1.258\times10^{-8}$
1600	0.023374020219	$2.008\times10^{-8}$	$8.795\times10^{-10}$	$1.052\times10^{-13}$	$1.097\times10^{-8}$

The method-to-method difference falls by about four when $N$ doubles, as expected for the centered second derivative. With $V=0$, the exact answer is $[\pi/(40-(-20))]^2=0.002741556778$; the last free-box error is below $10^{-9}$. For $\ell=2$ the potential factorizes as

$$V_{\mathrm{RW}}=\frac{6(r-2M)(r-M)}{r^4}\gt0 \qquad (r\gt2M).$$

Thus the continuum Dirichlet quadratic form is positive independently of either numerical method. The sampled minimum $1.253\times10^{-5}$ on this finite interval checks that the implemented tortoise-coordinate inversion preserves that sign.

The underlying proof of positivity is $\omega^2\int\lvert\Psi\rvert^2dx =\int(\lvert\Psi'\rvert^2+V\lvert\Psi\rvert^2)dx$, with the boundary term vanishing at the two Dirichlet walls. The measure is $dx$, and the operator domain is $H^2\cap H_0^1$ on this finite interval. A small matrix residual only shows that the discrete eigenpair was solved accurately; it is not this continuum argument.

Transferring the workflow to a non-GR branch

First compute a background with quantified field-equation and boundary residuals. Then derive the perturbation system from the same action, eliminate only genuine constraints, inspect the kinetic and gradient matrices before division by background-dependent coefficients, and impose regular horizon, center, surface, and asymptotic conditions in regular variables. Continue each eigenvalue along a named branch rather than sorting spectra anew at every point; a zero mode should meet a bifurcation or turning-point prediction when one exists.

Scalarized Einstein–Gauss–Bonnet black holes show why branch labels matter. Blázquez-Salcedo and collaborators derive a radial master equation and its horizon/infinity conditions in Blázquez-Salcedo et al. 2018, § 3, Eqs. (16)–(22), pp. 5–6 (PDF), then find that nodal branches are radially unstable while the fundamental branch depends on the coupling function Blázquez-Salcedo et al. 2018, § 4, pp. 6–15 (PDF). The conclusion is not “scalarized holes are stable” or “unstable”; it is a result conditional on branch, coupling, sector, and parameter range.

Beyond eigenvalues: strong-field viability

A positive sampled spectrum is meaningful only where the perturbation equations remain hyperbolic. Vanishing kinetic coefficients, a change of characteristic signature, or a singular field redefinition can invalidate the master equation before an eigenvalue crosses zero. The cited scalarized solutions include a small-mass region where the radial equation loses hyperbolicity; formal mode results there are not a well-posed evolution statement Blázquez-Salcedo et al. 2018, § 4 and Conclusions, pp. 12–16 (PDF).

Complete the analysis with axial and polar gravitational sectors, extra-field polarizations, radial fluid modes for stars, and rotational or superradiant sectors where applicable. Check kinetic-sign and gradient matrices, constraint propagation, asymptotic flux, and cutoff hierarchy. Repeat domain-size and boundary-condition studies: a boxed mode can drift as either wall moves, while a genuine localized instability should converge. Spurious roots often reveal themselves through unresolved oscillations, failed residual convergence, or sensitivity to the eliminated variable.

Exercise: interpret a negative mode

Suppose a modified self-adjoint boxed perturbation operator returns $\omega_0^2=-0.04/M^2$. What time dependence follows, and which claim is justified?

Solution

With perturbations proportional to $e^{-i\omega t}$, the roots are $\omega=\pm0.2i/M$. The root $+0.2i/M$ grows as $e^{0.2t/M}$, so its e-folding time is $5M$. If the eigenfunction satisfies the physical boundary conditions and converges under grid and domain refinement, this disproves mode stability in that sector on that branch. It does not identify a nonlinear endpoint, prove another branch stable, or test sectors omitted from the reduction. If the principal part is non-hyperbolic, even the exponential interpretation is not a valid well-posed evolution claim.

In physical units the time in the exercise is $5GM_{\rm phys}/c^3$. For the positive Schwarzschild potential and Dirichlet walls actually used above, a negative eigenvalue would contradict the quadratic-form proof and reveal a numerical or implementation error; the hypothetical value is relevant only after the operator or boundary conditions have been changed.

References

• Jose Luis Blázquez-Salcedo, Daniela D. Doneva, Jutta Kunz, and Stoytcho S. Yazadjiev, “Radial perturbations of the scalarized Einstein-Gauss-Bonnet black holes,” Physical Review D 98 (2018) 084011. DOI; Open PDF.
• Tullio Regge and John A. Wheeler, “Stability of a Schwarzschild Singularity,” Physical Review 108 (1957) 1063–1069. DOI.