Mergers of black hole binaries driven by misaligned circumbinary discs

With hydrodynamical simulations we examine the evolution of a highly misaligned circumbinary disc around a black hole binary including the effects of general relativity. We show that a disc mass of just a few percent of the binary mass can significantly increase the binary eccentricity through von-Zeipel--Kozai-Lidov (ZKL) like oscillations provided that the disc lifetime is longer than the ZKL oscillation timescale. The disc begins as a relatively narrow ring of material far from the binary and spreads radially. When the binary becomes highly eccentric, disc breaking forms an inner disc ring that quickly aligns to polar. The polar ring drives fast retrograde apsidal precession of the binary that weakens the ZKL effect. This allows the binary eccentricity to remain at a high level and may significantly shorten the black hole merger time. The mechanism requires the initial disc inclination relative to the binary to be closer to retrograde than to prograde.


INTRODUCTION
A body orbiting a black hole binary in an inclined orbit can cause oscillations of the binary eccentricity (Farago & Laskar 2010;Aly et al. 2015;Naoz 2016).These von-Zeipel-Kozai-Lidov oscillations (ZKL, von Zeipel 1910;Kozai 1962;Lidov 1962) can occur even in the case that the outer body is much less massive than the inner binary (Chen et al. 2019;Hamers 2021).Increased eccentricity drives the merger of the black holes more rapidly as a result of energy loss from gravitational radiation (Peters 1964;Antonini & Perets 2012;Silsbee & Tremaine 2017).For a very high mass third body, the maximum eccentricity is achieved when the third body is in a polar orbit relative to the binary orbit (Liu & Lai 2018;Liu et al. 2019a,b).In this configuration, the angular momentum vector of the third body orbit is aligned to the eccentricity vector of the binary.However, recently we found that for a lower mass third body, a mass of just a few percent of the binary mass, the largest binary eccentricity growth occurs for inclinations that are initially closer to retrograde than to prograde (Lepp et al. 2023a).
A black hole binary may accrete material in a chaotic fashion at random inclinations relative to the binary orbit.This is how supermassive black holes at the centres of galaxies are thought to grow (King & Pringle 2006;Nixon et al. 2012;King & Nixon 2015).Similarly chaotic accretion on to young binary stars is seen in simulations of star forming regions (Offner et al. 2010;Bate 2012Bate , 2018) ) and observations show that misaligned circumbinary discs are common (e.g.Chiang & Murray-Clay 2004;Köhler 2011;Brinch et al. 2016;Czekala et al. 2019).
The orbit of a massive body around a binary can undergo compli-★ E-mail: rebecca.martin@unlv.educated behaviour relative to the binary orbit (Verrier & Evans 2009;Farago & Laskar 2010;Doolin & Blundell 2011;Naoz 2016;Chen et al. 2019).A misaligned orbit undergoes nodal precession and the binary eccentricity oscillates on the same timescale.The nodal precession may be centered on the binary angular momentum vector for low initial inclination (circulating orbits) or centered on a stationary inclination for higher initial inclination (librating orbits).For a test particle, without the effects of general relativity (GR), the stationary inclination is at 90 • , a polar orbit in which the angular momentum vector of the outer body is aligned to the eccentricity vector of the inner binary.GR drives prograde apsidal precession of the binary (e.g.Ford et al. 2000;Naoz et al. 2017;Zanardi et al. 2018) that increases the stationary inclination (Lepp et al. 2022;Zanardi et al. 2023) while the mass of the particle lowers the stationary inclination (Chen et al. 2019;Martin & Lubow 2019).The maximum binary eccentricity is achieved when the third object is in a librating orbit (Chen et al. 2019).
At each radius in the circumbinary disc, the material feels the same torque as a point mass at that radius.However, the rings of a disc communicate with each other in a wave-like manner (if the disc aspect ratio, /, is larger than the Shakura & Sunyaev (1973) viscosity parameter, ) or through viscosity (if / ≪ ).If the communication is faster than the global precession timescale, then the disc can precess as a solid body (e.g.Lubow & Martin 2018) otherwise the disc can become warped or even break (Larwood et al. 1996;Nixon et al. 2013;Nixon & King 2016).Two parts of a broken disc can undergo nodal precession independently on different timescales.The evolution of the third body depends upon its mass and its orbital radius.Therefore, the evolution of a circumbinary disc depends upon the distribution of mass within the disc and breaking may lead to significantly different behaviour.
In this work, we show that a massive circumbinary disc can drive ZKL like oscillations of the binary eccentricity of an equal mass binary.This mechanism has been previously seen in the case of an extreme mass ratio binary composed of a planet and a star (Terquem & Ajmia 2010).In that case, the eccentricity of the planet orbit oscillates.In Section 2 we consider hydrodynamic simulations of a circumbinary disc around an equal mass black hole binary and body simulations for comparison.We show that a circumbinary disc with a mass of only a few percent of the binary mass can significantly increase the binary eccentricity and speed up the merger time.We find that disc breaking can lead to a highly eccentric binary for a long timescale rather than an oscillating eccentricity as is the case if the outer body is a particle rather than a disc (Lepp et al. 2023a).In Section 3 we consider the application to a range of black hole masses and we draw our conclusions in Section 4.

CIRCUMBINARY DISC AND PARTICLE COMPARISON
We consider a binary composed of two equal mass black holes with masses  1 =  2 = 30 M ⊙ and total mass  b =  1 +  2 .They are initially in an orbit with semi-major axis of  b = 10 R ⊙ (e.g. de Mink & King 2017;Martin et al. 2018) with eccentricity  b = 0.2.In this section we first describe the initial conditions for the circumbinary disc hydrodynamic simulation.We then show an equivalent massive particle -body simulation and then we compare this to the disc simulation.Finally we consider a four-body simulation for comparison to the circumbinary disc simulation in which the disc breaks into two components.

Circumbinary disc set-up
We use the smoothed particle hydrodynamics (SPH Bate et al. 1995) code Phantom (Price & Federrath 2010;Price et al. 2018) to model the evolution of a circumbinary disc.This code has been extensively tested for circumbinary discs (Nixon et al. 2012;Facchini et al. 2013;Aly et al. 2018;Cuello & Giuppone 2019;Abod et al. 2022;Smallwood et al. 2019).The binary is modelled with sink particles that accrete the mass and angular momentum of any SPH particles that move inside their sink radius (Bate et al. 1995).The sink radius for each black hole is 2.5 R ⊙ .We include a modification to the sinksink acceleration to model the effects of GR (Nelson 2000;Childs et al. 2023a) that drives prograde apsidal precession of the binary orbit (e.g.Zanardi et al. 2018).
The circumbinary disc is initially relatively far from the inner tidal truncation radius around the binary that is around 2 − 3  b for a coplanar disc (Artymowicz & Lubow 1994) but smaller for a polar disc (Franchini et al. 2019).The initial total disc mass is  d = 0.033  b and is composed of 500, 000 particles that are distributed in a relatively narrow ring with surface density profile Σ ∝  −3/2 between  in = 6.5  b and  out = 10  b .The disc is free to spread inwards and outwards.The disc is initially inclined to the binary orbit by  = 130 • with longitude of ascending node relative to the binary eccentricity vector of Ω = 90 • .We also considered larger initial tilts, of 170 • and 150 • since these can drive even larger eccentricities for a particle (Lepp et al. 2023a).However, we did not find significant eccentricity growth with  = 170 • .With  = 150 • we find similar behaviour to  = 130 • with only slightly lower eccentricity growth.We discuss this further in Section 3.
The disc is locally isothermal with sound speed  s ∝  −3/4 and disc aspect ratio / = 0.01 at  = 6  b .This choice allows the Shakura & Sunyaev (1973) viscosity parameter  and the shellaveraged smoothing length per scale height ⟨ℎ⟩ / to be constant over the disc.We take the viscosity parameter to be  = 0.1, typical for fully ionised discs (e.g.Martin et al. 2019).The viscosity is implemented by adapting the SPH numerical viscosity with  AV = 1.35 and  AV = 2 (Lodato & Price 2010).The disc is resolved with ⟨ℎ⟩ / ≈ 0.7 initially.

Three-body simulation
For comparison to the disc evolution, we first consider the evolution of a massive particle orbiting the binary.We use the -body code Rebound (Rein 2012) to model the evolution of the three-body system.
The angular momentum of the disc described in the previous section is where the Keplerian angular velocity in the disc is (2) The angular momentum of a circumbinary particle of mass  3 orbiting at semi-major axis  3 is where We solve  d =  3 to find that a particle with the same mass and angular momentum of the initial disc has semi-major axis  3 = 7.8  b .In the particle simulation we incline the particle orbit by  = 130 • and take Ω = 90 • initially.The inclination of the particle relative to the binary is where lp is a unit vector parallel to the particle angular momentum vector and lb is a unit vector parallel to the binary angular momentum vector.The longitude of ascending node is where êb = (  ,   ,   ) is the binary eccentricity vector (e.g.Chen et al. 2020).The phase angle of the eccentricity vector of the binary in the plane of the initial binary orbit is calculated with The left side of Fig. 1 shows the evolution of the binary eccentricity, semi-major axis and   in time for the particle simulation in blue.The upper two panels of the right side of Fig. 1 show the inclination and nodal phase angle of the third body relative to the binary for the three-body simulation in blue.The eccentricity of the binary oscillates up to a maximum of 0.85 while the inclination oscillates to lower values as a result of the ZKL effect.The binary separation remains constant.

Circumbinary disc simulation
The yellow lines in the left panel of Fig. 1 show the binary eccentricity, semi-major axis and   in the circumbinary disc simulation.
Initially the binary evolves in a very similar way to that in the particle simulation.However, as the disc evolves from its initial surface density distribution, the behaviour begins to differ.
The upper two panels of the right panel of Fig. 1 shows the disc evolution at three different radii.The middle of the disc (green line) is initially very similar to the particle.However, the inner parts of the disc (yellow) show that the inner disc undergoes nodal libration and aligns to polar ( = Ω = 90 • ).The later evolution of the outer parts of the disc show a misaligned and circulating disc since the longitude of ascending node varies over 360 • .The binary and the inner polar disc precess together with the combined effect of GR and the massive inner disc.The binary semi-major axis decreases in the disc simulation as a result of the accretion of low angular-momentum material from the disc.The times of high accretion rate (where the disc mass decreases more rapidly in the lower right panel of Fig. 1) correspond to the times where the binary orbit shrinks more rapidly (middle left panel of Fig. 1).Fig. 2 shows the surface density, inclination and longitude of ascending node of the disc as a function of radius at different times.Fig. 3 shows the disc from three different viewing locations at the same times as the lines in Fig. 2. The top row in Fig. 3 shows the initial setup.In the second row, at a time of  = 1000  orb , the gas has flowed inwards towards the binary and is strongly warped close to the binary.At a time of  = 2000  orb , a narrow inner ring has broken off inside a radius of about 2.6  b and aligned to retrograde (green line in Fig. 2).This ring is accreted relatively quickly (see also Nixon et al. 2011a).At at time of  = 2300  orb the disc breaks again, this time farther out at about  = 6  b because the binary by now is much more eccentric.This inner ring aligns to polar and remains there for the remainder of the simulation.Between about  = 2500  orb and  = 3000  orb there are three rings in the disc with the inner two being polar.At later times, the two inner rings join together and there is just one radially wide inner polar ring.
In order for the outer disc to drive ZKL oscillations of the binary, the driving of the apsidal precession of the binary must be dominated by the outer disc, otherwise its effect on the binary eccentricity averages to zero (Holman et al. 1997).The lower left panel of Fig. 1 shows the apsidal precession of the binary in the frame of the initial binary orbit.In the absence of other effects, GR drives prograde apsidal precession of the binary on a timescale  GR = 25, 000  orb for  b = 0.2 and this decreases to  GR = 17, 000  orb for  b = 0.8 (see equation 27 in Naoz et al. 2017).The disc drives prograde apsidal precession for inclinations close to coplanar (e.g.Lepp et al. 2023b) but faster retrograde apsidal precession for polar inclinations (e.g.Zhang & Fabrycky 2019; Childs et al. 2023b).The apsidal precession observed in the simulation is initially slightly prograde.This is a combination of the small prograde GR plus a small prograde component from the inclined circumbinary disc (that is inclined by  ≳ 100 • for all radii).The retrograde inner ring drives prograde apsidal precession.Once the binary becomes highly eccentric and the polar inner ring forms, the apsidal precession is much faster and retrograde.This retrograde precession prevents ZKL oscillations from continuing and the binary eccentricity remains high.
A polar disc ring drives retrograde apsidal precession of the binary at a rate three times faster than the prograde apsidal precession driven by an equivalent prograde (or retrograde) ring (e.g.Zhang & Fabrycky 2019;Childs et al. 2023b).The higher the binary eccentricity, the more likely it is that a polar ring can form since the larger the range of inclinations for which the orbits are librating (e.g.Verrier & Evans 2009;Chen et al. 2019).This suggests that the locking of the binary eccentricity through the formation of an inner polar ring is much more likely for high binary eccentricity.
We have also considered simulations with a larger disc aspect ratio.With / = 0.02, similar behaviour is observed with slightly reduced eccentricity growth.The peak eccentricity is about 0.6 and after the disc breaks, it settles to about 0.4.With / = 0.05, there is no disc breaking, but because of the much shorter viscous timescale, there is not significant eccentricity growth before the disc has largely dissipated.In the case that the disc is in good radial communication and does not undergo disc breaking, then the system behaves more similarly to the case with a third body instead of a disc and the binary eccentricity oscillates.This could occur for a disc with a much smaller viscosity  parameter than used in this work.
The merger of the black hole binary is driven by gravitational radiation on an approximate timescale 1964).A decrease in the semi-major axis of the binary of about 10% (as seen in our simulation) leads to a merger speed-up factor of about 1.5.While the gas disc is present, the ZKL effect leads to increased binary eccentricity that can remain high even after the gas disc has dissipated.The increase in the binary eccentricity up to about 0.7 leads to merger speed-up factor of about 10. Therefore the increased eccentricity provides a larger merger time speed-up factor than the accretion of material that directly shrinks the binary orbit.We note that larger speed-up factors can be achieved for a larger disc mass and wider binaries that drive higher eccentricities (Lepp et al. 2023a).
While we did not observe significant binary eccentricity growth in the simulation with  = 170 • , it is useful to compare the amount of binary shrinkage in this simulation since retrograde circumbinary discs are already known to shrink the orbit of a binary as a result of the accretion of material with negative angular momentum (Nixon et al. 2011a,b).For  = 170 • , the binary semi-major axis decreases at a steady rate of that is less than half the rate observed in our  = 130 • case.There is no disc breaking for  = 170 • and the disc mass decreases about a factor of two more slowly.This suggests that the binary shrinkage observed in the  = 130 • simulation is a result of the accretion of low angular momentum material.

Four-body simulation
In order to understand the mechanism more clearly, we consider now some -body simulations in Fig. 4. The blue lines show the same three-body simulation as described in Section 2.2 and shown in Fig. 1.The green lines show a simulation in which we replace the outer particle with two particles to represent the inner and outer disc close to the time when the disc breaks,  = 2500  orb .In the disc simulation, at  = 2700  orb , the mass of the inner disc is 0.22  d0 and the mass of the outer disc is 0.68  d0 , where  d0 is the initial disc mass.The angular momentum ratio of the inner disc to the outer disc is 0.22.In the four-body simulation, the two particles have the same masses as the two parts of the disc and the same angular momentum ratio.The inner particle is at an orbital radius of 3.6  b and the outer is at 7.7  b .The inner particle begins at  = Ω = 90 • and the outer particle has the same properties as the three-body particle at time  = 2500  orb except the mass and semi-major axis are slightly smaller.Fig. 4 shows the binary evolution.The evolution of the binary eccentricity and   are remarkably similar to the disc simulation.The eccentricity oscillations become much smaller in amplitude with the polar inner particle.The apsidal precession also becomes much faster and in a retrograde direction.Finally, the yellow lines show a three-body simulation but with the outer particle mass set to 0.68  d0 and orbital radius of 7.7  b .This shows that the effect of the decreasing disc mass is not so significant as the formation of the inner polar ring.

EXTENSION OF THE PARAMETER SPACE
We now consider the applicability of the mechanism to varying disc mass, black hole masses and disc parameters.In order to explore a wide range of parameter space that is not feasible with circumbinary disc simulations, we use just the three-body simulations to explore the effect of changing the mass of the third body and varying black hole masses.

Mass of the companion object
Fig. 5 shows the maximum eccentricity of the binary in the three-body simulations in which we vary the mass, orbital radius and inclination of the third body.The largest eccentricities are reached for a third body that is inclined by 170 • (top panel), however, the radial range for the orbit is narrow.A disc is unable to drive such high eccentricity growth beginning at this inclination because it radially spreads on a shorter timescale than the ZKL oscillation can be driven.For lower initial disc inclination, the radial range of orbits that can drive ZKL of the binary increases.Therefore ZKL driving from a disc is easier at  = 130 • .This is also the most effective inclination since for lower inclinations the maximum eccentricity decreases.

Mass of the black hole binary
Fig. 5 shows that there is maximum radius for which the particle can drive binary eccentricity growth.Material inside of this radius can drive the growth.If there is any material outside of this radius, it is precessing too slowly to keep up with the binary and is found in circulating orbits (Childs et al. 2023a).We can estimate analytically the radius at which the material needs to be accreted for this mechanism to work.Lepp et al. (2023a) found the critical radius at which the stationary inclination reaches 180 • to be where the angular frequency of the binary is Ω b = √︃   b / 3 b .This is the radius outside of which there are no librating orbits.The apsidal precession of the binary driven by GR is and the precession rate for the longitude of the periastron of the binary caused by the small mass companion is given by A more general expression for  2 is given in Lepp et al. (2023a) but at an inclination of 180 • it simplifies to (Zanardi et al. 2018;Zhang & Fabrycky 2019;Lepp et al. 2023a).For the initial parameters of the disc simulation, we have  c1 / b = 10.3.However, during the simulation the binary eccentricity evolves, as does the surface density profile of the disc.We now examine how the critical radius scales with mass and length.For a fixed binary eccentricity, we have  GR ∝

√︃
b / 3 b .Since the event horizon of a black hole scales linearly with mass it is reasonable to consider the case where both mass and length share the same overall scaling, so  b ∝  b .In this case all of the terms scale as 1/ b and  c / b is scale free.It should be noted, the same is true for the other critical radii in Lepp et al. (2023a).So long as the masses and lengths are scaled by the same amount, the critical radii relative to  b remains unchanged.This was also seen in Childs et al. (2023a) for the test particle case in which the apsidal precession is driven only by GR.
Thus the results described in this work can be scaled to different mass black hole binaries by scaling the semi-major axis by the same factor.For example, if we considered a binary of two 3 × 10 6 M ⊙ black holes, this is 10 5 times the mass that we have considered.If the binary separation is increased by the same factor of 10 5 , so that  b = 10 6 R ⊙ = 0.02 parsec, then the critical radius relative to the binary separation would be the same.Although our analytical scaling only works at 90 • and 180 • we have run numerical tests that suggest that this scaling works for this system.

Disc lifetime
The considerations in the previous two subsections have not taken into account the lifetime of the disc.For this mechanism to operate requires the disc lifetime to be much longer than the ZKL oscillation timescale.The lifetime of the disc is approximately the viscous timescale that is given by (e.g.Pringle 1981).Note that   / orb is independent of the binary mass.An approximate ZKL timescale can be found with for a circular orbit outer body (Ford et al. 2000;Kiseleva et al. 1998), where  d is the disc mass.This suggests that for a disc mass that is a fixed fraction of the binary mass ( d ∝  b ), the ZKL timescale is independent of the binary mass.Note that we require   ≫  ZKL since as the disc mass decreases, the ZKL effect weakens.
For the parameters of our simulation, with  = 7.8  b , we have   / orb ≈ 3.5 × 10 5 , much longer than the ZKL oscillation timescale that is a few thousand binary orbits (see Fig. 1).For a disc with a larger disc aspect ratio, / = 0.05, the viscous timescale can become comparable to the ZKL oscillation timescale within a factor of a few.In this case, the disc is depleted before the binary eccentricity grows, as discussed in Section 2.3.In general, we expect the mechanism to work for a range of black hole masses provided that the disc mass is scaled accordingly and that the disc aspect ratio is sufficiently small.

CONCLUSIONS
We have found that a highly misaligned disc around a binary black hole can significantly increase the eccentricity of a black hole binary through the ZKL mechanism.Disc breaking that occurs during the high eccentricity phase can lead to an inner polar ring.The polar ring drives fast retrograde apsidal precession of the binary orbit which suppresses the ZKL effect and the binary eccentricity can remain high.Therefore, the capture of retrograde circumbinary material can speed up the merger more than the capture of a third body.This has significant implications for the merger timescale for black hole binaries that accrete material through chaotic accretion processes.

Figure 1 .Figure 2 .
Figure 1.Comparison of the circumbinary disc simulation to a three-body simulation.Left: The binary eccentricity (upper panel), binary semi-major axis scaled to the initial semi-major axis (middle panel) and phase angle of the binary eccentricity vector (lower panel) as a function of time for the -body simulation (blue) and the disc simulation (yellow).Right: The inclination (upper panel) and longitude of ascending node (middle panel) of the third body (blue) and disc (at radius  = 4  b in yellow,  = 7  b in green and  = 10  b in red) in the frame of the binary orbit.The lower panel shows the total mass of the disc.

Figure 3 .
Figure 3.Each row shows three views of the same disc: the disc in the  −  plane in which the binary orbits initially (left), the  −  plane (middle) and the  −  plane (right).The times shown from top to bottom are  = 0, 1000, 2000, 2300, 2700 and 4200  orb .The red circles show the binary.

Figure 4 .
Figure 4. Binary eccentricity and   for the -body simulations.The blue lines here show the three-body simulation that is the same as the blue lines of Fig. 1.This simulation has  1 =  2 = 30 M ⊙ ,  3 = 0.033  b ,  3 = 7.8  b ,  = 130 • and Ω = 90 • initially.The yellow lines show a simulation in which the mass of the third body decreases at time  = 2500  orb to 0.68  3 and its semi-major axis moves to 7.7  b .The green lines show the four-body simulation in which the outer particle is replaced by two particles to represent the two parts of a broken disc, the inner polar disc and the outer misaligned disc.The particles have masses of 0.22  3 and 0.68  3 and semi-major axes of  3 = 3.6  b and  4 = 7.7  b , respectively.
The colours show the maximum eccentricity of the inner binary for a three-body system for varying companion mass,  3 , and orbital radius, , for four different initial inclinations.The inner binary is equal mass with  b = 60 M ⊙ ,  b = 10 R ⊙ and  b = 0.2 initially.