Stellar-Dynamical Hardening

A massive black-hole binary hardens by giving orbital energy and angular momentum to stars that pass close enough for three-body slingshots. In the ideal full-loss-cone limit, the inverse semimajor axis grows nearly linearly with time. Whether a real nucleus sustains that supply depends on geometry, relaxation, centrophilic orbits, and the binary’s own depletion of phase space. The central calculation is therefore both local and global: calibrate one encounter, then demonstrate how the host refills encounters.

Helpful background. Black-Hole Geometry supplies relativistic length scales; Relativistic Matter and Compact Objects supplies the nuclear components; Black-Hole Astrophysics supplies host and population context; Two-Body Dynamics supplies Keplerian energy and angular momentum; Gravitational-Wave Astronomy supplies the later inspiral.

The hard-binary scale

Let $M_1\ge M_2$, $M=M_1+M_2$, and $\mu=M_1M_2/M$. The binary binding energy is

$$E_b=-\frac{GM_1M_2}{2a}.$$

A common convention defines the hard-binary separation by

$$a_h=\frac{G\mu}{4\sigma^2},$$

where $\sigma$ is the surrounding one-dimensional stellar velocity dispersion. Other conventions replace $\mu$ by $M_2$ or change the numerical coefficient; any quoted $a_h$ must state its definition. For $a\lesssim a_h$, an interacting star typically leaves with positive energy and the binary becomes more tightly bound.

Scattering experiments compress many encounters into the dimensionless hardening coefficient $H$:

$$s\equiv\frac{d}{dt}\left(\frac1a\right) \simeq H\frac{G\rho}{\sigma}.$$

For a hard, roughly equal-mass circular binary in an isotropic full loss cone, $H$ is commonly of order $10$–$20$, but it depends on mass ratio, eccentricity, stellar speed, and how encounters are sampled. Quinlan 1996, §§ 2–4 derives these scattering coefficients and separates energy hardening from eccentricity evolution.

The associated energy transfer is explicit:

$$\frac{d|E_b|}{dt} =\frac{GM_1M_2}{2}\,s.$$

Stars carry this energy outward or escape, creating a mass deficit and changing the density and anisotropy that set later encounters. The background cannot remain an infinite fixed reservoir.

Loss-cone supply

At a given stellar energy, the loss cone contains angular momenta whose periapses intersect the binary. In a spherical collisionless nucleus, those orbits are depleted and refilled only by relaxation; the measured hardening rate can then fall as particle number increases. In a triaxial remnant, centrophilic orbit families can transport stars collisionlessly. Rotation, a third black hole, massive perturbers, gas, and repeated galaxy mergers provide other supply routes.

This distinction explains why a constant $H$ is not a universal galaxy-scale solution. It is the full-loss-cone kernel. A global calculation must measure the actual flux into interaction phase space, the energy and angular momentum carried by ejected stars, and the change in the host distribution. Merritt and Milosavljević 2005, §§ 3–5 connects scattering coefficients to depletion, core scouring, and the final-parsec question.

Eccentricity is not a passive parameter. Individual encounters can increase or decrease it depending on mass ratio and stellar distribution. Because gravitational radiation is much faster at high $e$, a hardening model that evolves only $a$ can misestimate the transition to radiation reaction by orders of magnitude.

Numerical scattering experiment

A reproducible experiment has two layers. First, integrate single stars past a binary with declared $q$, $a$, $e$, incoming speed, impact-parameter distribution, orbital phase, and orientation. Record the change in binary energy and angular momentum and the star’s asymptotic state. Sample until $H$ and the eccentricity coefficient have stable Monte Carlo errors.

Second, embed the binary in a self-consistent nucleus. Pin the density profile, velocity anisotropy, shape and rotation, black-hole masses, particle mass, force law, softening, regularized region, timestep, and escape boundary. The public ketju–GADGET-4 implementation is one suitable environment because black-hole interactions inside a regularized region are unsoftened. Mannerkoski et al. 2023, §§ 2–3 reports hardening convergence once the black-hole-to-stellar-particle mass ratio is sufficiently large and quantifies integrator-tolerance tests.

Refine particle number, stellar softening, regularized-region radius, timestep tolerance, encounter impact-parameter range, and integration duration. Compare a deliberately replenished full-loss-cone benchmark with an isolated spherical nucleus and a triaxial remnant. Accept the astrophysical result only if $1/a(t)$, $e(t)$, stellar energy gain, loss-cone flux, and global energy residual converge over a common interval.

The principal failure modes are artificial two-body relaxation, Brownian wandering of an over-light binary, discontinuous softened/unsoftened forces, an ejection boundary that removes bound stars, confusing empty and full loss cones, and measuring a slope before the initial transient ends.

Reduced hardening benchmark

For constant $s$, the exact solution is

$$a(t)=\frac{a_0}{1+s a_0t}.$$

The teaching script below integrates $da/dt=-sa^2$ with explicit Euler steps for $a_0=1$, $s=0.8$, and $t=2$. The exact answer is $a=0.384615$. The decreasing error verifies the reduced integrator, not a stellar nucleus.

Step count $N$	Semimajor-axis result	Absolute error
50	0.380042	0.004573
100	0.382341	0.002274
200	0.383482	0.001134
400	0.384049	0.000566

Use Python 3.12.5 and NumPy==2.1.3. Save the following program as stellar_hardening.py. The script writes stellar_hardening.csv.

from pathlib import Path
import numpy as np

a0, hardening_rate, final_time = 1.0, 0.8, 2.0
exact = a0 / (1.0 + hardening_rate * a0 * final_time)
rows = []

for steps in (50, 100, 200, 400):
    semimajor = a0
    dt = final_time / steps
    for _ in range(steps):
        semimajor -= dt * hardening_rate * semimajor**2
    error = abs(semimajor - exact)
    rows.append((steps, semimajor, error))
    print(f"{steps:4d}  a={semimajor:.6f}  error={error:.6f}")

np.savetxt(
    Path("stellar_hardening.csv"),
    np.asarray(rows, dtype=float),
    delimiter=",",
    header="steps,semimajor_axis,absolute_error",
    comments="",
    fmt=["%d", "%.6f", "%.6f"],
)

Create an isolated environment and rerun with:

python3.12 -m venv .venv-stellar-hardening
.venv-stellar-hardening/bin/python -m pip install numpy==2.1.3
.venv-stellar-hardening/bin/python stellar_hardening.py

The script produces the table above. To analyze an N-body run, replace the constant-rate evolution with a fit to $1/a(t)$ over a preregistered post-transient window, use block resampling for correlated fluctuations, and close the slope against the stars’ measured energy gain.

Interpretation and boundaries

A nearly constant $s$ in a simulation does not show that every galaxy avoids stalling. It establishes a rate for the modeled density, geometry, and loss-cone supply. Conversely, declining hardening in a low-$N$ spherical model may reflect artificial relaxation or binary wandering rather than the collisionless limit. The convergence trend with $N$ and the orbital structure of the remnant decide which interpretation is credible.

Here the result is $a(t)$, $e(t)$, and a stellar budget for the assumed environment. Two-Body Dynamics develops the general encounter mechanics, and Gravitational-Wave Astronomy explains inference from the subsequent gravitational-wave signal.

Estimate a hardening interval

Take $M_1=10^8M_\odot$, $M_2=3\times10^7M_\odot$, $\sigma=200\,{\rm km\,s^{-1}}$, $\rho=10^4M_\odot\,{\rm pc^{-3}}$, and $H=15$. Compute $a_h$ with the convention above. Assuming the full-loss-cone rate stays constant, how long does the binary take to shrink from $a_h$ to $0.10\,{\rm pc}$?

Solution

The reduced mass is

$$\mu=\frac{M_1M_2}{M_1+M_2} \simeq2.31\times10^7M_\odot.$$

Using $G=4.3009\times10^{-3}\,{\rm pc}\,({\rm km\,s^{-1}})^2M_\odot^{-1}$ gives

$$a_h=\frac{G\mu}{4\sigma^2} \simeq0.620\,{\rm pc}.$$

The hardening rate is

$$s=H\frac{G\rho}{\sigma} \simeq3.23\,\frac{{\rm km\,s^{-1}}}{{\rm pc}^2} \simeq3.30\,{\rm pc^{-1}\,Myr^{-1}}.$$

Therefore

$$\Delta t =\frac{1/a_f-1/a_h}{s} =\frac{10.0-1.61}{3.30}\,{\rm Myr} \simeq2.54\,{\rm Myr}.$$

This short value is conditional on a dense, permanently full loss cone. Density depletion, anisotropy, shape evolution, or an empty loss cone changes the result.

References

• Mannerkoski, Matias, et al. “ketju—Resolving Small-scale Supermassive Black Hole Dynamics in GADGET-4.” Monthly Notices of the Royal Astronomical Society 524 (2023): 4062–4082. DOI.
• Merritt, David, and Miloš Milosavljević. “Massive Black Hole Binary Evolution.” Living Reviews in Relativity 8 (2005): 8. DOI.
• Quinlan, Gerald D. “The Dynamical Evolution of Massive Black Hole Binaries. I. Hardening in a Fixed Stellar Background.” New Astronomy 1 (1996): 35–56. DOI.