A Multi-Resolution Method for Modelling Galaxy and Massive Black Hole Mergers

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.


INTRODUCTION
The detection of Gravitational Waves (GWs) from two merging black holes by the LIGO scientific collaboration (Abbott et al. 2016) has opened a new window on the Universe, providing fresh insights into the astrophysical sources of measurable GWs.The signal has proven the existence of black holes and confirmed the prediction that binaries of black holes coalesce.Since this first detection some 90 events have been observed to date, including both black hole and neutron star mergers (The LIGO Scientific Collaboration et al. 2021), providing direct measurements of black hole masses and spins.These have provided new impetus to ongoing efforts to detect GWs from binaries of supermassive black holes (BHBs), which represent the loudest sources of GWs in the Universe (Peters 1964).
Supermassive black holes (MBHs; 10 5 ≲  ≲ 10 9 ) are found in the centres of all but the smallest galaxies (Ferrarese & Ford 2005;Kormendy & Ho 2013), and their masses correlate tightly with the properties of their host galaxies through scaling relations such as the MBH mass-bulge mass relation and the MBH massstellar velocity dispersion relation (Kormendy & Ho 2013;Reines & Volonteri 2015).MBHs grow over time through gas/star accretion, and via mergers with other MBHs following galactic mergers (Begelman et al. 1980).The evolution of BHBs proceeds over three ★ E-mail: k.attard@surrey.ac.uk major phases: dynamical friction, hardening and gravitational wave inspiral and coalescence.First, dynamical friction causes MBHs to sink to the centre of the merger remnant until they form a bound pair.The pair then hardens due to three-body interactions with stars during the rapid gravitational slingshot phase (Quinlan 1996;Sesana et al. 2008).Following these strong encounters, stars are ejected with velocities comparable to the binary's orbital velocity and a core is carved in the stellar distribution (e.g.Merritt & Milosavljević 2005;Gualandris & Merritt 2012).This 'core scouring' by a BHB is the leading mechanism for explaining the observation of large stellar cores in massive elliptical galaxies (Merritt 2006).Once all stars initially on losscone orbits have been ejected and the binary is formally bound and hard (at roughly parsec scale separations for typical 10 8 M ⊙ MBHs), further hardening relies on efficient repopulation of the losscone through angular momentum diffusion, a mechanism which has been proven efficient in non-spherical models like merger remnants (Vasiliev et al. 2014(Vasiliev et al. , 2015;;Gualandris et al. 2017;Frigo et al. 2021).Such 'collisionless losscone refilling' is able to drive the BHB to separations (of order milliparsec or less) where emission of GWs commences.Inspiral then proceeds quickly in this final GW phase, leading the MBHs to coalescence.
Mergers of the most massive MBHs ( > 10 7 M ⊙ ) are detectable by Pulsar Timing Arrays (PTAs) in the frequency range (10 −9 −10 −6 ) Hz (Desvignes et al. 2016;Reardon et al. 2016;Perera et al. 2019), while MBHs in the mass range (10 4 − 10 7 ) M ⊙ will be the main target of the Laser Space Interferometer (LISA) in the frequency range (10 −4 − 1) Hz (Amaro-Seoane et al. 2017;Barack et al. 2019).LISA detections will provide accurate estimates of masses, spins and orbital parameters up to a redshift of  = 20.Unlike electromagnetic radiation, GWs travel essentially unobstructed across space and allow us to probe the early evolution of MBHs and their host galaxies at high redshift (Colpi et al. 2019).
Recently, the European Pulsar Timing Array (EPTA), the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) and the Parks Pulsar Timing Array (PPTA) reported evidence for a spatially correlated stochastic red noise signal consistent with expectations for a stochastic Gravitational Wave Background (GWB) Antoniadis et al. (2023a); Agazie et al. (2023); Reardon et al. (2023).This represents the incoherent superposition of unresolved binary mergers, which manifests itself as a long timescale, low frequency common signal across the pulsars in the array.The astrophysical interpretation of the signal will provide information on the large-scale structure of the Universe, the occupation fraction of MBHs in galaxies, the merger rate of galaxies, scaling relations between MBHs and host galaxies and the efficiency of pairing and merging of BHBs (Antoniadis et al. 2023b).
In this context, it is crucial to model GW sources at low frequencies and estimate detection rates for both current and upcoming instruments.These predictions help constrain expected GW event numbers and aid in the interpretation of any signal and its astrophysical implications (Amaro-Seoane et al. 2022).Current estimates come mainly from semi-analytical models and cosmological simulations, and are affected by significant uncertainties.
Semi-analytical models have been used to follow the evolution of MBHs and their host galaxies, including dark matter (DM) halos, within the context of the ΛCDM framework (e.g.Somerville et al. 2008;Volonteri & Natarajan 2009;Sesana et al. 2014;Lacey et al. 2016).Starting from a DM merger tree, structures are evolved using semi-analytic prescriptions.Recent models include prescriptions for MBH growth and feedback, and have been shown to reproduce observed properties such as morphology, colour, star formation rates and luminosity functions at low redshift (e.g.Ricarte & Natarajan 2018;Barausse et al. 2020).Being computationally inexpensive relative to -body simulations, they can be used to explore a wide parameter space (Lagos et al. 2018;Izquierdo-Villalba et al. 2022).However, they rely on simplified assumptions and prescriptions regarding the pairing and evolution of BHBs, with key phases like the binary hardening due to encounters with surrounding stars being modelled based on the results of isolated three-body scattering experiments (Sesana et al. 2006;Sesana & Khan 2015;Izquierdo-Villalba et al. 2022).Both low and high-mass BH seeds must be considered for the MBH population at high redshift, as seed masses are currently largely unconstrained (Volonteri 2010;Amaro-Seoane et al. 2022).
By contrast, cosmological -body and hydrodynamic simulations model baryonic and dark matter on large scales from early times, following subsequent galactic mergers and gas accretion simultaneously.Galactic merger mass ratios and encounter geometries are consistently modelled in subsequent events, with realistic dynamical friction times.Massive halos can be seeded with black holes, which then grow over time due to gas accretion and mergers.The inclusion of feedback effects results in gas expulsion, and the quenching of both star formation and MBH growth (Di Matteo et al. 2005;Pontzen et al. 2017;Ricarte et al. 2019;Nelson et al. 2019a).Formation of BHBs during mergers is then followed down to the resolution limit, generally set by the choice of the softening length, which is of order 1 kpc in Illustris (Genel et al. 2014) and of order 300 pc in the recent Illustris TNG50-1 (Nelson et al. 2019b;Pillepich et al. 2019).In the context of BHB evolution, MBHs have only just begun pairing at these distances, and are still unbound to each other.BHBs are then assumed to merge promptly, without any modelling of the hardening phase.A trade-off between simulation volume and resolution is also present, limiting the statistics of BHBs in higher resolution simulations.
Predictions for GW detectors have been obtained a posteriori using the properties of merging BHBs in cosmological simulations as input to semi-analytical modelling (Kelley et al. 2017;Katz et al. 2020;Barausse et al. 2020).These generally include an externally imposed delay to account for the time spent by the binary in the dynamical friction and hardening phase, resulting in significantly fewer mergers, especially for models with high mass seeds, with the lowest LISA rates reaching only a few mergers per 4 yr nominal mission duration (Barausse et al. 2020).Typical LISA merger rates vary significantly for different MBH assumptions, from ∼ (20 − 110) yr −1 (Sesana et al. 2011) to ≲ 1 yr −1 (Katz et al. 2020), but appear higher when MBH triple interactions are taken into account (Bonetti et al. 2019).Current predictions give GWB amplitudes in the PTA frequency range, including characteristic strains of 6.9 × 10 −16 from the MassiveBlack-II cosmological simulation (Sykes et al. 2022), of 1.4 × 10 −16 − 1.1 × 10 −15 from the SHARK semi-analytic model (Curyło & Bulik 2022) and 1.2 × 10 −15 from the L-Galaxies model (Izquierdo-Villalba et al. 2022).
An alternative approach to ensure higher resolution in the central regions of galaxies where MBHs bind into binaries is provided by high accuracy simulations of individual mergers by means of direct summation simulations (e.g.Gualandris et al. 2022), hybrid collisionless/collisional simulations (e.g.Nasim et al. 2021), and hybrid collisionless/regularised simulations (e.g.Rantala et al. 2017).Direct summation simulations achieve low force errors by brute force calculations of all pairwise gravitational forces, but are limited to particle numbers of order a million by their O ( 2 ) computational scaling (Gualandris & Merritt 2012;Gualandris et al. 2022;Khan et al. 2020).At these low resolutions, stochastic effects lead to spurious Brownian motion of the BHB (Bortolas et al. 2016) as well as a significant scatter in the eccentricity of the binary (Nasim et al. 2020).Increasing resolution while maintaining a bound on force errors is challenging, but has been achieved in the Fast Multiple Method (FMM) code griffin with a sophisticated monitoring of errors and adaptive choice of parameters (Dehnen 2014), which we employ in this work.Here, the MBH-MBH and MBH-star interactions are modelled with direct summation while the star-star interactions are modelled with the more efficient, yet still accurate, FMM.An alternative hybrid approach is implemented in the KETJU code (Rantala et al. 2017).This uses a standard tree code for most force calculations (building upon the GADGET-3 collisionless code), but an accurate regularised treatment for the BHB and BHB-star encounters, including relativistic corrections through the Post Newtonian formalism (Mannerkoski et al. 2019(Mannerkoski et al. , 2022)).
Modelling galaxies and galactic mergers at the large resolutions  ≳ 10 7 required to avoid stochastic effects (Nasim et al. 2020) is challenging even with more efficient codes like griffin due to the large range of spatial and temporal scales involved.In this paper, we present a numerical scheme to increase resolution in the central regions of galaxies based on mass refinement.We present results of both isolated galaxy models and merging galaxies and show that mass refinement schemes are effective at increasing resolution in galactic centres, without any adverse effects on stability or dynamical evolution.In particular, we show that a mass refined  = 10 6 galaxy model behaves similarly to a  = 10 7 reference galaxy model.
This paper is organised as follows: we first describe the multi-mass In Sec. 5 we extrapolate the evolution to late times and estimate the total coalescence timescales.We conclude with a discussion and summary in Sec. 6.

A MULTI-MASS RESOLUTION SCHEME
We implement a multi-mass refinement scheme to increase the resolution in the central regions of galaxy models, with the aim to accurately model the evolution of BHBs formed in galactic mergers at a reduced computational cost.In the case of multi-component models featuring a central MBH, a stellar bulge, and a DM halo, the scheme is applied independently to the bulge and the halo.
Starting from a desired total particle number for each component, which we denote  b and  h , stars are first oversampled by a factor , termed the multiplication factor, and then coarse-grained based on radial position.A number of radial zones can be defined such that stars in the innermost zone are left unchanged, while only a given fraction of stars is retained in the outer zones.This fraction decreases moving away from the galaxy centre.The mass of each particle in the outer zones is then increased by the same factor  to ensure that the total mass in the galaxy and the mass density profile are unaffected.We adopt a larger number of zones for higher multiplication factors in order to reduce the increase in particle masses at the outermost radii.In this way, we increase the resolution in the central region of the galaxy whilst avoiding large jumps in particle masses between the radial zones, which would otherwise induce spurious mass segregation.A commonly adopted mass factor to avoid mass segregation is of order 10 (Alexander & Hopman 2009), and this is consistent with the results shown in Fig. 13.
We implement two different schemes, labelled 'a' and 'b', with parameters as listed in Table 1.Scheme 'a' is our standard implementation, characterised by  = 3 zones, while scheme 'b' is a more aggressively refined scheme with  = 4 zones.The latter scheme can be considered physically motivated, as its innermost radial bin is equal to twice the radius of the sphere of influence of the central MBH, one of the fundamental scales in the system (Peebles 1972) where  bh represents the MBH mass and  the local 1D stellar velocity dispersion.The radial boundaries and mass multiplication factors are given in Table 1.In scheme 'a', particles are initially over-seeded by a factor 5, with the first (innermost) zone retaining all particles, and its radius set by the requirement of containing 50% of all particles (  and  ℎ , respectively, for bulge and halo).The particle number is then reduced by a factor of 5 for each successive radial zone, and the multiplication factor for the remaining particle masses increases correspondingly.This multiplication factor determines the placement of the boundary between the outer radial zones.Scheme 'b' is similar, with an additional zone being added to prevent an excessively large increase in mass in the outermost radial zone.In this scheme, the innermost 50% of the particles is divided between zone 1 (with 1%) and zone 2 (with 49%), while the outermost 50% is divided between zone 3 (with 37.5%) and zone 2 (with 12.5%).

Galaxy models
We consider a set of multi-component models representative of massive elliptical galaxies consisting of a stellar bulge, a dark matter (DM) halo, and a central MBH.The stellar bulge follows a Sérsic profile (Sérsic 1963;Sersic 1968) in projection where  0 is the surface intensity,   is the scale radius,  is the Sérsic index, and   is a function of , where   ≈ 2 − 1 3 + 4 405 (Ciotti & Bertin 1999).We set  = 4, appropriate for elliptical galaxies, as in the de Vaucouleurs profile (de Vaucouleurs 1948).The DM halo follows an NFW profile (Navarro et al. 1996) where  is the scale radius and  0 is a normalisation factor.We extract merger trees from the IllustrisTNG-300-1 cosmological simulation (Pillepich et al. 2019;Springel et al. 2018;Marinacci et al. 2018;Nelson et al. 2019a;Naiman et al. 2018) and select a major merger at low redshift whose central MBH is an expected PTA source.The MBH has a mass  bh = 7.14 × 10 8 M ⊙ .We then set the parameters for the bulge and halo components using observational scaling relations.The stellar bulge mass  b is derived from the MBH-stellar mass relation log (Kormendy & Ho 2013;Reines & Volonteri 2015) with parameters  = 8.95 ± 0.09 and  = 1.40 ± 0.21 taken from fits to observed elliptical galaxies (Reines & Volonteri 2015).This gives  b = 8.54× 10 10 M ⊙ .
The halo mass  h is taken from the halo-stellar mass scaling relation of Chae et al. (2014) for massive elliptical galaxies, with a typical value  h = 5.04 × 10 12 M ⊙ .The virial radius  200 is derived from the relation with the halo mass (Binney & Tremaine 2008), which gives  200 = 353.69kpc.The scale radius of the NFW profile is given by  =  200 /, where the concentration parameter  is derived from the  −  h relation in Dutton & Macciò (2014).Lastly, we set the half-mass radius of the stellar component using the relation  1/2 = 0.015 200 (Kravtsov 2013), giving  1/2 = 5.30 kpc.
We utilise Agama, an action-based galaxy modelling software library (Vasiliev 2019), to sample the total gravitational potential of the multi-component galaxy system and construct equilibrium initial conditions for our simulations.We adopt units such that  = 1,  = 1 kpc and  =  b .We generate isolated models, denoted with identifier 'I', at two reference resolutions,  = 10 6 and  = 10 7 , and apply the mass refinement schemes described in section 2 for each resolution.We then set up equal mass mergers of two independent realisations of the same multi-component model, denoted with identifier 'M', at the same resolutions.Galaxies are placed on bound Keplerian orbits of eccentricity  = 0.7 and at an initial distance of  = 378.6 kpc.The initial eccentricity is relatively high, in agreement with typical orbits in large scale cosmological simulations (Khochfar & Burkert 2006), but not so high to result in a head-on galaxy collision.The initial distance is set to  ≃ 9   , where   represents the scale radius of the NFW model describing the halo of the galaxies, defined as the radius where the particle number drops to half of the total.The parameters of all models are listed in Table 2, including the total particle number , the number of halo and bulge particles, and the ratio of halo/bulge particle mass to the MBH mass.The final column gives the numerical multiplier  adopted in the mass refined schemes.
While the ratio between a halo particle mass and a bulge particle mass is equal for all models ( ℎ /  ∼ 6.7), the MBH to bulge particle mass ratio at the innermost radial shell decreases as the refinement scheme becomes more aggressive.In particular, I6b has the same  b / bh as I7, and the smallest ratio is obtained for I7b.For all the merger models with  = 10 6 , we produce 7 different random realisations, denoted i to vii, in order to measure the scatter in the orbital parameters.

Numerical simulations
We evolve all models with griffin (Dehnen 2014), a Fast Multiple Method (FMM) code that monitors force errors and utilises adaptive parameters to ensure a distribution of force errors similar to that in a direct summation code, while retaining the O ()) scaling of the FMM technique.The code has been shown to perform similarly to a direct summation code in modelling the evolution of BHBs formed in galactic mergers (Nasim et al. 2020).
The simulations reported here use a multipole expansion order  = 5 and a relative total force error of ≲ 1×10 −3 for the  = 10 6 models and ≲ 5 × 10 −5 for the  = 10 7 models.The softening parameter is initially set to  = 30 pc for star-star interactions and to  bh = 10 pc for MBH-MBH and MBH-star interactions for all models.We then reduce the softening length  bh in the merger simulations to 5 pc after a time ∼ 2 Gyr corresponding to just prior to the end of the dynamical friction phase.This setup allows to reduce computational time in the early stages of the merger, driven by dynamical friction, and yet to accurately model the evolution of the BHB through the rapid hardening phase dominated by encounters with background stars.The softening kernel is a near-Plummer variation for smooth source particles with a density  ∝ ( 2 +  2 ) −3.5 .All simulations are evolved until the binary reaches a separation of order the softening length  bh , below which the dynamical evolution is no longer reliable.We ensure that this separation is always smaller than the hard-binary separation.

EVOLUTION OF THE MULTI-MASS MODELS
Mass refinement schemes of the type presented in this work have been adopted in the literature with the aim of increasing resolution in the central regions of a single galaxy model, with Cole et al. (2012) allowing a reduction in the total particle number by half at the same central resolution without any adverse effects.Here we test the viability of our multi-zone multi-mass schemes, for both isolated and merger models.Galactic mergers are violent events that affect the distribution of stars at both large and small scales, and as such A significant increase in the particle number < 10 kpc is observed for the models with the multi-mass scheme applied, and no significant expansion occurs over time for these models.it is not guaranteed that the refinement will work over time in a merger simulation.However, the 'mixing theorem' shows that body simulations should preserve central phase space density, even through mergers (Dehnen 2005).As such, we can anticipate that the high resolution regions at the centre of each merging galaxy should overlap once the merger is complete.
The distribution of bulge particles in the isolated models, with and without mass refinement, is shown in Fig. 1 at an early time and a late time in the evolution.This illustrates the efficiency of the refinement schemes, with increased population in the central regions of models 'a' and 'b'.The mass refined models retain their higher central resolution to later times, in agreement with earlier studies (Cole et al. 2012).Crucially, this holds true in the merger models as well, as shown in Fig. 2, despite the violent mixing happening during the merger process.Models M6a and M6b achieve a population within the central ∼ 100 pc similar to model I7.
The radial density profiles shown in Fig. 3 and 4 support the  conclusions from the particle distributions.The schemes are effective at populating the central regions of the models, and the distributions are stable over time.The flattening observed in the inner density profiles of the merger models can be attributed to core scouring during the hardening phase of the black hole binary.
The shape and kinematic properties of the models are also preserved by the refinement schemes.The shape of a galaxy or merger remnant can be characterised by the axis ratios of an ellipsoid fit to the stellar distribution, defined as where  >  >  > 0 are the axes of symmetry of the ellipsoid.The ratios / and / are equal to unity for spherical models, and decrease for axisymmetric (/ < 1) or triaxial (/ < 1 and / < 1) models (de Zeeuw 1985).The triaxiality parameter provides a measure of the departure from sphericity, computed as with values  > 0.5 for prolate spheroids,  < 0.5 for oblate spheroids, and  = 0.5 for maximum triaxiality.We calculate the axis ratios using Pynbody, a Python package for analysis of astrophysical simulations (Pontzen et al. 2013).The radial dependence of the axis ratios for the stellar component is shown in Figure 5 for the isolated models and in Figure 6 for the merger models, including also the triaxiality parameter .In all isolated models, the sphericity is maintained over the simulation time, save for a noticeable departure at the innermost radii in model I6.This is likely due to low resolution in the innermost radial bins, which is resolved by the mass refinement in models I6a/b.The shape is preserved by the schemes also for the merger models, which display significant triaxiality, here in the form of a prolate shape ( > 0.5), in agreement with results from Gualandris et al. (2017); Bortolas et al. (2018).
We present 3D stellar velocity dispersion profiles in Fig. 7 for the isolated models and Fig. 8 for the merger models, showing a small degree of deviation between the low and high resolution models.In particular, models I6 and M6 show a higher velocity dispersion at all radii.The mass refined models, on the other hand, are consistent at both resolutions, and consistent with M7.
A similar behaviour can be seen in the anisotropy profiles (see The spherical models maintain their shape over time even when the mass refinement schemes are applied. where   () is the tangential velocity dispersion and   () is the radial velocity dispersion.The isolated models are isotropic at all radii, with the exception of model I6 which displays a small radial bias over time.The merger models display the expected anisotropy with an inner tangential bias, imparted from core scouring as stellar scatterings during the hardening phase of binary evolution preferentially eject stars on radial orbits (Quinlan & Hernquist 1997;Thomas et al. 2014).Interestingly, there is a modest discrepancy at low resolution even when the schemes are applied, suggesting that the binary hardening phase is sensitive to the specific sequence and properties of encounters with stars experienced over time.No discrepancy is observed at the larger resolution, with or without mass refinement.
The Lagrangian radii shown in Fig. 11 and Fig. 12 for isolated and merger models, respectively, reveal a small expansion in the bulge particles, and predominantly those at lower radii.This owes to relaxation effects at small radii and at low resolution, as indicated by the    A modest discrepancy between the models is observed at the lower resolution, while at the higher resolution the models are indistinguishable.
Figure 11.Lagrangian radii for bulge (bottom four curves) and halo (top four curves) particles for all isolated models with  = 10 6 (left) and  = 10 7 (right) across the simulation time.The curves correspond to radii containing mass fractions of 3%, 10%, 13%, and 20% of the total bulge mass.I6 shows a significant expansion at the lowest bulge radius, which is notably reduced in the multi-mass models.The halo particles remain stable in all models.
dependence of the effect on particle number and radius.The expansion lessens with increasing radius, due to an increase in relaxation time with radius, and is almost negligible in both the refined  = 10 6 models and all the  = 10 7 models.Furthermore, it does not affect the halo particles.The strong variation at  = 2.0 − 2.5 Gyr in the merger models is a signature of the merger process, and can be taken as a measure of the merger time.
The results presented in Fig. 3 and 4 are in agreement with expectations based on a simple calculation of the dynamical friction timescale for point mass particles.A particle of mass  in a system of lighter particles will sink from a radius   on a timescale given by Chandrasekhar's dynamical friction formula (Binney & Tremaine 2008) where  is the local stellar velocity dispersion.The timescale  df of the halo particles is plotted in Fig. 13 for both isolated and merger models as a function of radial distance.The chosen radial positions correspond to the centres of the scheme zones.The horizontal line indicates the total simulation time of 8 Gyr for the isolated models, and 4 Gyr for the merger models.Halo particles in the I6 models are affected by dynamical friction against the sea of bulge particles at distances  ≲ 2 kpc while mass refined models I6a and I6b are only affected within  ≲ 1 kpc.At higher resolution, sinking is expected in model I7 only for particles within the central ≲ 500 pc, which is reduced to the central ≲ 200 pc in the mass refined models I7a and I7b.We notice that the multi-mass models experience a downturn in the calculated  df at large radii  ≳ 10 kpc, due to the increase in particle masses introduced by the schemes to compensate for the reduced particle number.Even so, the schemes appear beneficial out to ∼ 40 kpc.The merger models exhibit a very similar behaviour to the isolated models.We note that the total wall-clock time of the refined and unrefined models is comparable, and the mass refinement schemes do not produce any appreciable slow-down in our Griffin runs.

EVOLUTION OF THE MASSIVE BLACK HOLE BINARIES
We follow the evolution of the MBHs through the merger, the formation of a bound binary and the collisionless losscone refilling phase and employ a semi-analytical model to compute the merger timescale for the different models due to GW emission.We examine the effect of resolution and of the multi-mass schemes on the orbital elements of the binary as well as the coalescence timescale.Using the 7 random realisations of the  = 10 6 merger models, we quantify the scatter in the orbital elements of the binary and the total merger time.
We adopt the dispersion in semi-major axis, eccentricity and merger timescales as our criterion for convergence.
The three characteristic phases of binary evolution can clearly be seen in Fig. 14, which shows the distance between the MBHs as a function of time.The galaxies first evolve under the effects of dynamical friction, which brings them together during the merger and, at the same time, causes the MBHs to sink to the centre of the merger remnant, where they form a pair and eventually a bound binary.At the end of this phase, the MBHs are approximately at a separation  f , where the enclosed stellar mass is equal to twice the mass of the secondary MBH, as Three-body interactions between the binary and background stars become increasingly important until they dominate the evolution, as the influence of dynamical friction wanes once the binary begins to reach thermal equilibrium with surrounding stars (Antonini & Merritt 2012;Kelley et al. 2017).Such interactions are known as slingshot interactions as the stars remove energy and angular momentum from the binary and are ejected to large distances.As a result, the MBHs spiral closer together (Begelman et al. 1980;Merritt 2013) until all stars initially belonging to the losscone are ejected; this occurs roughly at the hard-binary separation  h , which is of order a parsec for black holes of  ∼ 10 8 M ⊙ .There are different definitions of the hard-binary separation in the literature, here we adopt the formal one given by the separation where the relative velocity of the binary surpasses the velocity dispersion of local stars (Merritt 2013) where  =  1  2 /( 1 +  2 ) is the reduced mass of the binary.
Beyond  h , binary hardening relies on collisionless repopulation of the losscone, until GW emission becomes important (at roughly milliparsec scales) and drives the MBHs to coalescence.Our numerical simulations are terminated when the separation between the MBHs reaches  BH , the BH softening length, which we have set to be smaller than  h to ensure that the evolution is followed self-consistently well into the losscone refilling phase.Fig. 14 shows that all models reproduce the same evolution for the separation between the MBHs all the way down to the BH softening length.The same is not true for the evolution of the orbital elements of the binary (see Fig. 15), which are characterised by a significant spread among the different random realisations at the lower resolution.The spread decreases with increasing resolution (from left to right), though much more slowly for the eccentricity.This appears to be more sensitive to perturbations and stochastic effects, as may be expected based on the results of Nasim et al. (2020); Gualandris et al. (2022); Rawlings et al. (2023).In particular, a few random realisations show binaries forming with very large eccentricity, despite a moderate initial orbital eccentricity for the galactic merger.At the highest resolution, model M7 shows a non negligible difference in eccentricity with respect to models M7a/b, implying that mass refinement is required even at  = 10 7 to reduce the scatter in eccentricity.Because the coalescence timescale of BHBs due to GW emission is strongly dependent on the orbital elements, and in particular on the eccentricity, we expect that these results will affect the time to coalescence.
We follow the late evolution of the binaries with a semi-analytical model that solves the coupled differential equations for the orbital elements under the effects of both stellar hardening and GW emission (Gualandris et al. 2022) where the first term represents the contribution from three-body stellar interactions and the second term models GW emission.
The contribution from the stellar interactions can be written as where  represents the time-dependent hardening rate of the binary while  is the eccentricity growth rate as defined by Quinlan (1996).
The contribution from GW emission can be derived for two point masses directly from the Einstein field equations in the non relativistic limit (Peters 1964) where   =  1 +  2 is the total mass of the binary.We compute the time-dependent numerical hardening rate in the hardening phase simulated with griffin and fit an exponential decay over time.The resulting fits are shown in Fig. 16.We observe some variation in hardening rate in the M6 models, which is significantly reduced in the M6b models.This can explain the differences in the evolution of the semi-major axis in the M6 models.The hardening rate is fully converged at the larger resolution, with no discernible difference past 3 Gyr.
The eccentricity growth rate is approximately constant over the simulated hardening phase in all cases.Therefore, for each model, we adopt its average over the hardening phase to be used in Eqs. 13 and 14, together with the fitted hardening rate.
Figure 14.The evolution of the distance  between the MBHs for all models at the lower (left) and higher (right) resolution.Overall, all models consistently reproduce the evolution of  over time.However the M6 and M6a models show small discrepancies at late times, which disappear in M6b and at  = 10 7 .The dashed horizontal line marks the softening length for the BH-BH and BH-particle interactions.The dotted lines in the rightmost panel mark the critical separations of the binary:   , approximately where the dynamical friction phase ends, and  ℎ , the hard binary separation.T (Gyr) Figure 15.The evolution of the semi-major axis  (top) and the eccentricity  (bottom) of the BHB with time, computed from the griffin simulations, for all models.The lower resolution models show a larger scatter in semi-major axis and especially in eccentricity with respect to the higher resolution models.At  = 10 6 , the refinement schemes prove to be effective at reducing the spread in the orbital elements.Model M6b is the most consistent with the higher resolution models.In addition, we find a clear trend of reduced scatter in the semi-major axis from left to right, with models M7/M7a/M7b showing excellent agreement.The same trend is seen in the eccentricity, though clearly with a larger scatter, due to stronger stochasticity.
We start the semi-analytical integration from time 2.6 Gyr for the M6 models and from time 2.5 Gyr for the M7 models, to ensure consistency with the -body integrations, and calculate the orbital elements until coalescence, as determined by the numerical solver.In order to reduce noise, we take an average of the semi-major axis and eccentricity values within ±50 Myr of the starting time.
The resulting evolution of the binary orbital elements is shown in Fig. 17 for all models.The predicted trajectories, shown by the solid lines, are consistent with the prior evolution from the griffin integrations, shown by the points.We note that stochasticity in the body evolution significantly affects the extrapolations to late times, as both three-body encounters and GW emission depend strongly on the semi-major axis and, especially, the eccentricity.Given that the hardening rate and eccentricity growth rate are very similar in all models except the M6s, the merger timescales are primarily set by the value of eccentricity at binary formation.For example, among the M6a models, the binary in model M6avi is the least eccentric and the one with the longest merger timescale.Overall, however, the refinement schemes are effective at reducing the scatter in eccentricity, with models M6b showing the smallest scatterand merger timescales more consistent with the M7 models.
Table 3 gives the time spent in each phase of evolution, namely dynamical friction, stellar hardening and GW emission.The latter two processes coexist in the hardening phase, we therefore define an arbitrary boundary between the two when the rate of change of semi-major axis due to each process is equal.All times are calculated as total evolutionary times since the onset of the galactic merger.The lowest resolution models give timescales in the range 3−13 Gyr, with the longest being ∼ 13.3 Gyr for model M6avi.All high resolution models give timescales under 11 Gyr, a desirable result in terms of detection rates for PTA.
We calculate the dispersion in semi-major axis, eccentricity (both  taken at time  = 2.6 Gyr) and merger timescales.These are shown in Fig. 18 for all models, as a function of the number of particles within 25 infl of the binary, as a measure of the effective resolution achieved in each model.We find that the dispersion in all quantities decreases with increasing resolution (from models M6, to M6a, M6b and M7), meaning that the spread is reduced by the schemes and by increasing particle number.However, the dispersion in eccentricity is reduced much more slowly than the one in semi-major axis, due to its intrinsically stochastic nature.Because of the strong dependence of the merger timescale on initial binary eccentricity, convergence in the total timescale is not achieved even in models M7a/b, and would require even higher resolution and/or more random realisations to verify.
Finally, we have explored the effect of the choice of a fixed softening on convergence by re-running one realisation each of the M6, M6a, M6b models with a mass dependent softening , according to  ℎ ∝  1/3 .We find that the measured dispersion in eccentricity and merger timescale is largely unchanged.

DISCUSSION AND CONCLUSIONS
We have presented a multi-mass refinement scheme to improve resolution at the centres of multi-component galaxy models, without impacting their stability and dynamical evolution.The initial model is over-seeded by a set amount, with both halo and bulge particles divided into radial zones.The defined central zone remains unchanged at the higher resolution produced by the over-seeding, whilst in the other radial zones we remove particles and correspondingly increase the masses of the remaining particles to compensate.In this way, the resulting galaxy model has a higher central resolution, necessary to Table 3. Characteristic timescales for the binaries throughout the evolution.From left to right: the simulation identifier, the time to reach   (an estimate of the time spent in the dynamical friction phase), the time spent in stellar hardening, the time spent in GW emission, the total time from the start of the simulation to BHB merger, the dispersion of the eccentricity at the binary formation time, the dispersion of the semi-major axis and the dispersion of the binary timescale.Figure 18.Dispersion of the BHB eccentricity (solid line), semi-major axis (dotted line) and merger timescale (dashed line) calculated from the random realisations of the  = 10 6 models and the three M7 models, shown as a function of the number of particles within 25 infl of the binary (i.e for models M6, M6a, M6b, and M7s from left to right).The dispersion in all quantities decreases with increasing resolution, proving that the mass refinement schemes are effective.However, convergence in semi-major axis is significantly faster than in eccentricity, resulting in a large dispersion in merger timescales even at the highest resolution.
accurately model MBH dynamics, at the same particle number, total stellar mass, and mass density profile.We introduced two schemes with varying degrees of over-seeding, the standard model 'a' and the more aggressive model 'b', with 'b' models having one more radial zone to compensate for the increased particle masses at outer radii.These implementations are applied to a set of isolated and merging multi-component galaxy models, each comprised of a dark matter halo, a stellar bulge, and a central MBH.We evolved a set of isolated and equal-mass merger models using the FMM code griffin, at the reference resolutions of  = 10 6 and  = 10 7 , including a set of 7 random realisations for each of the  = 10 6 resolution models.
The isolated models show that the multi-mass schemes are effective at increasing resolution in the central regions, and this is maintained over time.Both schemes are also effective at reducing the relaxation-driven expansion of the bulge particles observed in the lower resolution models, and expansion is nearly eliminated in the M6b implementation.We show that the introduction of mass refinement does not affect the stability of the models nor the subsequent evolution.Density profiles, merger remnant shape and anisotropy profiles remain unchanged.Binary formation occurs at ∼ 2.5 Gyr for all models.However, while all models are able to reproduce the same behaviour for the separation between the MBHs over time, the orbital elements are sensitive to resolution.
The M6 models show a large scatter in semi-major axis compared to mass refined or higher resolution models.Models M6a and M6b are consistent with each other and with the M7 models.The dispersion in semi-major axis calculated from the different random realisations shows convergence at a resolution of  = 10 6 , if a refinement scheme is applied.
The binary eccentricity, however, shows more variation than the semi-major axis, as it is more sensitive to perturbations and low- effects.A significant scatter is present at the lowest resolution of  = 10 6 and some variation still exists at the highest resolution of  = 10 7 .Models M7a and M7b, however, have very similar eccentricity evolution.
We have computed the time to coalescence for the BHBs in all models by extending the -body simulations with a semi-analytical calculation of the evolution of the orbital elements under the combined effects of stellar hardening and GW emission.We have determined the time-dependent hardening rate directly from the -body integrations and used an exponential fit in the models.We find that merger timescales can vary by a several Gyrs in the M6 models and by a few Gyr in the M7 models.Model M6b, with the most aggressive scheme, gives the most similar evolution and merger timescale to the higher resolution models.By contrast, the  = 10 7 models show a smaller spread in predictions, with models M7a and M7b coming within ∼ 1 Gyr of each other.
Variations in merger timescales have also been reported in Gualandris et al. ( 2022) at particle numbers of order  = 10 6 , while Nasim et al. (2021) show that a resolution in excess of  = 10 7 within the half-light radius is required to reduce the scatter in the merger timescale to ∼ 10%.We calculate that this corresponds to about  = 4 × 10 6 particles1 within 5  infl , which is a factor 3 larger than what we achieve in models M7a/b.Taking the results of Nasim et al. (2021) at face value, this implies that the effective resolution in models M7a/b results in a 20 − 30% error in the predicted merger timescale.Achieving a 10% error is hence possible through a combination of  (< 5 infl ) = 4 × 10 6 and application of the mass refinement scheme.
We conclude that mass-refinement schemes of the type described in this work are extremely effective at increasing resolution in the central regions, and should be combined with an appropriate particle number to reduce scatter in key quantities: a particle number of  (< 5 infl ) = 6 × 10 4 is sufficient (e.g.model M6b) to reach convergence in the evolution of the semi-major axis, while a particle number of  (< 5 infl ) = 10 6 is required to approach convergence in the evolution of the eccentricity and in the time to coalescence.In order to achieve a 10% error on the merger timescale, a resolution  (< 5 infl ) = 4 × 10 6 is required.This is due to the strong dependence of both stellar hardening and GW emission on eccentricity.
We note that our galaxies are all set on initial orbits of  = 0.7, i.e. moderately eccentric.Gualandris et al. (2022) have shown that binaries with lower eccentricities exhibit a greater spread in merger timescales.We therefore expect somewhat lower variations for more eccentric orbits ( ∼ 0.9) (e.g.Nasim et al. 2020Nasim et al. , 2021;;Mannerkoski et al. 2022;Rawlings et al. 2023).
Merger timescales can span a range from a few to several Gyr, making it likely that a third intruder will interact with the merger remnant before the BHB has reached coalescence.Subsequent mergers and triple galaxy encounters need to be modelled to make reliable predictions of merger timescales for PTA.A study of the merger tree selected from Illustris-TNG-300-1 presented in this work is underway.A refinement scheme of the type presented here is particularly beneficial for such studies given the need for high accuracy.While in this work we have only tested the mass refinement schemes on equal mass mergers of spherical galaxies, preliminary results show that they can be reliably applied to different morphologies (axisymmetric or triaxial models) and unequal mass mergers.This opens up the possibility of adopting mass refinement schemes as a best practice in any -body simulation where large resolution is key, including cosmological simulations.

Figure 1 .
Figure1.Radial distribution of the bulge particles for all isolated models with  = 10 6 particles (top panels) and  = 10 7 (bottom panels) at two characteristic times:  = 0 Gyr (left) and  = 2.8 Gyr (right), corresponding to a late time in the evolution.The grey dotted line represents the 5 infl radius.A significant increase in the particle number < 10 kpc is observed for the models with the multi-mass scheme applied, and no significant expansion occurs over time for these models.

Figure 2 .
Figure 2. Radial distribution of the bulge particles for all merger models with  = 10 6 particles (left) and  = 10 7 (right), both at the later time  = 2.8 Gyr past the end of the merger process.The grey dotted line represents the 5 infl radius.A higher central resolution is maintained for the models with the multimass scheme applied post-merger, most significantly in the M6 models.

Figure 3 .
Figure 3. Spatial density profiles of the isolated models with  = 10 6 particles (top panels) and  = 10 7 particles (bottom panels) at time  = 0.0 Gyr (left) and at the later time  = 2.8 Gyr (right).The dashed vertical line represents the star-star softening length.The profiles show the same trends as the bulge particle distributions, with all the multi-mass models showing stability at the innermost radii.

Figure 4 .
Figure 4. Spatial density profiles of the merger models with  = 10 6 particles (left) and  = 10 7 (right) at the later time  = 2.8 Gyr.The dashed vertical line represents the star-star softening length.

Figure 5 .
Figure 5. Axis ratios of the isolated models as a function of radius for  = 10 6 particles (top) and  = 10 7 particles (bottom) at time  = 0.0 Gyr (left) and at the later time  = 2.8 Gyr (right): / (solid lines) and / (dashed lines).The spherical models maintain their shape over time even when the mass refinement schemes are applied.

Figure 9 .
Figure 9. Stellar velocity anisotropy profile of the isolated models with  = 10 6 particles (top panels) and  = 10 7 particles (bottom panels) at time  = 0.0 Gyr (left) and at later time  = 2.8 Gyr (right).All models maintain an isotropic profile over time, except for I6 which shows a moderate tangential anisotropy at late times, likely due to low resolution.

Figure 10 .
Figure 10.Stellar velocity anisotropy profile of the merger models with  = 10 6 particles (left) and  = 10 7 particles (right) at the later time  = 2.8 Gyr.A modest discrepancy between the models is observed at the lower resolution, while at the higher resolution the models are indistinguishable.

Figure 12 .
Figure12.Lagrangian radii for bulge (bottom four curves) and halo (top four curves) particles for all merger models with  = 10 6 (left) and  = 10 7 (right) across the simulation time.The curves correspond to radii containing mass fractions of 3%, 10%, 13%, and 20% of the total bulge mass.The strong disturbance at  ∼ 2.5 Gyr marks the completion of the merger process.The M6 model shows a significant expansion at the smallest radii, similarly to I6.

Figure 13 .
Figure 13.Dynamical friction timescale of the halo particles for all models at time  = 2.8 Gyr, for the isolated models (left) and the post-merger models (right).The dashed horizontal lines mark the duration of the simulations (black).The points mark the centres of each radial bin in the multi-mass schemes.The particles are affected by dynamical friction at  ≲ 2 kpc in I6, and at  ≲ 1 kpc in the multi-mass models.In model I7, dynamical friction only affects particles within the central ≲ 500 pc, which is reduced to the central ≲ 200 pc in the multi-mass models.The multi-mass models show a downturn in  df at the outer radii, where the particle masses are most increased.The merger models show similar trends to the isolated models.

Figure 16 .
Figure16.Evolution of the hardening rate, computed numerically as in Eq. 15, for all models.The scatter decreases with increasing resolution, from left to right.

Figure 17 .
Figure17.The evolution of the semi-major axis  (top) and the eccentricity  (bottom) of the binary with time, for all models.The evolution simulated with griffin is shown with points, whilst the subsequent evolution modelled with the semi-analytic model including GWs is shown with solid lines.

Table 1 .
The parameters of the refinement scheme: scheme identifier, zone number  (for both bulge and halo), percentage of total particle number within zone , mass multiplication factor , and start and end radius ( 0 - 1 , kpc) of each zone, for bulge and halo.

Table 2 .
Properties of the adopted galaxy models: model identifier, type of model (isolated or merger), total particle number, number of halo particles, number of bulge particles, number of particles within the 5 infl in the initial models, minimum halo star to MBH mass ratio, minimum bulge star to MBH mass ratio, and central resolution multiplication factor for the models with mass refinement.