Constructing Alternative-Theory Solution Branches

A compact-object solution is not specified by differential equations alone. One must also choose a regular center or horizon, admissible data at infinity, a gauge, a matter closure when a star is present, and a rule for following one connected branch rather than jumping to another. This page builds that boundary-value problem and tests its essential pieces on a horizon-to-infinity scalar–Gauss–Bonnet calculation. It establishes a construction method, not the stability or empirical viability of the resulting branch.

Helpful background. Einstein Equations provide the field equations and constraint identities used below; Gravitational Actions and Boundary Terms distinguish fixed boundary data from derived conditions; Formulations of General Relativity organize gauge fixing and elliptic reduction; Black-Hole Geometry and Mechanics supplies horizon regularity and invariant charges; Black-Hole Perturbation Theory supplies mode and asymptotic boundary conditions; and Relativistic Matter and Compact Objects supplies stellar equations of state and surface matching.

A solution branch is a boundary-value problem

Let a diffeomorphism-invariant theory have fields $\Psi=(g_{ab},\phi^A,A_a^I,\ldots)$ and equations $E[\Psi;\lambda]=0$, where $\lambda$ denotes couplings. After a stationary symmetry reduction and gauge choice, a compact solution is a root of

$$E[u;\lambda]=0, \qquad B_{\rm in}[u;\lambda]=0, \qquad B_{\infty}[u;\lambda]=0.$$

Here $u$ is the set of reduced fields. At a nonextremal Killing horizon, smoothness is most transparently imposed in coordinates regular across the future horizon, or by a Euclidean polar expansion when a bifurcation surface exists. At a stellar center, parity gives even scalar and fluid variables and an enclosed mass $m(r)=O(r^3)$. At infinity, the admissible falloffs depend on whether the solution is asymptotically flat, Kaluza–Klein, de Sitter, or anti-de Sitter; those alternatives are different problems, not interchangeable numerical regulators.

For a second-order system, only the defining boundary data should be imposed independently. Remaining coefficients in the local series follow from the equations. Over-prescribing both a coefficient and the relation that determines it can make a consistent solution appear nonexistent. The distinction between defining and derived conditions, together with horizon and asymptotic examples, is developed in Dias, Santos, and Way 2016, § V, pp. 51–66.

Inner boundary	Defining data	Derived checks	Outer data
Nonextremal horizon	horizon scale, angular velocity or potentials, regular field branch	finite invariants, common surface gravity, vanishing forbidden flux	asymptotic frame, sources, normalizable modes
Extremal horizon	near-horizon throat data and admissible end behavior	compatibility with the degenerate equations	charges and asymptotic moduli
Stellar center	central pressure or density and regular field amplitudes	$m=O(r^3)$ and vanishing odd radial derivatives	zero pressure at the surface, exterior matching, asymptotic fields

Gauge fixing, charges, and branch labels

The reduced equations must have an invertible principal part after gauge freedom is removed. In the Einstein–DeTurck method one solves

$$R_{ab}-\nabla_{(a}\xi_{b)}=8\pi G\left(T_{ab}-\frac{1}{2}Tg_{ab}\right), \qquad \xi^a=g^{bc}\left(\Gamma^a{}_{bc}-\bar\Gamma^a{}_{bc}\right),$$

and subsequently verifies $\xi^a=0$; a small Newton update is not that verification. The reference metric must share the intended causal structure and boundary conditions. Ellipticity, the possibility of Ricci solitons, and the required check of the DeTurck vector are stated in Dias, Santos, and Way 2016, § IV.B.2, pp. 41–42.

A branch label should contain enough invariant information to survive a coordinate change: horizon topology, node number, conserved charges, asymptotic moduli, and a continuation parameter such as horizon radius or central density. Pseudo-arclength continuation is preferable near a fold, where $dQ/d\lambda$ can diverge even though the curve itself is smooth. The general-relativistic branch is recovered by taking the declared coupling to zero while holding physical boundary data fixed. A singular rescaling that changes the cutoff or the matter metric is not a GR limit.

A horizon-to-infinity verification problem

For this benchmark use geometrized units $G=c=1$: $M=GM_{\rm phys}$ is a length, $\phi$ is dimensionless, and $\alpha$ has dimensions of length squared. Consider the decoupling-limit scalar equation on Schwarzschild with $M\gt0$,

$$\square\phi+\alpha\mathcal G=0, \qquad f(r)=1-\frac{2M}{r}, \qquad \mathcal G=\frac{48M^2}{r^6}.$$

For a static spherical field this becomes

$$\frac{d}{dr}\left[r^2f(r)\phi'(r)\right] =-\frac{48\alpha M^2}{r^4}.$$

Writing $p=r^2f\phi'$, regularity at $r_h=2M$ requires $p(r_h)=0$; otherwise $\phi$ contains a logarithm at the horizon. With $\phi(\infty)=0$, direct integration gives

$$\phi_{\rm exact}(r)=\frac{2\alpha}{M} \left(\frac{1}{r}+\frac{M}{r^2}+\frac{4M^2}{3r^3}\right), \qquad D=\lim_{r\to\infty}r\phi=\frac{2\alpha}{M}.$$

This is a stringent but bounded test: the horizon is a regular singular point, infinity must be truncated and controlled, the scalar charge has two extractions, and $\alpha=0$ must return $\phi=0$. The full backreacted problem has additional metric equations and a finite-size restriction; near-horizon regularity fixes the scalar charge rather than leaving it independent in the explicit branch of Sotiriou and Zhou 2014, § III.B–D, pp. 5–7.

Reproducible collocation calculation

Restoring $G$ gives $D=2\alpha/(GM_{\rm phys})$. The driver below uses $\alpha/M^2=1$ to test a linear equation at a convenient amplitude; it does not claim a controlled backreacted expansion at that coupling. A physical perturbative prediction additionally requires $\lvert\alpha\rvert/M^2\ll1$ and the metric correction at the retained order.

The retained calculation uses CPython 3.9.6, NumPy 2.0.2, and SciPy 1.13.1, IEEE-754 binary64 arithmetic, solve_bvp, and tolerances shown in the driver. It starts at $r_h+10^{-5}M$, imposes the exact regular flux there, and uses the leading asymptotic Robin condition at $R$. Save the Python fence as branch_bvp.py. The script must produce the refinement table and exit nonzero if charge convergence or the GR-branch benchmark fails.

import sys
import numpy as np
from scipy.integrate import solve_bvp

M, ALPHA = 1.0, 1.0
RH, EPS = 2.0 * M, 1.0e-5
EXACT_CHARGE = 2.0 * ALPHA / M

def solve_case(n, rmax, tol, alpha=ALPHA):
    r0 = RH + EPS
    grid = np.linspace(r0, rmax, n)
    p0 = 16.0 * alpha * M**2 * (r0**-3 - RH**-3)

    def equations(r, y):
        phi, flux = y
        return np.vstack((flux / (r * (r - RH)),
                          -48.0 * alpha * M**2 / r**4))

    def boundary(ya, yb):
        return np.array((ya[1] - p0,
                         yb[0] + yb[1] / (rmax - RH)))

    guess = np.zeros((2, grid.size))
    sol = solve_bvp(equations, boundary, grid, guess,
                    tol=tol, max_nodes=32*n)
    if not sol.success:
        raise RuntimeError(sol.message)

    sample = np.linspace(0.65*rmax, rmax, 200)
    phi, flux = sol.sol(sample)
    design = np.column_stack((np.ones_like(sample), 1/sample,
                              1/sample**2, 1/sample**3))
    charge_fit = np.linalg.lstsq(design, phi, rcond=None)[0][1]
    charge_flux = -sol.y[1, -1]
    mismatch = np.max(np.abs(boundary(sol.y[:, 0], sol.y[:, -1])))
    residual = np.max(sol.rms_residuals)
    return (sol.x.size, charge_flux, charge_fit, residual,
            mismatch, np.max(np.abs(sol.y[0])))

def main():
    cases = ((64, 20.0, 1e-6), (128, 40.0, 1e-8),
             (256, 80.0, 1e-10))
    print("N R charge_flux charge_fit fit_error residual mismatch")
    last_error = None
    for n, rmax, tol in cases:
        nodes, qf, qa, residual, mismatch, _ = solve_case(n, rmax, tol)
        error = abs(qa - EXACT_CHARGE)
        print(f"{nodes:d} {rmax:.0f} {qf:.12f} {qa:.12f} "
              f"{error:.3e} {residual:.3e} {mismatch:.3e}")
        if last_error is not None and error >= last_error:
            raise SystemExit("asymptotic charge did not converge")
        last_error = error
    _, _, _, _, _, gr_amp = solve_case(64, 20.0, 1e-8, alpha=0.0)
    print(f"GR_BRANCH_MAX_ABS_PHI={gr_amp:.3e}")
    if gr_amp > 1e-12 or last_error > 1e-7:
        raise SystemExit("independent benchmark failed")

if __name__ == "__main__":
    main()

Create an isolated environment and rerun exactly as follows:

python3 -m venv .venv-branch
. .venv-branch/bin/activate
python3 -m pip install numpy==2.0.2 scipy==1.13.1
python3 branch_bvp.py

Convergence and independent checks

The following output was reproduced with the pinned environment. N is the final adaptive node count. The flux charge retains the expected $O(R^{-3})$ finite-boundary error, whereas the four-term asymptotic fit isolates the $1/r$ coefficient.

N	Outer radius $R/M$	Charge from flux	Charge from fit	Fit error	ODE residual	Boundary mismatch
142	20	1.998000031971	1.999999310526	$6.895\times10^{-7}$	$9.702\times10^{-7}$	$0.000$
621	40	1.999750000182	1.999999969452	$3.055\times10^{-8}$	$9.996\times10^{-9}$	$6.939\times10^{-18}$
2892	80	1.999968750001	1.999999997829	$2.171\times10^{-9}$	$9.999\times10^{-11}$	$0.000$

Three independent facts agree: direct quadrature gives $D=2$, asymptotic fitting approaches that value, and setting $\alpha=0$ returns GR_BRANCH_MAX_ABS_PHI=0.000e+00. This separates equation residual, boundary mismatch, and observable error; none can substitute for the other.

Interpretation and failure modes

The fit includes a constant because the finite-radius Robin condition shifts the field’s zero by a small amount. The exact solution of this equation is a constant plus the displayed cubic polynomial in $1/r$, so there is no omitted $1/r^4$ term for the fit to remove. The remaining fitted-charge error is numerical. Radius, initial mesh, and tolerance all change between rows; their joint improvement is a consistency check, not a measurement of one method’s convergence order.

The benchmark verifies the mechanics of a regular horizon-to-infinity solve, not a nonlinear hairy black hole. A production branch calculation must additionally check all metric and matter equations, constraints, curvature invariants, charges from independent surfaces, a Smarr or first-law relation when applicable, and spectral convergence after domain decomposition. It must also show that curvature and field gradients remain below the effective-theory cutoff.

Branch jumping. Newton iteration finds a nearby root, not every root. A node-count change, discontinuous charge, or failed tangent predictor signals that continuation may have crossed branches.

Singular inner data. A finite coordinate field can still produce a divergent freely falling invariant. Expand in a horizon-regular frame and reject a root that violates the regularity discriminant.

False convergence. A small update norm can coexist with a large equation residual or DeTurck vector. Refine the domain, tolerance, and outer boundary separately, and compare invariant outputs rather than coordinate samples alone.

Changed physical problem. Holding a dimensionless coupling fixed while rescaling the horizon can change the dimensionful theory. Record which coupling, asymptotic modulus, equation of state, and cutoff are held fixed along the branch.

Empirical likelihoods and numerical bounds are beyond this construction problem and belong to Tests of Gravity. A mathematically regular branch is only a theory prediction until stability, formation, environmental, waveform, and detector-response questions are supplied.

Exercises

For the last two rows of the table, compute the fitted-error reduction factor. Explain why it is not an independently measured convergence order, why the flux error and fit error differ, and whether the finest row proves a new nonlinear branch.

Solution

The fitted error decreases by about $14.1$, whose base-two logarithm is $3.81$. This is not a convergence order in $R$: the collocation tolerance and mesh also changed, and the exact scalar already lies in the fit’s $\{1,1/r,1/r^2,1/r^3\}$ basis. By contrast, the exact finite-radius flux is $-p(R)=2-16/R^3$ at $M=\alpha=1$, so its extraction error is analytically $16/R^3$. To measure numerical convergence, vary one control while keeping the others sufficiently accurate. Neither result proves a nonlinear branch: the metric is fixed to Schwarzschild. A branch claim requires the coupled equations, invariant horizon checks, continuation, and stability.

References

• Óscar J. C. Dias, Jorge E. Santos, and Benson Way, “Numerical Methods for Finding Stationary Gravitational Solutions,” Classical and Quantum Gravity 33 (2016) 133001. DOI. Open PDF.
• Thomas P. Sotiriou and Shuang-Yong Zhou, “Black Hole Hair in Generalized Scalar-Tensor Gravity: An Explicit Example,” Physical Review D 90 (2014) 124063. DOI. Open PDF.