Continuous-Time Modelling of Black Hole Binary Evolution with Neural ODEs

TL;DR

The study tackles the challenge of predicting the long-term evolution and merger times of supermassive black hole binaries formed in galaxy mergers, which is crucial for interpreting the gravitational-wave background observed by pulsar timing arrays. It introduces a parameterised neural ordinary differential equation (PNODE) trained on an ensemble of high-fidelity N-body simulations, explicitly conditioned on simulation parameters to model the BH pair’s secular evolution in terms of the orbital energy and angular momentum. The PNODE reliably reproduces the evolution of these descriptors and the derived Keplerian elements across the ensemble, with modest extrapolation to higher resolutions, and when coupled to a semi-analytical model yields merger times consistent with direct N-body results within current uncertainties. This approach provides a computationally efficient, resolution-aware surrogate that enables broader exploration of merger timescales and GWB constraints without sacrificing fidelity in the relevant dynamical regime.

Abstract

Pulsar timing arrays (PTAs) can detect the low-frequency stochastic gravitational-wave background (GWB) generated by an ensemble of supermassive black hole binaries (BHBs). Accurate determination of BHB merger timescales is essential for interpreting GWBs and constraining key astrophysical quantities such as black hole (BH) occupation fractions and galaxy coalescence rates. High-accuracy $N$-body codes such as \texttt{Griffin} can resolve sub-pc BHB dynamics but are too costly to explore a wide range of initial conditions, motivating the need for surrogate models that emulate their long-term evolution at much lower computational cost. We investigate neural ordinary differential equations (NODEs) as surrogates for the secular orbital evolution of BHBs. Our primary contribution is a parameterised NODE (PNODE) trained on an ensemble of $N$-body simulations of galaxy mergers spanning a two-dimensional parameter space defined by the initial orbital eccentricity and particle resolution $(e_i, N)$, with the learned vector field explicitly conditioned on these parameters. A single PNODE thereby learns a simulation-parameter-conditioned dynamical model for the coupled evolution of the BH pair's orbital state across the ensemble, yielding smooth trajectories from which stable hardening and eccentricity growth rates can be extracted. The PNODE accurately reproduces the secular evolution of the specific orbital energy and angular momentum, and the corresponding Keplerian orbital elements, for held-out trajectories, with modest generalisation to a partially unseen high-resolution case. Combining PNODE predictions with semi-analytical prescriptions for stellar hardening and gravitational-wave emission yields BHB merger timescales consistent with those obtained from direct $N$-body inputs within current theoretical uncertainties.

Figures (5)

• Figure 1: NODE computational flow: the initial state $\mathbf{h}(0)$ is evolved forward in time via a learnt vector field $f_\theta$, integrated using an ODE solver to produce the final state $\mathbf{h}(T)$. Here, $t$ denotes a continuous variable that may represent either an auxiliary depth parameter or the real physical time, depending on the context.
• Figure 2: Predictions of PNODE for $\varepsilon$ and $h$ (dashed lines) against truth (simulation) for selected BH trajectories in test set for initial galactic orbital eccentricities $e_{i}=0.9, 0.97, 0.99$ at resolutions (32M, 64M, 128M). All quantities are shown in the dimensionless code units of the simulations.
• Figure 3: Predictions of PNODE for $e$ and $a$ (dashed lines) against truth (simulation) for selected BH trajectories in test set for initial galactic orbital eccentricities $e_{i}=0.9, 0.97, 0.99$ at resolutions (32M, 64M, 128M). All quantities are shown in the dimensionless code units of the simulations.
• Figure 4: Evolution of the orbital elements $a$ and $e$ of the BHB as a function of time, from the start of the $N$-body simulation to coalescence of the MBHs. The lines represent the elements calculated via the semi-analytical model, with hardening rates $s$ and $K$ from a direct numerical differentiation and fit of both the $N$-body data (solid lines) and the predictions from PNODE (dashed lines). The circles represent the original $N$-body data, , from the time of binary formation to the end of the simulation, while the stars represent the orbital elements calculated at early times, when the BHs are still unbound, from the apocentre and pericentre values. The total evolutionary time from the beginning of the simulation to coalescence provides an estimate of the total merger timescale of the BHBs, as predicted through extrapolation.
• Figure 5: Forecasts of the Keplerian orbital elements $(\ln(a), e)$ for a single $N$-body simulation under several classical and neural models (FFNN, BNN, LSTM, XGBoost, VARIMA, NODE) at different chronological train/test splits. The blue dashed line indicates the start of the forecasting regime, and the semi-major axis $a$ is expressed in kpc (we plot $\ln(a)$). Some models (e.g. XGBoost, VARIMA) exhibit visibly noisier trajectories due to the absence of inherent temporal smoothing; these are shown for completeness and comparison, and across all methods the forecast errors over this approximately linear interval are broadly comparable, with no clear performance winner and no systematic advantage for more flexible models such as the NODE or less flexible but strongly regularised models such as the BNN.