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 ∼ 5 × 105, 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 ∼ 107 within the stellar half-light radius suffices to reduce the scatter in the merger timescale to ∼10%. Our results suggest that at least some of the uncertainty in low-frequency GW rates owes to insufficient numerical resolution.