Abstract Binaries of supermassive black holes (MBHBs) represent the primary sources of the gravitational wave background (GWB) detectable by Pulsar Timing Arrays (PTAs). The eccentricity with which binaries form in galactic mergers is the key parameter determining their evolutionary timescale from pairing to coalescence. However, accurately determining the binary eccentricity at formation is difficult in simulations due to stochastic effects. We present a numerical study of the formation and evolution of MBHBs that are potential PTA sources. We simulate mergers of equal-mass galaxies on different initial orbits and follow the dynamics of the MBHBs through the hardening phase. We find that low-resolution simulations are affected by stochasticity due to torques from the stellar distribution acting at pericentre passages. The dispersion in binary eccentricity decreases with increasing central resolution, as expected for a Poisson process. We provide a fitting formula for the resolution requirement of an N-body simulation of MBHB formation and evolution as a function of the initial eccentricity of the merger, e0, and the required accuracy in the binary eccentricity, eb. We find that binaries experience a torque at first pericentre that is approximately independent of initial eccentricity, producing a general trend in which the binary eccentricity decreases above sufficiently large initial orbital eccentricities. While this behaviour is generic, the precise cross-over eccentricity (e0 ∼ 0.97 in our models) and the sharpness of the drop-off depend on the galaxy initial conditions. We provide a fitting formula for eb(e0) that can be used in semi-analytical models to determine the merger timescales of MBHBs as well as the amplitude and slope of the GWB.
Gualandris et al. (2026) studied this question.