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 similar to 5 x 10(5), typical for contemporary studies, the eccentricity of the BHBs can vary significantly for different random realizations of the same initial condition, resulting in a substantial variation of the merger time-scale. This scatter owes to the stochasticity of stellar encounters with the BHB and decreases with increasing N. We estimate that N similar to 10(7) within the stellar half-light radius suffices to reduce the scatter in the merger time-scale to similar to 10 per cent. Our results suggest that at least some of the uncertainty in low-frequency GW rates owes to insufficient numerical resolution.