Numerical Spectra with Nonflat Boundaries

Numerical spectra with a cosmological horizon or timelike AdS boundary are defined as much by their endpoint domains as by the radial differential equation. A converged eigenvalue with the wrong cosmological-horizon orientation or an accidentally asymptotically flat outer condition is a precise answer to the wrong problem. This page builds and checks a horizon-regular Schwarzschild–AdS benchmark. Its de Sitter branch is only an analytic formula smoke check: it does not discretize a de Sitter operator or test either horizon orientation.

Helpful background. Black-Hole Geometry and Mechanics supplies horizon parameters; Partial Differential Equations and Green Functions supplies generalized eigenproblems; Formulations of General Relativity supplies field-specific master variables; Gravitational Actions and Boundary Terms supplies flux tests; Causality and Global Structure supplies nonflat endpoints; Black-Hole Perturbation Theory supplies the $e^{-i\omega t}$ convention.

A horizon-regular Schwarzschild–AdS problem

Set $G=c=L=1$ and consider a massless scalar on four-dimensional Schwarzschild–AdS,

$$f(r)=1+r^2-\frac{r_+(1+r_+^2)}{r}.$$

Here $r_+$ is the event-horizon radius, not the mass parameter: $2GM=r_+(1+r_+^2)$ in these units. With $dr_*/dr=f^{-1}$ and $\phi=e^{-i\omega t}Y_{\ell m}\,u(r)/r$, the tortoise equation is

$$\frac{d^2u}{dr_*^2}+[\omega^2-V_\ell(r)]u=0, \qquad V_\ell=f\left[\frac{\ell(\ell+1)}{r^2}+\frac{f'}r\right].$$

Factor the future-horizon phase by writing $u=e^{-i\omega r_*}y(r)$. The equation becomes linear in $\omega$,

$$f y''+(f'-2i\omega)y' -\left[\frac{f'}r+\frac{\ell(\ell+1)}{r^2}\right]y=0.$$

Primes here denote $d/dr$. In ingoing coordinates $v=t+r_*$, the full field is $e^{-i\omega v}Y_{\ell m}y/r$. The physical ingoing branch is analytic in $r-r_+$ after phase extraction. Merely requiring a finite $y$ is weaker, as we will check explicitly.

At infinity the scalar branches are a constant source and a fast $r^{-3}$ falloff. This page selects only the latter, so $y\sim r^{-2}$; it does not impose an asymptotically flat outgoing wave at the timelike boundary. These conditions and the regular-coordinate equation are given by Horowitz and Hubeny 2000, arXiv v2, Eqs. (2.11)–(2.17), pp. 6–7, and § 3, pp. 9–11, PDF. Compactify with $x=r_+/r\in[0,1]$ and set $y=x^2F(x)$. The AdS boundary is $x=0$ and the horizon is $x=1$.

Since $d/dr=-(x^2/r_+)d/dx$, substitution gives the following equation after division by a common $x^3$. Primes on $F$ now denote $d/dx$:

$$a_2F''+(a_1+\omega b_1)F'+(a_0+\omega b_0)F=0,$$

with

$$\begin{aligned} a_2&=x\left[1+\frac{x^2}{r_+^2}-x^3\left(1+\frac1{r_+^2}\right)\right],\\ a_1&=4-7x^3+\frac{x^2(6-7x)}{r_+^2}, &b_1&=\frac{2ix}{r_+},\\ a_0&=-9x^2+\frac{x[6-\ell(\ell+1)-9x]}{r_+^2}, &b_0&=\frac{4i}{r_+}. \end{aligned}$$

Because $a_2$ vanishes at both endpoints, a smooth solution satisfies first-order endpoint relations. At the AdS boundary, $F'(0)=-i\omega F(0)/r_+$. Keeping $F$ bounded excludes the divergent source branch. At the horizon the endpoint row is necessary but does not alone select analytic ingoing behavior.

To see why, put $z=r-r_+$ and insert $y\sim z^\rho$ into the radial equation. Its most singular coefficient gives

$$\rho\bigl(f'(r_+)\rho-2i\omega\bigr)=0.$$

For our $r_+=1$, $\ell=0$ benchmark, $f'(r_+)=4$ and the unwanted branch is

$$y_{\mathrm{out}}=z^\rho(1+O(z)),\qquad \rho=\frac{i\omega}{2}\simeq1.335603+1.399112i.$$

Both $y_{\mathrm{out}}$ and its first derivative tend to zero, so this branch also satisfies the limiting first-order row. Its second derivative is singular; it is not analytic at the future horizon. The polynomial collocation ansatz approximates the smooth branch, but a small endpoint residual by itself cannot certify that domain. This is precisely the outgoing-branch distinction in Horowitz and Hubeny 2000, arXiv v2, Eq. (2.17) and discussion, p. 7, PDF.

Collocation, environment, and mode selection

Use Chebyshev–Lobatto points on $x\in[0,1]$. Replacing derivatives by matrices gives

$$(A+\omega B)\mathbf F=0,$$

which is solved as the generalized eigenproblem $A\mathbf F=-\omega B\mathbf F$. The benchmark environment is CPython 3.9.6, NumPy 2.0.2, and SciPy 1.13.1 using scipy.linalg.eig in complex double precision. Inputs are $r_+=L=1$, $\ell=0$, the fast falloff, and resolutions $N=16$ through $40$.

The target branch is identified by the published value $\omega L\simeq2.7982-2.6712i$: the four-dimensional, $\ell=0$, $r_+=L$ row of Horowitz and Hubeny 2000, arXiv v2, Table 1, p. 12, PDF. The code uses a longer approximate reference to select the same nearby root, not as independent evidence for those extra digits. Proximity is a benchmark label, not a discovery algorithm or proof that no slower mode exists. A parameter scan should continue eigenvectors by overlap and check conditioning, not select “the least damped number” from an unfiltered pencil. The analytic horizon-series method is reviewed by Berti, Cardoso, and Starinets 2009, arXiv v2, § 4.4, pp. 28–29, PDF; their § 6.1, pp. 43–48 discusses distinct black-hole size and angular-momentum regimes.

Convergence and independent checks

The 2026-09-12 rerun of the unchanged collocation program reproduced:

$N$	Tracked $\omega L$	Change from previous row	Pencil residual
16	$2.7982405124015-2.6712350024504i$	—	$1.3\times10^{-14}$
20	$2.7982230527723-2.6712057079494i$	$3.4\times10^{-5}$	$5.4\times10^{-14}$
24	$2.7982232441512-2.6712058260310i$	$2.2\times10^{-7}$	$5.4\times10^{-14}$
28	$2.7982232428215-2.6712058258154i$	$1.3\times10^{-9}$	$6.6\times10^{-14}$
32	$2.7982232428290-2.6712058257850i$	$3.1\times10^{-11}$	$1.6\times10^{-13}$
40	$2.7982232428626-2.6712058258080i$	$4.1\times10^{-11}$	$1.4\times10^{-13}$

The residual is

$$\epsilon_{\mathrm{pencil}}= \frac{\lVert A\mathbf F+\omega B\mathbf F\rVert_\infty} {\lVert A\mathbf F\rVert_\infty+ \lvert\omega\rvert\lVert B\mathbf F\rVert_\infty},$$

including both endpoint rows. It measures satisfaction of the finite matrix problem, not the continuum eigenvalue error. The plateau after $N\approx32$ is consistent with double-precision roundoff; this sweep alone does not isolate its cause. No higher-precision or mapping sweep was performed. The long digits record numerical outputs, not a certified error bound. The independent integration below tests the same frequency while enforcing the analytic horizon domain differently. Negative imaginary part is also a necessary consistency check for this passive scalar benchmark, not a substitute for those boundary tests.

Dimensional consistency provides another inexpensive check. Denote the scalar field mass by $\mu_\phi$ to distinguish it from the black-hole length $\mu=GM$ used elsewhere in this volume. The Schwarzschild–AdS benchmark above sets $\mu_\phi=0$; the pure-de Sitter smoke check below displays the general $\mu_\phi L$ formula and then evaluates its massless case. Restoring $L$, the inputs enter only through $r_+/L$, $\mu_\phi L$, and angular integers, while the output is $\omega L$. Rescaling every length by the same factor cannot change the displayed dimensionless eigenvalue. The endpoint equations also have a useful asymmetry: the horizon row is necessary for the chosen analytic future-ingoing branch, whereas the boundary row encodes the chosen AdS falloff. Swapping those rows or reversing the compact coordinate should leave the physical result invariant only after the derivative orientation is transformed consistently.

Independent shooting and the horizon domain

An independent calculation integrates the original radial equation in $r$, not the collocation matrices or their compactified coefficients. Fix $r_+=L=1$, $\ell=0$. At a large finite radius $R_b$, initialize the normalizable branch from

$$y=r^{-2}\sum_{n=0}^{N_b}b_n r^{-n},\qquad b_0=1.$$

Substitution into the original equation gives

$$\begin{aligned} 0={}&n(n+3)b_n+2i\omega(n+1)b_{n-1}\\ &+n(n+1)b_{n-2}-2n^2b_{n-3}, \end{aligned}$$

with negative-index coefficients zero. Thus $b_1=-i\omega$ and $b_2=-3(\omega^2+1)/5$. The fixed $b_0$ is an arbitrary nonzero normalization, not a boundary source.

Integrate $(y, dy/dr)$ inward using DOP853. At a positive offset $z_m=r-1$, construct two independent local series,

$$y_{\mathrm{in}}=\sum_{n=0}^{N_h}a_n z^n,\qquad y_{\mathrm{out}}=z^{i\omega/2} \sum_{n=0}^{N_h}\widetilde a_n z^n,$$

each with leading coefficient one. Powers use the real logarithm of $z\gt0$. Multiplying the radial equation by $r^3$ gives a polynomial operator $P y''+Q y'+U y=0$, where

$$\begin{aligned} P&=r^5+r^3-2r^2,\\ Q&=2r^4+2r-2i\omega r^3,\qquad U=-2r^3-2. \end{aligned}$$

Expand these polynomials about $r=1$ and equate successive powers. For exponent $\rho$, the coefficient multiplying the new $a_n$ is $(n+\rho)[4(n+\rho)-2i\omega]$; lower coefficients determine its right-hand side. In particular, $a_1^{\mathrm{in}}=2/(2-i\omega)$. The code below implements this recurrence for both $\rho=0$ and $\rho=i\omega/2$. Resonant denominators require a different local construction and are excluded from this benchmark.

At the matching point solve

$$\begin{pmatrix} y_{\mathrm{in}}&y_{\mathrm{out}}\\ y'_{\mathrm{in}}&y'_{\mathrm{out}} \end{pmatrix} \begin{pmatrix}A\\B\end{pmatrix} =\begin{pmatrix}y\\y'\end{pmatrix}.$$

Then find a zero of the complex ratio $B/A$, starting from $\omega=2.8-2.7i$. This explicitly removes the unwanted outgoing branch. It does not identify a mode by testing whether the total field happens to be small at the horizon. A nonzero $A$ and a reasonably conditioned local basis are essential.

The 2026-09-12 calculation varied boundary radius, horizon offset, both series orders and ODE tolerances in eleven fixed cases. Four representative rows are shown below; the full sequence is in the executable. The orders denote $(N_b,N_h)$ and tolerances denote (relative, absolute).

Case	$R_b$	$z_m$	Orders	Tolerances	Computed $\omega L$
Coarse series	8	0.05	(6,4)	$(10^{-9},10^{-12})$	$2.7981660684-2.6711670595i$
Boundary refinement	10	0.02	(12,6)	$(10^{-11},10^{-14})$	$2.7982232425-2.6712058256i$
Horizon refinement	20	0.005	(16,8)	$(10^{-11},10^{-14})$	$2.7982232428-2.6712058258i$
Tighter integration	20	0.005	(16,8)	$(10^{-12},10^{-15})$	$2.7982232428-2.6712058258i$

All eleven cases agree with the collocation target to four significant figures, $\omega L=2.798-2.671i$. Even the coarse series differs by only $6.91\times10^{-5}$ in complex absolute frequency; refined rows agree much more closely. The two-column basis condition numbers range from about 2.29 to 13.40 in the stated $(y,y')$ scaling. These are finite convergence and conditioning observations, not a rigorous eigenvalue enclosure. A tiny $B/A$ at one chosen truncation is insufficient—the coarse row also has a tiny root residual. The series and endpoint variations test the error that root finding alone cannot detect.

The de Sitter smoke check

For a Schwarzschild–de Sitter implementation, compactify the interval between $r_+$ and $r_c$, factor the future-event-horizon and future-cosmological-horizon powers separately, and linearize the resulting quadratic pencil. A separate limiting check uses pure de Sitter space: there the inner condition is regularity at the center, not an event-horizon condition. With outgoing behavior at the cosmological horizon, the scalar frequencies are

$$\omega L =-i(2n+\ell+h_\pm), \qquad h_\pm=\frac32\pm\sqrt{\frac94-\mu_\phi^2L^2},$$

with $n=0,1,\ldots$ in the appropriate static-patch convention. The script below exposes this analytic ds-limit case only as a formula and Fourier-sign smoke check. Because it substitutes numbers directly into the closed-form expression, it cannot test a de Sitter differential operator, compactification, Frobenius factor, signed current, or cosmological-horizon orientation. It also does not pretend that those values are Schwarzschild–de Sitter modes.

Current limitation. No de Sitter eigenvalue computed by the code below is accepted on this page. A genuine Schwarzschild–de Sitter benchmark still requires the two-horizon boundary-value problem described in the preceding paragraph, a resolution and precision sweep, reconstruction in regular $u$ and $v$ coordinates, and comparison with an independent method.

Reproduce the calculations

Use Python 3.9.6 with NumPy 2.0.2 and SciPy 1.13.1 to reproduce the recorded environment. The calculations use complex double precision; different numerical libraries or platforms may change the final digits. The local error tolerances are not final frequency error bars SciPy 1.13.1, solve_ivp, DOP853 and tolerance parameters.

Collocation program

Save the code as nonflat_qnm.py. Run python3 nonflat_qnm.py --case ads for the collocation sweep or python3 nonflat_qnm.py --case ds-limit for the pure-de Sitter analytic formula smoke check. The AdS branch emits eigenvalues, changes, and pencil residuals; the de Sitter branch explicitly reports that no operator or orientation was tested.

import argparse, json, platform
import numpy as np
import scipy
from scipy.linalg import eig

REFERENCE = 2.7982232428 - 2.6712058258j

def chebyshev(n):
    j = np.arange(n + 1)
    t = np.cos(np.pi * j / n)
    c = np.ones(n + 1)
    c[[0, -1]] = 2
    c *= (-1.0) ** j
    T = np.tile(t, (n + 1, 1))
    D = np.outer(c, 1 / c) / (T.T - T + np.eye(n + 1))
    D -= np.diag(D.sum(axis=1))
    return (1 + t) / 2, 2 * D

def ads_mode(n, rh=1.0, ell=0):
    x, D = chebyshev(n)
    D2 = D @ D
    a2 = x * (1 + x*x/rh**2 - x**3 * (1 + 1/rh**2))
    a1 = 4 - 7*x**3 + x*x*(6 - 7*x)/rh**2
    a0 = -9*x*x + x*(6 - ell*(ell + 1) - 9*x)/rh**2
    A = np.diag(a2) @ D2 + np.diag(a1) @ D + np.diag(a0)
    B = np.diag(2j*x/rh) @ D + np.diag(np.full(n + 1, 4j/rh))
    values, vectors = eig(A, -B)
    k = int(np.nanargmin(np.abs(values - REFERENCE)))
    w, v = values[k], vectors[:, k]
    numerator = np.linalg.norm(A @ v + w * (B @ v), np.inf)
    scale = np.linalg.norm(A @ v, np.inf) + abs(w) * np.linalg.norm(B @ v, np.inf)
    return w, float(numerator / scale)

def main():
    p = argparse.ArgumentParser()
    p.add_argument("--case", choices=["ads", "ds-limit"], required=True)
    args = p.parse_args()
    meta = {"python": platform.python_version(), "numpy": np.__version__,
            "scipy": scipy.__version__, "case": args.case}
    if args.case == "ds-limit":
        ell, muL, n = 1, 0.0, 0
        hp = 1.5 + np.sqrt(2.25 - muL*muL)
        hm = 1.5 - np.sqrt(2.25 - muL*muL)
        meta.update({"ell": ell, "muL": muL, "n": n,
                     "test_kind": "analytic-formula-smoke-check",
                     "operator_solved": False,
                     "horizon_orientation_checked": False,
                     "omegaL": [str(-1j*(2*n + ell + hm)),
                                str(-1j*(2*n + ell + hp))]})
    else:
        rows, previous = [], None
        for n in [16, 20, 24, 28, 32, 40]:
            w, residual = ads_mode(n)
            rows.append({"N": n, "omega": [w.real, w.imag],
                         "change": None if previous is None else abs(w - previous),
                         "residual": residual})
            previous = w
        meta.update({"rh_over_L": 1.0, "ell": 0,
                     "boundary": "fast scalar falloff", "rows": rows})
    print(json.dumps(meta, indent=2))

if __name__ == "__main__":
    main()

Independent shooting program

Save this second block as sads_shooting.py. The default command python3 sads_shooting.py runs the refined case; add —sweep for all eleven cases, or —output followed by a new JSON filename to retain results. The program refuses to overwrite that file. It requires NumPy and SciPy, not the collocation program.

"""Massless ell=0 SAdS shooting, rh=L=1; e^{-i omega t}.

Requires NumPy and SciPy. Default: one refined run. --sweep: eleven finite
convergence cases. --output creates a NEW JSON file, never overwriting one.
y=e^{i omega r*}u is normalizable at infinity and ANALYTIC at the horizon.
The endpoint row alone also admits a C1 outgoing branch at this frequency.
"""
import argparse
import json
from pathlib import Path
import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import root

SEED = 2.8-2.7j
# name, outer radius, horizon offset, boundary/horizon orders, rtol, atol
CASES = [
    ('coarse local expansions',8.,.05,6,4,1e-9,1e-12),
    ('boundary order 8',10.,.02,8,6,1e-11,1e-14),
    ('boundary order 12',10.,.02,12,6,1e-11,1e-14),
    ('boundary order 16',10.,.02,16,6,1e-11,1e-14),
    ('boundary radius 20',20.,.02,16,6,1e-11,1e-14),
    ('boundary radius 40',40.,.02,16,6,1e-11,1e-14),
    ('horizon offset 0.005',20.,.005,16,8,1e-11,1e-14),
    ('horizon offset 0.001',20.,.001,16,8,1e-11,1e-14),
    ('horizon offset 0.0002',20.,.0002,16,8,1e-11,1e-14),
    ('looser integration',20.,.005,16,8,1e-9,1e-12),
    ('tighter integration',20.,.005,16,8,1e-12,1e-15),
]


def boundary_state(w,r,order):
    # y=r^-2 sum b_n r^-n, b0=1, b1=-i omega.
    b = [1.+0j]
    for n in range(1,order+1):
        value = 2j*w*(n+1)*b[n-1]
        if n>=2: value += n*(n+1)*b[n-2]
        if n>=3: value -= 2*n*n*b[n-3]
        b.append(-value/(n*(n+3)))
    powers = -2-np.arange(order+1)
    terms = np.asarray(b)*r**powers
    return np.array([np.sum(terms),np.sum(powers*terms)/r])


def horizon_state(w,z,rho,order):
    # r^3 times the radial ODE: P y''+Q y'+R y=0, z=r-1.
    P = [0,4,11,11,5,1]
    Q = [4-2j*w,10-6j*w,12-6j*w,8-2j*w,2]
    R = [-4,-6,-6,-2]
    a = [1.+0j]
    for n in range(1,order+1):
        value = 0j
        for m,am in enumerate(a):
            k = n+1-m
            if k<len(P): value += P[k]*(rho+m)*(rho+m-1)*am
            k = n-m
            if k<len(Q): value += Q[k]*(rho+m)*am
            k = n-1-m
            if k<len(R): value += R[k]*am
        denominator = (rho+n)*(4*(rho+n)-2j*w)
        if abs(denominator)<1e-12: raise ValueError('resonant local exponent')
        a.append(-value/denominator)
    powers = rho+np.arange(order+1)
    terms = np.asarray(a)*np.exp(powers*np.log(z))  # z>0, real logarithm
    return np.array([np.sum(terms),np.sum(powers*terms)/z])


def projection(w,R,z,nb,nh,rtol,atol):
    def rhs(r,state):
        f,fp = 1+r*r-2/r,2*r+2/r**2
        y,dy = state
        return [dy,((fp/r)*y-(fp-2j*w)*dy)/f]
    result = solve_ivp(rhs,(R,1+z),boundary_state(w,R,nb),
                       method='DOP853',rtol=rtol,atol=atol)
    if not result.success: raise RuntimeError(result.message)
    basis = np.column_stack((horizon_state(w,z,0,nh),
                             horizon_state(w,z,.5j*w,nh)))
    A,B = np.linalg.solve(basis,result.y[:,-1])
    if abs(A)<1e-12: raise RuntimeError('ingoing amplitude too small')
    return B/A,float(np.linalg.cond(basis))


def shoot(case):
    label,*parameters = case
    def objective(value):
        ratio,_ = projection(complex(*value),*parameters)
        return [ratio.real,ratio.imag]
    result = root(objective,[SEED.real,SEED.imag],method='hybr',
                  options={'xtol':1e-9})
    w = complex(*result.x)
    ratio,condition = projection(w,*parameters)
    if not result.success or abs(ratio)>1e-7 or abs(w-SEED)>.5:
        raise RuntimeError('targeted shooting root did not converge')
    return {'case':label,'parameters':parameters,'omega':[w.real,w.imag],
            'outgoing_to_ingoing_abs':float(abs(ratio)),
            'basis_condition_2':condition}


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--sweep',action='store_true')
    parser.add_argument('--output',type=Path)
    args = parser.parse_args()
    if args.output and (args.output.exists() or not args.output.parent.is_dir()):
        parser.error('--output must be a new file in an existing directory')
    data = {'rh_over_L':1,'ell':0,'massless':True,
            'horizon_domain':'analytic ingoing, extracted outgoing coefficient zero',
            'rows':[shoot(case) for case in (CASES if args.sweep else CASES[-1:])]}
    payload = json.dumps(data,indent=2,sort_keys=True,allow_nan=False)+'\n'
    if args.output:
        with args.output.open('x',encoding='utf-8') as stream: stream.write(payload)
    else: print(payload,end='')


if __name__=='__main__':
    main()

Failure modes and acceptance tests

Endpoint falloff mixing is exposed by projecting the numerical solution onto both asymptotic branches. A wrong cosmological-horizon orientation is exposed by reconstructing the solution in regular $u$ or $v$ coordinates and checking the signed current. Branch swapping is exposed by eigenvector overlap and parameter continuation. Spectral pollution is exposed when a root fails to stabilize under $N$, precision, or mapping changes. Near extremality requires a map resolving the shrinking surface-gravity scale; simply increasing a global polynomial order can give deceptively small algebraic residuals.

For a new spectrum or parameter survey, seek frequency and eigenfunction convergence, differential and endpoint residuals, an explicit boundary domain, independent-method agreement, continuous mode labels and appropriate limit tests. These requirements are not all accomplished by this one example. Here the performed checks are the six-resolution eigenvalue sweep, complete pencil residual, original-radial-equation shooting with eleven endpoint/order/ tolerance variants, and comparison with the published four-decimal benchmark. The derivative matrix and full coordinate transformation were also checked independently. Eigenfunction convergence, higher precision, mapping changes, parameter continuation and the pure-AdS limit were not computed.

The de Sitter formula output supplies none of the differential-operator or boundary-domain checks and is not an accepted de Sitter spectrum. For a new production table, vary precision, map and extraction treatment and retain machine-readable results rather than copying displayed digits. The AdS boundary-condition page explains why changing the boundary domain changes the physical problem.

Exercises

Distinguish convergence from a small residual

The table has a residual of $1.3\times10^{-14}$ already at $N=16$, while the frequency changes by $3.4\times10^{-5}$ at the next resolution. Explain why the first number does not certify the eigenvalue to thirteen digits, and give two tests that do measure the relevant error.

Solution

The residual says only that the computed vector solves the truncated matrix pencil accurately. It does not bound the difference between that pencil and the differential operator, nor rule out a well-resolved spurious branch. Resolution variation directly samples truncation error: changes fall from $10^{-5}$ to about $10^{-11}$. Agreement with an independently formulated shooting or series method tests the operator and endpoint implementation. Eigenfunction-tail decay, mapping variation, and higher precision provide additional discriminants.

Reject a false horizon test

Suppose a numerical solution contains an unwanted local contribution $Bz^\rho$, where $B\ne0$ and $\rho=1.335603+1.399112i$. Does checking that its field and first derivative vanish as $z\to0^+$ exclude an outgoing wave? What additional check is needed?

Solution

No. The magnitudes of this contribution and its first derivative scale as $z^{1.335603}$ and $z^{0.335603}$, so both vanish. Its second derivative has magnitude proportional to $z^{-0.664397}$ and is singular. Consequently the branch can pass a first-order endpoint-row test while failing the analytic future-ingoing domain. Resolve the two local Frobenius branches at a finite positive offset and check that the coefficient $B$ vanishes relative to the nonzero ingoing coefficient, with offset and series-order convergence. Repeatedly reducing the endpoint residual without this branch test does not establish the physical boundary condition.

References

• Berti, E., Cardoso, V., and Starinets, A. O. (2009). “Quasinormal modes of black holes and black branes.” Classical and Quantum Gravity 26, 163001. DOI. Open PDF, arXiv v2. Page locators refer to this edition.
• Horowitz, G. T., and Hubeny, V. E. (2000). “Quasinormal modes of AdS black holes and the approach to thermal equilibrium.” Physical Review D 62, 024027. DOI. Open PDF, arXiv v2. Page locators refer to this edition.
• SciPy community (2024). SciPy 1.13.1 Reference Guide: solve_ivp.