Resolving massive black hole binaries evolution via adaptive particle-splitting

The study of the interaction of a massive black hole binary with its gaseous environment is crucial in order to be able to predict merger rates and possible electromagnetic counterparts of gravitational wave signals. The evolution of the binary semi-major axis resulting from this interaction has been recently debated, and a clear consensus is still missing, also because of several numerical limitations, i.e. fixed orbit binaries or lack of resolution inside the cavity carved by the binary in its circumbinary disc. Using on-the-fly particle-splitting in the 3D meshless code gizmo, we achieve hyper-Lagrangian resolution, which allows us to properly resolve the dynamics inside the cavity, and in particular for the first time the discs that form around the two components of a live binary surrounded by a locally isothermal gaseous circumbinary disc. We show that the binary orbit decays with time for very cold and very warm discs and that the result of the interaction in the intermediate regime is strongly in uenced by the disc viscosity as this essentially regulates the fraction of mass contained in the discs around the binary components as well as the fraction that is accreted by the binary. We find the balance between these two quantities to determine whether the binary semi-major axis decreases with time.


INTRODUCTION
Massive black holes (MBHs) are expected to reside in the nuclei of (virtually all) massive galaxies (e.g. Kormendy & Ho 2013). When two massive galaxies merge the MBHs in their centers migrate towards the center of the newly formed galaxy owing to the effect of dynamical friction (Chandrasekhar 1943) against the background of stars and gas. The two MBHs are expected to form a bound massive black hole binary (MBHB) at parsec scales. At these scales dynamical friction becomes inefficient and further evolution of the binary orbit towards coalescence requires a mechanism to extract angular momentum and energy from the binary. The main mechanisms proposed in the literature are three-body scattering of stars intersecting the binary orbit (e.g., Quinlan 1996;Sesana et al. 2007) or the interaction with a circumbinary gaseous disc (Escala et al. 2005;Dotti et al. 2007;Cuadra et al. 2009) depending on the binary sep-aration scale (see Bortolas et al. 2021 for details on the competition between stellar and gaseous interaction).
Each of the two massive galaxies brings with it a significant amount of gas that sinks to the center forming a circumbinary accretion disc. This gaseous disc might facilitate the MBHB merger and give rise to electromagnetic counterparts of the gravitational wave (GW) emission (Begelman et al. 1980;Armitage & Natarajan 2002;Lodato et al. 2009). A coplanar accretion disc can extend down to a few times the binary separation (Artymowicz & Lubow 1994). The disc orbits resonate with the binary orbit at discrete locations (outer Lindblad resonances), leading to the exchange of angular momentum between the disc and the binary (Lynden-Bell & Kalnajs 1972;Lin & Papaloizou 1986). The magnitude of the resonant torques depends on the binary potential, i.e. its mass ratio and eccentricity, and is proportional to the disc surface density at the resonance locations (Goldreich & Tremaine 1979). Therefore, the amount of angular momentum transferred from the binary to the disc at the resonances depends on the disc properties as well.
Very early numerical simulations (Artymowicz & Lubow 1994, 1996 investigated the interaction of a bi-nary with its gaseous circumbinary disc finding that only a small amount of material is able to enter the cavity carved by the binary and to accrete onto the binary components. The main finding of these works is that the binary semi-major axis decreases with time owing to the interaction with the disc. This picture has been recently challenged by a few works (Muñoz et al. 2019;Duffell et al. 2020;Muñoz et al. 2020) employing 2D (and one 3D simulation, see Moody et al. 2019) static or moving-mesh grid numerical simulations with fixed binary orbits. In particular, these studies found that the secular angular momentum transfer onto the binary is strongly positive within the range of binary and disc parameters explored. More recently Tiede et al. (2020) found, using the same numerical techniques, that the sign of the torque exerted by the disc onto the binary depends on the disc temperature, i.e. on its aspect ratio H/R, for locally isothermal discs. In particular, they found discs with H/R 0.04 to shrink the binary. Using 3D smoothed particle hydrodynamics (SPH) simulations of locally isothermal discs, Heath & Nixon (2020) found instead the threshold value for binary expansion to be H/R 0.2. Other works that employed SPH simulations in the regime where the disc self-gravity cannot be neglected, and the disc temperature changes with time, found binary shrinking as a result of the interaction with massive discs regardless of the disc temperature (Cuadra et al. 2009;Roedig et al. 2012;Franchini et al. 2021). The discrepancy in the results inferred using different numerical techniques has been argued to originate from the lack of numerical resolution in SPH simulations which, in constrast with grid-based ones, were unable to properly resolve the dynamics of the gas streams entering the cavity, artificially suppressing the positive torque associated to such gas that forces the binary to expand (Muñoz et al. 2019).
In this letter, we employ the public version of the code gizmo (Hopkins 2015) in its mesh-less finite mass (MFM) mode, coupled with adaptive particle-splitting for numerical refinement of the gas dynamics inside the disc cavity, in the aim at accurately measuring the gravitational and accretion torques that originate from the discs that form around the binary components (also called mini-discs) onto the binary itself. 1 In Section 2 we describe the details of the numerical method and present the results of our simulations in Section 3. We discuss the relevance of the parameters explored in Section 4 and finally draw our conclusions in Section 5.

NUMERICAL SETUP
The initial conditions for this work consist of a live binary surrounded by a circumbinary gaseous disc. The 3D distribution of the 10 6 equal mass gas particles sampling the disc and the initial orbit of the two sink particles of the binary are generated using the SPH code phantom (Price et al. 2017a). The equal mass binary has an initial mass M = M 1 + M 2 = 1 and a separation a = 1. The circumbinary disc initially extends from R in = 2a to R out = 10a, with a fixed aspect ratio H/R = 0.1. The gas is described by a locally isothermal equation of state with the sound speed c s defined by Eq. (4) in Farris et al. (2014), and a surface density profile scaling as Σ ∝ R −3/2 , normalised to get a total mass M disc = 0.1M .
The simulations are performed with the MFM method in gizmo, keeping the binary live, i.e. its orbit is allowed to change during the evolution. The MFM is a mesh-free Lagrangian approach which encapsulates the advantages of both grid-based and particle-based codes. In particular, by partitioning the domain with discrete tracers (particles/cells), and solving the integral form of the fluid dynamics equations via a finite-volume Godunov method, this numerical technique allows to obtain intrinsic adaptivity in resolution while minimizing the numerical advection of angular momentum through the 'cells' and also ensures shock capturing without the need of the extra artificial viscosity terms commonly used in SPH codes. We here also include the effect of gas viscosity entering the Navier-Stokes fluid equations as described in Hopkins (2016), assuming a shear viscosity in the disc ν = αc s H parametrised using a viscosity parameter α = 0.1 (Shakura & Sunyaev 1973), and no bulk viscosity.

Gas accretion and gravitational forces
Every time a gas particle approaches one of the sinks, entering its sink radius r sink , we flag it as eligible to be accreted. During the accretion event, conservation of mass, and linear and angular momentum are ensured, in the same way it is done in the phantom code (Bate et al. 1995). 2 In this work, we neglect the disc self- Figure 1. Column density plots of the circumbinary disc and the discs around the binary components (shown by the white circles). The view is of the x-y plane (i.e. the binary orbital plane) and the density has been integrated through z. The color scale spans about three orders of magnitude in density and is the same for all the plots. The green dashed circle corresponds to the initial binary orbit, and the cyan dashed circle to the location of the strongest Lindblad resonance, i.e. 2.08a. gravity, and only include the gravitational interaction between i) the two sink particles, and ii) sink and gas particles (see Franchini et al. 2021, for a discussion about the role of self-gravity). In order to get a high accuracy in the dynamical evolution, we suitably modify the code as follows. First, we include the gravitational acceleration exerted by the gas particles onto the two black holes via direct summation, which guarantees the most accurate estimation of the acceleration on both sinks. Secondly, we introduce a new timestep criterion for the sink particles that forces both of them to evolve on the shortest common timestep. Finally, we re-set the centre of mass and the centre of mass velocity of the entire system to the origin at each course timestep, in order to avoid the build-up of small conservation errors in the linear momentum that, over the long integration times we consider, might displace the binary from the center. Note that this does not influence the dynamics of the system, since it simply corresponds to a rigid shift of the positions and velocity of the centre of mass of all sink and gas particles to the origin.

Particle splitting
In order to increase the resolution when and where necessary, without globally slowing down our simulations, we employ an on-the-fly adaptive particle-splitting approach, which is similar in spirit to the adaptive mesh refinement of grid-based codes. Such technique represents a natural generalisation of the refinement/derefinement scheme in arepo (Springel 2010) and gizmo (Hopkins 2015) to maintain an almost constant mass per cell during the simulations when a finite-volume scheme is employed. In this work, we exploit the particle splitting algorithm in gizmo, splitting gas particles that enter a sphere of radius r ref centred on the centre of mass of the binary-disc system. In particular, particles at radii r < r ref are split when their mass m χ ref m max , with m max a scale mass and χ ref the refinement factor, a suitably defined function, which we choose of the form Here, ξ = r/r ref , ξ min = r min /r ref with r min the radius at which the maximum resolution is reached, χ min ref = 1/32 and the coefficients are defined as In particular, we choose r min = 0.75a, that guarantees that the flow around the two components of the binary is always at the maximum resolution, and r ref = 4a, that allows to properly resolve the edge of the cavity and the dynamics of the gas streaming inwards.
When a particle is flagged for refinement, two children are created, each with mass 0.5 m, and located at a distance dr = min(0.25 h, 0.35 d nbg ), where h is the particle kernel size and d nbg is the distance to the nearest neighbor, from the position of the parent and on opposite sides along a random direction. The numerical factors are necessary to minimize perturbations to hydrodynamic quantities and avoid overlapping of fluid elements (Anglés-Alcázar et al. 2020). All the remaining properties are instead directly inherited by the children. This approach allows the almost perfect conservation of momentum, angular momentum and energy throughout the simulation. The density is re-calculated immediately after the splitting using the standard approach. In addition to splitting, particles are also allowed to merge when their mass is significantly below the resolution requirements, i.e. m < χ ref m min , with m min the scale mass for merger. In our numerical simulations, we set m max equal to the initial mass of the particles in the disc, and m min to 1/1000th of it, which translates in neglecting merging, although we also explored different values and the inclusion of mergers, finding negligible differences.

RESULTS
We now describe our main results, exploring how different parameters and resolution requirements affect our conclusions.

Lagrangian resolution
As a first check of our numerical and physical setup, we performed a simulation using standard Lagrangian resolution, i.e. switching off particle splitting, and a large sink radius r sink = 0.2, in order to obtain a fair comparison with the results presented in Heath & Nixon (2020). We find a general agreement between our simulation with H/R = 0.1 and that represented by the red dashed curve in Fig. 4 of Heath & Nixon (2020), with the binary shrinking over time. We also performed the same test with a thicker disc, i.e. H/R = 0.2, in this case finding binary expansion, as already found by Heath & Nixon (2020).
We note that the rate at which the binary shrinks (or expands) is in general slightly higher in our runs than in Heath & Nixon (2020), likely because of the different treatment of hydrodynamics (e.g. different numerical treatment of shocks) in the two codes and the slightly different equation of state employed for the gas. However, a precise comparison between the codes is beyond the scope of this work.

Hyper-Lagrangian resolution
The use of pure Lagrangian resolution has been argued to possibly lead to a less accurate treatment of the interaction between a binary and its circumbinary disc, because of the small amount of particles entering the cavity, resulting in a poor resolution around the two MBHs. In order to overcome this potential limitation, we employ here on-the-fly particle splitting which allows us to significantly increase the resolution of our simulations in the cavity at a moderate computational cost.
We here show and discuss the result of our fiducial case with particle splitting (implemented according to Eq. 1), disc aspect ratio H/R = 0.1, constant throughout the disc, and viscosity parameter α = 0.1. Compared to the purely Lagrangian tests, the sink size is now reduced to r sink = 0.05, to ensure that the gas dynamics around each MBH is properly resolved. Figure 1 shows the column density plots of the circumbinary disc together with the resolved discs around the two binary components (white circles) obtained with particle splitting at t = 500, 1500, 2000, 2500 P b . We can clearly resolve the dynamics inside the cavity carved by the live binary using our hyper-Lagrangian approach. The cavity does not remain circular, as already observed in previous studies (Ragusa et al. 2016;Muñoz et al. 2019;Heath & Nixon 2020), and the formation of shocks at the cavity edge is also expected due to the part of the streams that gets flown back to the disc inner edge after entering the cavity.
We test the conservation of angular momentum in the simulations with and without particle splitting, following the procedure outlined in Franchini et al. (2021) (see also Muñoz et al. 2019), in order to ensure that the particle splitting algorithm does not significantly affect the conservation. We divide the contribution of the disc to the change in the binary angular momentum in a gravitational and an accretion part, the latter including also the fraction of the accreted angular momentum that is converted into the spin of the two MBHs.
The left panel of Figure 2 shows the contributions to the binary angular momentum change as computed directly from our fiducial simulation. In particular, from the comparison between the black straight line (i.e. sum of gravitational and accretion contribution) and the red line that shows the live binary angular momentum change, it is evident that the angular momentum is conserved in our fiducial run almost exactly. The small discrepancy is due to our approximate calculation of the gravitational torque between two subsequent snapshots instead of calculating it at each timestep. As a further consistency check, we also compared the same quantities in a simulation without particle splitting (keeping r sink = 0.05), finding the impact of the refinement scheme on the conservation to be negligible.
Since the simulation conserves angular momentum, we can write the evolution of the binary semi-major axis aṡ where L z = µ GM a(1 − e 2 ) = L z,grav + L z,acc is the z component of the binary angular momentum, µ is the reduced mass and e is the binary eccentricity. The first term is always positive and represents the accretion of angular momentum onto the binary (dotted line in Fig.  2). The second term is given by the sum of the positive contribution of the discs around the two components and the negative contribution of the circumbinary disc and corresponds to the dashed line in Fig. 2. The terms due to the accretion of mass (i.e. third and fourth terms) are additional negative contributions to the semi-major axis evolution. The last term comes from the binary eccentricity evolution and we find it to be negligible in our runs. Therefore we essentially have two positive terms whose effect is to drive the binary apart and two negative contribution that remove angular momentum from the binary driving it towards merger. In more physical terms, the angular momentum change due to the gravitational and accretion torques translates into a change of the different elements of the system, i.e. the masses and orbital elements. Since the eccentricity contribution is negligible in our runs, the evolution of the binary separation depends on the fraction of the angular momentum exchanged with the disc goes into the mass accretion terms. The right panel of Figure 2 shows the evolution of the live binary semi-major axis with time, directly calculated from our fiducial run using energy and angular momentum conservation. The binary semi-major axis decreases within the first ∼100 binary orbits, due to a transient phase as the disc adjusts to an equilibrium configuration. After the first 100 orbits, the shape (and size) of the cavity changes significantly as the streams are flown back to the circumbinary disc. This causes the amount of material at the 2:1 Lindblad resonance (represented by the cyan dashed circle in Fig. 1) to decrease by a factor of ∼5-10 weakening the strength of the resonance. As a consequence, the positive torque provided by the discs around the binary components, which is able to overcome the negative contribution of the circumbinary disc, leads to the increase of the binary separation. 3 We find the expansion phase to be very mild and the semi-major axis to increase at a progressively slower pace. This might be due to the progressive decrease of the disc mass owing to accretion, which drops to M d ≈ 0.04M after 4000 orbits. We note that this phase appears to be in broad agreement with the findings of grid code simulations (Muñoz et al. 2019;Duffell et al. 2020) indicating that indeed resolving the dynamics inside the cavity is important in order to properly investigate the binary-disc interaction.

The impact of spatial/mass resolution and of live/fixed binary orbits
Although the extremely high resolution achieved owing to the hyper-Lagrangian refinement enabled us to better resolve the gas flowing into the cavity, it might be argued that the resolution achieved is insufficient to properly measure torques in the cavity. We have therefore estimated the average inter-particle separation of the gas in the discs surrounding the binary components in our fiducial simulation, finding ∆r 0.015a. This corresponds to, e.g., the mid-resolution run in Tiede et al. (2020), already confirming the good accuracy of our measured torques. Nonetheless, we have further investigated the impact of resolution inside the cavity by decreasing the minimum refinement factor down to χ min ref = 1/64, i.e. by a factor 2. In order to avoid the initial transient in the orbital evolution, we have restarted our fiducial run from the 500th orbit, carefully emptying the cavity to prevent numerical instabilities arising from a sudden increase in resolution near the binary components. Note that the amount of mass removed by this process is negligible compared to the disc mass, and the removal does not affect the subsequent evolu-tion, since the discs around the two components reform in less than one orbit. We find the structure and density of these discs to remain essentially unaltered and unaffected by the increase in resolution. The blue lines in the left panel inset of Figure 2 show the gravitational and accretion torque measured in this case. We find the accretion torque to remain the same while the gravitational torque is slightly lower in this high resolution run. This reflects on the evolution of the binary semi-major axis (see the blue line in the right panel inset in Fig. 2) which increases at a slightly slower pace.
The majority of previous studies on circumbinary discs is based on simulations in which the binary is kept on a fixed orbit. While this assumption might be reasonable when the change in binary properties is small, a proper analysis of its impact is in order. For this reason, we have also performed a simulation identical to our fiducial run, in which we kept the binary orbit fixed, as in previous studies (Muñoz et al. 2019;Duffell et al. 2020;Moody et al. 2019). The comparison between the fixed and live binary orbit run, in terms of gravitational and accretion torques, is shown by the green lines in the left panel inset of Figure 2. We started the simulation from the 1000th binary orbit of our fiducial live binary run since, by this time, the disc is in a quasi-steady state. We find that the accretion and gravitational torques do not vary significantly compared to our fiducial live binary simulation. We note however that the surface density of the discs that form around the two binary components is slightly lower in the simulation with fixed binary orbit. This does not cause the torques to be significantly different over a few hundred orbits but it is an effect that will be further analysed in a future work.

Disc temperature
In order to determine whether there is a critical condition that discriminates between binary shrinking and expansion, we also explored different disc temperatures, keeping the small sink radius r sink = 0.05 and using particle splitting.
First, we decreased the disc aspect ratio to H/R = 0.03, i.e. within the regime where both grid-code (Tiede et al. 2020) and SPH (Heath & Nixon 2020) simulations find binary shrinking, and we found our simulations with particle splitting to be in agreement with previous studies. We note that, given the low disc temperature, the mass in the discs around the binary components is very low regardless of the refinement and this reduces the positive contribution to the gravitational torque. Heath & Nixon (2020) found, using pure SPH simulations without particle splitting, that circumbinary discs with aspect ratios H/R ≥ 0.2 cause the binary orbit to expand owing to the large amount of material that is able to enter the cavity and to accrete onto the binary. For completeness, and a further check for our fiducial run, we also explored the case of a thicker disc with H/ R=0.2. 4 The left panel of Figure 3 shows the column density plots of the circumbinary disc together with the resolved discs around the two binary components (white circles) obtained with particle splitting for a thick accretion disc. The density scale is the same as in Fig. 1. We find that the material that enters the cavity is significantly larger compared to the thinner disc simulation.
We can see from the right panel of Figure 3 that the binary semi-major axis decreases with time, after undergoing a very brief and mild expansion within the first 20 orbits, even though the disc thickness is in the regime Heath & Nixon (2020) found for expansion. Despite the gravitational and accretion torques being positive throughout the duration of the simulation, we find the negative contribution of the accretion of mass (and not angular momentum) to dominate the evolution of the binary.

DISCUSSION
Interestingly we find the binary to shrink for thin and very thick discs and that there is a regime of intermediate disc thickness where the negative and positive torques balance out, leading to a phase of possible binary expansion.
We compared our simulation with H/R = 0.1 with and without particle splitting, and found that having a small sink radius without resolving the discs that form around the binary components might result in an overestimate of the (positive) gravitational torque and this in turn leads to a much faster binary expansion phase. This lasts until the accretion of mass onto the binary components becomes large enough to overcome the positive torque and is very sensitive to the resolution inside the cavity. In our fiducial case the expansion is significantly milder and slower, however it seems to last for much longer due to the positive torque from the discs around the two components being strong enough to balance the negative torque from the accretion of mass. In our low resolution simulation we find the binary semimajor axis to decrease again after ∼ 800 orbits. We further explored the effect of resolution by increasing the resolution inside the cavity by a factor 2 with re- spect to our fiducial run. We find that our results are overall not affected by the numerical resolution.
The balance between the positive and negative torques is essentially regulated by the disc viscosity, which determines the amount of material that enters the cavity, forms the discs and then accretes onto the binary components. Therefore different treatments and different α values are very likely to change the binary evolution. For instance, decreasing the value of α by a factor two in our fiducial run with aspect ratio H/R = 0.1 leads to the binary semi-major axis decreasing with time without transiting to the expansion phase.
Since the MFM method has been proven to give more accurate results on standard tests compared to the pure SPH approach (Hopkins 2015), even with a small number of neighbours, we reran our fiducial simulation with particle splitting, but employing only 32 neighbours instead of the 58 used in the rest of this work. We find the expansion phase to be significantly milder.
We also note that, by 3500 P b , the binary in our fiducial run has already accreted ∼ 60% of the initial disc mass and the semi-major axis is still increasing, even though at a progressively slower pace. In order to understand whether this phase is just a transient, in a future work we will implement gas particle injection to provide a constant source of gas into the circumbinary disc to feed the discs inside the cavity.
Finally, we do not find binary expansion in the simulation presented in Section 3.4 with H/R = 0.2 using particle splitting with a small accretion radius. The reason for this change resides in the gas distribution in the cavity. The gas distribution inside the discs around the binary components in the thicker disc case is uniform and therefore removing the inner regions by enlarging the accretion radius leads to a higher specific angular momentum being transferred to the binary during accretion, which in turn causes the binary to expand. Although whether there exist a threshold for binary expansion in terms of the disc aspect ratio still remains to be determined, such large disc thicknesses are well beyond the typical aspect ratio we may expect for circumbinary discs surrounding massive black hole binaries, if their structure resembles that of AGN discs (Collin-Souffrin & Dumont 1990).

CONCLUSIONS
In this work, we use hyper-Lagrangian refinement to resolve for the first time the gas dynamics inside the cavity carved by a live binary in its circumbinary disc, using 3D hydrodynamical simulations employing the mesh-less finite mass numerical algorithm of the code gizmo. In particular, resolving the discs that form around the binary components allows us to better measure the torques they exert on the binary. We find the gravitational torques by these discs to be positive, in broad agreement with the results of previous grid-based simulations (Muñoz et al. 2019;Duffell et al. 2020). However, we find these not to be strong enough to overcome the negative circumbinary disc torque and the negative angular momentum variation due to the accretion of mass both in the case of thin and very thick discs.
We therefore conclude that the binary semi-major axis typically decreases with time as a result of the interaction with the gaseous accretion disc, if the disc is thin, but even if the disc is very geometrically thick. Interestingly, we find an intermediate regime where the positive torque exerted by the discs around the binary components is able to counter balance the negative torques, leading to a possibly temporary phase of expansion. However, we caution that this is very sensitive to the viscosity treatment, and therefore worth investigating more in detail in a follow-up paper.