# Dynamics and spin alignment in massive, gravito-turbulent circumbinary discs
around supermassive black hole binaries

Martin A. Bourne1,2, Davide Fiacconi1,2, Debora Sijacki1,2, Joanna M.
Piotrowska2,3 and Sophie Koudmani1,2,4,5
1 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge,
CB3 0HA, UK
2 Kavli Institute for Cosmology, University of Cambridge, Madingley Road,
Cambridge, CB3 0HA, UK
3 Cahill Center for Astronomy and Astrophysics, California Institute of
Technology, Pasadena, CA, USA
4 Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue,
New York, NY 10010, USA
5 St Catharine’s College, University of Cambridge, Trumpington Street,
Cambridge CB2 1RL, UK
E-mail<EMAIL_ADDRESS>
###### Abstract

Parsec-scale separation supermassive black hole binaries in the centre of gas-
rich galaxy merger remnants could be surrounded by massive circumbinary discs
(CBDs). Black hole mass and spin evolution during the gas-rich binary inspiral
are crucial in determining the direction and power of relativistic jets that
radio observations with LOFAR and SKAO will probe, and for predicting
gravitational wave (GW) emission that IPTA and LISA will measure. We present
3D hydrodynamic simulations capturing gas-rich, self-gravitating CBDs around a
$2\times 10^{6}$ M⊙ supermassive black hole binary, that probe different mass
ratios, eccentricities and inclinations. We employ a sub-grid Shakura-Sunyaev
accretion disc to self-consistently model black hole mass and spin evolution
together with super-Lagrangian refinement techniques to resolve gas flows,
streams and mini-discs within the cavity, which play a fundamental role in
torquing and feeding the binary. We find that higher mass ratio and eccentric
binaries result in larger cavities, while retrograde binaries result in
smaller cavities. All of the simulated binaries are expected to shrink with
net gravitational torques being negative. Unlike previous simulations, we do
not find preferential accretion onto the secondary black hole. This implies
smaller chirp masses at coalescence and hence a weaker GW background.
Critically this means that spin-alignment is faster than the binary inspiral
timescale even for low mass ratios. However, we find that mini-disc and hence
spin alignment is not guaranteed in initially misaligned systems, potentially
leading to a significant fraction of recoiled remnants displaced from their
host galaxies if chaotic accretion is the dominant feeding channel.

###### keywords:

black hole physics – accretion, accretion discs – quasars: supermassive black
holes – methods: numerical – black hole mergers – gravitational waves

††pagerange: Dynamics and spin alignment in massive, gravito-turbulent
circumbinary discs around supermassive black hole binaries–Dynamics and spin
alignment in massive, gravito-turbulent circumbinary discs around supermassive
black hole binaries††pubyear: 2023

##  1 Introduction

Within the standard $\Lambda$-cold dark matter ($\Lambda$CDM) cosmological
model of our Universe (White & Rees, 1978; Peebles, 1982; Blumenthal et al.,
1984; Davis et al., 1985), cosmic structure formation proceeds hierarchically
in a bottom-up fashion. Small mass structures form first and then grow in mass
via “smooth” inflows and mergers to form larger mass objects as cosmic time
progresses (Navarro et al., 1995; Mo & White, 2002; Fakhouri & Ma, 2008;
Fakhouri et al., 2010; Guo & White, 2008; Genel et al., 2009; O’Leary et al.,
2021). Mergers of cosmic structures have been studied extensively (Toomre &
Toomre, 1972; Barnes, 1992; Boylan-Kolchin et al., 2005; Cox et al., 2006;
Naab & Trujillo, 2006; Conselice, 2014; Fiacconi et al., 2015; Talbot et al.,
2023), and are a crucial mechanism in galaxy evolution. Not only does this
lead to the mass assembly of dark matter halos and galaxies within but, in
some cases significant amounts of gas can be funnelled to the centre of the
remnant galaxy that can drive both star bursts (Barnes & Hernquist, 1991,
1996; Mihos & Hernquist, 1996; Di Matteo et al., 2005; Hopkins et al., 2006;
Cox et al., 2008; Hopkins et al., 2013), and also potentially black hole
growth and activity (Di Matteo et al., 2005; Hopkins et al., 2006, 2007;
Johansson et al., 2009; Barai et al., 2014; Gao et al., 2020; Pierce et al.,
2023).

While the exact nature of gas delivery to supermassive black holes in the
aftermath of a galaxy merger is not yet fully understood, supermassive black
holes are likely residing in the centre of the majority of massive galaxies
(see e.g., Ferrarese & Merritt, 2000; Häring & Rix, 2004; Gültekin et al.,
2009; Kormendy & Ho, 2013; McConnell & Ma, 2013; Sahu et al., 2019) and these
galaxies are expected to undergo multiple major mergers with similarly sized
galaxies over their cosmic lifetime (see e.g., O’Leary et al., 2021; Huško et
al., 2022a, and references therein). Hence it is natural to expect that during
the merging process, due to the dynamical friction, supermassive black holes
present in parent galaxies will sink towards the centre of the merger remnant
system and form a $\sim$ parsec scale supermassive black hole binary (Begelman
et al., 1980; Mayer et al., 2007; Mayer, 2013; Khan et al., 2016, 2018; Colpi,
2014). Depending on the baryon content of the local environment, the hardening
processes of such a supermassive black hole binary prior to gravitational wave
(GW) emission may proceed through several different channels, such as three-
body encounters with stars (Quinlan, 1996; Milosavljević & Merritt, 2001;
Berczik et al., 2006; Sesana et al., 2006; Rantala et al., 2017) or a third
black hole (Hoffman & Loeb, 2007; Bonetti et al., 2018; Mannerkoski et al.,
2021), as well as interactions with a gaseous circumbinary disc (CBD, e.g.,
Cuadra et al., 2009; Nixon et al., 2011a; Nixon et al., 2013; Roedig et al.,
2012; D’Orazio et al., 2013; Farris et al., 2014; Moody et al., 2019; Muñoz et
al., 2019; Muñoz et al., 2020; D’Orazio & Duffell, 2021; Franchini et al.,
2021, 2022). If both supermassive black holes are accreting then they should
appear as dual AGN, with a number of such objects observed on $\gtrsim$ kpc
scales (see Fig. 1 in Chen et al., 2022). On scales of $\lesssim 1-10$ pc only
a handful of binary candidates have been observed, for example using VLBA
(Rodriguez et al., 2006; Bansal et al., 2017; Kharb et al., 2017), while on
even smaller sub-parsec scales candidate binaries have been inferred, for
example, from light curves (Sillanpaa et al., 1988; Graham et al., 2015; Chen
et al., 2020; Millon et al., 2022). Notwithstanding the difficulties of these
observational efforts and possible systematic biases, the low number of
observed binaries tantalisingly suggests that the hardening process may be
efficient with the majority of binaries merging into a larger supermassive
black hole within the age of the Universe.

The Advanced LiGO and VIRGO collaboration opened the GW window onto our
Universe when they made the first detection of merging stellar mass black
holes (Abbott et al., 2016), and has been followed by the observation of an
intermediate-mass black hole (Abbott et al., 2020). More recently, the first
tentative detections of a nHz GW background, expected to be produced by a
cosmic population of merging supermassive black holes, were made by members of
the International Pulsar Timing Array (IPTA Verbiest et al., 2016; The
International Pulsar Timing Array Collaboration et al., 2023); NANOGrav
(Agazie et al., 2023a), the Parkes PTA (Reardon et al., 2023) and joint
results from the European and Indian PTAs (EPTA Collaboration et al., 2023).
An upper limit on the strain amplitude of individual supermassive black hole
binaries was also presented (Agazie et al., 2023b). Future space-based
missions, such as LISA, will provide a fantastic opportunity to observe $\sim
0.01-100$ mHz frequency GWs emitted by individual supermassive black hole
merger events over a wide range of masses (peaking between $\sim
10^{4}-10^{7}$ M⊙) and out to very high redshifts (Amaro-Seoane et al., 2017,
2023). It is hence of great importance to provide detailed theoretical
predictions on the merger rate of supermassive black holes, binary parameters,
and also on the likely black hole spin magnitudes and orientation before
merger. The concerted effort of theoretical predictions and future
observations will be a unique way not only to probe spin evolution and hence
black hole radiative efficiencies but also the amount of gravitational recoil
merger remnants will experience (Campanelli et al., 2007; González et al.,
2007; Lousto & Zlochower, 2011, 2013; Lousto & Healy, 2019; Sperhake et al.,
2020). This will shed light on the occupation fraction of supermassive black
holes in galaxies, which impacts predicted GW detection rates (e.g. Sesana,
2007; Sesana et al., 2009a; Blecha & Loeb, 2008), the likelihood of offset
holes and the amount and spatial delivery of feedback from AGN. As well as
determining individual GW events, merger rates and binary properties
additionally directly impact the expected form of the stochastic GW background
(e.g., Rajagopal & Romani, 1995; Phinney, 2001; Enoki et al., 2004; Sesana et
al., 2004; Sesana et al., 2009b; Kelley et al., 2017; Siwek et al., 2020).

A very interesting regime, especially for multi-messenger science, comprises
of merging supermassive black hole binaries in gas-rich environments. Here it
is expected that both supermassive black hole feeding and feedback processes
could be of relevance, which may lead to significant spin evolution as well.
Simulations of merging gas-rich galaxies suggest that a non-negligible
fraction of gas loses angular momentum and is hence efficiently funnelled
towards the central region of the merger remnant (Barnes & Hernquist, 1991,
1996; Mihos & Hernquist, 1996; Di Matteo et al., 2005; Hopkins et al., 2006;
Mayer et al., 2007, 2010; Cox et al., 2008; Johansson et al., 2009; Hopkins et
al., 2013). Such a gas-rich environment provides a prime opportunity for the
formation of massive CBDs around the supermassive black hole binary, for
example via the circularisation of infalling gas clouds (Dunhill et al., 2014;
Goicovic et al., 2016, 2017, 2018).

In prograde binaries, gravitational torques from the CBD have traditionally
been invoked to facilitate the extraction of angular momentum to drive the
shrinking of the binary before it enters the GW regime (see e.g., Pringle,
1991; Artymowicz & Lubow, 1994, 1996; Ivanov et al., 1999; Gould & Rix, 2000;
Armitage & Natarajan, 2002, 2005; Lodato et al., 2009). However, gas that
resides inside the cavity of a CBD, such as streams or mini-discs, can also
affect the binary orbital evolution by supplying material for accretion and
acting as an additional source of gravitational torques (see e.g., MacFadyen &
Milosavljević, 2008; Cuadra et al., 2009; Shi et al., 2012; Roedig et al.,
2011, 2012; D’Orazio et al., 2013; Farris et al., 2014; Miranda et al., 2017;
Muñoz et al., 2019; Tiede et al., 2020; Franchini et al., 2022; Siwek et al.,
2023b). In particular, simulations that resolve the mini-discs and gas flows
within the binary orbit ($R<a$) find that net torque from this region is
generally positive while the net torque from larger radii is generally
negative (Roedig et al., 2012; Muñoz et al., 2019; Tang et al., 2017; Tiede et
al., 2020; Franchini et al., 2022; Siwek et al., 2023b). With this in mind,
several recent studies have found that the total net torque acting on a binary
can, in some scenarios, be positive and lead to binary expansion (e.g.,
Miranda et al., 2017; Tang et al., 2017; Moody et al., 2019; Muñoz et al.,
2019; Muñoz et al., 2020; Duffell et al., 2020; Heath & Nixon, 2020; Tiede et
al., 2020; D’Orazio & Duffell, 2021; Siwek et al., 2023b; Wang et al., 2023).
These simulations typically invoke a (locally-)isothermal equation of state
and neglect self-gravity. There is only a limited number of simulations that
instead consider self-gravitating, massive CBDs (Cuadra et al., 2009; Roedig
et al., 2011, 2012; Roedig & Sesana, 2014; Franchini et al., 2021). These
simulations, which additionally employ an adiabatic equation of state with
$\beta-$cooling (see Section 2) to regulate the disc stability, universally
find that binaries shrink. A further complication that has only been addressed
by a handful of works is the effect of misaligned (Nixon et al., 2013; Dunhill
et al., 2014; Aly et al., 2015; Moody et al., 2019) and retrograde (Nixon et
al., 2011a; Dunhill et al., 2014; Roedig & Sesana, 2014; Heath & Nixon, 2020;
Tiede & D’Orazio, 2023) CBDs, with the latter behaving quite differently to
the prograde case due to lack of resonant torquing by the CBD and providing a
promising channel to shrink binaries.

The main aim of this paper is to build on existing literature and model for
the first time both the mass and spin evolution of black holes in massive
gravito-turbulence CBDs, as well as analyse the general CBD and binary
evolution. To achieve this we perform a suite of high-resolution simulations
of gas-rich CBDs around live binaries whose properties are expected to place
them in the observing range of LISA (Amaro-Seoane et al., 2017, 2023). We
consider different configurations by varying the binary mass ratio,
eccentricity and inclination angle (including fully retrograde binaries). In
order to resolve streams and mini-discs within the CBD cavity we employ a
super-Lagrangian refinement technique (Curtis & Sijacki, 2015) that vastly
improves the resolution around each black hole (see also, Franchini et al.,
2022). Additionally, given the important implications of black hole spin on
both the GW emission and black hole recoil velocities, we employ the sub-grid
accretion disc model of Fiacconi et al. (2018) to track both the mass and spin
evolution of the black holes. The paper is organised as follows; in Section 2
we present the numerical code adopted, physical modules included and our setup
of the initial conditions for the simulations performed. In Section 3 we
outline the main results, including the evolution of the CBD and binary, an
analysis of the expected torques experienced by the binary and the black hole
growth and spin evolution. In Section 4 we put these results into the context
of the existing literature and discuss their implications while in Section 5
we summarise our main conclusions.

##  2 Numerical simulations

###  2.1 Initial conditions

We consider supermassive black hole binaries with initial total mass $M_{\rm
bin}=M_{\bullet,1}+M_{\bullet,2}=2\times 10^{6}$ M☉ and mass ratio $0.1\leq
q\equiv M_{\bullet,2}/M_{\bullet,1}\leq 1$. We chose these parameters to
bracket the binary parameter space that would be of the largest interest for
the Laser Interferometer Space Antenna (LISA; Amaro-Seoane et al. 2013). The
binary is initially on a Keplerian orbit with eccentricity of either $e=0$ or
$e=0.5$ and initial semi-major axis $a=2$ pc, corresponding to a period of
$t_{\rm bin}=0.187$ Myr. The peri- and apo-centre of the binary are aligned
along the $x$ axis.

A gaseous CBD surrounds the binary in the $xy$ plane and initially extends
between $R_{\rm in}=2a=4$ pc and $R_{\rm out}=7a=14$ pc from its centre of
mass. The binary and the disc are not necessarily coplanar. For misaligned
systems, the orbital plane of the binary crosses the disc midplane only
through the $y$ axis with an inclination angle $i$; for reference, $i=0\degr$
and $i=180\degr$ correspond to prograde and retrograde coplanar orbits,
respectively. The disc follows the surface density profile

$\Sigma(R)=\Sigma\left(\frac{R}{R_{\rm in}}\right)^{-\alpha},$ (1)

where $\alpha=2$. The normalisation $\Sigma$ is chosen to enforce $M_{\rm
cbd}/M_{\rm bin}=0.1$, where $M_{\rm cbd}$ is the total mass of the CBD.

The vertical structure of the CBD is constructed by assuming that the aspect
ratio $H/R\approx c_{\rm s}(R)/V_{\rm K}(R)$ is constant, where $c_{\rm s}(R)$
and $V_{\rm K}(R)$ are the local isothermal sound speed and the Keplerian
velocity for $M_{\rm bin}$, respectively. Therefore, the temperature profile
of the CBD, $c_{\rm s}^{2}(R)=(H/R)^{2}V_{\rm K}^{2}(R)\propto R^{-1}$, can be
normalised after choosing the value of the aspect ratio. Since we include
self-gravity in the CBD, we use the Toomre parameter $Q$ to impose the
normalisation. Specifically, we choose $Q\approx 3(H/R)(M_{\rm bin}/M_{\rm
cbd})=1.5$, which corresponds to an initial aspect ratio $H/R\approx 0.5M_{\rm
cbd}/M_{\rm bin}$. We do not include any additional explicit source of
viscosity beyond the intrinsic (low) numerical viscosity, instead rellying on
the angular momentum transport due to self-gravity (e.g., Cuadra et al.,
2009). The vertical structure of the CBD is locally sampled from a Gaussian
distribution with a standard deviation given by $H(R)$. The CBD is initially
in nearly Keplerian rotation with azimuthal velocity

$V_{\phi}^{2}(R)=\frac{GM_{\rm
bin}}{R}\left[1+\frac{3}{4}\frac{q}{(1+q)^{2}}\left(\frac{R}{a}\right)^{-2}\right]-3c_{\rm
s}^{2}(R),$ (2)

where the term in the square brackets represents the correction for the tidal
field of the binary and the last term accounts for the gas pressure. We do not
add any additional velocity dispersion to the gas.

Table 1 summarises the properties of the supermassive black hole binaries
explored in this work. We explore 1:1, 1:3 and 1:10 mass ratio binaries
denoted by q01, q03, and q10, respectively, and consider circular ($e=0$),
eccentric ($e=0.5$), co-planar prograde ($i=0\degr$), misaligned
($i=45\degr$), and retrograde ($i=180\degr$) binaries. For example, the 1:1
mass ratio, circular, co-planar prograde binary simulation is labelled
q01e00i00. The binary and CBD, as expected, torque each other from the outset
of the simulation, i.e. also during the relaxation phases (see Sections 2.2
and 2.3). To consider the $45\degr$ inclined binary evolution in a pre-relaxed
system, we have an additional simulation, q03e00i45mod, in which we take the
q03e00i00 run after the two initial relaxation phases, rotate the binary and
cavity region about the $y-$axis in order to be inclined with respect to the
CBD by $45\degr$ and track its subsequent evolution.

Table 1: Summary of the characteristics of the performed runs. Note that: (i) $M_{\rm bin}=2\times 10^{6}$ M☉; (ii) $q\equiv M_{\bullet,2}/M_{\bullet,1}\leq 1$. From left to right: simulation label, initial mass ratio $q$, initial eccentricity $e$, initial inclination angle $i$ between the binary orbital plane and the CBD midplane, initial Eddington ratio $f_{\rm Edd,1}$ for $M_{\bullet,1}$, initial Eddington ratio $f_{\rm Edd,2}$ for $M_{\bullet,2}$, initial $\alpha$ accretion disc and black hole angular momentum ratio for $M_{\bullet,1}$, initial $\alpha$ accretion disc and black hole angular momentum ratio for $M_{\bullet,2}$, misalignment angle between $\mn@boldsymbol{J}_{\bullet}$ and $\mn@boldsymbol{J}_{\rm d}$ for $M_{\bullet,1}$, misalignment angle between $\mn@boldsymbol{J}_{\bullet}$ and $\mn@boldsymbol{J}_{\rm d}$ for $M_{\bullet,2}$. Label | $q$ | $e$ | $i$ | $f_{\rm Edd,1}$ | $f_{\rm Edd,2}$ | $(J_{\rm d}/J_{\bullet})_{1}$ | $(J_{\rm d}/J_{\bullet})_{2}$ | $\theta_{1}$ | $\theta_{2}$
---|---|---|---|---|---|---|---|---|---
|  |  | (°) |  |  |  |  | (°) | (°)
q01e00i00 | $1$ | 0 | 0 | $3\times 10^{-3}$ | $3\times 10^{-3}$ | 0.71 | 0.71 | 76.6 | 64.7
q03e00i00 | $1/3$ | 0 | 0 | $2\times 10^{-3}$ | $6\times 10^{-3}$ | 0.65 | 0.81 | 64.7 | 76.6
q10e00i00 | $1/10$ | 0 | 0 | $1.5\times 10^{-3}$ | $1.5\times 10^{-2}$ | 0.65 | 1.02 | 76.6 | 64.7
q01e05i00 | $1$ | 0.5 | 0 | $3\times 10^{-3}$ | $3\times 10^{-3}$ | 0.71 | 0.71 | 64.7 | 76.6
q03e05i00 | $1/3$ | 0.5 | 0 | $2\times 10^{-3}$ | $6\times 10^{-3}$ | 0.65 | 0.81 | 64.7 | 76.6
q10e05i00 | $1/10$ | 0.5 | 0 | $4\times 10^{-3}$ | $1\times 10^{-2}$ | 0.49 | 1.15 | 64.7 | 76.6
q01e00i45 | $1$ | 0 | 45 | $3\times 10^{-3}$ | $3\times 10^{-3}$ | 0.71 | 0.48 | 86.1 | 96.3
q03e00i45 | $1/3$ | 0 | 45 | $2\times 10^{-3}$ | $6\times 10^{-3}$ | 0.65 | 0.55 | 55.0 | 104.4
q10e00i45 | $1/10$ | 0 | 45 | $1.5\times 10^{-3}$ | $1.5\times 10^{-2}$ | 0.44 | 1.02 | 113.3 | 14.6
q03e00i45mod | $1/3$ | 0 | 45 | $2\times 10^{-3}$ | $6\times 10^{-3}$ | 0.44 | 0.81 | 110.0 | 57.2
q03e00i180 | $1/3$ | 0 | 180 | $1\times 10^{-3}$ | $3\times 10^{-3}$ | 0.79 | 0.99 | 64.7 | 76.7
q10e00i180 | $1/10$ | 0 | 180 | $1\times 10^{-3}$ | $1\times 10^{-2}$ | 0.72 | 1.15 | 64.7 | 76.6

###  2.2 Numerical code and simulation setup

The simulations presented in Table 1 are performed using the moving-mesh
hydrodynamical code arepo (Springel, 2010; Pakmor et al., 2016), which solves
the Euler equations on an unstructured Voronoi tessellation that moves with
the flow. The fluid dynamics are solved using a finite-volume, second-order
reconstruction scheme coupled with an exact Riemann solver. The unstructured
mesh adapts by locally refining and de-refining in a quasi-Lagrangian way to
maintain a typical gas mass $m_{\rm g}^{\rm target}$ per cell. The initial
conditions are sampled with $10^{6}$ initial mesh elements in the CBD
corresponding to a target mass $m_{\rm g}^{\rm target}=0.2$ M☉. We employ the
hierarchical time integration scheme described in Springel et al. (2021),
which splits the Hamiltonian into “fast” and “slow” components, i.e. those on
short versus long time-steps, respectively. With this approach, interactions
between particles/cells on the smallest timesteps (up to $N_{\rm direct}^{\rm
max}=2048$) are solved using direct summation, while a standard octree
approach is used for all other interactions. This is primarily done for
efficiency, to avoid unnecessary tree constructions for small numbers of
particles, but has the added benefit of more accurately calculating the
gravitational forces between the black holes and their nearest neighbours,
which are on the shortest timesteps, via direct summation. We discuss the
implications of using this approach, compared to using solely direct summation
or a tree to solve for the gravity in Section 3.2.2. The black holes have
gravitational softening lengths $\epsilon_{\bullet}=0.01$ pc, while gas cell
softenings adapt with the cell size down to a minimum of $\epsilon_{\rm
g}=0.005$ pc.

We evolve the system assuming an equation of state for an ideal fluid,
$P=(\gamma-1)\rho u$, where $\rho$ is the gas density, $u$ is the specific
internal energy, and $\gamma=5/3$ is the adiabatic index, appropriate for a
monoatomic ideal gas. We also allow the gas to cool radiatively through a
simple $\beta$-cooling prescription, with the rate of cooling given by

$\left(\frac{{\rm d}u}{{\rm d}t}\right)_{\rm cool}=-\frac{u}{t_{\rm
cool}(R)}\,,$ (3)

where the local cooling timescale $t_{\rm cool}(R)$ is proportional to the
local orbital time $\Omega_{K}^{-1}(R)$, namely $t_{\rm cool}(R)=\beta(GM_{\rm
bin}R^{-3})^{-1/2}$. We set $\beta=10$; this choice allows the gas to cool
slowly without fragmenting (Gammie, 2001; Rice et al., 2005; Cuadra et al.,
2009; Franchini et al., 2021), and eventually makes the CBD settle into a
marginally stable configuration with $Q\sim 1$. This choice implicitly assumes
a source of heating that prevents rapid cooling and fragmentation of the CBD,
which could, for example, come from stellar heating (Thompson et al., 2005;
Lodato et al., 2009).

The initial conditions are evolved for 50 binary orbits corresponding to about
$1.5t_{\rm cool}$ at the outer edge of the CBD (note this is the first of two
relaxation periods, see also Section 2.3). During this relaxation period, the
binary perturbs the inner edge of the disc, inducing shocks and spiral waves
that propagate from the cavity through the CBD. The shocks increase the
temperature in the disc, in particular close to the inner rim of the cavity,
where the aspect ratio of the CBD, $H/R$, can reach $\approx 0.1$ and the
Toomre parameter grows to a few before cooling can efficiently radiate the
excess energy. After a few binary orbits, streams of low-density gas from the
cavity edge start to develop and bring fresh material close to the two black
holes.

###  2.3 Black hole refinement

Figure 1: Effect of the super-Lagrangian refinement on the evolution of the
cavity. Panel (a) and (b) compare the surface density of the gas (top row) and
the mesh structure (bottom row) after 100 binary orbits without and with
super-Lagrangian refinement switched on during the last 50 binary orbits,
respectively. Panel (c) and (d) show the evolution of the binary semi-major
axis (normalised to the initial value) and eccentricity, respectively, for the
cases without (orange solid curve) and with (light blue dashed curve) this
extra refinement, respectively. The light blue dot marks when the super-
Lagrangian refinement is activated.

Inflow from the CBD is set by the interplay between the binary gravitational
torque on the inner edge of the cavity and the mechanism of angular momentum
transport in the CBD, i.e. gravito-turbulence in our configuration. In
addition the binary dynamics are determined by a combination of physical
processes including gravitational torques from the CBD, streams and mini-discs
within the cavity and accretion. As such, proper modelling of mass transport
through the cavity is critically important and relies on the capability of
resolving relatively low-density features, such as streams funnelling gas
towards the binary. This is intrinsically hard to achieve in
(nearly-)Lagrangian approaches that would naturally maintain low resolution in
under-dense regions. We can overcome this limitation in arepo by using the
super-Lagrangian refinement scheme presented by Curtis & Sijacki (2015, 2016)
(see also the recent simulations of Franchini et al., 2022, for a similar
approach). Super-Lagrangian refinement allows us to focus both mass and
spatial resolution locally around each black hole particle in the simulation.
Specifically, the scheme splits and merges cells within $R_{\rm ref}$ around
each black hole to force the maximum cell size of the mesh elements to
decrease linearly from $R_{\rm max}^{\rm cell}$ at $R_{\rm ref}$ to $R_{\rm
min}^{\rm cell}$ at the centre. Cell radii are allowed to remain within a
factor 2 of the local maximum. We set $R_{\rm max}^{\rm cell}=0.2$ pc and
$R_{\rm min}^{\rm cell}=0.01$ pc; $R_{\rm max}^{\rm cell}$ is chosen to match
the typical size of the mesh cell just inside the CBD cavity, while $R_{\rm
min}^{\rm cell}$ is equal to the gravitational softening of the black holes.
We have chosen an $e$-dependent value for $R_{\rm ref}$, namely $R_{\rm
ref}=3$ pc and $R_{\rm ref}=5$ pc for $e=0$ and $e=0.5$, respectively;
moreover, we set $R_{\rm ref}=2$ pc for runs with $i=180\degr$ because the
cavity edge forms at closer distances from the binary (see Section 3.1). These
choices guarantee that the refinement region slightly crosses the edge of the
cavity into the CBD, enforcing continuity in the mesh structure. The
refinement may lead to the creation of very low mass cells; such occurrences
are prevented by the lower limit $m_{\rm min}=10^{-5}m_{\rm g}^{\rm target}$
for cells within the refinement region.

Fig. 1 shows the effect of the super-Lagrangian refinement on the evolution of
the cavity for run q01e00i00 as a representative example. We compare the
evolution after $100$ binary orbits without any refinement and with the
refinement activated during the last $50$ orbits to allow for the system to
relax. Gaseous streams connecting the cavity edge to the black holes that can
be barely seen without refinement are significantly better resolved thanks to
the super-Lagrangian refinement. The streams effectively bring material close
to the black holes and lead to the formation of gaseous “mini-discs”,
otherwise absent without refinement. The development of these features also
has an impact on the dynamical evolution of the binary. As already shown by
previous work (Roedig et al., 2011, 2012; Tang et al., 2017; Muñoz et al.,
2019; Tiede et al., 2020; Siwek et al., 2023b; Franchini et al., 2022) and
further discussed in Section 3.2.2, the gravitational torques caused by the
gas streams inside the cavity fundamentally contribute to the binary
evolution. As the super-Lagrangian refinement allows for the streams
penetrating the cavity to be better resolved, the gravitational torques are
more prominent. While in the run without refinement the shrinking of the
binary semi-major axis appears to slow after $\sim 50~{}t_{\rm bin}$, the
binary in the run with refinement continues to shrink after this time, though
we do not see appreciable differences in the evolution of the eccentricity.
After the $\beta$-cooling relaxation, all simulated systems are furthermore
relaxed with the super-Lagrangian refinement for $50$ additional binary
orbits. Thus, the overall relaxation phase lasts $\approx 18.7$ Myr,
corresponding to $\approx 3t_{\rm cool}(R_{\rm out})$, after which the CBD
attains a nearly steady state and we use the resulting snapshot as initial
conditions to be further evolved for $500~{}t_{\rm bin}$ including the
accretion and spin evolution model originally presented by Fiacconi et al.
(2018) and discussed in the next section.

###  2.4 Black hole spin model

The accretion disc model of Fiacconi et al. (2018) assumes that each black
hole is surrounded by a sub-grid accretion disc following the analytical thin
$\alpha$-disc solution111In the following, we will refer to the large-scale,
hydrodynamically-resolved gaseous disc surrounding the binary as the “CBD”, to
the compact, hydrodynamically-resolved gaseous discs around each black hole as
the “mini-discs”, and to the sub-grid, analytic $\alpha$-discs as “accretion
discs”. (Shakura & Sunyaev, 1973). The black holes and accretion discs are
described by their masses and angular momenta, self-consistently coupled
through mass accretion and the Bardeen-Petterson effect (Bardeen & Petterson,
1975; Scheuer & Feiler, 1996; King et al., 2005; Martin et al., 2007; Perego
et al., 2009; Dotti et al., 2013). Moreover, the system composed of the black
hole and the accretion disc is subjected to an inflow of matter and angular
momentum as dictated by the hydrodynamical simulations. We refer the reader to
Fiacconi et al. (2018) for a thorough description of the model.

The accretion model requires the initial black hole mass and spin, the initial
mass and accretion rate of the accretion disc, and the directions of the
angular momenta to be specified (see Table 1). The spin of a black hole
depends on its merger and accretion history, with the properties of the
accretion flow (i.e. how chaotic or coherent it is) determining whether or not
a black hole is effectively spun up (Volonteri et al., 2005, 2007; King &
Pringle, 2006, 2007; Dotti et al., 2013; Sesana et al., 2014). We assume that
both black holes are initially highly spinning as a result of previous
evolution and we assign black hole spins of $a_{\bullet,1,2}=0.9$, consistent
with current observational constraints for $\sim 10^{6}$ M⊙ black holes
(Reynolds, 2019). The timescale, $\tau_{\rm align}$, on which black hole spin
aligns with respect to the accretion disc angular momentum scales as
$a_{\bullet}^{5/7}$ (see Equation (12)) and so in choosing high spin values we
are able to apply something of an upper limit on values of $\tau_{\rm align}$.
The accretion disc is initiated with a mass of $M_{\rm
d}=10^{-3}M_{\bullet,1,2}$, and angular momentum that is aligned with the
angular momentum of the surrounding mini-disc. On the other hand, the black
hole spins have been given arbitrary initial orientations with respect to the
accretion disc angular momentum (shown in Table 1, noting that some of the
accretion disc angular momenta are counter-rotating with respect to the black
hole spin in some of the misaligned set-ups). Given that the initial spins of
individual black holes will be determined by their accretion and merger
histories prior to forming the binary, they do not necessarily need to be
aligned, neither with each other nor their initial mini-discs, at the outset
of our simulations. The initial accretion rate through the sub-grid accretion
disc is set to be about the same for both black holes and close to the initial
mass inflow rate measured inside the mini-discs. Finally, the initial
accretion disc - black hole angular momentum ratio, $J_{\rm d}/J_{\bullet}$,
is set to ensure consistency between the initial $f_{\rm Edd}$ and other black
hole and accretion disc properties (see Equation (8)).

Finally, we note an important modification we made to the original accretion
model. In Fiacconi et al. (2018), the mass inflow rate $\dot{M}_{\rm inflow}$
and the specific angular momentum of the inflowing material
$\mn@boldsymbol{L}_{\rm inflow}$ were calculated with a kernel-weighted
interpolation from a volume encompassing a mass of $32m_{\rm g}^{\rm target}$,
typically slightly larger than the size of the mini-discs. While we employ the
same approach to estimate $\dot{M}_{\rm inflow}$, we have adjusted how we
measure $\mn@boldsymbol{L}_{\rm inflow}$ as despite the central weighting this
approach tends to overestimate its magnitude. Instead, we only take the
direction of $\mn@boldsymbol{L}_{\rm inflow}$ and we assign $L_{\rm
inflow}=L_{\rm d}$, where $L_{\rm d}$ is the instantaneous specific angular
momentum of the accretion disc. This means that the condition originally
imposed in Fiacconi et al. (2018) requiring $L_{\rm inflow}\leq L_{\rm d}$ is
automatically satisfied and we assume that all inflowing gas would be able to
circularise onto the sub-grid accretion disc. In practice we typically find
reasonable agreement between $L_{\rm d}$ and the specific angular momentum of
gas on scales of the sub-grid accretion disc.

##  3 Results

Figure 2: Gas surface density maps projected on the $xy$ plane along the $z$
axis of all the simulations in Table 1 at time $t=100~{}t_{\rm bin}$. For each
pair of images, the one on left-hand side shows the CBD, while the one on the
right-hand side shows a zoom in the inner cavity region. The blue and red
filled circles mark the position of the primary and secondary black hole,
respectively. The red curves for the first three simulations show the Roche
lobes projected on the $xy$ plane. The red arrows in the last two simulations
highlight the location of the dynamical friction wake. Both the appearance and
properties of the CBD and the mini-discs are affected by the binary mass ratio
and orbital parameters. Figure 3: Time-evolution of radial profiles of key CBD
properties for our simulation suite. For each row, the top panels show the gas
surface density, with dashed curves highlighting the initial $\propto R^{-2}$
profile; the middle panels show the aspect ratio $H/R$ estimated as $c_{\rm
s}/V_{\phi}$; and the bottom panels show the Toomre parameter $Q$. Light-grey
and dark-grey shaded areas highlight the marginal stability and instability
regions for $Q$, respectively, in the limit of a razor-thin disc. The vertical
dotted lines in all panels indicate the original inner edge of the CBD at
$2a=4$ pc, while for comparison the vertical grey lines indicate the location
of the CBD density peak, illustrating that for prograde binaries the CBD moves
outwards (especially in the case of eccentric binaries), while for retrograde
binaries it moves inwards.

###  3.1 Properties of the circumbinary discs

After the relaxation phases, we run the system for an additional $500~{}t_{\rm
bin}$ during which the CBD outside the cavity evolves in a nearly steady
state. Fig. 2 compares density projections for all runs at time
$t=100~{}t_{\rm bin}$. Each run is presented by two panels, one showing the
whole disc (orange/purple colourmap) and the other zoomed in to the cavity
region (blue/green colourmap). As expected, all simulations show rather
similar CBDs at large distances from the cavity edge, except for the cases of
misaligned binaries with $i=45\degr$, where the CBD becomes warped owing to
the gravitational torque of the binary (see further discussion in Section
3.2.3). Since we have aligned the image in the binary orbital plane, the CBD
appears distorted because of projection effects (i.e. it does not lie in the
binary orbital plane), while otherwise, it resembles the companion runs with
$i=0\degr$. Prominent spiral arms across the CBD piling up close to the cavity
edge represent a common feature across different simulations (including
previous works that model self-gravitating CBDs, e.g., Cuadra et al., 2009;
Roedig et al., 2012; Franchini et al., 2021). The spiral arms are maintained
by the interplay of cooling and local shock heating induced by self-gravity in
the CBD, preventing fragmentation into self-bound clumps and maintaining a
gravito-turbulent state in the disc. Gravito-turbulence is the main process
responsible for mass and angular momentum transport in the CBD. According to
the analysis of Gammie (2001), gravito-turbulence may be described through the
effective $\alpha$-parameter $\alpha_{\rm
SG}=(4/9)\left[\gamma(\gamma-1)\beta\right]^{-1}\approx 0.04$, where the last
approximate equality holds for our parameter choices. For typical values of
gas surface density and sound speed as measured from the simulations close to
the cavity edge, we expect mass transport rates $\dot{M}=3\pi\nu\Sigma\approx
5\times 10^{-3}$ M☉ yr-1, which, when considering the fact that the CBD is
torqued by the central binary, is indeed qualitatively in agreement with, if a
little higher than, time-averaged inflow rates $\gtrsim 10^{-4}$ M☉ yr-1
measured directly from the simulations, with fluctuations of factor $\sim$a
few between different runs. This suggests that at least qualitatively and on
average over time, the CBD behaves like an $\alpha$-disc with
$\alpha=\alpha_{\rm SG}$.

It is interesting to note that the strength and spatial distribution of spiral
arms and the overall shape of the CBD is affected by the binary mass ratio and
our initial choice of orbital parameters. However, most of the morphological
diversity between various runs appears close to and inside the cavity. The
cavity edge is not circular as its shape is continuously stirred by the
gravitational perturbation from the central binary. This in turn triggers the
formation of streams of gas penetrating the inner region of the cavity and
providing fresh gas to the mini-discs (see also, Roedig et al., 2012; Shi et
al., 2012; D’Orazio et al., 2013; Farris et al., 2014; Ragusa et al., 2016;
Miranda et al., 2017; Muñoz et al., 2019; Tiede et al., 2020; Franchini et
al., 2022, noting, however, that not all of these simulations capture the
mini-discs themselves.).

Except for the $q=1$ case, where the whole system is symmetric and both black
holes alternatively induce the formation of the streams depending on surface
density fluctuations close to the cavity edge, the secondary black hole is
typically associated with the most prominent stream. The size of the cavity
decreases for lower values of the mass ratio $q$, with a similar trend being
seen in figure 3 of Siwek et al. (2023a). This indicates that the torques
exerted by the binaries that move the inner edge of the CBD outwards become
less effective at doing so for low q values. Indeed, we find that the cavity
becomes larger for all prograde binaries, except for the q10e00i00 case, which
stays roughly constant albeit with a less well-defined edge.

On the other hand, the retrograde ($i=180\degr$) cases all exhibit shrinking
cavities as the CBD inner edge moves inwards. This is expected due to
differences between the binary-disc interaction in the prograde vs. retrograde
cases. In prograde discs, gravitational torques act to reduce the binary
separation while transporting angular momentum and the inner edge of the CBD
outwards (e.g., Pringle, 1991; Cuadra et al., 2009; Lodato et al., 2009). On
the other hand, due to a lack of resonances in retrograde CBDs the net
gravitational torques are significantly weaker, which allows continued inflow
of material and a smaller cavity (Nixon et al., 2011a; Roedig & Sesana, 2014;
Bankert et al., 2015).

Typically, the sharp inner edge at $2a$ is lost rapidly during relaxation. For
the eccentric cases in particular the inner edges of the CBDs are disturbed
and move outwards. In the high eccentricity runs, the cavity edge is less
sharp than in the cases with nearly circular orbits, with cavities also being
more extended and containing more diffuse material. The diffuse material is
produced by the hydrodynamical interaction between the streams and the central
binary. At the apocentre, the cavity edge is maximally perturbed and the
streams start to develop; by the time the black holes reach the pericentre,
their relative velocity with respect to the streams increases and the faster
motion of the black holes makes streams wind up. The streams are eventually
mixed and dispersed just after the pericentric passage when the separation
between the two black holes increases again.

The binary parameters also have a large impact on the sizes of the mini-discs.
As already pointed out by Artymowicz & Lubow (1994), the radial extent of the
mini-discs changes with the mass ratio of the binary and the eccentricity of
the orbit. For circular orbits, the mini-discs surrounding the primary and the
secondary black holes grow and decrease, respectively, as the binary mass
ratio decreases, approximately following the shape of the Roche lobes as shown
in Fig. 2. When highly eccentric orbits are considered, the sizes of the mini-
discs are reduced by the tidal truncation induced at the pericentre. In all
cases, the mini-discs show spiral arm patterns induced by the interaction with
and the circularisation of the infalling streams, as well as bridges of
material between the two mini-discs typically following the primary.

Barring the mini-disc around the primary black hole in the q10e00i180 case,
mini-discs appear significantly less developed in retrograde binaries compared
to the prograde cases. In the $q=1/3$ retrograde run, a tiny circum-primary
mini-disc forms, while the secondary black hole does not develop a clear mini-
disc and accretion mostly proceeds through episodic and very thin streams of
gas plunging toward the black hole. Similarly, no mini-disc forms around the
secondary black hole in the retrograde $q=1/10$ run either. The retrograde
binary orbits in a cavity that is less depleted of material when compared with
the corresponding prograde cases. No obvious gas streams seem to penetrate the
cavity. The binaries develop tiny perturbations in the background material
that resemble the wakes caused by dynamical friction, in particular for the
$q=1/10$ simulations, where more material diffuses in from the CBD inner edge,
and the whole cavity resembles a Type II gap between the mini-disc around the
primary and the CBD, as observed in protoplanetary disc calculations (Rafikov,
2002; Nelson & Papaloizou, 2003; Crida & Morbidelli, 2007) and has previously
been suggested as a mechanism for bringing black hole binaries towards merger
(e.g., Goldreich & Tremaine, 1980; Lin & Papaloizou, 1986; Gould & Rix, 2000;
Armitage & Natarajan, 2002).

Figure 4: Time evolution of the semi-major axis (top panels), and eccentricity
(bottom panels), of all the simulated binaries. From left to right: comparison
between circular and eccentric binaries; comparison between co-rotating co-
planar and co-rotating misaligned binaries; comparison between co-planar and
counter-rotating binaries. Blue, purple and yellow curves refer to mass ratios
$q=1$, $q=1/3$, and $q=1/10$, respectively. The dark and light grey shaded
areas indicate the initial relaxation phases without and with the super-
Lagrangian refinement, respectively. The middle panels also include the
q03e00i45mod binary evolution. Overall we find that the interaction of the
binaries with their CBD leads to shrinking of the semi-major axis and
eccentricity growth, and expect the binaries to be driven towards merger.

Fig. 3 shows the time evolution of the CBD surface density, aspect ratio, and
Toomre parameter $Q$ profiles for our different simulations. Despite
fluctuations due to the formation of spiral arms, the CBD roughly maintains
the original $\propto R^{-2}$ profile over time. The surface density profiles
peak near the cavity edge, which we indicate by the vertical grey lines in
each panel. This can be compared to the vertical black dotted line, which
indicates the original location of the CBD inner edge and illustrates whether
or not the CBD moves outwards (as for prograde binaries) or inwards (as for
retrograde binaries). The surface density filling the cavity (i.e. to the left
of the grey vertical line) of binaries with mass ratios $q\geq 1/3$ and
$i=0\degr$ tends to decrease over time, particularly for circular binaries.
Similar but weaker behaviour is seen for $q=1/10$ binaries (i.e. run q10e00i00
and q10e05i00); however, the cavity surface density is typically larger,
therefore confirming the tendency of the $q=1/10$ cases to have more diffuse
material inside the cavity.

The evolution of the cavity gas content is determined by the balance between
the inflow/outflow of material from/to the CBD driven by the interaction with
the binary and accretion onto the black holes that slowly depletes the mini-
discs that are built up during the second relaxation phase of the simulation.
On the other hand, we see an opposite trend for counter-rotating binaries,
where, particularly in the case of the q10e00i180 binary, the surface density
inside the cavity slightly grows over time, as the torque exerted by the
binary is less effective in maintaining the cavity than in the prograde cases.
For the eccentric cases, the CBD edge is significantly further out and is not
so sharp due to interaction between the secondary black hole and the CBD.

For systems in which the binary and CBD are initially coplanar ($i=0\degr$ and
$i=180\degr$), the aspect ratio and Toomre parameter are generally constant
within the CBD for $R\lesssim 6-7a$. For these systems, the aspect ratio
typically ranges between $\approx 0.03$ and $\approx 0.06$ and the Toomre
parameter is $\approx 1.3-1.5$, with little dependence on the binary
properties. Such values are indeed suggestive of a gravito-turbulent steady
state in the CBDs, where spiral arms continuously form and disperse with no
fragmentation into self-gravitating bound clumps. However, deviations from
these values can be seen at early times for binaries that are not initially
coplanar with their CBD ($i=45\degr$), which undergo significant evolution as
the binary and warped CBD come into alignment (see Section 3.2.3). For all
binaries, inside the cavity, the aspect ratio grows to $\gtrsim 0.1$ because
of the presence of the mini-disc and the Toomre parameter $Q\gg 1$; this
indicates that the mini-discs are generally less dense and hence thicker than
the CBD as expected. However, the values of these quantities within about
$1.5a$ should be taken with a grain of salt as the geometry of the system
significantly departs from axisymmetry as assumed by the profiles. On the
other hand, outside $\sim 7a$, i.e. beyond the original extent of the CBD in
the initial conditions, all the profiles evolve significantly. We typically
see a progressive reduction of both the aspect ratio and the Toomre parameter
over time, corresponding to an increase in the surface density. This happens
because the very outer and under-dense regions of the CBD, produced by the
initial readjustment of the disc during the relaxation phase, are still
experiencing the slow action of radiative cooling. As the temperature
decreases, the outer regions of the discs become progressively thinner and
denser with time, explaining the evolution of the $\Sigma$, $H/R$ and $Q$
profiles.

###  3.2 Binary evolution and torques

####  3.2.1 General behaviour

While some studies of CBDs include analytic binaries, with the binary orbit
pre-determined and fixed (e.g., Muñoz et al., 2019; Muñoz et al., 2020; Siwek
et al., 2023b) in the simulations presented here we include a live binary (see
e.g., Cuadra et al., 2009; Roedig et al., 2011, 2012; Franchini et al., 2021)
whereby its evolution arises naturally through its interaction with the gas
via gravitational torques and accretion. Fig. 4 shows the time evolution of
every binary’s semi-major axis and eccentricity. The figure is further divided
into three columns that compare circular and eccentric binaries (left), co-
rotating co-planar and co-rotating misaligned binaries (middle), and co-planar
and counter-rotating binaries (right). As already seen in previous studies
that consider the evolution of a live binary in a self-gravitating CBD (Cuadra
et al., 2009; Roedig et al., 2012; Franchini et al., 2021), the general
behaviour is that the semi-major axis of the binary decreases and the
eccentricity increases over time. As a consequence, the angular momentum of
the binaries, defined as:

$J_{\rm bin}=\mu\sqrt{GMa\left(1-e^{2}\right)},$ (4)

where $\mu=M_{1}M_{2}/(M_{1}+M_{2})$ and $M=M_{1}+M_{2}$ are the binary
reduced and total mass, respectively222$M_{i}=M_{\bullet,i}+M_{{\rm d},i}$ is
the sum of the mass of the black hole and the accretion disc., decreases over
time (by $\sim$a few$-10\%$), even considering the mass evolution of each
component due to external accretion that may lead to an increase in the binary
angular momentum (Roedig et al., 2012). We should also caveat that, as we
discuss further in Section 3.2.2, the gravitational torques experienced by the
binary can be sensitive to numerical effects, in particular the choice of
gravity solver.

We do not observe any obvious trend in the evolution of $a$ as a function of
the mass ratio $q$; however, we notice that for eccentric, inclined and
retrograde $q=1/10$ binaries, their semi-major axis decreases most
significantly. As we discuss in Section 3.2.2, the evolution of the binary
separation and hence angular momentum is dominated by the gas gravitational
torques acting on the black holes. Different regions provide positive and
negative torques on the binaries and it is the balance of torques from gas in
the vicinity of the binary at $R<a$ (typically positive), and gas at $R>a$
that largely comprises of streams and the CBD (typically negative) that are
important in determining the net torque experienced by the binary. In all
cases where the $q=1/10$ binaries shrink more rapidly we find that they
experience strong negative torques (see Fig. 7), while the q10e00i00 binary
experiences a relatively weaker simulated torque. However, upon accurately
calculating torques in post-processing we find that the simulated torque is
likely an underestimate in this case and that the q10e00i00 binary should
likely shrink more rapidly (see Section 3.2.2 for further discussion).

Some recent studies that employ an exact analytic binary potential have been
used to explore binary eccentricity evolution, with the expectation being that
binaries with zero initial eccentricity should remain circular, while non-zero
initial eccentricities should either shrink or grow towards some critical
value (e.g., D’Orazio & Duffell, 2021; Zrake et al., 2021; Siwek et al.,
2023b). This being said studies that include live binaries find that
eccentricity can grow even when the binary is initially circular (e.g., Cuadra
et al., 2009; Franchini et al., 2021). We find mild eccentricity growth in all
of our binaries except for in the retrograde cases (see also, Tiede &
D’Orazio, 2023, for retrograde eccentricity evolution). The eccentricity
generally tends to grow slightly faster for the $q=1/3$ cases than for equal-
mass binaries, regardless of the other parameters. This tendency with the
binary mass ratio is inverted for the $q=1/10$ cases, where the eccentricity
growth is the slowest, except for the binaries starting with high initial
eccentricity $e=0.5$. Cuadra et al. (2009) found that eccentricity growth of
an initially circular binary only starts once the CBD develops asymmetries.
However, while our CBDs develop asymmetries fairly quickly, the eccentricity
for our initially circular binaries starts to grow from $t=0$. This being
said, recent work of Franchini et al. (2021), also finds mild early
eccentricity growth when simulating massive gravito-turbulent CBDs around
circular binaries.

Figure 5: Time-averaged surface density and torque maps for all the coplanar
runs between $200~{}t_{\rm bin}$ and $300~{}t_{\rm bin}$. Each simulation is
identified by its label from Table 1. For eccentric runs ($e=0.5$), both the
averages at the apocentre and at the pericentre are shown. The inner dashed
circle indicates the radius corresponding to the semi-major axis for runs with
$e=0$ and to the apo- and pericentre for the runs with $e=0.5$; the outer
dashed circle at $R=2a=4$ pc indicates the initial cavity rim location. It can
be seen that there are both positive and negative torques acting on the
binary, in the mini-discs, the central cavity, including the gas streams and
further out in the CBD. The spatial location and strength of these torques
depends on the binary mass ratio and orbital parameters.

A very different behaviour is shown by the counter-rotating binaries. Their
eccentricity does not significantly grow over time, as has been seen for
circular retrograde binaries in the recent work of Tiede & D’Orazio (2023).
However, Nixon et al. (2011a) discuss that the eccentricity growth depends on
accretion onto the secondary, showing that if the secondary is able to accrete
efficiently at apocentre (but not pericentre) then an initially circular orbit
can gain eccentricity. Given the secondary black holes in our retrograde
simulation accrete a negligible amount of mass (see Section 3.3.1), this
processes is not occurring in our simulations. The overdense wake lagging
behind the counter-rotating binary and the low eccentricity are features
reminiscent of the action of dynamical friction on a massive body embedded in
a rotating background (Dotti et al., 2007, 2009). Additionally, as the cavity
is more permeated by diffuse material, the counter-rotating black holes
experience a headwind gas motion. As we show in Section 3.2.2, the retrograde
binaries experience strong negative torques, with the continuous refilling of
the cavity leading to the steep and continuous shrinking. This is in contrast
to the prograde cases in which the cavity is slowly cleared of gas, the CBDs
expand and the shrinking of the binary appears to slow down.

####  3.2.2 Torque Analysis

Gravitational torques from the CBD and gas within the cavity, including the
mini-discs and streams, are critical in determining the binary evolution
(e.g., Cuadra et al., 2009; Roedig et al., 2012; Muñoz et al., 2019). In this
section, we calculate the gravitational torques exerted by the gas on the
binaries, and decompose them by their spatial contributions. This means we are
able to analyse the dominant origins of such torques and their net influence
on the binary separation. As well as allowing us to predict whether or not the
binaries should grow or shrink, comparisons between the simulated binary
evolution and post-processed torques from both direct-summation and tree based
force calculations let us assess the appropriateness of the numerical
techniques we employ.

Fig. 5 shows time-averaged maps of the gas surface density and torque surface
density for every coplanar run. Panels on the left-hand side show initially
circular binaries and are split into prograde (top) and retrograde (bottom)
runs, while panels on the right-hand side are for eccentric binaries and show
maps at both pericentre (top) and apocentre (bottom). For each run, the top
panel shows gas surface density, while the bottom panel shows the surface
density of the torque component directed along the binary angular momentum
vector. The maps are produced by averaging over $\sim 100$ snapshots between
$t=200~{}t_{\rm bin}$ and $300~{}t_{\rm bin}$. Torque surface density maps
show both positive (red) and negative (blue) torques with respect to the
binary angular momentum vector, which act to change the binary angular
momentum. The surface density maps further illustrate many of the features
discussed in Section 3.1 and provide a useful reference when considering which
gas primarily torques the binaries.

If we first consider the prograde, circular binaries, we find that all the
binaries experience torques from gas in the dense inner edge of the CBD and
just inside the cavity. In the case of $q=1/10$ binary, there are clear
torques from the gas that has moved into the cavity, with the stream feeding
the secondary clearly visible. We also note that there are significant torque
contributions from the gas in and around the black holes. For the retrograde
binaries, the inner edge of the CBD has moved inwards and there is more
diffuse gas within the cavity that torques the binary. For the $q=1/3$ case
only small, tenuous mini-discs form around the black holes, although gas
around the black holes and in the streams feeding them still contributes to
the torques. In the $q=1/10$ case a clear mini-disc forms around the primary
and contributions to the torques from gas around the primary and streams
feeding the secondary are seen. Finally, in the case of the eccentric
binaries, the torquing of the binary varies depending on the stage of the
orbit – at pericentre, when the black holes are furthest from the CBD inner
edge, the dominant contribution to the torques comes from gas around the black
holes, while at apocentre there are additional important contributions from
material in the CBD and from the streams and diffuse gas that exists in the
cavity.

Figure 6: Average radial profiles of the gas torque density, ${\rm
d}{\boldsymbol{\tau}}/{\rm d}R$, acting along $\hat{\bf J}_{\rm bin}$, for all
in-plane binaries. The left and right most columns show the circular pro-grade
and retrograde binaries, respectively, while the central two columns show the
initially eccentric binaries at pericentre (left) and apocentre (right). Each
row shows binaries with mass ratios of $q=1$, $1/3$ and $1/10$, from top to
bottom, respectively. Each profile is produced by averaging over $100$
orbits${}^{\ref{foot:orbit_ave}}$, while the line colours indicate the time
range as shown in the legend. Profiles are split into two radial regions
representing gas inside and outside the binary region ($R=a=2$ pc for circular
orbits or $R=a(1\pm e)=2\pm 1$ pc for eccentric runs) with solid/dashed lines
indicating whether the integrated net torque in each region is
positive/negative. Vertical dashed lines indicate the peak CBD density from
Fig. 3. In general we find that material close to the binary, i.e. comprising
the mini-discs, acts to torque up the binary, while gas at larger radii
provides a net negative torque. Overall, the net torque acting on the binaries
tends to be negative, driving them to shrink.

Taking a more quantitative approach, similar to Roedig et al. (2012), we
calculate average profiles of the gas torque density, ${\rm
d}{\boldsymbol{\tau}}/{\rm d}R$, acting on the binary, along the direction of
the binary angular momentum. From left to right, the columns in Fig. 6 show
circular prograde binaries, eccentric binaries at pericentre, eccentric
binaries at apocentre and circular retrograde binaries, respectively. From top
to bottom, the panels are for binary mass ratios of $q=1$, $1/3$ and $1/10$,
respectively (misaligned binaries are discussed in Section 3.2.3). For each
binary, we create five profiles, each of which is generated by averaging over
$100$ orbits333To ensure consistency with the eccentric binary profiles, for
which the apocentre/pericentre snapshot from each orbit is used, for the non-
eccentric binaries we randomly select one snapshot per orbit when calculating
the average profiles. The same approach is additionally taken for the profiles
in Fig. 10.. We additionally split each profile into the gas torques from two
radial regions; namely close to the binary (within $R=a=2$ pc for circular
orbits or $R=a(1\pm e)=2\pm 1$ pc for eccentric runs) indicated by the grey
shaded region, and at larger radii without shading. For each region, we
integrate the profile to find the net torque and use solid lines if the net
torque in the region is positive (i.e. would increase the binary angular
momentum) and dashed lines if the net torque in the region is negative (i.e.
would decrease the binary angular momentum). The vertical dashed line in each
panel corresponds to the CBD density peak from Fig. 3 and indicates the
location at which the CBD transitions to the cavity region.

All of the runs show spikes in the torque profiles at small radii
(corresponding to gas in and around the mini-discs), followed by smoother
variations in the cavity and CBD. Specifically, there is typically a negative
torque within the cavity followed by varying positive and negative peaks at
the outer edge of the cavity/inner edge of the CBD. First considering the
circular, prograde binaries, as highlighted by Roedig et al. (2012), while
peaks in the torque profiles at large radii (such that they are within the
main body of the CBD) could be associated with outer Lindblad resonances
(OLRs, see also, Artymowicz & Lubow, 1994), one should be careful applying
such an interpretation to peaks that appear within the cavity region (but
beyond the mini-discs). These are associated with gas streams and are expected
to have a kinetic origin. Gas within $R<a$ provides net positive torques on
the binary, while outside of this region, the gas provides a net negative
torque, dominated by the gas in the cavity and inner edge of the CBD. As we
discuss later (see Fig. 7 and Table 2), on average the net torque acting on
the binary is negative and acts to shrink it. This is due to the magnitude of
the net negative torque at $R\gtrsim a$ generally being larger than the
magnitude of the net positive torque at $R\lesssim a$, however, as we discuss
in Section 4 the exact balance is likely sensitive to the gas thermodynamics
(see e.g., Roedig et al., 2012; Franchini et al., 2022; Wang et al., 2023).

In the case of the eccentric binaries, we see that, as discussed in relation
to Fig. 5, at pericentre the torques are dominated by gas close to the black
holes, with its net torque being either positive or negative. At apocentre,
the profile features look somewhat more akin to the circular case, with net
torques from gas close to the black hole being positive but those from gas in
the cavity and CBD being negative. Again, on average the overall net torque
acting on the binary is negative.

While the rough location of peaks and troughs in the prograde cases remain the
same, their magnitude reduces over time, particularly in the central region.
This is likely a reflection of the declining average surface density in the
cavity (see Fig. 3), which arises due to a combination of black hole accretion
and the CBD moving outwards due to being torqued by the binary (see discussion
in Section 3.3.1).

In Table 2 we provide the average torques experienced by each binary (directed
along its own angular momentum vector) split into separate time bins of $\sim
100$ orbits each and decompose the contributions into being from gas at
$R<a_{\rm sep}$, $R>a_{\rm sep}$, and all gas. Specifically, $a_{\rm sep}$ is
the measured binary separation at the time torques are calculated, i.e. this
would be $a$ for circular binaries and vary between $a(1-e)$ and $a(1+e)$ for
eccentric binaries moving between apo- and pericenter, taking into account the
evolution of $a$ and $e$ over the course of the simulation (see Fig. 4). The
colour coding scales linearly from $0$ to $\pm 5$, i.e., brighter colours
represent stronger torques, with red and blue representing positive and
negative torques, respectively. We find that in the majority of cases, even
though the absolute average torque exerted from inner and outer regions
typically decrease over time, the net gravitational torque from all gas only
shows relatively mild evolution (although short-term fluctuations do occur,
see Fig. 7).

Finally, the retrograde circular binary exhibits many more sharp variations in
the torque profiles, indicative of the more complex distribution of gas that
fills the cavity region. As discussed in Nixon et al. (2011a), a retrograde
circular binary does not excite resonances within the CBD, and so peaks and
troughs seen here should not be interpreted as such. The fact that the radial
location of peaks and troughs within the CBD region itself vary further
supports this by showing that the torques are not produced from consistent
well-defined resonant radii. Overall, while throughout the simulation some
regions switch from providing net positive and net negative torques to the
binary, on average the binary experiences a strong negative torque, as shown
quantitatively in Table 2. Specifically, we find that the torques within the
cavity and CBD roughly cancel each other out in the retrograde run and it is
the gas within the vicinity of the black holes that dominates the net torquing
of the binary.

Figure 7: Average values for the torque component along the binary angular momentum vector experienced by each binary compared to the estimated gravitational and accretion contributions. Coloured squares show the average torques over $500$ binary orbits while error bars show the one-$\sigma$ variation on an orbit by orbit basis. Specifically we show the actual torque experienced by the binary based the evolution of its angular momentum (yellow, $\tau_{\rm bin}^{\rm sim}$ as well as two estimates for the gravitational torque using either a pure tree approach (blue, $\tau_{\rm grav}^{\rm tree}$) or pure direct summation (purple, $\tau_{\rm grav}^{\rm direct}$) and an estimated accretion torque (pink, $\tau_{\rm accretion}^{\rm est}$). Shaded squares show the estimated total torque experienced by the binary when adding $\tau_{\rm accretion}^{\rm est}$ to either $\tau_{\rm grav}^{\rm tree}$ or $\tau_{\rm grav}^{\rm direct}$. The average gravitational torques are always negative, while accretion torques are always positive. All total average torque experienced by binaries in the simulation is negative, with most values being generally consistent with those found by post processing methods. Some binaries, particularly misaligned binaries can show larger differences between actual and estimated torques, however, only the $q01e00i45$ binary exhibits a qualitative disagreement between the method in terms of the sign of the overall torque. Table 2: Averages of the gravitational torque component along the binary angular momentum vector (in units of $10^{-6}GM_{\rm bin}^{2}a^{-1}$) experienced by the different binary systems. We show the total torque exerted by all gas, as well as individual contributions from gas at $R<a_{\rm sep}$ and $R>a_{\rm sep}$, where $a_{\rm sep}$ is the measured binary separation at each time torques are calculated ($a$ for circular binaries and varying between $a(1-e)$ and $a(1+e)$ for eccentric binaries). Averages are calculated in five time bins of $100$ orbits each. From left to right the columns show: simulation label, region over which torque is calculated, and five columns for separate time bins. The colour coding scales linearly from $0$ to $\pm 5$, i.e., whereby brighter colours represent stronger torques, with red and blue representing positive and negative torques, respectively. In the majority (although not all) of cases the absolute torques exerted by the inner and outer regions tend to decrease, although only mild evolution is seen in the overall average gravitational torque acting on the binary. Label | Region | $0-100\times t_{\rm bin}$ | $100-200\times t_{\rm bin}$ | $200-300\times t_{\rm bin}$ | $300-400\times t_{\rm bin}$ | $400-500\times t_{\rm bin}$
---|---|---|---|---|---|---
Circular prograde binaries
q01e00i00 | all gas | $-1.4$ | $-1.173$ | $-1.312$ | $-1.335$ | $-1.314$
| $R<a_{\rm sep}$ | $1.893$ | $2.331$ | $1.586$ | $1.378$ | $1.051$
| $R>a_{\rm sep}$ | $-3.293$ | $-3.504$ | $-2.898$ | $-2.712$ | $-2.364$
q03e00i00 | all gas | $-0.975$ | $-1.238$ | $-1.821$ | $-1.797$ | $-1.922$
| $R<a_{\rm sep}$ | $3.279$ | $2.569$ | $1.895$ | $1.179$ | $0.858$
| $R>a_{\rm sep}$ | $-4.254$ | $-3.806$ | $-3.715$ | $-2.976$ | $-2.779$
q10e00i00 | all gas | $-1.279$ | $-1.224$ | $-1.108$ | $-1.109$ | $-0.941$
| $R<a_{\rm sep}$ | $2.839$ | $2.55$ | $2.108$ | $1.965$ | $1.949$
| $R>a_{\rm sep}$ | $-4.118$ | $-3.774$ | $-3.216$ | $-3.074$ | $-2.89$
Eccentric binaries
q01e05i00 | all gas | $-0.13$ | $-0.556$ | $-1.077$ | $-1.164$ | $-0.414$
| $R<a_{\rm sep}$ | $4.805$ | $2.797$ | $1.69$ | $1.357$ | $1.644$
| $R>a_{\rm sep}$ | $-4.929$ | $-3.348$ | $-2.763$ | $-2.519$ | $-2.055$
q03e05i00 | all gas | $-0.809$ | $-1.127$ | $-1.041$ | $-1.235$ | $-1.226$
| $R<a_{\rm sep}$ | $3.475$ | $2.697$ | $2.144$ | $1.665$ | $1.35$
| $R>a_{\rm sep}$ | $-4.266$ | $-3.814$ | $-3.181$ | $-2.898$ | $-2.574$
q10e05i00 | all gas | $-3.225$ | $-3.216$ | $-2.731$ | $-2.604$ | $-2.858$
| $R<a_{\rm sep}$ | $1.191$ | $1.075$ | $0.89$ | $0.822$ | $0.602$
| $R>a_{\rm sep}$ | $-4.41$ | $-4.29$ | $-3.62$ | $-3.427$ | $-3.46$
Misaligned binaries
q01e00i45 | all gas | $-0.913$ | $0.355$ | $-0.149$ | $-0.708$ | $-0.826$
| $R<a_{\rm sep}$ | $4.792$ | $4.102$ | $2.719$ | $1.798$ | $1.797$
| $R>a_{\rm sep}$ | $-5.706$ | $-3.747$ | $-2.867$ | $-2.506$ | $-2.623$
q03e00i45 | all gas | $-1.045$ | $-0.347$ | $-0.135$ | $-0.041$ | $-0.834$
| $R<a_{\rm sep}$ | $4.705$ | $3.463$ | $2.963$ | $3.367$ | $2.65$
| $R>a_{\rm sep}$ | $-5.753$ | $-3.809$ | $-3.097$ | $-3.409$ | $-3.484$
q03e00i45mod | all gas | $-0.423$ | $-0.477$ | $-0.913$ | $-1.14$ | $-0.512$
| $R<a_{\rm sep}$ | $4.045$ | $4.893$ | $2.018$ | $0.615$ | $1.181$
| $R>a_{\rm sep}$ | $-4.467$ | $-5.37$ | $-2.931$ | $-1.755$ | $-1.692$
q10e00i45 | all gas | $-2.069$ | $-1.902$ | $-2.124$ | $-1.499$ | $-0.961$
| $R<a_{\rm sep}$ | $2.531$ | $3.149$ | $3.177$ | $4.706$ | $5.182$
| $R>a_{\rm sep}$ | $-4.597$ | $-5.05$ | $-5.3$ | $-6.204$ | $-6.143$
Retrograde binaries
q03e00i180 | all gas | $-2.146$ | $-2.74$ | $-2.861$ | $-2.688$ | $-2.534$
| $R<a_{\rm sep}$ | $-2.107$ | $-2.842$ | $-2.963$ | $-2.775$ | $-2.6$
| $R>a_{\rm sep}$ | $-0.033$ | $0.108$ | $0.108$ | $0.092$ | $0.071$
q10e00i180 | all gas | $-2.689$ | $-3.804$ | $-3.676$ | $-3.004$ | $-2.007$
| $R<a_{\rm sep}$ | $-2.187$ | $-3.08$ | $-2.99$ | $-2.433$ | $-1.652$
| $R>a_{\rm sep}$ | $-0.502$ | $-0.724$ | $-0.685$ | $-0.57$ | $-0.355$

To get a sense of the overall importance of the gravitational torques on the
evolution of the binary we can compare the total gravitational torque exerted
on the binary by the gas to the accretion torque (the effective torque
experienced by the binary due to linear momentum conservation when accreting
gas) and to the simulated rate of change of the binary angular momentum. The
total gravitational torque acting on the binary is given by

${\boldsymbol{\tau}}_{\rm grav}=M_{1}({\bf r}_{1}\times{\bf a}_{1})+M_{2}({\bf
r}_{2}\times{\bf a}_{2})\,,$ (5)

where $M_{i}$, ${\bf r}_{i}$ and ${\bf a}_{i}$ are the mass, position and
gravitational acceleration444In practice, ${\bf a}_{i}$ can be either the
total gravitational acceleration or that only due to the gas, because for each
black hole the acceleration due to the other black hole provides no
contribution to the total torque in Equation (5). of each black hole particle
(which includes the sub-grid $\alpha$ accretion disc) and $i=1$ or $2$ refer
to the primary and secondary black hole, respectively. Additionally, the
linear momentum of accreted material is added to each black hole during each
accretion event, which results in a change in the $i^{th}$ black hole’s
angular momentum of

$\Delta{\bf J}_{{\rm acc},i}=\Sigma_{j}\left({\bf r}_{i}\times m_{j}{\bf
v}_{j}\right)\,,$ (6)

where $m_{j}$ and ${\bf v}_{j}$ are, respectively, the mass accreted from and
velocity of the the $j^{th}$ cell and the sum is performed over all cells from
which gas is accreted. In practise we do not log every accretion event
explicitly and for the purpose of analysing the torques experienced by the
binary we estimate it from the simulation output. To do this, we calculate the
expected specific angular momentum contribution of accreted gas for each black
hole, i.e., ${{\bf L}}_{{\rm gas},i}=\Delta{\bf J}_{{\rm
acc},i}/\Sigma_{j}m_{j}$ from each snapshot and estimate the total change in
angular momentum due to accretion as

$\Delta{\bf J}_{\rm acc}^{\rm est}=\Delta M_{1}\bar{{\bf L}}_{{\rm
gas},1}+\Delta M_{2}\bar{{\bf L}}_{{\rm gas},2}\,,$ (7)

where $\Delta M_{i}$ and ${\bar{\bf L}}_{\rm gas,i}$ are, respectively, the
change in mass of the $i^{th}$ black hole particle and the specific angular
momentum of accreted gas averaged over consecutive snapshots. We then define
an effective accretion torque as ${\boldsymbol{\tau}}_{\rm accretion}^{\rm
est}=\Delta{\bf J}_{\rm acc}^{\rm est}/\Delta t$, where $\Delta t$ is the time
between snapshots.

Torques can act to change both the magnitude and direction of ${\bf J}_{\rm
bin}$. To isolate the torque that acts only to change the magnitude of the
binary angular momentum, which determines whether the binary shrinks or
expands, we take only the component along the binary angular momentum vector.
For in-plane binaries this is essentially the total torque, however, this is
not the case for misaligned binaries (see Section 3.2.3). We calculate
averages for the relevant torque component over 500 binary orbits and show
values found for each simulation by the coloured squares in Fig. 7. For the
gravitational torques, we make two estimates, one in which the black hole
accelerations are calculated via a pure tree method ($\tau_{\rm grav}^{\rm
tree}$, shown in blue) and one in which the accelerations are calculated via
direct summation over all gas cells ($\tau_{\rm grav}^{\rm direct}$, purple).
This allows us to compare these methods with each other and to the simulated
binary angular momentum evolution – which in general shows reasonably
consistent results. Estimated accretion torques, $\tau_{\rm accretion}^{\rm
est}$, are shown by the pink squares, which we additionally add to the
gravitational torque estimates and show as faded squares to illustrate the
total expected torque in each case. Finally, by directly measuring the change
in $|{\bf J}_{\rm bin}|$ from the simulations we can calculate the average
total torque, $\tau_{\rm bin}^{\rm sim}$, experienced by the binary along its
angular momentum vector, which is shown by the yellow squares. The vertical
bars represent the one sigma scatter in orbit averaged values, i.e. how much
the torques vary on an orbit-by-orbit basis.

In the majority of cases, the gravitational torques measured using the tree
versus direct summation are in reasonable agreement with each other and, at
least for coplanar binaries, with torques derived from the binary evolution
when taking estimated accretion torques into account. However, non-negligible
differences between the tree and direct summation method are evident in some
runs, particularly for misaligned binaries. In all cases, the mean torque
experienced by the binary is closer to the value calculated using the tree
approach. This suggests that for these runs the hybrid force calculation we
employ does not perform so well at recovering the torque one would expect from
a pure direct summation force calculation and would indicate that we may
require a larger value of $N_{\rm direct}^{\rm max}$. This being said, it does
result in less variation in the torques experienced by the binary that are in
general more consistent with variations expected from direct summation
calculations. For all binaries, we find that the gravitational torque
estimates result in the same qualitative behaviour in providing negative
torquing of the binary.

Referring back to Fig. 4 and the discussion in Section 3.2, apart from in the
circular prograde systems, the $q=1/10$ mass ratio binaries shrink more
rapidly compared to their higher mass ratio counterparts. From Fig. 7 we can
now see that this is because the $q=1/10$ binaries typically experience strong
negative net torques leading to rapid extraction of their angular momentum. In
the case of the q10e00i00 binary, the average simulated torque is likely
under-estimated by a factor of $\sim 2$, and so like the other $q=1/10$
binaries, we would expect it to show more rapid shrinking than its higher mass
ratio counter-parts (although still likely less rapidly than other $q=1/10$
systems). While we reserve a discussion of torques in the misaligned binaries
to Section 3.2.3, we note here that the faster shrinking of inclined binaries
compared to their in-plane counterparts exhibited in Fig. 4 is potentially due
to over-estimating the gravitational torques acting on the binary, with the
predicted behaviour of the q01e00i45 binary from this torque analysis to be
that it gains angular momentum opposed to losing it. These small
inconsistencies highlight the need to accurately calculate gravity forces on
the fly for full self-consistent evolution of a live binary. We also note,
however, that our estimates for the accretion torques may be over-smoothed as
we only estimate their values between snapshots and may be missing episodes of
low or high angular momentum gas accretion at intermediate times. Finally, for
retrograde binaries, we find a negligibly small accretion torque. This is
because the accretion prescription we employ, which requires net gas inflow
onto the black hole, results in very little accretion onto the secondary (see
Section 3.3.1). If we were to use less conservative accretion criteria for
which more accretion occurs onto the secondary black hole, we would expect it
to accrete material with angular momentum opposite to that of the binary and
hence lead to more rapid shrinking of the binary (e.g., Nixon et al., 2011a).

####  3.2.3 Misaligned binary evolution

Figure 8: Surface density maps for the q03e00i45mod binary illustrating the
evolution of the CBD and black hole binary. The $x-y$ and $x-z$ projections
are shown in the top and bottom rows, respectively with the time of the image
shown in the top right of each panel in the bottom row. The primary and
secondary black holes are shown by the blue and red points, while the green
and white arrows show the angular momentum vectors of the binary and CBD,
respectively. As the system evolves, the binary and CBD torque each other such
that their angular momenta move to align, resulting in the production of a
warp in the CBD that propagates outwards with time.

So far our analysis has focused on binaries for which the angular momentum
vector is aligned or anti-aligned with that of the CBD. In this section, we
consider binaries that are initially inclined by $45\degr$ with respect to the
CBD. Fig. 8 shows evolving density projections for the q03e00i45mod run in the
$x-y$ (top) and $x-z$ (bottom) planes. As is most clearly seen in the bottom
row, the binary and gas torque each other in such a way as to align the
angular momentum of the CBD and the binary, as shown by the white and green
arrows, respectively. This proceeds by driving a warp through the CBD such
that gas at the inner edge of the disc aligns first and gas at progressively
larger radii aligns as the system evolves. This process has been explored in
several previous works (e.g., Nixon et al., 2013; Dunhill et al., 2014; Aly et
al., 2015; Moody et al., 2019), and is analogous to the Bardeen-Petterson
effect that acts to align black holes and their accretion discs (Bardeen &
Petterson, 1975; King et al., 2005; Nixon et al., 2011b).

The evolution of the torques acting on misaligned binaries is explored in Fig.
9, which shows the evolution of the $x$, $y$ and $z$ components of the gas
torque acting on the binary. We calculate both the torque experienced by the
binary from the rate of change of the binary angular momentum ($\tau_{\rm
bin}^{\rm sim}$, solid) and that calculated in post-processing using a direct
summation method ($\tau_{\rm grav}$, dashed), as well as adding the estimated
accretion torque ($\tau_{\rm grav}+\tau_{\rm acc}$, dotted, although this is
essentially indistinguishable from the $\tau_{\rm grav}$ curve). The lines
represent the torque averaged over $10$ orbits while the shaded regions
indicate the one-$\sigma$ scatter. Each panel represents a different
simulation, with all four initially misaligned binaries shown, as labelled.
There are only very small differences in $\tau_{\rm bin}^{\rm sim}$ and
$\tau_{\rm grav}$, showing that the instantaneous evolution of the binary
angular momentum is well captured. Additionally, there is a negligible impact
from accretion torques, showing that the global binary angular momentum
evolution is dominated by the gravitational interaction with the gas. However,
as we discuss below this is not necessarily the case for the magnitude of the
angular momentum. The binaries initially experience large torques, with the
amplitude of each component oscillating and damping as the binary and CBD
align. As discussed in Section 2.1, the runs in which the binary is initially
misaligned undergo some realignment during the relaxation phase and to verify
this has a limited impact on our results we perform an additional run,
q03e00i45mod, in which we artificially misalign the cavity region of the
q03e00i00 run to be at $45\degr$ to the CBD after the initial relaxation
phase. Comparing this run to the q03e00i45 run, we see that they show similar
behaviour albeit with a time offset. The offset roughly corresponds to $100$
orbits, i.e., the relaxation time, although this is not exact.

Isolating the torques that solely change the magnitude of the binary angular
momentum, Fig. 10 shows, as in Fig. 6, average radial profiles of the gas
torque density, ${\rm d}{\boldsymbol{\tau}}/{\rm d}R$, acting along the
direction of the binary angular momentum vector for the inclined runs. Each
panel represents a different binary and colours illustrate different time
ranges with each profile being produced by averaging over $100$
orbits${}^{\ref{foot:orbit_ave}}$. The profiles show broadly similar trends to
their in-plane circular counterparts with several sharp features from gas
within the binary orbit ($R<a=2$ pc) and smoother peaks and troughs at larger
radii. Additionally, as in the in-plane binaries, the gas at smaller radii
provides a net positive torque to the binary, while gas at larger radii
provides a net negative torque, which is also shown quantitatively in Table 2.
Comparing to the equivalent coplanar binaries, we expect from the direct
summation calculations, that the misaligned $q=1$ and $q=1/3$ systems should
experience stronger positive torques in the $R<a_{\rm sep}$ region and
generally less negative gravitational torques overall. The misaligned binaries
typically also experience mildly stronger positive accretion torques. For the
q01e00i45 system in particular, combining the estimated accretion torque with
the torque estimated from direct summation would lead to a net positive torque
on the binary. This suggests that accretion torques could play a particularly
important role in determining the binary evolution for inclined systems.
However, as we outline in Section 4, even though the total torque acting on
the q01e00i45 binary that we estimate in post-processing is positive, we still
expect the binary to shrink.

We finally note that comparing the q03e00i45 and q03e00i45mod runs we find
that the ‘modified’ run experiences a stronger negative torque leading to its
more rapid shrinking at later times as seen in Fig. 4. However, it is worth
noting that the post-processed torques estimated from direct summation are
much more similar between these two runs, and differences in the accretion
torques that may lead to different binary evolution are likely an artefact of
how the modified binary was created, as at least initially, it continues to
experience somewhat coherent accretion from the mini-discs that were
artificially realigned with the binary after relaxation, while the q03e00i45
run experiences a less coherent flow throughout its evolution (see Section
3.3.2 and Fig. 14).

Figure 9: Evolution of the $x$ (blue), $y$ (pink) and $z$ (yellow) components
of the torque experienced by the binary in misaligned cases, as shown by the
solid line. For comparison the post processed gravitational torque experienced
by the binary estimated with a direct summation method is shown by the dash
lines, while the estimated gravitational plus accretion torque is shown by the
dotted lines. The actual and post-processed torques match very well and show
that the torquing is dominated by the gravitational component. Initially large
torques oscillate and damp as the binary and CBD angular momentum realign. The
q03e00i45mod run shows a $\sim 100$ $t_{\rm bin}$ offset compared to the other
misaligned binaries due to their misalignment being initialised after the
relaxation. Figure 10: As in Fig. 6 but showing the average radial profiles of
the gas torque density, ${\rm d}{\boldsymbol{\tau}}/{\rm d}R$, acting along
$\hat{\bf J}_{\rm bin}$ for all initially misaligned binaries. Each panel is
labelled with the binary run name, the time period over which each profile is
calculated is indicated by the colours given in the legend and as in Fig. 6,
solid and dashed lines represent regions that exert positive and negative net
torques on the binary, respectively. Vertical dashed lines indicate the peak
CBD density from Fig. 3. The profiles generally show similar features to the
in-plane binaries, with material close to the binary typically exerting net
positive torques, while material at larger radii exert net negative torques.
Figure 11: Evolution of inflow and accretion rates onto each black hole.
Specifically the primary and secondary black holes are shown in blue and red,
respectively, with shaded regions indicating the range of inflow rates
calculated onto the $\alpha$ accretion discs, while the solid (red and blue)
lines show the accretion rate onto each black hole from the sub-grid $\alpha$
accretion disc. Additionally, inflow/outflow rates through the cavity are
shown by the green solid/dashed lines, respectively. In general, primary black
holes experience higher inflow rates and accretion rates. Inflow rates onto
secondary black holes, as well as onto black holes in high eccentricity
binaries, exhibit more fluctuations than inflows onto primary black holes and
for circular binaries. Figure 12: Time sequence of gas surface density maps
focusing on the cavity region for one representative orbit of an initially
circular (q10e00i00, top row) and eccentric (q10e05i00, bottom row) binary.
The sequences show more erratic dynamics in the eccentric case, with
triggering, development and destruction of gas streams reaching the secondary
black hole first and then partially feeding the mini-disc surrounding the
primary black hole.

###  3.3 Black hole and accretion disc evolution

Having considered the evolution of the binary and its interaction with the
CBD, in this section, we delve deeper into the evolution of the black holes
and their sub-grid accretion discs. Specifically, we look at gas accretion
onto the black holes and hence their mass as well as spin evolution. This is
important as both the binary mass ratio, the magnitude and the orientation of
the black hole spins have important consequences for gravitational wave
emission and the remnant black hole properties, including its recoil velocity.

####  3.3.1 Accretion

Considering accretion onto the black holes, Fig. 11 shows mass inflow rates
onto each black hole particle by the blue (primary) and red (secondary) shaded
regions. Specifically, these regions show the $5^{\rm th}$ to $95^{\rm th}$
percentile range of the inflow rates within time bins of length $5~{}t_{\rm
bin}$. The corresponding solid blue and red lines show the accretion rate
through the sub-grid accretion disc onto each black hole. While there are
fluctuations in the inflow rates on to the accretion disc-black hole system,
the rate at which material ultimately reaches a black hole is mediated by the
sub-grid accretion disc (see Equation (8) below for analytic form of Eddington
rate through the accretion disc). This results in the black hole accretion
rates being smoother, and at times lower, than the instantaneous inflow rates.
The green line shows mass flow rates through the cavity, calculated at
$R=2a=4$ pc for circular binaries and $R=3a=6$ pc for eccentric binaries.
Solid lines show net inflow and dashed lines show net outflow, both averaged
over $\sim 10$ binary orbits. Mass flows through the inner edge of the CBD and
cavity fluctuate between being inflowing and outflowing, with the net rate
being roughly consistent with, although generally slightly lower than, the
binary accretion rate. This is likely an effect of the relaxation period we
impose during which accretion is prohibited and the mini-discs build up
reservoirs from which the black holes subsequently accrete. While in our
simulations the CBD provides a roughly constant supply of material to the
cavity, leading to a somewhat steady-state, in reality black hole and/or
stellar feedback, or a more dynamic environment, may change this scenario.

As discussed in Section 2.4, the sub-grid accretion disc model requires black
hole accretion rates to be initialised and are set to the values provided in
Table 1 such that $\dot{M}_{\bullet,i}=f_{\rm Edd,i}\times M_{\bullet,i}$.
Subsequent evolution then depends on the mass and angular momentum evolution
of the sub-grid accretion disc, which itself depends on the resolved gas
inflows feeding the disc. In all cases, we find inflow rates onto the primary
black hole to exceed those onto the secondary. Inflow rates onto the secondary
black hole in a system often exhibit much larger fluctuations when compared to
inflow rates onto the primary, which are typically smoother. Counter to our
findings, previous works have found that due to the secondary black hole
typically inducing a stronger gas stream from the cavity it experiences
preferential accretion over the primary (e.g., Roedig et al., 2012; Farris et
al., 2014; Duffell et al., 2020; Muñoz et al., 2020; Dittmann & Ryan, 2021;
Siwek et al., 2023a). We expect that differences in our findings are related
to the modelling of the gas thermodynamics within the cavity, with many of the
mentioned works employing a (locally-)isothermal equation of state, and the
exact modelling of the black hole accretion, with previous works making use of
sink particle prescriptions opposed to the inflow rate calculation that we
employ. We discuss these points and their implications for our results further
in Section 4.3. In the most extreme case, inflow rates onto the retrograde
binaries, particularly onto the secondaries can be greatly inhibited. This
leads to declining black hole accretion rates as the sub-grid accretion disc
mass is depleted. This is likely a result of our conservative accretion
criterion of requiring net inflow onto the black hole and may have important
implications for the evolution of retrograde binaries (see e.g., Nixon et al.,
2011a).

Comparing the top two rows, inflow rates onto eccentric binaries show
increased scatter due to the more complicated gas dynamics within the cavity
compared to the circular binaries. This is illustrated in Fig. 12, which shows
a face-on gas density projection for the q10e00i00 and q10e05i00 binaries in
the top and bottom rows, respectively. In the circular orbit case, the streams
feeding both black holes and the primary mini-disc persist throughout the
orbit. However, for the eccentric binary, the streams form and disperse, and
although a mini-disc remains around the primary it is much less extended and
receives new gas from streams more sporadically than in the circular binary
case. The misaligned binaries also give rise to more complex gas dynamics in
the cavity, however, as shown in the third row of Fig. 11, instead of leading
to a general increase the scatter of inflow rates (compared to the in-plane
circular binaries), exhibit more long time variations, particularly onto the
secondary black holes. One feature to note is for the q10e00i45 run, in which
the sub-grid accretion rate onto the primary experiences a slight spike
followed by a sudden drop at $\sim 120$ orbits. This is simply down to the
nature of the equation determining the black hole accretion rate,
$\dot{M}_{\bullet}=f_{\rm Edd}M_{\bullet}/\eta\tau_{\rm Salp}$, where from
Fiacconi et al. (2018) the dimensionless Eddington rate is given as

$f_{\rm Edd}\approx 0.76\\\
\times\left(\frac{\eta}{0.1}\right)\left(\frac{M_{\rm
d}}{10^{4}M_{\odot}}\right)^{5}\left(\frac{M_{\bullet}}{10^{6}M_{\odot}}\right)^{-47/7}\left(\frac{a_{\bullet}J_{\rm
d}/J_{\bullet}}{3}\right)^{-25/7}\,,$ (8)

where $\eta$ is the spin-dependant radiative efficiency and $\tau_{\rm Salp}$
is the Salpeter time. Variations in $f_{\rm Edd}$ are dominated by changes in
the disc mass, $M_{\rm d}$ and angular momentum, $J_{\rm d}$, which both
decrease during the first $\sim 120$ orbits as $\dot{M}_{\bullet}>\dot{M}_{\rm
in}$ and almost balance with each other. However, after a brief increase in
the accretion rate the $M_{\rm d}$ decreases sufficiently to result in a sharp
decline in the black hole accretion rate. After this
$\dot{M}_{\bullet}<\dot{M}_{\rm in}$ and the accretion rate stabilises for the
remainder of the simulation. A similar although less dramatic process is seen
after $\sim 140$ orbits in the q03e00i45mod run.

####  3.3.2 Spin evolution

Figure 13: Spin alignment and binary inspiral timescales extracted from the
simulations and based on the analytical model of Gerosa et al. (2015).
Analytic predictions are illustrated by the shaded blue and red regions, which
show spin-alignment timescales as a function of binary mass ratio for the
primary and secondary black hole, respectively, and the dotted green line,
which shows the binary inspiral time. Simulated black hole alignment time-
scales (with associated scatter) are shown by the blue and red points for the
primary and secondary black hole, respectively, with different runs indicated
by the lines styles. The range of binary inspiral times found in the
simulations are shown by the green hatched region. The top and bottom panels
differ by the assumed accretion rate ratio between the primary and secondary
black hole, $f_{\rm Edd,2}/f_{\rm Edd,1}=q^{-\beta}$, with $\beta=2.0$ and
$1.0$, respectively. In general we find that black hole spins align with their
accretion discs on timescales of $\lesssim$ the binary inspiral time and in
the majority of cases much more quickly. The lack of preferential accretion
onto secondary black holes means that spin alignment can also efficiently
occur for the primary black hole in low mass ratio systems.

As well as modelling black hole mass growth, the sub-grid black hole model
allows us to study the evolution of the black hole spin and angular momentum
of the sub-grid accretion disc. In particular, an understanding of black hole
spin alignment is important to determine the recoil velocities of black hole
merger remnants (Campanelli et al., 2007; González et al., 2007; Lousto &
Zlochower, 2011, 2013; Lousto & Healy, 2019; Sperhake et al., 2020).
Specifically, relativistic numerical simulations predict that very high recoil
velocities, up to $\sim 5000$ km s-1, can be achieved if the black hole spins
are large and misaligned (Lousto & Zlochower, 2011, 2013; Lousto & Healy,
2019). Such velocities would be able to eject black holes even from very
massive galaxies, remove them from the pool of future black hole mergers and
hence impact the GW background (Rajagopal & Romani, 1995; Phinney, 2001; Enoki
et al., 2004; Sesana et al., 2004). On the other hand, if the black hole spins
are both aligned with the binary angular momentum, much smaller recoil
velocities are expected. Furthermore, black hole spin magnitude and
orientation is fundamental if we want to study black hole-driven jet (and
associated wind) feedback (see e.g., Tchekhovskoy et al., 2011; Yuan &
Narayan, 2014; Liska et al., 2017; Talbot et al., 2021, 2022).

Gerosa et al. (2015) developed an analytic model of black hole binary
evolution in a CBD assuming the system is co-planar. Based on this model, the
timescale on which the spin of each black hole would align with its accretion
disc’s angular momentum is given as

$\tau_{\rm align}=3.4\frac{\tau_{\rm Salp}}{f_{\rm CBD}}\frac{\alpha_{\rm
AD}}{F(q)}\left(\alpha_{2}\frac{H}{R}\right)^{2/3}\\\
\times\begin{cases}a_{\bullet,1}^{2/3}\left(1+q^{\beta-1}\right)/\left(q^{\beta-1}+q^{\beta}\right)&\text{Myr
[Primary]}\\\
a_{\bullet,2}^{2/3}\left(q+q^{\beta}\right)/\left(1+q\right)&\text{Myr
[Secondary]}\,,\end{cases}$ (9)

where $f_{\rm CBD}=\dot{M}_{\rm CBD}\tau_{\rm Salp}/M_{\rm bin}$ is the
dimensionless Eddington rate of material feeding the cavity from the CBD,
$\alpha_{\rm AD}=0.1$ is the assumed $\alpha-$viscocity of the accretion disc,
$\alpha_{2}$ is the vertical $\alpha-$viscosity parameter given as (Pringle,
1992; Ogilvie, 1999)

$\alpha_{2}=\frac{1}{\alpha}\frac{2\left(1+7\alpha^{2}\right)}{4+\alpha^{2}},$
(10)

$H/R$ is the accretion disc aspect ratio, which for this analysis we set to
$0.01$, and $F(q)$ is the fraction of material flowing into the cavity that is
accreted by the binary components and takes the form
$F(q)=0.8054+0.984\log(q)+0.3818\log^{2}(q)$, where the numerical coefficients
were found by Gerosa et al. (2015) as best fits to the simulation results of
D’Orazio et al. (2013). Finally, $\beta$ determines how material accreted by
the binary is distributed between the primary and secondary black holes such
that $f_{\rm Edd,2}/f_{\rm Edd,1}=q^{-\beta}$, whereby values of $\beta>1$
imply preferential accretion onto the secondary black hole. Gerosa et al.
(2015) compared the alignment timescales to the predicted binary inspiral time
given by

$\tau_{\rm in}=\frac{3}{4}\frac{\left(1+q\right)M_{\rm cbd}/M_{\rm
bin}+1}{\left(1+q\right)^{2}}\frac{\tau_{\rm Salp}}{f_{\rm CBD}}\text{ Myr}.$
(11)

The full derivation of this equation is given in Gerosa et al. (2015), but in
essence is an interpolation between the analytical predictions for binary
inspiral times in three different mass regimes, either $M_{\bullet,2}<<M_{\rm
cbd}<<M_{\bullet,1}$ (Lin & Papaloizou, 1979; Artymowicz & Lubow, 1994;
Armitage & Natarajan, 2002), $M_{\bullet,2}\sim M_{\rm cbd}<<M_{\bullet,1}$
(Syer & Clarke, 1995; Ivanov et al., 1999; Lodato et al., 2009; Baruteau &
Masset, 2013) and $M_{\bullet,2}\lesssim M_{\bullet,1}$ (Rafikov, 2013).

Applying the Gerosa et al. (2015) model to our setup parameters, in Fig. 13,
we plot the predicted spin alignment timescales as a function of binary mass
ratio, $q$, for the primary and secondary black hole in the blue and red
shaded regions, respectively, while the green dotted line shows the predicted
inspiral time. The alignment timescales are represented by shaded regions to
illustrate the $\sim 2.5$ uncertainty in warp propagation theory. The top
panel shows the results of the analytic modelling assuming preferential
accretion onto the seconday black hole, i.e., $\beta=2$, as originally assumed
by Gerosa et al. (2015). This model leads to inhibited accretion onto the
primary for low $q$ binaries and results in spin alignment timescales for the
primary black hole that exceed the inspiral time and hence spin alignment is
not expected for these systems. This issue can be alleviated if a larger
fraction of material feeding the binary ends up accreting on to the primary
black hole. To illustrate this in the bottom panel we set $\beta=1$, implying
equal accretion rates onto each black hole. In this case, alignment timescales
reduce for the primary and increase for the secondary black hole, with all
times being less than the inspiral time down to $q=0.1$.

We stress that behaviour discussed in the above paragraph arises from the
analytic modelling of Gerosa et al. (2015). We additionally plot alignment
times taken directly from our sub-grid accretion disc model to illustrate
expected behaviour from our simulated systems. Spin-alignment timescales from
our simulations are given by (Fiacconi et al., 2018)

$\tau_{\rm
align}\approx\frac{0.17}{\cos\left(\pi/7\right)}\left(\frac{M_{\bullet}}{10^{6}M_{\odot}}\right)^{-2/35}\left(\frac{f_{\rm
Edd}}{\eta_{0.1}}\right)^{-32/35}a_{\bullet}^{5/7}{\rm Myr}\,.$ (12)

Median alignment timescales for primary and secondary black holes are shown by
the blue circles and red squares, respectively, while different line styles
indicate different binary eccentricities and inclination angles, as shown on
the legend. Vertical bars on each point illustrate the $5^{th}$ to $95^{th}$
percentile range. We only include snapshots after the first $100$ orbits to
limit any possible influence from the initial sub-grid parameter choices. The
green hatched regions indicates the range of inspiral times calculated from
the simulation data over the same time period. The upper bound of the green
hatched region follows the analytical expectation for inspiral times from
Equation (11), but in general, the simulations exhibit faster rates of binary
shrinking notwithstanding that our CBD is in a quasi steady state.

Several interesting features are clear from our simulations, firstly, while
the alignment timescales for secondary black holes are close to the region
predicted by the analytical model, the majority of primary black holes also
have alignment timescales in this region and do not follow the blue track
predicated by Gerosa et al. (2015). This behaviour largely stems from
differences in the black hole accretion rates in our simulated binaries and
those assumed by Gerosa et al. (2015), with the fact that our simulations do
not exhibit preferential accretion onto the secondary black hole playing a key
factor. We find that $\beta$ is always less than one and often close to zero,
and a larger fraction of material entering the cavity accretes onto the
primary black hole (and onto the binary in general) in low $q$ systems. These
conditions result in more efficient spin-alignment for primary black holes
even at low $q$ values. The only outliers from this general behaviour occur
for the retrograde binaries, which, apart from the $q=1/10$ primary black
hole, experience much less accretion than in other binary setups. As such
these points are almost entirely determined by the initial parameters of the
sub-grid accretion disc and their points could be seen as something of an
upper bound. Had these black holes experienced similar levels of accretion as
that observed in the other binaries then they would show much shorter
alignment timescales. Barring these outliers, spin alignment occurs on
timescales shorter than the inspiral timescale, often at rates $\lesssim
0.2\times\tau_{\rm in}$, meaning that for the parameter space studied here,
and for our assumed accretion disc model (see Section 4.5 for a discussion of
caveats) we would expect spins to be aligned by the time black holes reach
merger.

Figure 14: Polar projections of the evolution of the black hole spins (thick
solid lines initiated by circles) and sub-grid accretion disc angular momenta
(thin solid lines initiated by squares), with quantities associated with the
primary and secondary black holes are shown by blue and red, respectively. The
angular momenta of mini-disc material is shown by the blue and pink points
(around the primary and secondary black hole, respectively), with their
transparency level scaling with the gas inflow rate onto the corresponding
black hole. For misaligned runs, the binary and CBD angular momenta are
additionally shown by thick dashed green and purple lines initiated by
diamonds, respectively. In general, coplanar runs tend towards alignment
between black hole spins and the global angular momentum vector, however,
mini-discs in low $q$ misaligned systems do not align rapidly which inhibits
spin alignment in such systems.

There is, however, a caveat to the Gerosa et al. (2015) model and analysis we
perform here – namely that for this line of argument to work one has to assume
that the whole system is coplanar. If this is not the case, while black hole
spins may align with their accretion flow, there is no requirement for the
accretion flows and hence the black hole spins to be aligned with each other.
To explore this further, in Fig. 14 we show the angular momentum evolution of
different components of the system as polar projections. Quantities associated
with the primary and secondary black holes are shown by blue and red,
respectively, with spins shown by thick lines initiated by circles and
accretion disc angular momenta by thin lines initiated by squares. We also
show the angular momenta of accreted material, which we refer to as the mini-
disc angular momenta, although note that in some instances a well defined
mini-disc may not exist and instead this is simply representative of the gas
local to the black hole. We only show mini-disc angular momenta at times when
there is inflow onto the black hole as illustrated by faded blue and pink
points for the primary and secondary black hole, respectively. The points have
transparency levels that scale with the instantaneous inflow rate. The binary
and CBD angular momenta are shown for misaligned binary runs by thick dashed
lines initiated by diamonds. Each projection represents a different binary
setup as labelled.

First, for all in-plane binary cases the accretion disc angular momenta
generally remain aligned with the $z-$axis (i.e. at the pole), albeit with
some small deviation and possible precession. The black hole spins, which are
initially not aligned with this axis, progressively re-orientate towards the
$z-$axis. Concerning the mini-disc angular momentum, there are a few
behaviours to note. In the circular binary case, the mini-disc angular momenta
are on average aligned with the $z-$axis, although the distributions of
orientations in the case of the secondary black holes become broader for low
mass ratios. Larger deviations from $z-$axis alignment generally correspond to
times in which inflows are weaker, which fits with a picture in which the
lower mass mini-discs around secondary black holes in low $q$ systems are more
susceptible to disruption. Somewhat larger variations in the mini-disc angular
momenta are seen in the eccentric binaries, for both the primary and secondary
black holes, reflecting the dynamics of the cavity gas and episodic nature of
the mini-discs as the binary moves between apo- and pericentre. In the
retrograde binaries, there are clear inflows onto the primary black holes that
are on average aligned with the $z-$axis, with more scatter in the $q=1/3$
case compared to the $q=1/10$ case, reflecting that the latter exhibits a
long-lived coherent mini-disc. On the other hand, the $q=1/10$ binary shows no
inflows onto the secondary black hole, and while in the $q=1/3$ binary there
is some inflow, its angular momentum can vary over a wide range of angles,
with no coherent or long-lived mini-disc or accretion flow forming. Barring
the primary black hole in the $q=1/10$ case, the lack of coherent long-lived
inflows in the retrograde binaries results in limited evolution of the black
hole spin orientations.

In all of the coplanar cases, the nature of the systems means that the
assumptions used for the analysis presented in Fig. 13 are valid. Next, we
consider cases in which the binary and CBD are initially misaligned. In these
cases, we additionally show the binary and CBD angular momenta by the green
and purple dashed lines, which evolve towards aligning with each other. In the
$q=1$ case, the mini-disc angular momenta also tend to reorient towards a
common axis with little scatter, with the accretion disc angular momentum
approaching the same direction and the black hole spins also tending to align
by the end of the simulation. This suggests that in the $q=1$ case, the
analysis presented in Fig. 13 should still hold. However, for lower mass ratio
binaries, while the binary and CBD angular momentum align along a common axis,
the mini-disc angular momentum varies significantly (see also, Moody et al.,
2019). This leads to the accretion disc angular momenta also varying and while
for the secondary black hole, the accretion disc angular momenta and black
hole spin appear to move towards the common axis of the CBD/binary angular
momenta, the accretion disc around the primary exhibits a precession-like
motion and the black hole spin does not appear to align with the CBD/binary
angular momentum during the timescales explored in the simulation. In this
scenario, the assumptions used in the analysis of Fig. 13 potentially break
down and longer simulations would be required to understand the late time
evolution of the black hole spin orientation in low-$q$, misaligned systems.
We additionally note that for two such systems (q10e00i45 and q03e00i45mod),
the black hole spin and accretion disc are initially highly misaligned and in
a retrograde configuration. The extent of any initial misalignment of the
black hole spin could be an important factor in the alignment process that
requires further exploration.

##  4 Discussion

The simulation of binaries and their interaction with surrounding CBDs is a
rich area of research (see e.g., Lai & Muñoz, 2023, for a recent review),
however, the vast parameter space in terms of both the binary and CBD
properties make developing a comprehensive understanding of this topic a
formidable task. Many simulations of CBDs (e.g., D’Orazio et al., 2013;
Miranda et al., 2017; Moody et al., 2019; Muñoz et al., 2019; Muñoz et al.,
2020; D’Orazio & Duffell, 2021; Siwek et al., 2023b) make two key assumptions:
i) the gas can be treated as isothermal (or locally isothermal), meaning that
heating due to compression or shocks is assumed to rapidly radiate away, and
ii) the CBD self-gravity can be neglected. We instead consider the self-
gravitating regime, which has only been studied in a handful of previous works
(e.g., Cuadra et al., 2009; Roedig et al., 2011, 2012; Franchini et al.,
2021), and similarly employ an adiabatic equation of state with $\beta$
-cooling such that the CBD remains marginally Toomre stable. The CBD, which is
prevented from fragmenting, forms spiral structures that act to transport
angular momentum (e.g., Lodato & Rice, 2004, 2005; Franchini et al., 2021). We
have performed a parameter study with respect to the binary mass ratio, $q$,
eccentricity, $e$ and inclination angle, $i$, while fixing the initial total
binary mass ($M_{\rm bin}=2\times 10^{6}$ M☉) and properties of the CBD. In
this section, we discuss our results further and put them into the context of
the wider literature.

###  4.1 Evolution of prograde binaries

One of the key questions to address is how the interaction between the binary
and CBD alters the separation of the binary and whether or not CBDs can drive
binaries from pc to mpc scales on which GW emission can bring the binary to
coalescence. Specifically, both gravitational torques (Section 3.2.2) and
accretion onto the black holes (Section 3.3.1) can modify the binary angular
momentum. Gravitational torquing by a CBD has traditionally been invoked as a
mechanism to extract angular momentum from a binary and drive it towards
merger (e.g. Pringle, 1991; Artymowicz & Lubow, 1994, 1996; Ivanov et al.,
1999; Gould & Rix, 2000; Armitage & Natarajan, 2002, 2005; Lodato et al.,
2009). In particular, OLRs within prograde CBDs provide a mechanism through
which the binary and CBD can exchange angular momentum (Lynden-Bell & Kalnajs,
1972; Goldreich & Tremaine, 1979, 1980; Lin & Papaloizou, 1986). However, as
we have shown, gas within the cavity, including streams and mini-discs, also
torque the binary and provide material for accretion, which offers additional
channels for binary evolution (see also, MacFadyen & Milosavljević, 2008;
Cuadra et al., 2009; Shi et al., 2012; Roedig et al., 2011, 2012; D’Orazio et
al., 2013; Farris et al., 2014; Miranda et al., 2017; Muñoz et al., 2019;
Tiede et al., 2020; Franchini et al., 2022; Siwek et al., 2023b, for previous
examples of such processes).

Roedig et al. (2012) utilised high-resolution simulations to analyse the
gravitational torques acting on the binary. Similar to our results, they found
that gas within the cavity makes a telling contribution; gas at $R<a$ provides
a net positive torque and gas at larger radii provides net negative torques.
Similar results have also been found in multiple simulations that resolve
mini-discs (e.g. Muñoz et al., 2019; Tang et al., 2017; Tiede et al., 2020;
Franchini et al., 2022; Siwek et al., 2023b). Roedig et al. (2012) attribute
positive torque within $R<a$ to super-Keplerian inflows that bend in front of
the black holes, while at $R>a$, there is a combination of torques facilitated
by OLRs in the CBD and torques within the cavity that are purely kinematic in
origin. This picture is consistent with the torques identified in our prograde
binary systems, for example in Fig. 5 and Fig. 6.

Accretion onto the binary can additionally provide effective torquing due to
both changing the black hole masses and velocities (e.g., Roedig et al., 2012;
Franchini et al., 2021, 2022; Shi et al., 2012; Muñoz et al., 2019; Tiede et
al., 2020; Siwek et al., 2023b). The anisotropic component of the accretion
torque (arising from differences in velocity between the black hole and
accreting material), is generally found to be small compared to the
gravitational torque (Tang et al., 2017; Tiede et al., 2020; Duffell et al.,
2020; D’Orazio & Duffell, 2021), although not necessarily negligible (Muñoz et
al., 2019; Siwek et al., 2023b). Accretion torques presented in this paper
refer to the total change in binary angular momentum due to accretion, which
is always positive and in the majority of cases subdominant compared to the
gravitational torques given the low accretion rates. The overall evolution of
the semi-major axis is given as (Roedig et al., 2012)

$\frac{\dot{a}}{a}=\frac{2\tau_{\rm bin}}{J_{\rm
bin}}+\frac{2e}{1-e^{2}}\dot{e}-\frac{\dot{M}_{\rm bin}}{M_{\rm
bin}}-\frac{2\dot{\mu}}{\mu}\,,$ (13)

which is found by differentiating Equation (4), rearranging and setting
$\dot{J}_{\rm bin}=\tau_{\rm bin}$. Neglecting the eccentricity term and
rearranging Equation (13), we can define the maximum torque under which an
equal mass, circular binary shrinks,

$\tau_{\rm bin}<\frac{3}{8}\dot{M}_{\rm bin}\sqrt{GM_{\rm bin}a},\\\ \simeq
1.12\times\left(\frac{\dot{M}_{\rm bin}}{2\times 10^{-4}~{}{\rm
M_{\odot}yr}^{-1}}\right)10^{-6}GM_{\rm bin}^{2}a^{-1}\,,$ (14)

where we use $M_{1}=M_{2}=M_{\rm bin}/2$ and assume
$\dot{M}_{1}=\dot{M}_{2}=\dot{M}_{\rm bin}/2$. As such, binaries experiencing
a positive torque can still shrink provided $\dot{M}_{\rm bin}$ is
sufficiently large (see e.g., Roedig et al., 2012). From our post-processed
direct summation torque estimates, the q01e00i45 binary should experience a
positive torque (see Section 3.2.2 and Fig. 7), that would increase $J_{\rm
bin}$. However, given the binary’s average accretion rate of $\sim 2\times
10^{-4}$ M⊙yr-1, Equation (14) suggests that the limiting torque is roughly an
order of magnitude greater than the average torque experienced by the
q01e00i45 binary. Therefore we expect this binary to shrink.

In all bar one case (the q01e00i45 binary) we expect the average total torque
(i.e. gravitational plus accretion) experienced by binaries in our simulations
to be negative, and that every binary should shrink. Previous works simulating
high-mass, self-gravitating CBDs have similarly found binary shrinking (e.g.,
Cuadra et al., 2009; Roedig et al., 2012; Roedig & Sesana, 2014; Franchini et
al., 2021). However, this behaviour is not universal. Recent works, which
typically neglect self-gravity and assume a (locally-)isothermal equation of
state, have found that net torques acting on a binary are sufficiently
positive to drive expansion (e.g., Miranda et al., 2017; Tang et al., 2017;
Moody et al., 2019; Muñoz et al., 2019; Muñoz et al., 2020; Duffell et al.,
2020; Heath & Nixon, 2020; Tiede et al., 2020; D’Orazio & Duffell, 2021; Siwek
et al., 2023b; Wang et al., 2023). As highlighted in the review of Lai & Muñoz
(2023), circular binaries above a mass ratio of $q\sim 0.1-0.3$ expand, while
for mass ratios below $q\sim 0.05-0.1$ they shrink (e.g., Duffell et al.,
2020; Muñoz et al., 2020; Siwek et al., 2023b). Eccentric binaries have also
been observed to expand for particular combinations of $e$ and $q$ (e.g. Muñoz
et al., 2019; D’Orazio & Duffell, 2021; Siwek et al., 2023b), although
D’Orazio & Duffell (2021) notes that initially eccentric binaries ($e\gtrsim
0.1$) tend towards an equilibrium eccentricity ($e\simeq 0.4$) where the
binaries are likely to shrink. Clearly, to elucidate this issue in future
works it will be of crucial importance to explore both the roles of gas self-
gravity and realistic thermodynamic modelling both within CBDs and critically
the mini-discs.

Simulations of massive CBDs have traditionally employed Lagrangian simulations
leading to restricted resolution in the cavity. We address this issue using a
super-Lagrangian refinement scheme. Similarly, Franchini et al. (2022)
compared traditional “Lagrangian” simulations of a massive ($M_{\rm
cbd}/M_{\rm bin}=0.1$) CBD to simulations employing “hyper-Lagrangian”
refinement. They find that while positive torques from gas at $R\lesssim a$
can lead to mild binary expansion, this is sensitive to choices of disc
temperature and thickness. Like many other simulations that find binary
expansion, Franchini et al. (2022) assumed an isothermal equation of state and
did not account for self-gravity. Varying the level of $\beta-$cooling, from
$\beta=0$ (locally-isothermal) to $\beta=4$, Wang et al. (2023) found binary
expansion for the majority of values explored except for
$0.7\lesssim\beta\lesssim 0.9$. The simulations presented in our work, as well
as others that regularly find binary shrinking (e.g. Cuadra et al., 2009;
Roedig et al., 2012; Roedig & Sesana, 2014; Franchini et al., 2021) instead
employ higher values of $\beta\simeq 10$. More massive CBDs, such as those
considered by these simulations, potentially drive more rapid initial
shrinking of the binary due to stronger gravitational interactions (Franchini
et al., 2021). We also note that Roedig et al. (2012) found that using an
isothermal equation of state in the cavity leads to significantly weaker
gravitational torques compared to the adiabatic case, leading to slower
evolution of the binary semi-major axis.

The viscous torque, which is proportional to $\alpha\left(H/R\right)^{2}$,
determines how rapidly material is transported to OLR locations, with the
expectation that thicker and more viscous discs promote resonant torquing and
hence more effective shrinking of the binary (see e.g. discussions in Heath &
Nixon, 2020; Franchini et al., 2021). Despite our simulations having low
intrinsic viscosity and equilibrium thickness ($H/R\sim 0.05$), we find that
binaries shrink. Other recent works have found that thicker discs and/or
higher viscosities may lead to binary expansion. For example, Heath & Nixon
(2020) found that for isothermal CBDs, with fixed $\alpha-$viscosity, thick
discs led to binary expansion due to accretion of high angular momentum
material (see also, Franchini et al., 2022). While by varying $H/R$ at fixed
kinematic viscosity, Tiede et al. (2020) found expansion of binaries in thick
discs and shrinking in thin discs driven by gravitational torques. Other works
have instead varied $\alpha$ for a fixed $H/R$, finding only negligible
differences in the accretion rate normalised torques (Muñoz et al., 2020;
Duffell et al., 2020), although Franchini et al. (2022) found that the level
of viscosity determined whether or not their binary entered a phase of
expansion.

###  4.2 Retrograde and misaligned binaries

While many simulations model coplanar prograde systems, the relative alignment
between a CBD and binary could be arbitrary. Black hole growth via chaotic
accretion would occur via many small randomly orientated accretion events
(King & Pringle, 2006, 2007), during which a binary may interact with many
misaligned accretion events and CBDs before eventually merging (Nixon et al.,
2013). For misaligned binaries ($i=45\degr$), we find, as expected from
analytical theory (Nixon et al., 2011b), that the binary and CBD torque each
other and eventually align via the outward propagation of a warp (see Fig. 8).
The radial gravitational torque profiles (Fig. 10) are qualitatively similar
to those of in-plane counterparts, although the $q=1$ and $q=1/3$ misaligned
binaries experience overall weaker negative gravitational torques, which is a
consequence of stronger positive net torques from gas at $R<a$ (see Table 2),
and in general experience stronger positive accretion torques. Simulations of
a misaligned $q=1$ binary by Moody et al. (2019) exhibit similar general
behaviours to ours in terms of the warping and binary-CBD alignment, although
for both misaligned and in-plane systems they find binary separation growth.
Similar to our work, they find that the mini-disc orientations are variable
due to CBD alignment occurring on timescales shorter than mini-disc alignment.
As we discuss in Sections 3.3.2 and 4.4, this has potentially important
consequences for spin alignment, remnant recoil velocities and GW detection
rates.

While in our simulations, and those of Moody et al. (2019), the disc warps,
under certain conditions the disc may instead break or tear (Nixon et al.,
2013; Aly et al., 2015). Nixon et al. (2013) found that discs can break at a
wide range of inclination angles. They provide analytical limits for the
radius within which disc breaking is expected for both viscous and inviscid
CBDs. Based on their viscous limit we might naively expect to see disc
breaking in our simulations. However, there are several differences between
our work and that of Nixon et al. (2013), such as the gas equation of state
(isothermal vs adiabatic with $\beta-$cooling), the inclusion of self-gravity
that induces stabilising spiral arms, and the treatment of viscosity. The
inviscid regime whereby the warp propagates as a wave (Papaloizou & Pringle,
1983), and may be more applicable to our low viscosity simulations, suggests a
breaking radius within the cavity and hence we would not expect disc breaking.
While all CBDs in our simulations end up aligning with the binary, Nixon et
al. (2011b) show that counter-alignment is also possible if $\cos i<-J_{\rm
cbd}/2J_{\rm bin}$. Aly et al. (2015) further found that alignment to stable
polar orbits is possible around highly eccentric binaries. Both Nixon et al.
(2013) and Aly et al. (2015) find that disc tearing can lead to angular
momentum cancellation and accretion of low angular momentum gas by the binary
that could aid binary shrinking.

Dunhill et al. (2014) studied the formation of CBDs from infalling turbulent
gas clouds on both prograde and retrograde orbits, finding that misaligned
discs can align or counter-align with the binary, with the retrograde CBDs
driving more significant binary evolution. Retrograde CBDs, therefore, provide
a promising avenue to drive binaries towards coalescence. Resonant torques,
facilitated by OLRs in prograde CBDs, do not occur in retrograde systems
(Nixon et al., 2011a). Therefore, while prograde CBDs are expected to move
outwards, retrograde CBDs are not, with gas more readily filling the cavities
(e.g., Nixon et al., 2011a; Roedig & Sesana, 2014; Bankert et al., 2015). This
behaviour is seen in our simulations (e.g., in Figs. 2, 3 and 5), and while
shrinking of $a$ appears to slow for prograde binaries, it continues to
decrease quite rapidly for retrograde cases (see Fig. 4). We find that the
evolution of retrograde binaries is dominated by gravitational torques
originating from gas at $R<a$ (see Table 2), which, unlike the prograde case,
provide net negative torques with respect to the binary angular momentum.
Similar results were found by Tiede & D’Orazio (2023), who additionally found
that while initially, circular binaries remain circular (as we also find, see
Fig. 4), eccentric binaries become ever more eccentric. Tiede & D’Orazio
(2023) additionally find that for $e\gtrsim 0.55$, retrograde Lindblad
resonances can be driven in the CBD (see also, Nixon & Lubow, 2015). An
alternative mechanism to shrinking retrograde binaries via gravitational
torques was shown by Nixon et al. (2011a) who instead find that accretion can
drive binaries to coalescence. If accretion is efficient at apocentre but not
pericentre, eccentricity growth can occur even for an initially circular
binary. Negligible growth of the secondary black holes in our simulated
retrograde binaries is likely due to a combination of differences in the gas
thermodynamic modelling, the lack of a robust mini-disc around these black
holes, and the requirement for net inflow in order to accrete gas. With
respect to the latter, a less restrictive accretion criteria may lead to
eccentricity growth similar to that seen in Nixon et al. (2011a).

###  4.3 Evolution of black hole properties

The varied population of supermassive black holes that reside in the centres
of galaxies (Ferrarese & Merritt, 2000; Häring & Rix, 2004; Gültekin et al.,
2009; Kormendy & Ho, 2013; McConnell & Ma, 2013; Sahu et al., 2019), including
dwarfs (Reines et al., 2013; Reines et al., 2020; Mezcua et al., 2016; Kaviraj
et al., 2019; Birchall et al., 2020), gain their mass through a combination of
accretion and merger events over cosmic time. Numerous simulations have been
performed to study accretion onto binaries (e.g., Shi et al., 2012; MacFadyen
& Milosavljević, 2008; D’Orazio et al., 2013; Farris et al., 2014; Miranda et
al., 2017; Duffell et al., 2020; Muñoz et al., 2020; Wang et al., 2023; Siwek
et al., 2023a). As highlighted by Roedig et al. (2012), simulations that
excise the binary region commonly assume a limiting case whereby all material
inflowing at $R=a$ is accreted, however, they find a reduction in the
accretion rate of $\sim 25\%$ as material is stirred and slingshotted back
towards the CBD. While total black hole accretion rates in our simulations are
roughly equal to the mass flow through the cavity, many exceed it slightly,
which is likely a result of the second relaxation phase we employ, during
which accretion is prohibited and the well-defined mini-discs form.

We show inflow rates onto black holes in Fig. 11 and similarly to previous
works we find that they vary on times scales of $\sim t_{\rm bin}$ (e.g.,
MacFadyen & Milosavljević, 2008; Cuadra et al., 2009; Roedig et al., 2012;
Muñoz & Lai, 2016; Lai & Muñoz, 2023; Wang et al., 2023). While these
fluctuations, driven by the binary orbit and induced accretion streams, are
common among simulations, some additionally show long-term variability on a
timescale of $\sim 5t_{\rm bin}$ associated with the formation of a lump at
the inner edge of the CBD (Shi et al., 2012; MacFadyen & Milosavljević, 2008;
D’Orazio et al., 2013; Muñoz & Lai, 2016; Miranda et al., 2017; Wang et al.,
2023; Lai & Muñoz, 2023), although this is only expected when gas can cool
rapidly Wang et al. (2023).