Massive Black-Hole Binary Formation

A common route to a massive black-hole binary follows a hierarchy of merging structures. The host galaxies first interpenetrate; each black hole sinks with its surviving nuclear envelope; the two nuclei dissolve into a common potential; and mutual binding eventually defines a binary. This page follows that galaxy-merger route. Its durable calculation is a scale-matched orbital-energy budget. Treating two unresolved sink particles as “merged” replaces that calculation with a numerical convention.

Helpful background. Black-Hole Geometry supplies compact-object scales; Relativistic Matter and Compact Objects supplies the remnant and gas components; Black-Hole Astrophysics supplies population and host context; Two-Body Dynamics supplies bound-orbit mechanics; Gravitational-Wave Astronomy supplies the later detector connection.

From galaxy merger to bound pair

Let the black holes have masses $M_1$ and $M_2$. Before their separation is small compared with the remnant’s density scale, each moves mainly in the galaxy potential and loses orbital energy to stars, dark matter, and gas. In an isotropic Maxwellian background, Chandrasekhar’s drag on a perturber of effective mass $M_p$ is

$$\frac{d\boldsymbol v}{dt} =-4\pi G^2M_p\rho\ln\Lambda\, \frac{\mathcal F(X)}{v^3}\boldsymbol v,$$ $$\mathcal F(X)= \operatorname{erf}(X) -\frac{2X}{\sqrt{\pi}}e^{-X^2}, \qquad X=\frac{v}{\sqrt2\sigma}.$$

Chandrasekhar 1943, pp. 255–262 derives the local homogeneous-background result. A galaxy merger violates several of its assumptions: the background is inhomogeneous and time dependent, the orbit is eccentric, the satellite is stripped, and a rotating or gaseous wake can alter the force. The formula is therefore a benchmark and reduced model, not a subgrid truth.

For a circular orbit in a singular isothermal sphere, $\rho=v_c^2/(4\pi Gr^2)$ and $X=1$. With $\mathcal F(1)\simeq0.428$,

$$\frac{dr}{dt} =-\frac{\mathcal F(1)GM_p\ln\Lambda}{v_cr},$$

so the idealized time to reach the center from $r_i$ is

$$t_{\rm df} \simeq\frac{1.17\,v_cr_i^2} {GM_p\ln\Lambda}.$$

The $M_p$ here includes whatever stellar or dark nuclear envelope remains bound to the black hole. Tidal stripping lowers $M_p$ and can lengthen the remaining decay dramatically.

Once the two holes share a nucleus, a necessary binding test is

$$\epsilon_{\rm rel} =\frac{v_{\rm rel}^2}{2} -\frac{G(M_1+M_2)}{r}<0.$$

It is not sufficient in a strongly varying external potential. Also require that the pair remain bound for multiple relative orbits and that its mutual force dominate the differential background force. A useful transition occurs near the radius enclosing stellar mass comparable to $M_1+M_2$. Merritt and Milosavljević 2005, §§ 2–3 separates galaxy-scale sinking, bound-pair formation, and subsequent stellar hardening.

What a merger calculation must resolve

A controlled initial model declares the halo, bulge, nuclear cluster, disk, and gas density profiles; black-hole masses; galaxy mass ratio; orbital energy and angular momentum; spin orientations if relevant; and physical versus comoving coordinates. Redshift enters through galaxy structure and available cosmic time, not as a decorative label. Star formation, feedback, cooling, and black-hole accretion are additional closures whose source terms must be archived.

The central numerical problem is dynamic range. Kiloparsec-scale tidal evolution can be treated with a tree, mesh, or fast-multipole method and softened particle gravity. Parsecs around the holes require sufficiently small particle masses and either much smaller softening or regularized integration. The public ketju module couples GADGET-4 to a regularized close-region integrator so black-hole–star and black-hole–black-hole forces need not be softened. Mannerkoski et al. 2023, §§ 2–3 documents the algorithm and verifies dynamical friction, hardening convergence, and tolerance control.

For a production rerun, pin the exact ketju and GADGET-4 commits, compiler, MPI and HDF5 libraries, compile flags, random seed, parameter file, initial-condition generator, and unit system. Archive particle IDs and snapshots, not only a separation curve. Measure:

• each hole’s orbital energy and angular momentum relative to the evolving remnant;
• its bound stellar-envelope mass and the local density and velocity distribution;
• the separation, relative energy, enclosed background mass, and binding persistence;
• global energy and momentum residuals, including gas and feedback source terms;
• the resolved radius compared with softening, regularized-region size, and influence radius.

Increase particle number while reducing particle mass; refine softening, timestep, and regularized-integrator tolerance independently; and vary the nuclear density profile. The measured pair-formation time should converge before it is passed to a population model. A sink merger at a fixed cell length is a failure mode, not evidence of binary formation.

Reduced decay benchmark

The following teaching calculation verifies the isothermal-sphere decay law before any galaxy simulation is interpreted. It integrates $dr/dt=-A/r$ to a known final radius with explicit Euler steps. The analytic solution is $r^2=r_0^2-2At$. The manufactured dimensionless setup uses $A=\mathcal F(1)$, $r_0=1$, and evaluates the radius at the exact time when $r=0.2$.

Step count $N$	Final-radius result	Absolute error
100	0.217292	0.017292
200	0.209105	0.009105
400	0.204684	0.004684
800	0.202377	0.002377

Use Python 3.12.5 and NumPy==2.1.3. Save the program as massive_binary_decay.py; the script writes massive_binary_decay.csv.

from pathlib import Path
import math
import numpy as np

factor = math.erf(1.0) - 2.0 * math.exp(-1.0) / math.sqrt(math.pi)
r0, target = 1.0, 0.2
exact_time = (r0**2 - target**2) / (2.0 * factor)
rows = []

for steps in (100, 200, 400, 800):
    radius = r0
    dt = exact_time / steps
    for _ in range(steps):
        radius -= dt * factor / radius
    error = abs(radius - target)
    rows.append((steps, radius, error))
    print(f"{steps:4d}  radius={radius:.6f}  error={error:.6f}")

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

Create an isolated environment and rerun with:

python3.12 -m venv .venv-binary-decay
.venv-binary-decay/bin/python -m pip install numpy==2.1.3
.venv-binary-decay/bin/python massive_binary_decay.py

The script produces the table above. It verifies the integrator and analytic normalization only; it does not represent a galaxy. A production comparison replaces the right-hand side with forces measured from the evolving host and must separately demonstrate particle, force, and profile convergence.

Interpretation and failure modes

Gas can accelerate sinking through a dense wake or disk torque, but clumps can scatter a hole and feedback can evacuate the center. A dry remnant can form a pair efficiently if it retains a steep, triaxial nucleus. There is therefore no universal delay assigned solely by black-hole mass ratio.

Common numerical failures include artificially large cores made by softening, spurious two-body heating from overly massive particles, black-hole wandering, a discontinuous change of integration scheme at the edge of a regularized region, premature sink merging, and a truncated outer halo that changes the wake. Vary each numerical scale relative to the physical influence and tidal radii.

Canonical two-body dynamics and detector likelihoods lie outside this page’s boundary. The former is developed in Two-Body Dynamics; the later waveform and selection problem remains in Gravitational-Wave Astronomy. Here the endpoint is a demonstrably bound massive pair with a converged formation time and environmental state.

Estimate sinking and binding

A nucleus has $v_c=200\,{\rm km\,s^{-1}}$, $\ln\Lambda=10$, and a black hole plus bound envelope of $M_p=10^8M_\odot$ at $r_i=1\,{\rm kpc}$. Estimate $t_{\rm df}$. How does complete stripping to $10^7M_\odot$ change it? At $r=20\,{\rm pc}$, two holes have total mass $2\times10^8M_\odot$ and relative speed $150\,{\rm km\,s^{-1}}$; is their isolated relative energy negative?

Solution

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

$$t_{\rm df} =\frac{1.17(200)(1)^2} {(4.3009\times10^{-6})(10^8)(10)} \frac{{\rm kpc}}{{\rm km\,s^{-1}}} \simeq0.0544\frac{{\rm kpc}}{{\rm km\,s^{-1}}} \simeq53\,{\rm Myr}.$$

Reducing the effective perturber mass by ten increases the idealized time to about $530\,{\rm Myr}$. For the pair,

$$v_{\rm esc} =\sqrt{\frac{2G(M_1+M_2)}{r}} \simeq293\,{\rm km\,s^{-1}},$$

so $v_{\rm rel}<v_{\rm esc}$ and the isolated relative energy is negative. A simulation must still check persistence and the differential background potential before declaring a resolved binary.

References

• Chandrasekhar, Subrahmanyan. “Dynamical Friction. I. General Considerations: The Coefficient of Dynamical Friction.” Astrophysical Journal 97 (1943): 255–262. DOI.
• 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.