  The coalescence of the most massive black hole (MBH) binaries releases
gravitational waves (GWs) within the detectable frequency range of Pulsar
Timing Arrays (PTAs) $(10^{-9} - 10^{-6})$ Hz. The incoherent superposition of
GWs from MBH mergers, the stochastic Gravitational Wave Background (GWB), can
provide unique information on MBH parameters and the large-scale structure of
the Universe. The recent evidence for a GWB reported by the PTAs opens an
exciting new window onto MBHs and their host galaxies. However, the
astrophysical interpretation of the GWB requires accurate estimations of MBH
merger timescales for a statistically representative sample of galaxy mergers.
This is numerically challenging; a high numerical resolution is required to
avoid spurious relaxation and stochastic effects whilst a large number of
simulations is needed to sample a cosmologically representative volume. Here,
we present a new multi-mass modelling method to increase the central resolution
of a galaxy model at a fixed particle number. We follow mergers of galaxies
hosting central MBHs with the Fast Multiple Method code Griffin at two
reference resolutions and with two refinement schemes. We show that both
refinement schemes are effective at increasing central resolution, reducing
spurious relaxation and stochastic effects. A particle number of $N\geq 10^{6}$
within a radius of 5 times the sphere of influence of the MBHs is required to
reduce numerical scatter in the binary eccentricity and the coalescence
timescale to <30$\%$; a resolution that can only be reached at present with the
mass refinement scheme.
