Abstract
The coalescence of massive black hole binaries (BHBs) in galactic mergers is
the primary source of gravitational waves (GWs) at low frequencies. Current
estimates of GW detection rates for the Laser Interferometer Space Antenna and
the Pulsar Timing Array vary by three orders of magnitude. To understand this
variation, we simulate the merger of equal-mass, eccentric, galaxy pairs with
central massive black holes and shallow inner density cusps. We model the
formation and hardening of a central BHB using the Fast Multiple Method as a
force solver, which features a $O(N)$ scaling with the number $N$ of particles
and obtains results equivalent to direct-summation simulations. At $N \sim
5\times 10^5$, typical for contemporary studies, the eccentricity of the BHBs
can vary significantly for different random realisations of the same initial
condition, resulting in a substantial variation of the merger timescale. This
scatter owes to the stochasticity of stellar encounters with the BHB and
decreases with increasing $N$. We estimate that $N \sim 10^7$ within the
stellar half-light radius suffices to reduce the scatter in the merger
timescale to $\sim 10$\%. Our results suggest that at least some of the
uncertainty in low-frequency GW rates owes to insufficient numerical
resolution.