Scattering Amplification Factors

An amplification factor is a signed energy-flux ratio, not a ratio of raw Teukolsky amplitudes. This page gives a reproducible scalar Kerr benchmark: integrate a future-ingoing Frobenius solution, fit incident and outgoing waves at infinity, and compare the amplitude gain with the independently conserved Wronskian flux. The scalar calculation exercises the threshold and Schwarzschild checks without importing the spin-dependent normalization needed for electromagnetic or gravitational fields.

Current computational scope. Only the neutral massless scalar sector is numerically evaluated here. Electromagnetic and gravitational amplification are discussed as extensions, not reported as reproduced benchmarks.

Helpful background. Black-Hole Geometry and Mechanics supplies Kerr horizon data; PDEs and Green Functions for Gravitation supplies ODE and Wronskian methods; Formulations of General Relativity supplies the separated field equation; Relativistic Matter and Compact Objects supplies field content; Black-Hole Perturbation Theory supplies endpoint normalizations.

Scalar benchmark and flux formula

The run uses $G=c=M=1$ internally, so $a=\chi_{\mathrm{BH}}=0.9$, and selects the neutral massless scalar mode $s=0$, $\ell=m=1$. The horizon radii and angular velocity are

$$r_\pm=1\pm\sqrt{1-a^2},\qquad \Omega_H=\frac{a}{r_+^2+a^2}.$$

All frequencies below are real and positive. With

$$\Delta=(r-r_+)(r-r_-),\qquad K=(r^2+a^2)\omega-am,$$

the radial equation is

$$\frac{d}{dr}\left(\Delta\frac{dR}{dr}\right) +\left(\frac{K^2}{\Delta}-\lambda_{\ell m}\right)R=0,$$

where $\lambda_{\ell m}=A_{\ell m}+a^2\omega^2-2am\omega$ and $A_{\ell m}$ is the scalar oblate-spheroidal eigenvalue. Near the future horizon,

$$R=x^p\sum_{n=0}^{N_H}c_nx^n,\qquad x=r-r_+,\qquad p=-i\frac{(r_+^2+a^2)(\omega-m\Omega_H)}{r_+-r_-}.$$

The initializer derives the coefficients recursively from the radial equation; it does not differentiate a leading plane wave at a finite endpoint. It fixes $c_0=1$, an amplitude convention rather than unit incident energy. The tortoise coordinate obeys $dr_*/dr=(r^2+a^2)/\Delta$; its additive constant changes phases but not gain. At infinity the fit uses

$$rR= e^{+i\omega r_*}\sum_{n=0}^{N_\infty} a_n\left(\frac{r_f}{r}\right)^n +e^{-i\omega r_*}\sum_{n=0}^{N_\infty} b_n\left(\frac{r_f}{r}\right)^n.$$

Here $r_f$ is the beginning of the fitting interval; $r_f/r$ is dimensionless and at most one. The powers model finite-radius corrections to both wave branches. Thus $A_{\mathrm{out}}=a_0$ and $A_{\mathrm{in}}=b_0$. The two gain estimates are

$$Z_{\mathrm{amp}} =\left\lvert\frac{A_{\mathrm{out}}}{A_{\mathrm{in}}}\right\rvert^2-1,$$

and, for unit horizon-series coefficient,

$$Z_H =-\frac{(r_+^2+a^2)(\omega-m\Omega_H)} {\omega\lvert A_{\mathrm{in}}\rvert^2}.$$

The conserved scalar Wronskian current is $J=\Delta\,\operatorname{Im}(R^*R')$, with prime denoting $d/dr$. Multiplying the real-frequency radial equation by $R^*$ and taking its imaginary part gives $dJ/dr=0$. The two endpoint limits are

$$\begin{aligned} J_\infty&=\omega\bigl(|A_{\mathrm{out}}|^2-|A_{\mathrm{in}}|^2\bigr),\\ J_H&=-(r_+^2+a^2)(\omega-m\Omega_H). \end{aligned}$$

Their equality gives $Z_{\mathrm{amp}}=Z_H$. The physical outward Killing-energy current has an additional common positive factor proportional to $\omega$, which cancels in these ratios.

The canonical scalar scattering field is $\Psi=\sqrt{r^2+a^2}\,R$. Its current $\operatorname{Im}(\Psi^*\,d\Psi/dr_*)$ is exactly $J$: differentiating the real rescaling contributes no imaginary part. At infinity $\Psi\sim rR$, whereas unit raw horizon amplitude corresponds to canonical squared amplitude $r_+^2+a^2$. This is the origin of the horizon factor, not an extra gain or a factor that can be dropped at both endpoints. The scalar conversion appears in Teukolsky and Press 1974, § V, Eqs. (5.3) and (5.6), p. 452; gravitational fluxes require the separate normalization in Teukolsky and Press 1974, § IV, pp. 449–450. The relation between signed horizon flux, reflection coefficient, and superradiant threshold is reviewed with explicit Kerr scattering normalizations in Brito, Cardoso, and Pani 2020, §4.7, pp. 57–68.

Environment and code

The 2026-09-12 rerun used macOS 26.6.2 arm64, Python 3.9.6, NumPy 2.0.2, and SciPy 1.13.1. The code below is complete. SciPy’s obl_cv supplies the scalar oblate eigenvalue $A_{\ell m}$; the shifted radial constant is $\lambda_{\ell m}$ as defined above. Its angular equation is the $s=0$ case on Spin-Weighted Spheroidal Harmonics. Changing the angular convention without translating the radial shift is a detectable failure SciPy 1.13.1, obl_cv reference.

For the complex state $(R,R')$, DOP853 controls a local error estimate using the scale atol plus rtol times the magnitude of each state component. Neither tolerance bounds the final gain: initialization, propagation, fitting and subtraction contribute separate errors SciPy 1.13.1, solve_ivp, method and tolerance parameters.

import argparse
import csv
import numpy as np
from scipy.integrate import solve_ivp
from scipy.special import obl_cv


def horizon_data(a, ell, m, w, order):
    rp = 1.0 + np.sqrt(1.0 - a * a)
    rm = 1.0 - np.sqrt(1.0 - a * a)
    d = rp - rm
    omega_h = a / (rp * rp + a * a)
    k_h = w - m * omega_h
    k0 = (rp * rp + a * a) * k_h
    k1 = 2.0 * rp * w
    angular = float(obl_cv(abs(m), ell, a * w))
    lam = angular + a * a * w * w - 2.0 * a * m * w
    p = -1j * k0 / d

    a2 = {2: d * d, 3: 2.0 * d, 4: 1.0}
    a1 = {1: d * d, 2: 3.0 * d, 3: 2.0}
    a0 = {
        0: k0 * k0,
        1: 2.0 * k0 * k1 - lam * d,
        2: k1 * k1 + 2.0 * k0 * w - lam,
        3: 2.0 * k1 * w,
        4: w * w,
    }
    coeff = np.zeros(order + 1, dtype=complex)
    coeff[0] = 1.0
    for n in range(1, order + 1):
        known = 0j
        for j, value in a2.items():
            k = n - j + 2
            if 0 <= k < n:
                known += value * (p + k) * (p + k - 1) * coeff[k]
        for j, value in a1.items():
            k = n - j + 1
            if 0 <= k < n:
                known += value * (p + k) * coeff[k]
        for j, value in a0.items():
            k = n - j
            if 0 <= k < n:
                known += value * coeff[k]
        denominator = d * d * (p + n) ** 2 + k0 * k0
        coeff[n] = -known / denominator
    return rp, rm, omega_h, k_h, angular, lam, p, coeff


def gain(w, a=0.9, ell=1, m=1, horizon_offset=1e-5,
         horizon_order=8, r_max=400.0, fit_min=180.0,
         infinity_order=4, rtol=2e-11, atol=2e-13):
    rp, rm, omega_h, k_h, angular, lam, p, coeff = horizon_data(
        a, ell, m, w, horizon_order
    )
    x0 = horizon_offset
    powers = x0 ** np.arange(len(coeff))
    r0 = x0 ** p * np.sum(coeff * powers)
    dr0 = x0 ** (p - 1) * np.sum(
        (p + np.arange(len(coeff))) * coeff * powers
    )

    def rhs(r, y):
        delta = (r - rp) * (r - rm)
        delta_prime = 2.0 * r - 2.0
        radial_k = (r * r + a * a) * w - a * m
        return np.array([
            y[1],
            -(delta_prime * y[1] + (radial_k ** 2 / delta - lam) * y[0])
            / delta,
        ], dtype=complex)

    radii = np.linspace(fit_min, r_max, 1200)
    solution = solve_ivp(
        rhs, (rp + x0, r_max), [r0, dr0], t_eval=radii,
        method="DOP853", rtol=rtol, atol=atol
    )
    if not solution.success:
        raise RuntimeError(solution.message)
    radial, derivative = solution.y
    d = rp - rm
    r_star = (
        radii
        + (rp * rp + a * a) / d * np.log(radii - rp)
        - (rm * rm + a * a) / d * np.log(radii - rm)
    )
    scaled_inverse_r = fit_min / radii
    columns = []
    for sign in (1, -1):
        phase = np.exp(1j * sign * w * r_star)
        columns.extend(
            phase * scaled_inverse_r ** n
            for n in range(infinity_order + 1)
        )
    design = np.column_stack(columns)
    fitted, *_ = np.linalg.lstsq(design, radii * radial, rcond=None)
    a_out = fitted[0]
    a_in = fitted[infinity_order + 1]
    fit_residual = np.linalg.norm(
        design @ fitted - radii * radial
    ) / np.linalg.norm(radii * radial)

    gain_amplitude = abs(a_out / a_in) ** 2 - 1.0
    gain_horizon = -(
        (rp * rp + a * a) * k_h
    ) / (w * abs(a_in) ** 2)
    flux_residual = abs(
        w * (abs(a_out) ** 2 - abs(a_in) ** 2)
        + (rp * rp + a * a) * k_h
    ) / (
        w * (abs(a_out) ** 2 + abs(a_in) ** 2)
        + abs((rp * rp + a * a) * k_h)
    )
    delta = (radii - rp) * (radii - rm)
    current = delta * np.imag(np.conj(radial) * derivative)
    current_scale = max(np.max(np.abs(current)), 1e-30)
    current_drift = np.max(np.abs(current - current[0])) / current_scale
    return {
        "Momega": w,
        "OmegaH": omega_h,
        "angular_eigenvalue": angular,
        "gain_amplitude": gain_amplitude,
        "gain_horizon": gain_horizon,
        "flux_residual": flux_residual,
        "current_drift": current_drift,
        "fit_residual": fit_residual,
    }


parser = argparse.ArgumentParser()
parser.add_argument("--output", default="scalar_gain.csv")
args = parser.parse_args()
omega_h = 0.9 / (2.0 * (1.0 + np.sqrt(1.0 - 0.9 ** 2)))
rows = [gain(w) for w in (0.10, 0.20, 0.30, omega_h, 0.34)]
with open(args.output, "w", newline="", encoding="utf-8") as stream:
    writer = csv.DictWriter(stream, fieldnames=rows[0].keys())
    writer.writeheader()
    writer.writerows(rows)
for row in rows:
    print(row)

A clean environment and machine-readable rerun are produced by saving that block as scalar_gain.py and running with Python 3.9.6:

python3 -m venv .venv-superradiance
.venv-superradiance/bin/python -m pip install \
  numpy==2.0.2 scipy==1.13.1
.venv-superradiance/bin/python scalar_gain.py \
  --output scalar_gain.csv

Choose a new CSV filename if an earlier result must be preserved: this baseline program opens its output in write mode. It writes the five-frequency scan and diagnostics, not the separate variation suite or raw complex amplitudes. To reproduce a variation, call gain with the parameters in the convergence table in place of a baseline call.

Recorded scan and residuals

The following values are from the 2026-09-12 rerun of the unchanged program first recorded on 2026-08-12. The gain is dimensionless; positive values mean amplification and negative values mean net absorption.

$M\omega$	Regime	$Z_{\mathrm{amp}}$	Flux-balance residual	Infinity-fit residual
0.100000000	superradiant	$2.49575298\times10^{-4}$	$3.64\times10^{-12}$	$1.51\times10^{-11}$
0.200000000	superradiant	$1.34409395\times10^{-3}$	$1.33\times10^{-12}$	$1.60\times10^{-11}$
0.300000000	superradiant	$1.34134786\times10^{-3}$	$1.84\times10^{-13}$	$1.91\times10^{-11}$
0.313394503	threshold	$2.74\times10^{-13}$	$1.37\times10^{-13}$	$1.96\times10^{-11}$
0.340000000	absorptive	$-7.82231806\times10^{-3}$	$9.70\times10^{-13}$	$2.18\times10^{-11}$

For this spin, $M\Omega_H=0.313394503137$. The threshold result is numerical zero at the precision supported by the asymptotic fit. The nonthreshold relative current drift was at most $1.91\times10^{-10}$ over the 1200 sampled points in each fitting interval, not over the whole integration domain. At threshold the analytic current and all sampled currents were zero in this run. The code returns zero after applying its denominator floor, but a relative error with respect to an exact zero current is undefined; absolute current error or error normalized by incident flux is the meaningful test. The sign change agrees with $\omega(\omega-m\Omega_H)<0$ below threshold and with ordinary absorption above it. This is a controlled scalar benchmark, not a claim about the larger spin-1 or spin-2 maxima tabulated by Teukolsky and Press.

The Schwarzschild benchmark is structural: setting $a=0$ gives $\Omega_H=0$, so every positive-frequency row must have $Z\leq0$. The code also recovers $Z=0$ at the Kerr threshold without replacing the signed horizon factor by its absolute value.

Convergence and independent checks

At $M\omega=0.2$, the baseline result is $Z_{\mathrm{amp}}=0.00134409395454$. The following variations were rerun; the ODE column specifies relative tolerance, with absolute tolerance $2\times10^{-13}$ held fixed in all five rows:

Variation	Horizon offset/order	Far endpoint/window/order	ODE tolerance	$Z_{\mathrm{amp}}$
Coarse	$10^{-4}/6$	$300/[120,300]/3$	$10^{-10}$	0.00134409394446
Baseline	$10^{-5}/8$	$400/[180,400]/4$	$2\times10^{-11}$	0.00134409395454
Tight	$3\times10^{-6}/10$	$500/[220,500]/5$	$5\times10^{-12}$	0.00134409395476
Later fit window	$10^{-5}/8$	$500/[250,500]/4$	$2\times10^{-11}$	0.00134409395483
Higher infinity order	$10^{-5}/8$	$400/[180,400]/5$	$2\times10^{-11}$	0.00134409395635

The full spread is $1.189\times10^{-11}$, approximately $1.2\times10^{-11}$ in absolute gain. It shows stable gain near $1.34409\times10^{-3}$; the longer digits facilitate rerun comparisons, not a certified error bound. Last digits differ slightly from the older recording; identical solver versions do not guarantee bitwise-identical numerical linear algebra. These parameter refinements were performed only at $M\omega=0.2$.

The flux-balance residual divides the current imbalance by the sum of the positive infinity-current magnitudes and the magnitude of the horizon current. Calling it $\epsilon_{\mathrm{flux}}$, its relation to gain disagreement is

$$\lvert Z_{\mathrm{amp}}-Z_H\rvert =\epsilon_{\mathrm{flux}}\left(2+Z_{\mathrm{amp}}+\lvert Z_H\rvert\right).$$

For the baseline $\omega=0.2$ row, the absolute disagreement is $2.6651\times10^{-12}$, while $\epsilon_{\mathrm{flux}}$ is $1.3307\times10^{-12}$. These are different diagnostics, not inconsistent answers. The infinity-fit residual is a relative Euclidean norm of the sampled fit error; it is not an error bound on either leading coefficient. The signed horizon current is analytically known and independent of the outgoing infinity fit, although $Z_H$ still uses the fitted incident amplitude. Current conservation, fitting residual and parameter refinement test complementary aspects of the same ODE calculation; they are not a second independent numerical method.

Failure modes and higher-spin work not performed here

Near $\omega=m\Omega_H$, subtracting two nearly equal infinity fluxes loses relative precision. The signed horizon expression is the stable cross-check. An infinity window that begins too early contaminates the fit with omitted inverse-radius terms; a window that is too short cannot separate the two phases. A horizon offset that is too large exposes Frobenius truncation, while one too small can amplify floating-point phase error.

For electromagnetic and gravitational Teukolsky variables, the same workflow applies but the raw $R$ amplitudes are not scalar energy amplitudes. One must insert the Teukolsky–Starobinsky constants, use a horizon-regular tetrad normalization, and compare the correct spin-dependent flux expressions. Applying the scalar formula to $s=-1$ or $s=-2$ would generate plausible but wrong gains. High precision is also needed close to extremality, where endpoint scales separate strongly.

No spin-1 or spin-2 number is produced by this page’s code. A valid extension must solve the appropriate Teukolsky equation, convert its raw asymptotic coefficients with the Teukolsky–Starobinsky identities, evaluate the spin-dependent energy flux at both endpoints, and recover selected published gains before claiming a three-spin benchmark. Teukolsky and Press give the primary normalization and amplification curves for electromagnetic and gravitational radiation Teukolsky and Press 1974, full pp. 443–461. Citing those curves is not a rerun of them.

The frequency condition is derived on The Superradiant Frequency Condition. Repeated gain under confinement is treated in Black-Hole Bombs and Reflecting Boundaries.

Exercises

For the baseline $M\omega=0.2$ row, $\lvert A_{\mathrm{in}}\rvert=34.80500955$, $r_+^2+a^2=2.87177979$, and $\omega-m\Omega_H=-0.1133945031$. Compute $Z_H$ and compare it with the table.

Solution

The horizon formula gives

$$Z_H =-\frac{(2.87177979)(-0.1133945031)} {(0.2)(34.80500955)^2} =0.00134409395.$$

It agrees with $Z_{\mathrm{amp}}=0.00134409395$ at the displayed precision. The positive sign follows from the negative horizon-frame frequency. Replacing that frequency by its magnitude would predict negative gain and fail both the threshold scan and the conserved-current balance.

References

• Brito, R., Cardoso, V., and Pani, P. (2020). Superradiance: New Frontiers in Black Hole Physics. Lecture Notes in Physics 971. Springer. DOI.
• SciPy community (2024). SciPy 1.13.1 Reference Guide: obl_cv and solve_ivp.
• Teukolsky, S. A., and Press, W. H. (1974). “Perturbations of a rotating black hole. III. Interaction of the hole with gravitational and electromagnetic radiation.” Astrophysical Journal 193, 443–461. DOI.