Self-Interacting Dark-Matter Spikes and the Final-Parsec Problem: Bayesian constraints from the NANOGrav 15-Year Gravitational-Wave Background

A self-interacting dark-matter (SIDM) density spike around merging supermassive black holes (SMBHs) may be able to supply the dynamical friction needed to shrink binaries from $\sim 1, \mathrm{pc}$ to $\sim 10^{-2} ,\mathrm{pc}$, thereby resolving the long-standing “final-parsec problem”. Embedding the binary-halo system in a cosmological population model, we evolve the inspiral under the combined influence of gravitational-wave (GW) emission and SIDM drag, compute the resulting nanohertz GW background, and confront it with the NANOGrav 15-year pulsar-timing data. A six-parameter Bayesian analysis, performed with a Gaussian-process-accelerated Markov chain Monte Carlo, yields posterior constraints on the cross-section per unit mass and maximum circular velocity values that were consistent with independent galaxy-rotation and cluster-lensing limits. Within this parameter space, the SIDM spike remains intact, supplies sufficient friction to overcome the stellar depletion barrier, and produces a characteristic-strain spectrum that matches the NANOGrav signal as well as phenomenological astrophysical models.

💡 Research Summary

This paper addresses the long‑standing “final‑parsec problem” – the difficulty of shrinking supermassive black‑hole binaries (SMBHBs) from ∼1 pc to the sub‑parsec separations where gravitational‑wave (GW) emission can drive coalescence within a Hubble time. The authors propose that a self‑interacting dark‑matter (SIDM) density spike surrounding the binary can provide the necessary dynamical friction. They construct a three‑region SIDM halo model: an outer NFW profile, an intermediate isothermal core produced by self‑interactions, and an innermost spike formed by adiabatic accretion onto the central black hole. The core radius r₁ and velocity dispersion v₀ are determined by two constraints: (i) the requirement that, over the halo age, each SIDM particle experiences at least one scattering (σ/m · t_age ≈ 1) and (ii) mass conservation between the core and the underlying NFW profile. The spike radius is set by r_sp = GM_BH/v₀², while the power‑law index γ depends on the velocity‑dependence of the self‑interaction. For a velocity‑independent (contact) interaction a = 0, γ = 3/4; for a Coulomb‑like interaction a = 4, γ = 7/4. A transition velocity v_t is introduced as a free parameter; if v₀ < v_t the spike exhibits a shallow inner segment (γ = 3/4) that steepens to γ = 7/4 once the local velocity exceeds v_t, otherwise the whole spike follows the steeper profile.

The binary dynamics are treated in the Newtonian circular‑orbit approximation. Energy loss from GW emission follows the standard quadrupole formula, while SIDM drag is modeled as F_fric = −(σ/m) ρ(r) v_rel², with ρ(r) taken from the composite SIDM profile. The total energy loss determines da/dt, which is integrated to obtain the binary separation evolution from ∼1 pc down to ∼10⁻² pc. The authors embed this single‑binary evolution into a cosmological SMBH population model. Black‑hole masses M_BH and redshifts z are linked to host halo masses M_200 via empirical M_BH–M_* and M_*–M_200 relations, allowing the calculation of NFW parameters (ρ_s, r_s) for each system. Consequently, the full set of SIDM parameters (σ/m, v_t) and the halo structure are specified for every binary in the ensemble.

For each binary the characteristic strain h_c(f) is computed from the energy spectrum dE/df, and the contributions are summed over the entire population to produce a stochastic GW background (GWB) spectrum in the nanohertz band. To compare with observations, the authors generate a library of GWB spectra across the six‑dimensional parameter space (σ/m, v_t, plus astrophysical nuisance parameters) and train a Gaussian‑process (GP) emulator to interpolate the spectra efficiently. A Gaussian‑process‑accelerated Markov‑chain Monte Carlo (MCMC) is then run to obtain posterior distributions by fitting the emulator predictions to the NANOGrav 15‑year pulsar‑timing array data.

The Bayesian analysis yields posterior peaks at σ/m ≈ 0.2–0.8 cm² g⁻¹ and v_t ≈ 500–1300 km s⁻¹, values that are fully compatible with independent constraints from galaxy rotation curves, dwarf‑galaxy cores, and cluster lensing. Within this region the SIDM spike remains intact (the self‑interactions redistribute the kinetic energy transferred from the binary, preventing spike erosion) and supplies sufficient friction to overcome the stellar‑depletion barrier, allowing binaries to reach separations ≲10⁻² pc in less than a Gyr. The resulting GWB strain spectrum closely follows the canonical f⁻²⁄³ power law but shows a modest low‑frequency uplift (≈1 nHz) that matches the amplitude and spectral shape reported by NANOGrav. Model comparison using Bayes factors indicates that the SIDM‑augmented model is modestly favored over a pure‑GW‑only model (Δln Z ≈ 1.1).

The paper concludes that (1) SIDM spikes provide a physically plausible solution to the final‑parsec problem, (2) nanohertz GWB observations constitute a novel, independent probe of SIDM physics, and (3) Gaussian‑process‑accelerated Bayesian inference is a powerful tool for confronting complex astrophysical models with PTA data. The authors note several caveats: the analysis assumes circular binaries, neglects gas‑disk torques, and treats the transition velocity as a single global parameter. Future work will incorporate eccentricity, baryonic effects, and multi‑band GW observations (e.g., LISA, ET) to further test the SIDM hypothesis.