Higher-Curvature Black-Hole Branches

A corrected black-hole solution is physically connected to general relativity only when its fields, charges, horizon data, and residuals approach a controlled GR solution on the same branch. A nonlinear solver can also converge to a disconnected branch or step past a fold. This page combines the full backreacted construction strategy with an executable scalar–Gauss–Bonnet decoupling benchmark: the latter starts at Schwarzschild, follows scalar hair through a turn in the coupling, and checks the asymptotic charge independently. It validates the continuation logic but does not substitute for solving the backreacted metric equations.

Required background. Gravitational Actions and Boundary Terms supplies horizon and asymptotic variational data, while Curvature supplies the Gauss–Bonnet invariant evaluated below.

Helpful background. Einstein Equations supplies the GR seed; Formulations of General Relativity supplies boundary-value and constraint methods; Quantum Gravity supplies the cutoff interpretation of the higher-curvature coupling.

Branches, seeds, and boundary data

For a scalar–Gauss–Bonnet example, fix the normalization

$$S=\frac{1}{16\pi G}\int d^4x\sqrt{-g}\left[ R-\frac12(\nabla\phi)^2+\ell^2 f(\phi)\mathcal G \right].$$

A full static spherical calculation solves coupled metric and scalar equations, constructs a regular horizon series, imposes asymptotic flatness, and extracts the ADM mass and scalar charge. If $f_{,\phi}(0)=0$, Schwarzschild with $\phi=0$ can lose stability at a scalar zero mode; the mode locates a bifurcation but not the nonlinear branch by itself. The onset equation, horizon discriminant, and asymptotic data are exhibited in Silva et al. 2018, eqs. (1)–(15).

Continuation begins with a verified seed, not an arbitrary initial guess. At each step one should record the field-equation residual, constraint residual, horizon regularity condition, asymptotic charges, and—when the full equations are solved—a Smarr or first-law check. Near a turning point in the chosen coupling, amplitude or pseudo-arclength continuation replaces naive parameter stepping. The distinction between update norm and equation residual, together with spectral and Newton convergence tests, is developed in Dias, Santos, and Way 2016, § VI, pp. 66–71; § VII.A–B, pp. 71–73; Appendix A, pp. 101–111.

Executable decoupling benchmark

The retained reference calculation measures every length in the geometrized mass $M=GM_{\mathrm{phys}}$ of the Schwarzschild seed. Thus $r$ below means $r_{\mathrm{physical}}/M$, $Q$ means $Q_{\mathrm{physical}}/M$, and the dimensionless metric functions are

$$F(r)=1-\frac{2}{r}, \qquad \mathcal G_{\rm Schw}=\frac{48}{r^6},$$

and chooses $f(\phi)=\phi+\phi^4/4$. With $\lambda=\ell^2/M^2$, the static scalar equation is

$$\frac{d}{dr}\left[r(r-2)\phi'\right] +\frac{48\lambda}{r^4}(1+\phi^3)=0.$$

Writing $q=r(r-2)\phi'$, regularity and asymptotic flatness require

$$\phi'_h=-\frac32\lambda(1+\phi_h^3), \qquad \phi(r)=\frac{Q}{r}+O(r^{-2}).$$

The horizon amplitude $\phi_h$ is the continuation parameter. This crosses a turn in $\lambda(\phi_h)$ without a singular Newton step. The positive branch approaches $\phi=0$, $Q=0$, and Schwarzschild continuously as $\lambda\to0$. Because the metric is fixed, $r_h=2$, $M=1$, and the horizon area remain unchanged; the computed horizon datum and scalar charge are the nontrivial branch observables. Backreaction, Wald entropy, and a metric first law are deliberately not claimed.

The reference environment is CPython 3.9.6, NumPy 2.0.2, and SciPy 1.13.1, using IEEE-754 binary64 arithmetic. DOP853 integrates from $r=2+10^{-6}$ to $r=80$; its absolute tolerance is $10^{-2}$ times the stated relative tolerance. Amplitude continuation uses steps of $0.05$. The outer condition retains the asymptotic series through $r^{-4}$.

Rerun. Save the Python source below as /tmp/higher_curvature_branch.py. The script must emit the branch and refinement tables and exit nonzero if the turn, convergence trend, flux check, or independent weak-coupling benchmark fails.

import platform
import sys

import numpy as np
import scipy
from scipy.integrate import simpson, solve_ivp
from scipy.optimize import root_scalar

if platform.python_version() != "3.9.6":
    raise RuntimeError("reference run requires CPython 3.9.6")
if np.__version__ != "2.0.2" or scipy.__version__ != "1.13.1":
    raise RuntimeError("reference run requires NumPy 2.0.2 and SciPy 1.13.1")

RH, DELTA, RMAX = 2.0, 1.0e-6, 80.0


def coupling_prime(phi):
    return 1.0 + phi**3


def integrate(lam, amp, rtol, dense=False):
    phip_h = -1.5 * lam * coupling_prime(amp)
    r0 = RH + DELTA
    y0 = np.array([amp + phip_h * DELTA,
                   r0 * (r0 - RH) * phip_h])

    def rhs(r, y):
        phi, q = y
        return np.array([q / (r * (r - RH)),
                         -48.0 * lam * coupling_prime(phi) / r**4])

    sol = solve_ivp(rhs, (r0, RMAX), y0, method="DOP853",
                    rtol=rtol, atol=rtol * 1.0e-2,
                    dense_output=dense)
    if not sol.success:
        raise RuntimeError(sol.message)
    return sol


def outer_mismatch(lam, amp, rtol):
    phi, q = integrate(lam, amp, rtol).y[:, -1]
    charge = -q + 16.0 * lam / RMAX**3
    asymptotic_phi = (charge / RMAX + charge / RMAX**2
                      + 4.0 * charge / (3.0 * RMAX**3)
                      + (2.0 * charge - 4.0 * lam) / RMAX**4)
    return phi - asymptotic_phi


def solve_lambda(amp, guess, rtol):
    x0 = max(1.0e-10, 0.96 * guess)
    x1 = max(2.0e-10, 1.04 * guess + 1.0e-5)
    root = root_scalar(lambda z: outer_mismatch(z, amp, rtol),
                       x0=x0, x1=x1, method="secant",
                       xtol=max(1.0e-13, rtol * 1.0e-2), maxiter=40)
    if not root.converged or root.root <= 0.0:
        raise RuntimeError("shooting root failed")
    return root.root


def diagnostics(amp, lam, rtol, ncheck):
    sol = integrate(lam, amp, rtol, dense=True)
    u = np.linspace(0.0, 1.0, ncheck)
    r = RH + DELTA + (RMAX - RH - DELTA) * u**2
    phi, q = sol.sol(r)
    ode_residual = (np.gradient(q, r, edge_order=2)
                    + 48.0 * lam * coupling_prime(phi) / r**4)
    residual = np.max(np.abs(ode_residual[3:-3]))
    charge = -q[-1] + 16.0 * lam / RMAX**3
    asymptotic_phi = (charge / RMAX + charge / RMAX**2
                      + 4.0 * charge / (3.0 * RMAX**3)
                      + (2.0 * charge - 4.0 * lam) / RMAX**4)
    mismatch = abs(phi[-1] - asymptotic_phi)
    flux_error = abs(q[-1] - q[0]
                     + simpson(48.0 * lam * coupling_prime(phi) / r**4,
                               x=r))
    return charge, residual, mismatch, flux_error


def main():
    global RMAX, DELTA
    amplitudes = np.linspace(0.05, 1.50, 30)
    branch = []
    guess = 6.0 * amplitudes[0] / 11.0
    for amp in amplitudes:
        lam = solve_lambda(amp, guess, 1.0e-11)
        charge, residual, mismatch, flux = diagnostics(
            amp, lam, 1.0e-11, 2401)
        branch.append((amp, lam, charge, residual, mismatch, flux))
        guess = lam

    peak = max(range(len(branch)), key=lambda i: branch[i][1])
    if peak == 0 or peak == len(branch) - 1:
        raise RuntimeError("branch turn was not bracketed")

    amp_ref = 1.20
    refinements = []
    for rtol, ncheck in [(1.0e-6, 301), (1.0e-8, 601),
                         (1.0e-10, 1201), (1.0e-12, 2401)]:
        lam = solve_lambda(amp_ref, 0.35, rtol)
        diag = diagnostics(amp_ref, lam, rtol, ncheck)
        refinements.append((rtol, ncheck, lam, *diag))

    amp_weak = 1.0e-4
    lam_weak = solve_lambda(amp_weak, 6.0 * amp_weak / 11.0, 1.0e-12)
    q_weak, _, _, _ = diagnostics(amp_weak, lam_weak, 1.0e-12, 4001)
    weak_lambda_error = abs(lam_weak / amp_weak - 6.0 / 11.0)
    weak_charge_error = abs(q_weak / lam_weak - 2.0)

    # Separate integration error from finite-domain and horizon-start error.
    boundary_rows = []
    for outer_radius in (40.0, 80.0, 160.0, 320.0):
        RMAX = outer_radius
        lam = solve_lambda(amp_ref, 0.35, 1.0e-12)
        charge, _, _, _ = diagnostics(amp_ref, lam, 1.0e-12, 2401)
        boundary_rows.append((RMAX, DELTA, lam, charge))
    horizon_rows = []
    RMAX = 160.0
    for horizon_offset in (1.0e-4, 1.0e-5, 1.0e-6, 1.0e-7):
        DELTA = horizon_offset
        lam = solve_lambda(amp_ref, 0.35, 1.0e-12)
        charge, _, _, _ = diagnostics(amp_ref, lam, 1.0e-12, 2401)
        horizon_rows.append((RMAX, DELTA, lam, charge))
    RMAX, DELTA = 80.0, 1.0e-6

    print("BRANCH amp lambda charge residual mismatch flux_error")
    for index in [0, 7, 15, peak, 23, 29]:
        print(" ".join(["BRANCH"] + [f"{v:.10e}" for v in branch[index]]))
    print("REFINE rtol N lambda charge residual mismatch flux_error")
    for row in refinements:
        print("REFINE " + " ".join(f"{v:.10e}" if isinstance(v, float)
                                   else str(v) for v in row))
    print(f"WEAK lambda_over_A={lam_weak/amp_weak:.10e} "
          f"charge_over_lambda={q_weak/lam_weak:.10e} "
          f"errors={weak_lambda_error:.3e},{weak_charge_error:.3e}")
    for label, rows in (("OUTER", boundary_rows), ("HORIZON", horizon_rows)):
        print(label + " Rmax delta lambda charge")
        for row in rows:
            print(label + " " + " ".join(f"{value:.12e}" for value in row))

    passed = (
        branch[peak - 1][1] < branch[peak][1]
        and branch[peak + 1][1] < branch[peak][1]
        and abs(refinements[-1][2] - refinements[-2][2]) < 2.0e-8
        and refinements[-1][4] < refinements[0][4]
        and refinements[-1][5] < 1.0e-9
        and refinements[-1][6] < refinements[0][6]
        and weak_lambda_error < 2.0e-4
        and weak_charge_error < 3.0e-2
        and max(abs(boundary_rows[-1][j] - boundary_rows[-2][j])
                for j in (2, 3)) < 1.0e-9
        and max(abs(horizon_rows[-1][j] - horizon_rows[-2][j])
                for j in (2, 3)) < 1.0e-9
    )
    if not passed:
        print("FAIL: convergence, turn, or independent benchmark",
              file=sys.stderr)
        return 1
    print("PASS: branch turn, refinement, and weak-coupling benchmark")
    return 0


if __name__ == "__main__":
    raise SystemExit(main())

Run it without hidden inputs:

python3 --version
python3 -c 'import numpy, scipy; print(numpy.__version__, scipy.__version__)'
python3 /tmp/higher_curvature_branch.py

Branch turn and convergence

The following values were recomputed with the retained source. The residual is the maximum independent finite-difference value of $dq/dr+48\lambda(1+\phi^3)/r^4$, excluding three points at each boundary. The flux error independently integrates the same conservation law with Simpson quadrature.

Horizon amplitude $\phi_h$	Coupling $\lambda$	Scalar charge $Q$
0.050	0.02727079161	0.05454424276
0.400	0.2105893049	0.4315771924
0.800	0.3420492322	0.8109528642
1.000	0.3588881542	0.9652019443
1.200	0.3509981386	1.098553976
1.500	0.3145952717	1.273331690

The corresponding numerical checks are:

Horizon amplitude $\phi_h$	Field residual	Outer mismatch	Flux error
0.050	$4.95\times10^{-7}$	$1.46\times10^{-17}$	$8.09\times10^{-12}$
0.400	$4.27\times10^{-6}$	$4.91\times10^{-16}$	$7.12\times10^{-11}$
0.800	$1.33\times10^{-5}$	$9.59\times10^{-14}$	$2.27\times10^{-10}$
1.000	$2.25\times10^{-5}$	$6.53\times10^{-14}$	$3.82\times10^{-10}$
1.200	$3.59\times10^{-5}$	$7.74\times10^{-14}$	$6.01\times10^{-10}$
1.500	$6.34\times10^{-5}$	$5.10\times10^{-13}$	$1.05\times10^{-9}$

$\lambda$ rises through $\phi_h=1.0$ and then falls, so the sampled branch brackets a turn between neighboring amplitude steps. Parameter stepping in $\lambda$ alone would become double-valued there; amplitude continuation crosses it without changing branch.

At the post-turn point $\phi_h=1.2$, four refinements give:

Solver rtol	Check grid $N$	Coupling $\lambda$	Charge $Q$
$10^{-6}$	301	0.3509981206	1.098553946
$10^{-8}$	601	0.3509981382	1.098553975
$10^{-10}$	1201	0.3509981386	1.098553976
$10^{-12}$	2401	0.3509981386	1.098553976

The checks for these same four runs are:

Solver rtol	Field residual	Outer mismatch	Flux error
$10^{-6}$	$2.29\times10^{-3}$	$1.15\times10^{-13}$	$2.48\times10^{-6}$
$10^{-8}$	$5.74\times10^{-4}$	$1.82\times10^{-14}$	$1.54\times10^{-7}$
$10^{-10}$	$1.44\times10^{-4}$	$9.19\times10^{-17}$	$9.62\times10^{-9}$
$10^{-12}$	$3.59\times10^{-5}$	$1.39\times10^{-14}$	$6.01\times10^{-10}$

The coupling and charge stabilize while the independent residual falls by approximately four when $N$ doubles, consistent with the second-order finite-difference diagnostic. The outer mismatch is already root-solver limited; it need not decrease monotonically after reaching roundoff.

Solver tolerance alone does not control the truncation of the exterior or the near-horizon series. The additional loops in the program independently vary those two choices at $\phi_h=1.2$ and rtol $10^{-12}$:

Outer radius $R_{\max}$	Horizon offset $\delta$	Coupling $\lambda$	Charge $Q$
40	$10^{-6}$	0.3509981434221	1.098553985710
80	$10^{-6}$	0.3509981385942	1.098553975974
160	$10^{-6}$	0.3509981384468	1.098553975662
320	$10^{-6}$	0.3509981384423	1.098553975652
160	$10^{-4}$	0.3509981222548	1.098553748537
160	$10^{-5}$	0.3509981383576	1.098553972881
160	$10^{-7}$	0.3509981384468	1.098553975696

Compared with the nominal $R_{\max}=80$, $\delta=10^{-6}$ run, extending the outer boundary to $R_{\max}=320$ changes either quantity by less than $4\times10^{-10}$. The run with $R_{\max}=160$, $\delta=10^{-7}$ also agrees with the nominal result within that difference. This supports eight decimal places for this fixed-metric boundary problem; it is not an error estimate for omitted metric backreaction. The larger $\delta=10^{-4}$ run demonstrates why the horizon-start error cannot be inferred from the much smaller outer shooting mismatch.

Independent weak-coupling check

As $\phi_h\to0$, $f_{,\phi}=1+O(\phi_h^3)$. The radial equation can then be integrated exactly:

$$q(r)=\frac{16\lambda}{r^3}-2\lambda, \qquad \phi(r)=2\lambda\left( \frac1r+\frac1{r^2}+\frac{4}{3r^3} \right).$$

Therefore

$$\frac{\lambda}{\phi_h}\longrightarrow\frac6{11}, \qquad \frac{Q}{\lambda}\longrightarrow2.$$

At $\phi_h=10^{-4}$ the independent run gives $\lambda/\phi_h=0.5454545443$ and $Q/\lambda=2.0000000000$, with errors $1.14\times10^{-9}$ and $3.89\times10^{-11}$. This checks the horizon sign, charge normalization, outer expansion, and branch selection without reusing the nonlinear target value.

What the calculation does not establish

The decoupling solution can fail as an approximation when scalar stress-energy backreacts appreciably, the dimensionless curvature and scalar-gradient ratios in the full EFT cease to be small, gradients become unresolved, or the full system loses hyperbolicity. A complete branch must also test the metric horizon discriminant, ADM charge, Wald entropy, and first law. Further numerical failure modes include converging to another shooting root, starting from a singular horizon series, finite-$R_{\max}$ contamination, stiffness, solver stagnation, and sensitivity to the amplitude step. None is diagnosed by a small Newton update alone.

Empirical likelihoods and numerical bounds remain in Tests of Gravity, while the ultraviolet completion of the higher-curvature action remains a separate Quantum Gravity question.

Exercises

Use the refinement table to decide which digits of $\lambda$ and $Q$ are supported. Then explain why the observed turn is evidence for a double-valued decoupling branch but not for a regular backreacted black-hole branch.

Solution

The final two solver refinements change $\lambda$ below $10^{-10}$ and $Q$ near $10^{-9}$, while both independent residual diagnostics continue to fall. This tests integration accuracy at fixed endpoints. The separate outer-radius and horizon-offset tests are also needed before supporting eight decimal places for the stated boundary problem. The rise and subsequent decrease of $\lambda$ bracket a turn in the sampled continuation curve; an exact fold location and its nondegeneracy would require a tangent or finer continuation calculation. Because the metric was held Schwarzschild, the calculation does not test metric backreaction, horizon discriminants, entropy, hyperbolicity, or a gravitational first law; those require the coupled system.

References

• Dias, Óscar J. C., Jorge E. Santos, and Benson Way. 2016. “Numerical Methods for Finding Stationary Gravitational Solutions.” Classical and Quantum Gravity 33: 133001. DOI.
• Silva, Hector O., Jeremy Sakstein, Leonardo Gualtieri, Thomas P. Sotiriou, and Emanuele Berti. 2018. “Spontaneous Scalarization of Black Holes and Compact Stars from a Gauss–Bonnet Coupling.” Physical Review Letters 120: 131104. DOI.