Parsecscale separation supermassive black hole binaries in the centre of gasrich galaxy merger remnants could be surrounded by massive circumbinary discs CBDs. Black hole mass and spin evolution during the gasrich binary inspiral are crucial in determining the direction and power of relativistic jets that radio observations with LOFAR LowFrequency Array and SKAO Square Kilometer Array Observatory will probe and for predicting gravitational wave GW emission that the IPTA International Pulsar Timing Array and LISA Laser Interferometer Space Antenna will measure. We present 3D hydrodynamic simulations capturing gasrich selfgravitating CBDs around a inlineformulatexmath idTM0001 notationLaTeX2times 106texmathinlineformula Minlineformulatexmath idTM0002 notationLaTeXodot texmathinlineformula supermassive black hole binary that probe different mass ratios eccentricities and inclinations. We employ a subgrid ShakuraSunyaev accretion disc to selfconsistently model black hole mass and spin evolution together with superLagrangian refinement techniques to resolve gas flows streams and minidiscs 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 lowmass ratios. When considering initially misaligned systems the orientation of the minidiscs around each black hole can vary significantly. We discuss the implications of this behaviour for black hole spin alignment and highlight the need for broader parameter space studies of misaligned systems to understand the impact on black hole recoil velocities.
['Astrophysics - High Energy Astrophysical Phenomena', 'Astrophysics - Astrophysics of Galaxies']
Dynamics and spin alignment in massive gravitoturbulent circumbinary discs around supermassive black hole binaries
{'Dynamics and spin alignment in massive, gravito-turbulent circumbinary discs around supermassive black hole binaries': "Martin A. Bourne 1 , 2 , 3 ★ , Davide Fiacconi 1 , 2 , Debora Sijacki 1 , 2 , Joanna M. Piotrowska 2 , 4 and Sophie Koudmani 1 , 2 , 5 , 6  \n- 1 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK\n- 2 Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK\n- 3 Centre for Astrophysics Research, Department of Physics, Astronomy and Mathematics, University of Hertfordshire, College Lane, Hatfield AL10 9AB, UK\n- 4 Cahill Center for Astronomy and Astrophysics, California Institute of Technology, Pasadena, CA, USA\n- 5 Center for Computational Astrophysics, Flatiron Institute, 162 5 𝑡ℎ Avenue, New York, NY 10010, USA\n- 6 St Catharine's College, University of Cambridge, Trumpington Street, Cambridge CB2 1RL, UK  \n13 September 2024", '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 × 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. When considering initially misaligned systems, the orientation of the mini-discs around each black hole can vary significantly. We discuss the implications of this behaviour for black hole spin alignment and highlight the need for broader parameter space studies of misaligned systems to understand the impact on black hole recoil velocities.  \nKey words: black hole physics - accretion, accretion discs - quasars: supermassive black holes - methods: numerical - black hole mergers - gravitational waves', '1 INTRODUCTION': "Within the standard Λ -cold dark matter ( Λ 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. 2024), and are a crucial mechanism in galaxy evolution. Not only does this lead to the mass assembly of dark matter halos and  \ngalaxies 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).  \nWhile the exact nature of gas delivery to supermassive black holes in the aftermath of a galaxy merger is not yet fully understood, galaxies hosting supermassive black holes (see e.g., Ferrarese & Merritt 2000; Haring & Rix 2004; Gultekin et al. 2009; Kormendy & Ho 2013; McConnell & Ma 2013; Sahu et al. 2019) are expected to undergo multiple major mergers with similarly sized galaxies over their cosmic lifetime (see e.g., O'Leary et al. 2021; Huˇsko 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  \nthe merger remnant system and form a ∼ 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'c & Merritt 2001, 2003; 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, 2013; Roedig et al. 2012; D'Orazio et al. 2013; Farris et al. 2014; Moody et al. 2019; Mu˜noz et al. 2019, 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 ≳ kpc scales (see Fig. 1 in Chen et al. 2022). On scales of ≲ 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). We refer to the review of De Rosa et al. (2019) for further discussion of dual AGN in the multimessenger era. 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.  \nThe 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; Agazie et al. 2024); 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 ∼ 0 . 01 -100 mHz frequency GWs emitted by individual supermassive black hole merger events over a wide range of masses (peaking between ∼ 10 4 -10 7 M ⊙ ) and out to very high redshifts (AmaroSeoane 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'alez 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.,  \nRajagopal & Romani 1995; Phinney 2001; Enoki et al. 2004; Sesana et al. 2004, 2009b; Kelley et al. 2017; Siwek et al. 2020).  \nAvery 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 both via direct extraction of spin energy (Blandford & Znajek 1977; Tchekhovskoy et al. 2011; Talbot et al. 2021) and the impact on inflowing gas dynamics (Talbot et al. 2022, 2024). 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; Barnes 2002; 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; Capelo et al. 2015; Capelo & Dotti 2017). Such a gasrich 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) or the excavation of a cavity within a larger nuclear gas disc by the binary (Souza Lima et al. 2020).  \nIn 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'c 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˜noz et al. 2019; Tiede et al. 2020; Dittmann & Ryan 2022; Franchini et al. 2022; Siwek et al. 2023b). In particular, simulations that resolve the mini-discs and gas flows within the binary orbit ( 𝑅 < 𝑎 ) 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˜noz 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˜noz et al. 2019, 2020; Duffell et al. 2020; Heath & Nixon 2020; Tiede et al. 2020; D'Orazio & Duffell 2021; Dittmann & Ryan 2022, 2024; 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 selfgravitating, 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 𝛽 -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 2024) 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.  \nThe main aim of this paper is to build on existing literature and modelforthe first time both the mass and spin evolution of black holes  \nin 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.1 Initial conditions': 'The binary and CBD setup that we model is generally scale-free meaning that it is relevant to a wide range of scenarios and the majority of results we present can be re-scaled with respect to the system mass and length. However, for this work we employ the subgrid accretion disc model of Fiacconi et al. (2018) to track the black hole spin evolution, which requires us to define a physical mass and length scale for our simulations. Therefore, we consider supermassive black hole binaries with initial total mass 𝑀 bin = 𝑀 · , 1 + 𝑀 · , 2 = 2 × 10 6 M ⊙ and mass ratio 0 . 1 ≤ 𝑞 ≡ 𝑀 · , 2 / 𝑀 · , 1 ≤ 1. We chose these parameters to be in the binary parameter space that would be of the largest interest for the Laser Interferometer Space Antenna (LISA; Amaro-Seoaneetal.2013).Thebinaryisinitially on a Keplerian orbit with eccentricity of either 𝑒 = 0 or 𝑒 = 0 . 5 and initial semi-major axis 𝑎 = 2 pc, corresponding to a period of 𝑡 bin = 0 . 187 Myr. The peri- and apo-centre of the binary are aligned along the 𝑥 axis.  \nA gaseous CBD surrounds the binary in the 𝑥𝑦 plane and initially extends between 𝑅 in = 2 𝑎 = 4 pc and 𝑅 out = 7 𝑎 = 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 𝑦 axis with an inclination angle 𝑖 ; for reference, 𝑖 = 0 · and 𝑖 = 180 · correspond to prograde and retrograde coplanar orbits, respectively. The disc follows the surface density profile  \nΣ ( 𝑅 ) = Σ GLYPH<18> 𝑅 𝑅 in GLYPH<19> -𝛼 , (1)  \nwhere 𝛼 = 2. The normalisation Σ is chosen to enforce 𝑀 cbd / 𝑀 bin = 0 . 1, where 𝑀 cbd is the total mass of the CBD.  \nThe vertical structure of the CBD is constructed by assuming that the aspect ratio 𝐻 / 𝑅 ≈ 𝑐 s ( 𝑅 )/ 𝑉 K ( 𝑅 ) is constant, where 𝑐 s ( 𝑅 ) and 𝑉 K ( 𝑅 ) are the local isothermal sound speed and the Keplerian velocity for 𝑀 bin , respectively. Therefore, the temperature profile of the CBD, 𝑐 2 s ( 𝑅 ) = ( 𝐻 / 𝑅 ) 2 𝑉 2 K ( 𝑅 ) ∝ 𝑅 -1 , can be normalised after choosing the value of the aspect ratio. Since we include self-gravity  \nin the CBD, we use the Toomre parameter 𝑄 to impose the normalisation. Specifically, we choose 𝑄 ≈ 3 ( 𝐻 / 𝑅 )( 𝑀 bin / 𝑀 cbd ) = 1 . 5, which corresponds to an initial aspect ratio 𝐻 / 𝑅 ≈ 0 . 5 𝑀 cbd / 𝑀 bin . We do not include any additional explicit source of viscosity beyond the intrinsic (low) numerical viscosity, instead relying 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 𝐻 ( 𝑅 ) . The CBD is initially in nearly Keplerian rotation with azimuthal velocity  \n𝑉 2 𝜙 ( 𝑅 ) = 𝐺𝑀 bin 𝑅 " 1 + 3 4 𝑞 ( 1 + 𝑞 ) 2 GLYPH<18> 𝑅 𝑎 GLYPH<19> -2 # -3 𝑐 2 s ( 𝑅 ) , (2)  \nwhere 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 (see e.g., equation 2 in Mu˜noz et al. 2020). We do not include the contribution of the CBD potential in this equation, noting that despite this the CBD remains stable throughout the simulation where the CBD self-gravity is self-consistently taken into account. Note that we do not add any additional velocity dispersion to the gas. The black hole binary and CBD are the only components simulated in this work. While contributions from other components (e.g., stars and dark matter) could be important in dense galactic nucleii for lowmass/large separation binaries, we neglect them here to isolate the effects of the CBD-binary interaction and to ensure that the system remains scale-free.  \nTable 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 ( 𝑒 = 0), eccentric ( 𝑒 = 0 . 5), coplanar prograde ( 𝑖 = 0 · ), misaligned ( 𝑖 = 45 · ), and retrograde ( 𝑖 = 180 · ) binaries. For example, the 1:1 mass ratio, circular, coplanar 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 · 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 𝑦 -axis in order to be inclined with respect to the CBD by 45 · and track its subsequent evolution.', '2.2 Numerical code and simulation setup': "Thesimulations presented in Table 1 are performed using the movingmesh hydrodynamical code /a.pc/r.pc/e.pc/p.pc/o.pc (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 𝑚 target g per cell. The initial conditions are sampled with 10 6 initial mesh elements in the CBD corresponding to a target mass 𝑚 target g = 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 𝑁 max direct = 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 neigh-  \nTable 1. Summary of the characteristics of the performed runs. Note that: (i) 𝑀 bin = 2 × 10 6 M ⊙ ; (ii) 𝑞 ≡ 𝑀 · , 2 / 𝑀 · , 1 ≤ 1. From left to right: simulation label, initial mass ratio 𝑞 , initial eccentricity 𝑒 , initial inclination angle 𝑖 between the binary orbital plane and the CBD midplane, initial Eddington ratio 𝑓 Edd , 1 for 𝑀 · , 1 , initial Eddington ratio 𝑓 Edd , 2 for 𝑀 · , 2 , initial 𝛼 accretion disc and black hole angular momentum ratio for 𝑀 · , 1 , initial 𝛼 accretion disc and black hole angular momentum ratio for 𝑀 · , 2 , misalignment angle between 𝑱 · and 𝑱 d for 𝑀 · , 1 , misalignment angle between 𝑱 · and 𝑱 d for 𝑀 · , 2 .  \nh 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 𝜖 · = 0 . 01 pc, while gas cell softenings adapt with the cell size down to a minimum of 𝜖 g = 0 . 005 pc.  \nWe evolve the system assuming an equation of state for an ideal fluid, 𝑃 = ( 𝛾 -1 ) 𝜌𝑢 , where 𝜌 is the gas density, 𝑢 is the specific internal energy, and 𝛾 = 5 / 3 is the adiabatic index, appropriate for a monoatomic ideal gas. We also allow the gas to cool radiatively through a simple 𝛽 -cooling prescription, with the rate of cooling given by  \nGLYPH<18> d 𝑢 d 𝑡 GLYPH<19> cool = -𝑢 𝑡 cool ( 𝑅 ) , (3)  \nwhere the local cooling timescale 𝑡 cool ( 𝑅 ) is proportional to the local orbital time Ω -1 𝐾 ( 𝑅 ) , namely 𝑡 cool ( 𝑅 ) = 𝛽 ( 𝐺𝑀 bin 𝑅 -3 ) -1 / 2 , where, as stated above, 𝑀 bin is the total binary mass and 𝑅 is the distance to the binary centre of mass. We discuss the implications for the chosen cooling prescription and its implementation further in Section 4. We set 𝛽 = 10; this choice allows the gas to cool slowly without fragmenting 1 (Gammie2001; Rice et al. 2005; Cuadra et al. 2009; Franchini et al. 2021), and eventually makes the CBD settle into a marginally stable configuration with 𝑄 ∼ 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).  \nThe initial conditions are evolved for 50 binary orbits corresponding to about 1 . 5 𝑡 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, 𝐻 / 𝑅 , can reach ≈ 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': "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 /a.pc/r.pc/e.pc/p.pc/o.pc 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 𝑅 ref around each black hole to force the maximum cell size of the mesh elements to decrease linearly from 𝑅 cell max at 𝑅 ref to 𝑅 cell min at the centre. Cell radii are allowed to remain within a factor 2 of the local maximum. We set 𝑅 cell max = 0 . 2 pc and 𝑅 cell min = 0 . 01 pc; 𝑅 cell max is chosen to match the typical size of the mesh cell just inside the CBD cavity, while 𝑅 cell min is equal to the gravitational softening of the black holes. We have chosen an 𝑒 -dependent value for 𝑅 ref , namely 𝑅 ref = 3 pc and 𝑅 ref = 5 pc for 𝑒 = 0 and 𝑒 = 0 . 5, respectively; moreover, we set 𝑅 ref = 2 pc for runs with 𝑖 = 180 · 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 𝑚 min = 10 -5 𝑚 target g for cells within the refinement region.  \nFig. 1 shows the effect of the super-Lagrangian refinement on the evolution of the cavity for run q01e00i00 as a representative example. Wecompare 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 re-  \nFigure 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 superLagrangian refinement is activated.  \n<!-- image -->  \nfinement. 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˜noz 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 ∼ 50 𝑡 bin , the binary in the run with refinement  \ncontinues to shrink after this time, though we do not see appreciable differences in the evolution of the eccentricity. The lack of difference suggests that the eccentricity evolution is primarily due to the binary being torqued by the CBD, which remains unchanged between runs with and without refinement. This is consistent with the findings of Roedig et al. (2011) and Roedig et al. (2012), who claim that eccentricity evolution arises due to gravitational interactions between over-densities at the inner edge of the CBD and the binary. After the 𝛽 -cooling relaxation, all simulated systems are furthermore relaxed with the super-Lagrangian refinement for 50 additional binary orbits. Thus, the overall relaxation phase lasts ≈ 18 . 7 Myr, corresponding to ≈ 3 𝑡 cool ( 𝑅 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 𝑡 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 𝛼 -disc solution 2 (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.  \nThe 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 𝑎 · , 1 , 2 = 0 . 9, consistent with current observational constraints for ∼ 10 6 M ⊙ black holes (Reynolds 2019). The timescale, 𝜏 align , on which black hole spin aligns with respect to the accretion disc angular momentum scales as 𝑎 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 𝜏 align . The accretion disc is initiated with a mass of 𝑀 d = 10 -3 𝑀 · , 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  \n2 In the following, we will refer to the large-scale, hydrodynamicallyresolved 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 𝛼 -discs as 'accretion discs'.  \nsimulations. 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, 𝐽 d / 𝐽 · , is set to ensure consistency between the initial 𝑓 Edd and other black hole and accretion disc properties (see Equation (8)).  \nFinally, we note an important modification we made to the original accretion model. In Fiacconi et al. (2018), the mass inflow rate / 𝑀 inflow and the specific angular momentum of the inflowing material 𝑳 inflow were calculated with a kernel-weighted interpolation from a volume encompassing a mass of 32 𝑚 target g , typically slightly larger than the size of the mini-discs. While we employ the same approach to estimate / 𝑀 inflow , we have adjusted how we measure 𝑳 inflow as despite the central weighting this approach tends to overestimate its magnitude compared to that of the gas resolved on scales of the sub-grid accretion disc (typically by a factor of a few). Instead, we only take the direction of 𝑳 inflow and we assign the magnitude of the inflowing gas specific angular momentum to 𝐿 inflow = 𝐿 d , where 𝐿 d is the instantaneous specific angular momentum of the accretion disc. This means that the condition originally imposed in Fiacconi et al. (2018) requiring 𝐿 inflow ≤ 𝐿 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 𝐿 d and the specific angular momentum of gas on scales of the sub-grid accretion disc.", '3.1 Properties of the circumbinary discs': "After the relaxation phases, we run the system for an additional 500 𝑡 bin during which the CBD outside the cavity evolves in a nearly steady state. Fig. 2 compares density projections for all runs at time 𝑡 = 100 𝑡 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 𝑖 = 45 · , where the CBDbecomes 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 𝑖 = 0 · . 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 𝛼 -parameter 𝛼 SG = ( 4 / 9 ) [ 𝛾 ( 𝛾 -1 ) 𝛽 ] -1 ≈ 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 / 𝑀 = 3 𝜋𝜈 Σ ≈ 5 × 10 -3 M ⊙ yr -1 . However, the complex nature of the system, driven by binary torques and gas hydrodynamics in full 3D geometry, means that the mass flows through the disc vary with time and while instantaneous mass inflow rates similar to  \nthis estimate are achieved, gas can also move outwards leading to time-averaged net inflow rates ≳ 10 -4 M ⊙ yr -1 measured directly from the simulations, with fluctuations of factor ∼ a few between different runs (see Fig. 11). This suggests that at least qualitatively and on average over time, the CBD behaves like an 𝛼 -disc, albeit with a lower net inflow rate than predicted from idealised analytic estimates.  \nIt 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˜noz et al. 2019; Tiede et al. 2020; Franchini et al. 2022, noting, however, that not all of these simulations capture the mini-discs themselves.).  \nExcept for the 𝑞 = 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 𝑞 , with a similar trend being seen in Fig. 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.  \nOn the other hand, the retrograde ( 𝑖 = 180 · ) 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, in the retrograde CBDs the net gravitational torques are negative and due to a lack of resonances are significantly weaker, which promotes continued inflow of material and a smaller cavity (Nixon et al. 2011a; Roedig & Sesana 2014; Bankert et al. 2015).  \nTypically, the sharp inner edge at 2 𝑎 is lost rapidly during relaxation. For the eccentric cases in particular the inner edges of the CBDsaredisturbed 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.  \nThe 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  \nFigure 2. Gas surface density maps projected on the 𝑥𝑦 plane along the 𝑧 axis of all the simulations in Table 1 at time 𝑡 = 100 𝑡 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 𝑥𝑦 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.  \n<!-- image -->  \nFigure 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 ∝ 𝑅 -2 profile; the middle panels show the aspect ratio 𝐻 / 𝑅 estimated as 𝑐 s / 𝑉 𝜙 ; and the bottom panels show the Toomre parameter 𝑄 . Light-grey and dark-grey shaded areas highlight the marginal stability and instability regions for 𝑄 , 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 2 𝑎 = 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.  \n<!-- image -->  \nsizes 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.  \nBarring 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 𝑞 = 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 𝑞 = 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 𝑞 = 1 / 10 simulations, where more material diffuses in from the CBD inner edge. While the whole cavity resembles a Type II gap between the minidisc around the primary and the CBD, as observed in protoplanetary disc calculations (Rafikov 2002; Nelson & Papaloizou 2003; Crida  \n& Morbidelli 2007), we note that the processes controlling gap formation differ for prograde and retrograde systems (see e.g., Ivanov et al. 2015).  \nFig. 3 shows the time evolution of the CBD surface density, aspect ratio, and Toomre parameter 𝑄 profiles for our different simulations. Despite fluctuations due to the formation of spiral arms, the CBD roughly maintains the original ∝ 𝑅 -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 𝑞 ≥ 1 / 3 and 𝑖 = 0 · tends to decrease over time, particularly for circular binaries. Similar but weaker behaviour is seen for 𝑞 = 1 / 10 binaries (i.e. run q10e00i00 and q10e05i00); however, the cavity surface density is typically larger, therefore confirming the tendency of the 𝑞 = 1 / 10 cases to have more diffuse material inside the cavity.  \nTheevolution 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  \nFigure 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 coplanar and co-rotating misaligned binaries; comparison between coplanar and counter-rotating binaries. Blue, purple and yellow curves refer to mass ratios 𝑞 = 1, 𝑞 = 1 / 3, and 𝑞 = 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.  \n<!-- image -->  \nslowly 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.  \nFor systems in which the binary and CBD are initially coplanar ( 𝑖 = 0 · and 𝑖 = 180 · ), the aspect ratio and Toomre parameter are generally constant within the CBD for 𝑅 ≲ 6 -7 𝑎 . For these systems, the aspect ratio typically ranges between ≈ 0 . 03 and ≈ 0 . 06 and the Toomre parameter is ≈ 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 ( 𝑖 = 45 · ), 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 ≳ 0 . 1 because of the presence of the mini-disc and the Toomre parameter 𝑄 ≫ 1. The mini-discs are indeed generally less dense and hence thicker than the CBD , however, the exact values of these quantities within about 1 . 5 𝑎 should be taken with a pinch of salt due to the mini-discs and streams resulting in the geometry of the system departing significantly from axisymmetry as assumed by the profiles. On the other hand, outside ∼ 7 𝑎 , 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 Σ , 𝐻 / 𝑅 and 𝑄 profiles.", '3.2.1 General behaviour': "While some studies of CBDs include analytic binaries, with the binary orbit pre-determined and fixed (e.g., Mu˜noz et al. 2019, 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 coplanar and co-rotating misaligned binaries (middle), and coplanar 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:  \n𝐽 bin = 𝜇 √︃ 𝐺𝑀𝑎 GLYPH<0> 1 -𝑒 2 GLYPH<1> , (4)  \nwhere 𝜇 = 𝑀 1 𝑀 2 /( 𝑀 1 + 𝑀 2 ) and 𝑀 = 𝑀 1 + 𝑀 2 are the binary reduced and total mass, respectively 3 , decreases over time (by ∼ 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.  \nWe do not observe any obvious trend in the evolution of 𝑎 as a function of the mass ratio 𝑞 ; however, we notice that for eccentric, inclined and retrograde 𝑞 = 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  \nis 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 𝑅 < 𝑎 (typically positive), and gas at 𝑅 > 𝑎 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 𝑞 = 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).  \nSome 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 2024, for retrograde eccentricity evolution). The eccentricity generally tends to grow slightly faster for the 𝑞 = 1 / 3 cases than for equal-mass binaries, regardless of the other parameters. This tendency with the binary mass ratio is inverted for the 𝑞 = 1 / 10 cases, where the eccentricity growth is the slowest, except for the binaries starting with high initial eccentricity 𝑒 = 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 𝑡 = 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.  \nA 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 (2024). However, Nixon et al. (2011a) discuss that the eccentricity growth depends on accretion onto (or gas captured by) 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. Unlike in Nixon et al. (2011a), the secondary black holes in our retrograde simulations accrete a negligible amount of mass (see Section 3.3.1), as such, a similar channel of eccentricity growth due to accretion does not occur 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  \nevolution (e.g., Cuadra et al. 2009; Roedig et al. 2012; Mu˜noz 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.  \nFig. 5 shows time-averaged maps of the gas surface density and torque surface density for every coplanar run. Panels on the lefthand 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 ∼ 100 snapshots between 𝑡 = 200 𝑡 bin and 300 𝑡 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.  \nIf 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 𝑞 = 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 𝑞 = 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 𝑞 = 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.  \nTaking a more quantitative approach, similar to Roedig et al. (2012), we calculate average profiles of the gas torque density, d 𝝉 / d 𝑅 , 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 𝑞 = 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 orbits 4 . We additionally split each profile into  \n<!-- image -->  \nFigure 5. Time-averaged surface density and torque maps for all the coplanar runs between 200 𝑡 bin and 300 𝑡 bin . Each simulation is identified by its label from Table 1. For eccentric runs ( 𝑒 = 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 𝑒 = 0 and to the apo- and pericentre for the runs with 𝑒 = 0 . 5; the outer dashed circle at 𝑅 = 2 𝑎 = 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.  \n<!-- image -->  \nFigure 6. Average radial profiles of the gas torque density, d 𝝉 / d 𝑅 , acting along ˆ J 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 𝑞 = 1, 1 / 3 and 1 / 10, from top to bottom, respectively. Each profile is produced by averaging over 100 orbits 4 , 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 ( 𝑅 = 𝑎 = 2 pc for circular orbits or 𝑅 = 𝑎 ( 1 ± 𝑒 ) = 2 ± 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.  \n<!-- image -->  \nthe gas torques from two radial regions; namely close to the binary (within 𝑅 = 𝑎 = 2 pc for circular orbits or 𝑅 = 𝑎 ( 1 ± 𝑒 ) = 2 ± 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.  \nAll 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 𝑅 < 𝑎 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 𝑅 ≳ 𝑎 generally being larger than the magnitude of the net positive torque at 𝑅 ≲ 𝑎 , 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).  \nIn 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.  \nWhiletheroughlocationofpeaksandtroughsintheprogradecases 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).  \nIn Table 2 we provide the average torques experienced by each binary (directed along its own angular momentum vector) split into separate time bins of ∼ 100 orbits each and decompose the contributions into being from gas at 𝑅 < 𝑎 sep , 𝑅 > 𝑎 sep , and all gas. Specifically, 𝑎 sep is the measured binary separation at the time torques are calculated, i.e. this would be 𝑎 for circular binaries and vary between 𝑎 ( 1 -𝑒 ) and 𝑎 ( 1 + 𝑒 ) for eccentric binaries moving between apo- and pericenter, taking into account the evolution of 𝑎 and 𝑒 over the course of the simulation (see Fig. 4). The colour coding scales linearly from 0 to ± 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).  \nFinally, the retrograde circular binary exhibits many more sharp  \nFigure 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𝜎 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, 𝜏 sim bin ) as well as two estimates for the gravitational torque using either a pure tree approach (blue, 𝜏 tree grav ) or pure direct summation (purple, 𝜏 direct grav ) and an estimated accretion torque (pink, 𝜏 est accretion ). Shaded squares show the estimated total torque experienced by the binary when adding 𝜏 est accretion to either 𝜏 tree grav or 𝜏 direct grav . 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 𝑞 01 𝑒 00 𝑖 45 binary exhibits a qualitative disagreement between the method in terms of the sign of the overall torque.  \n<!-- image -->  \nvariations 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 welldefined 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.  \nTo 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  \n𝝉 grav = 𝑀 1 ( r 1 × a 1 ) + 𝑀 2 ( r 2 × a 2 ) , (5)  \nwhere 𝑀 𝑖 , r 𝑖 and a 𝑖 are the mass, position and gravitational acceleration 5 of each black hole particle (which includes the sub-grid 𝛼  \nTable 2. Averages of the gravitational torque component along the binary angular momentum vector (in units of 10 -6 𝐺𝑀 2 bin 𝑎 -1 ) experienced by the different binary systems. We show the total torque exerted by all gas, as well as individual contributions from gas at 𝑅 < 𝑎 sep and 𝑅 > 𝑎 sep , where 𝑎 sep is the measured binary separation at each time torques are calculated ( 𝑎 for circular binaries and varying between 𝑎 ( 1 -𝑒 ) and 𝑎 ( 1 + 𝑒 ) 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 ± 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.  \naccretion disc) and 𝑖 = 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 𝑖 𝑡 ℎ black hole's angular momentum of  \nΔ J acc ,𝑖 = Σ 𝑗 GLYPH<0> r 𝑖 × 𝑚 𝑗 v 𝑗 GLYPH<1> , (6)  \nwhere 𝑚 𝑗 and v 𝑗 are, respectively, the mass accreted from and velocity of the the 𝑗 𝑡 ℎ cell and the sum is performed over all cells from which gas is accreted. In practise we do not log every accretion event  \ndue 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).  \nexplicitly 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., L gas ,𝑖 = Δ J acc ,𝑖 / Σ 𝑗 𝑚 𝑗 from each snapshot and estimate the total change in angular momentum due to accretion as  \nΔ J est acc = Δ 𝑀 1 ¯ L gas , 1 + Δ 𝑀 2 ¯ L gas , 2 , (7)  \nwhere Δ 𝑀 𝑖 and ¯ L gas , i are, respectively, the change in mass of the 𝑖 𝑡 ℎ black hole particle and the specific angular momentum of accreted gas averaged over consecutive snapshots. We then define an effective accretion torque as 𝝉 est accretion = Δ J est acc / Δ 𝑡 , where Δ 𝑡 is the time between snapshots.  \nTorques can act to change both the magnitude and direction of J 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 ( 𝜏 tree grav , shown in blue) and one in which the accelerations are calculated via direct summation over all gas cells ( 𝜏 direct grav , 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, 𝜏 est accretion , 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 | J bin | from the simulations we can calculate the average total torque, 𝜏 sim bin , 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.  \nIn 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 𝑁 max direct . 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.  \nReferring back to Fig. 4 and the discussion in Section 3.2, apart from in the circular prograde systems, the 𝑞 = 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 𝑞 = 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 underestimated by a factor of ∼ 2, and so like the other 𝑞 = 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 𝑞 = 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  \nepisodes 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': "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 · with respect to the CBD. Fig. 8 shows evolving density projections for the q03e00i45mod run in the 𝑥 -𝑦 (top) and 𝑥 -𝑧 (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., Miller & Krolik 2013; 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).  \nThe evolution of the torques acting on misaligned binaries is explored in Fig. 9, which shows the evolution of the 𝑥 , 𝑦 and 𝑧 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 ( 𝜏 sim bin , solid) and that calculated in postprocessing using a direct summation method ( 𝜏 grav , dashed), as well as adding the estimated accretion torque ( 𝜏 grav + 𝜏 acc , dotted, although this is essentially indistinguishable from the 𝜏 grav curve). The lines represent the torque averaged over 10 orbits while the shaded regions indicate the one𝜎 scatter. Each panel represents a different simulation, with all four initially misaligned binaries shown, as labelled. There are only very small differences in 𝜏 sim bin and 𝜏 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 · 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.  \nIsolating 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, d 𝝉 / d 𝑅 , acting along the direction of the binary angular momentum vector for the inclined runs. Each panel represents a different binary and colours illustrate different  \nFigure 8. Surface density maps for the q03e00i45mod binary illustrating the evolution of the CBD and black hole binary. The 𝑥 -𝑦 and 𝑥 -𝑧 projections are shown in the top and bottom rows, respectively with the time of the image shown in the top left 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.  \n<!-- image -->  \ntime ranges with each profile being produced by averaging over 100 orbits 4 . The profiles show broadly similar trends to their in-plane circular counterparts with several sharp features from gas within the binary orbit ( 𝑅 < 𝑎 = 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 𝑞 = 1 and 𝑞 = 1 / 3 systems should experience stronger positive torques in the 𝑅 < 𝑎 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.  \nWe 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).", '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  \nFigure 9. Evolution of the 𝑥 (blue), 𝑦 (pink) and 𝑧 (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 ∼ 100 𝑡 bin offset compared to the other misaligned binaries due to their misalignment being initialised after the relaxation.  \n<!-- image -->  \nblack 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.  \nFigure 10. As in Fig. 6 but showing the average radial profiles of the gas torque density, d 𝝉 / d 𝑅 , acting along ˆ J 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.  \n<!-- image -->', '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 th to 95 th percentile range of the inflow rates within time bins of length 5 𝑡 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 𝑅 = 2 𝑎 = 4 pc for circular binaries and 𝑅 = 3 𝑎 = 6 pc for eccentric binaries. Solid lines show net inflow and dashed lines show net outflow, both averaged over ∼ 5 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 mentioned in Section 2.1, quantities on resolved scales of the simulation can generally be presented in a scale-free manner, however, the sub-grid accretion disc model requires us to set a physical scale. As such, physical values applicable to our fiducial model parameters are shown on the left-hand axis of Fig. 11 and can be applied to all quantities presented. To illustrate how quantities might change for different system parameters we additionally show scale-free quantities on the right-hand axis, however, strictly speaking these are only applicable to CBD inflow/outflow rates and the  \nsimulated inflow rates onto the black hole particles, but not the mass flow rates through the sub-grid accretion disc.  \nAs 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 / 𝑀 acc ,𝑖 = 𝑓 Edd , i 𝑀 · ,𝑖 / 𝜂 𝑖 𝜏 Salp , where 𝜏 Salp = 450 Myr is the Salpeter time, and the growth rate of the black hole is given as / 𝑀 · ,𝑖 = ( 1 -𝜂 ) / 𝑀 acc ,𝑖 . 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, a number of 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˜noz 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 viscosity and thermodynamics within the cavity. While we employ an adiabatic equation of state with 𝛽 -cooling many of the previous works mentioned above employ a (locally-)isothermal equation of state instead. In addition the exact modelling of the black hole accretion may have an impact, 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, noting that several other previous studies also find preferential accretion onto the primary black hole (e.g., Ochi et al. 2005; Hanawa et al. 2010; Young et al. 2015). 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).  \nComparing 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 ∼ 120 orbits. This is simply down to the nature of the equation determining the black hole accretion rate, / 𝑀 acc = 𝑓 Edd 𝑀 · / 𝜂𝜏 Salp , where from Fiacconi et al. (2018) the dimensionless Eddington rate is given as  \n𝑓 Edd ≈ 0 . 76  \n× GLYPH<16> 𝜂 0 . 1 GLYPH<17> GLYPH<18> 𝑀 d 10 4 𝑀 ⊙ GLYPH<19> 5 GLYPH<18> 𝑀 · 10 6 𝑀 ⊙ GLYPH<19> -47 / 7 GLYPH<18> 𝑎 · 𝐽 d / 𝐽 · 3 GLYPH<19> -25 / 7 , (8)  \nFigure 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 𝛼 accretion discs, while the solid (red and blue) lines show the accretion rate onto each black hole from the sub-grid 𝛼 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. Note that scale-free units provided on the right-hand side of the figure should only be applied to quantities within the simulation domain, i.e. inflow/outflow rates, and not the black hole accretion rates from the sub-grid accretion disc model.  \n<!-- image -->  \nwhere 𝜂 is the spin-dependant radiative efficiency and 𝜏 Salp is the Salpeter time. Variations in 𝑓 Edd are dominated by changes in the disc mass, 𝑀 d and angular momentum, 𝐽 d , which both decrease during the first ∼ 120 orbits as / 𝑀 acc > / 𝑀 in and almost balance with each other. However, after a brief increase in the accretion rate the 𝑀 d decreases sufficiently to result in a sharp decline in the black hole accretion rate. After this / 𝑀 acc < / 𝑀 in and the accretion rate stabilises for the remainder of the simulation. A similar although less dramatic process is seen after ∼ 140 orbits in the q03e00i45mod run.', '3.3.2 Spin evolution': "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,  \nan understanding of black hole spin alignment is important to determine the recoil velocities of black hole merger remnants (Campanelli et al. 2007; Gonz'alez 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 ∼ 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  \nFigure 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.  \n<!-- image -->  \nwind) feedback (see e.g., Tchekhovskoy et al. 2011; Yuan & Narayan 2014; Liska et al. 2018; Talbot et al. 2021, 2022).  \nGerosa et al. (2015) developed an analytic model of black hole binary evolution in a CBD assuming the system is coplanar. 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  \n𝜏 align = 3 . 4 𝜏 Salp 𝑓 CBD 𝛼 AD 𝐹 ( 𝑞 ) GLYPH<18> 1 𝛼 2 𝐻 𝑅 GLYPH<19> 2 / 3 ×       𝑎 2 / 3 · , 1 GLYPH<16> 1 + 𝑞 𝛽 -1 GLYPH<17> / GLYPH<16> 𝑞 𝛽 -1 + 𝑞 𝛽 GLYPH<17> [Primary] 𝑎 2 / 3 · , 2 GLYPH<16> 𝑞 + 𝑞 𝛽 GLYPH<17> /( 1 + 𝑞 ) [Secondary] , (9)  \nwhere 𝑓 CBD = / 𝑀 CBD 𝜏 Salp / 𝑀 bin is the dimensionless Eddington rate of material feeding the cavity from the CBD, 𝛼 AD = 0 . 1 is the assumed 𝛼 -viscocity of the accretion disc, 𝛼 2 is the vertical 𝛼 -viscosity parameter given as (Pringle 1992; Ogilvie 1999)  \n𝛼 2 = 1 𝛼 2 GLYPH<16> 1 + 7 𝛼 2 GLYPH<17> 4 + 𝛼 2 , (10)  \n𝐻 / 𝑅 is the accretion disc aspect ratio, which for this analysis we set to 0 . 01, and 𝐹 ( 𝑞 ) is the fraction of material flowing into the cavity that is accreted by the binary components and takes the form 𝐹 ( 𝑞 ) = 0 . 8054 + 0 . 984 log ( 𝑞 ) + 0 . 3818 log 2 ( 𝑞 ) , 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, 𝛽 determines how material accreted by the binary is distributed between the primary and secondary black holes such that 𝑓 Edd , 2 / 𝑓 Edd , 1 = 𝑞 -𝛽 , whereby values of 𝛽 > 1 imply preferential accretion onto the secondary black hole. Gerosa et al. (2015) compared the alignment timescales to the predicted binary inspiral time due to interaction with a CBD which is given by  \n𝜏 CBD in = 3 4 ( 1 + 𝑞 ) 𝑀 cbd / 𝑀 bin + 𝑞 ( 1 + 𝑞 ) 2 𝜏 Salp 𝑓 CBD . (11)  \nThefull 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 𝑀 · , 2 << 𝑀 cbd << 𝑀 · , 1 (Lin & Papaloizou 1979; Artymowicz & Lubow 1994; Armitage & Natarajan 2002), 𝑀 · , 2 ∼ 𝑀 cbd << 𝑀 · , 1 (Syer", '& Clarke 1995; Ivanov et al. 1999; Lodato et al. 2009; Baruteau & Masset 2013) and 𝑀 · , 2 ≲ 𝑀 · , 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, 𝑞 , 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 ∼ 2 . 5 uncertainty in warp propagation theory. Similarly to Fig. 11, we show physical values applicable to our fiducial setup on the left-hand axis, and scale-free values on the right-hand axis. The latter can be freely used for the inspiral timescales, but only approximately apply to the alignment timescales under certain simplifying assumptions, which are discussed later in this section. The top panel shows the results of the analytic modelling assuming preferential accretion onto the seconday black hole, i.e., 𝛽 = 2, as originally assumed by Gerosa et al. (2015). This model leads to inhibited accretion onto the primary for low 𝑞 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 𝛽 = 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 𝑞 = 1 / 10.  \nWe 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. The accretion disc model of Fiacconi et al. (2018) assumes that spin-alignment occurs on a timescale of  \n𝜏 align ≈ 0 . 17 cos ( 𝜋 / 7 ) GLYPH<18> 𝑀 · 10 6 𝑀 ⊙ GLYPH<19> -2 / 35 GLYPH<18> 𝑓 Edd 𝜂 0 . 1 GLYPH<19> -32 / 35 𝑎 5 / 7 · Myr . (12)  \nUsing this definition, 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 𝑡 ℎ to 95 𝑡 ℎ percentile range. We only include snapshots after the first 100 orbits to limit any possible  \nFigure 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 timescales (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, 𝑓 Edd , 2 / 𝑓 Edd , 1 = 𝑞 -𝛽 , with 𝛽 = 2 . 0 and 1 . 0, respectively. In general we find that black hole spins align with their accretion discs on timescales of ≲ 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. Note that the scale-free units provided on the right-hand side of the figure are exact for the binary inspiral time and at best provide only an approximate guide for spin alignment timescales (see the main text for the necessary assumptions).  \n<!-- image -->  \ninfluence from the initial sub-grid parameter choices. The evolution of 𝜏 align given by Equation (12) as well as the black hole spin-black hole spin, 𝜏 bh bh , and black hole spin - binary angular momentum, 𝜏 bin bh , alignment timescales measured directly from the simulations are shown in Fig. A1. We note here that in the case of initially coplanar binaries both the analytic timescale and the timescales measured from the simulations generally match very well. However, in initially misaligned systems, non-negligable differences can occur between these different alignment timescale estimates, which we discuss further below and in more detail in Appendix A. Finally, the green  \nhatched 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.  \nSeveral 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 𝛽 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 𝑞 systems. These conditions result in more efficient spin-alignment for primary black holes even at low 𝑞 values. The only outliers from this general behaviour occur for the retrograde binaries, which, apart from the 𝑞 = 1 / 10primaryblack 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 ≲ 0 . 2 × 𝜏 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.  \nIn this work, we have assumed that the binary lives in an environment where hardening due to stellar interactions can be neglected. However, if the binary were to reside in a massive stellar bulge or central nuclear star cluster, stellar scattering may provide an additional channel to shrink the binary (Quinlan 1996; Milosavljevi'c & Merritt 2001, 2003; Berczik et al. 2006; Sesana et al. 2006; Rantala et al. 2017). Combining results from Milosavljevi'c & Merritt (2003) with the 𝑀 -𝜎 relation of Ferrarese et al. (2006), Cuadra et al. (2009) provide an approximate estimate for the timescale of binary shrinking due to stellar interactions as a function of binary mass and separation. It suggests that for our binary mass ( 𝑀 bin = 2 × 10 6 M ⊙ ), on ∼ parsec scales stellar scattering could be more effective than a CBD at shrinking the binary, and would proceed on timescales shorter than spin alignment. However, this should be seen as an approximate guide, given the dependence on the chosen 𝑀 -𝜎 relation and more generally on the exact stellar properties of the system. In any case, for sub-parsec separations as well as scenarios in which the binary does not reside in a stellar bulge or nuclear star cluster we expect that the CBD would dominate binary shrinking. In addition, given that 𝜏 align used in our sub-grid model is approximately inversely proportional to 𝑓 Edd and only weakly depends on M · (see Equation (12)), one might expect the ratio of 𝜏 align / 𝜏 CBD in to be roughly constant upon re-scaling of the system, i.e., for spin alignment to always occur on timescales shorter than binary shrinking due to the CBD. However, this behaviour should be seen as a qualitative guide and only approximately applies under the assumption that the mass flow rate through the sub-grid accretion disc follows the mass inflow rate from larger scales. As evident from Equation (8), 𝑓 Edd depends on a range of factors that may change upon re-scaling the large scale simulation properties, not withstanding the fact that changes in accretion rates  \nthat would arise from such re-scaling could lead to different accretion disc structures that would alter the alignment timescales.  \nIn addition to the potential role of stellar scattering shrinking the binary, there is another 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. However, if this is not the case, while black hole spins may align with their accretion flow, there is no requirement, at least initially, 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, and the total change in the angle between the black hole spins, Δ 𝜃 bh bh , is shown in the bottom left of each panel.  \nFirst, for all in-plane binary cases the accretion disc angular momenta generally remain aligned with the 𝑧 -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 𝑧 -axis, and hence each other. 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 𝑧 -axis, although the distributions of orientations in the case of the secondary black holes become broader for low mass ratios. We find that larger deviations from 𝑧 -axis alignment generally correspond to times in which inflow rates onto the sub-grid accretion disc are weaker, which fits with a picture in which the lower mass mini-discs around secondary black holes in low 𝑞 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 𝑧 -axis, with more scatter in the 𝑞 = 1 / 3 case compared to the 𝑞 = 1 / 10 case, reflecting that the latter exhibits a long-lived coherent mini-disc. On the other hand, the 𝑞 = 1 / 10 binary shows no inflows onto the secondary black hole, and while in the 𝑞 = 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 𝑞 = 1 / 10 case, the lack of coherent long-lived inflows in the retrograde binaries results in limited evolution of the black hole spin orientations.  \nIn 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 𝑞 = 1 case, the mini-disc  \nangular momenta also tend to reorient towards a common axis with little scatter, with the accretion disc angular momentum approaching the same direction. The black hole spins also appear to move towards alignment, with any periods of growing misalignment between the black hole spins and binary angular momentum being short and occuring early in the evolution (see Fig. A1). This suggests that in the 𝑞 = 1 case, the analysis presented in Fig. 13 should still hold if the system remains coplanar in the long term.  \nFor lower mass ratio binaries, while the binary and CBD angular momentum move towards alignment along a common axis, the minidisc angular momentum varies significantly (see also, Moody et al. 2019), even after the misalignment between the binary and CBD has reduced significantly. This is especially evident for mini-disc material around primary black holes, and leads to the accretion disc angular momenta also varying; sometimes exhibiting a precession-like motion. Considering the evolution of the black hole spin, we find that the misaligned accretion can result in deviations with respect to the initially coplanar systems. This is particularly evident, for example, for the secondary black holes in the q03e00i45 and q10e00i45 runs, as well as the primary black hole in the latter run. This being said, the overall trend is for the relative angle between the black hole spins to decrease in all binaries considered (see Δ 𝜃 bh bh values on each panel of Fig. 14, as well as its evolution in Fig. A2). We do, however, note that when considering the alignment evolution between black hole spin and binary angular momentum (see Appendix A and Fig. A1), their relative motion leads to periods during which the spins become more misaligned, suggesting that, at least in the early stages, alignment between a black hole's spin and the binary can be inhibited. In the long term, if the system is able to effectively globally align to become coplanar and subsequently evolve unperturbed, spin alignment should be expected (Miller & Krolik 2013; Gerosa et al. 2015). We stress, however, that we have studied only a small parameter space in initial black hole spins, CBD inclination angles and eccentricities here. It is also important to highlight that we have not considered regimes in which disc breaking/tearing occurs, and that the sub-grid accretion disc angular momentum only couples to the resolvable scales via accretion and is not impacted by the external tidal field of the gas and companion black hole, which could also effect the evolution of both the accretion disc and black hole spin (Gerosa et al. 2020; Nealon et al. 2022). We discuss these points further in Sections 4.3 and 4.4 and note that we suggest that longer term simulations in a range of more realistic environments would be critical to fully understand the late time evolution of the black hole spin orientation in low𝑞 , misaligned systems and its implications for GW emission and black hole recoil velocities. 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˜noz 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˜noz et al. 2019, 2020; D'Orazio & Duffell 2021; Siwek et al. 2023b) make two key assumptions: i) the gas can be treated  \nFigure 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. The total change in the angle between the black hole spins, Δ 𝜃 bh bh , is shown in the bottom left of each panel. In general, coplanar runs clearly tend towards alignment between black hole spins and the global angular momentum vector, however, mini-discs in low 𝑞 misaligned systems do not align rapidly which may initially impact the global alignment in such systems.  \n<!-- image -->  \nas 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 selfgravitating 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 𝛽 -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, 𝑞 , eccentricity, 𝑒 and inclination angle, 𝑖 , while fixing the initial total binary mass ( 𝑀 bin = 2 × 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 GWemission 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'c 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˜noz et al. 2019; Tiede et al. 2020; Franchini et al. 2022; Dittmann & Ryan 2022, 2024; Siwek et al. 2023b, for previous examples of such processes).  \nRoedig 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 𝑅 < 𝑎 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˜noz 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 𝑅 < 𝑎 to super-Keplerian inflows that bend in front of the black holes, while at 𝑅 > 𝑎 , 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.  \nAccretion 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˜noz 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 & Duffell2021),althoughnotnecessarily negligible (Mu˜noz et al. 2019;  \nSiwek 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)  \n/ 𝑎 𝑎 = 2 𝜏 bin 𝐽 bin + 2 𝑒 1 -𝑒 2 / 𝑒 -/ 𝑀 bin 𝑀 bin -2 / 𝜇 𝜇 , (13)  \nwhich is found by differentiating Equation (4), rearranging and setting / 𝐽 bin = 𝜏 bin . Neglecting the eccentricity term and rearranging Equation (13), we can define the maximum torque under which an equal mass, circular binary shrinks,  \n𝜏 bin < 3 8 / 𝑀 bin √︁ 𝐺𝑀 bin 𝑎, ≃ 1 . 12 × GLYPH<18> / 𝑀 bin 2 × 10 -4 M ⊙ yr -1 GLYPH<19> 10 -6 𝐺𝑀 2 bin 𝑎 -1 , (14)  \nwhere we use 𝑀 1 = 𝑀 2 = 𝑀 bin / 2 and assume / 𝑀 1 = / 𝑀 2 = / 𝑀 bin / 2. As such, binaries experiencing a positive torque can still shrink provided / 𝑀 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 𝐽 bin . However, given the binary's average accretion rate of ∼ 2 × 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.  \nIn 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˜noz et al. 2019, 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˜noz (2023), circular binaries above a mass ratio of 𝑞 ∼ 0 . 1 -0 . 3 expand, while for mass ratios below 𝑞 ∼ 0 . 05 -0 . 1 they shrink (e.g., Duffell et al. 2020; Mu˜noz et al. 2020; Siwek et al. 2023b), although recent parameter studies suggest that such behaviours are sensitive to both the disc thickness and viscosity (Dittmann & Ryan 2022, 2024). Eccentric binaries have also been observed to expand for particular combinations of 𝑒 and 𝑞 (e.g. Mu˜noz et al. 2019; D'Orazio & Duffell 2021; Siwek et al. 2023b), although D'Orazio & Duffell (2021) notes that initially eccentric binaries ( 𝑒 ≳ 0 . 1) tend towards an equilibrium eccentricity ( 𝑒 ≃ 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.  \nSimulations 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 ( 𝑀 cbd / 𝑀 bin = 0 . 1) CBD to simulations employing 'hyper-Lagrangian' refinement. They find that while positive torques from gas at 𝑅 ≲ 𝑎 can lead to mild binary expansion, this  \nis 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 𝛽 -cooling, from 𝛽 = 0 (locallyisothermal) to 𝛽 = 4, Wang et al. (2023) found binary expansion for the majority of values explored except for 0 . 7 ≲ 𝛽 ≲ 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 𝛽 ≃ 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.  \nThe viscous torque, which is proportional to 𝛼 ( 𝐻 / 𝑅 ) 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 ( 𝐻 / 𝑅 ∼ 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 𝛼 -viscosity, thick discs led to binary expansion due to accretion of high angular momentum material (see also, Franchini et al. 2022). While by varying 𝐻 / 𝑅 (or Mach number) at fixed kinematic viscosity, 𝜈 , Tiede et al. (2020) found expansion of binaries in thick discs and shrinking in thin discs driven by gravitational torques, Dittmann & Ryan (2022) found that binary expansion is insensitive to the value of 𝜈 in thick discs; however, in thin discs binaries continue to expand for high values of 𝜈 , but shrink ever more rapidly for decreasing values of 𝜈 . Other works have instead varied 𝛼 for a fixed 𝐻 / 𝑅 , finding only negligible differences in the accretion rate normalised torques (Mu˜noz 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 initial relative alignment between a CBD and binary could be arbitrary. Black hole growth via randomly orientated accretion events, such as chaotic accretion via many small clouds (King & Pringle 2006, 2007), means that the binary may interact with many misaligned accretion events and CBDs before eventually merging (Nixon et al. 2013). For misaligned binaries ( 𝑖 = 45 · ), 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 𝑞 = 1 and 𝑞 = 1 / 3 misaligned binaries experience overall weaker negative gravitational torques, which is a consequence of stronger positive net torques from gas at 𝑅 < 𝑎 (see Table 2), and in general experience stronger positive accretion torques. Simulations of a misaligned 𝑞 = 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 also find that the mini-disc orientations are variable in misaligned systems.  \nWe discuss the implications of such variability in the context of our simulations on spin alignment, recoil velocities and GW detection rates in Sections 3.3.2 and 4.4.  \nWhile 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 𝛽 -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 𝑖 < -𝐽 cbd / 2 𝐽 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.  \nDunhill 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 𝑎 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 𝑅 < 𝑎 (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 (2024), 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 (2024) additionally find that for 𝑒 ≳ 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; Haring & Rix 2004; Gultekin et al. 2009; Kormendy & Ho 2013; McConnell & Ma 2013; Sahu et al. 2019), including dwarfs (Reines et al. 2013, 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'c 2008; D'Orazio et al. 2013; Farris et al. 2014; Miranda et al. 2017; Duffell et al. 2020; Mu˜noz 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 𝑅 = 𝑎 is accreted, however, they find a reduction in the accretion rate of ∼ 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.  \nWe show inflow rates onto black holes in Fig. 11 and similarly to previous works we find that they vary on times scales of ∼ 𝑡 bin (e.g., MacFadyen & Milosavljevi'c 2008; Cuadra et al. 2009; Roedig et al. 2012; Mu˜noz & Lai 2016; Lai & Mu˜noz 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 ∼ 5 𝑡 bin associated with the formation of a lump at the inner edge of the CBD (Shi et al. 2012; MacFadyen & Milosavljevi'c 2008; D'Orazio et al. 2013; Mu˜noz & Lai 2016; Miranda et al. 2017; Wang et al. 2023; Lai & Mu˜noz 2023), although this is only expected when gas can cool rapidly Wang et al. (2023).  \nMany simulations find that accretion occurs preferentially onto the secondary black hole and drives the binary mass ratio towards unity (e.g., Roedig et al. 2012; Farris et al. 2014; Duffell et al. 2020; Mu˜noz et al. 2020; Dittmann & Ryan 2021; Siwek et al. 2023a). In our simulations we do not observe this behaviour, instead finding that accretion occurs preferentially onto the primary black hole, which has important implications for spin alignment and GW emission. A key difference in our work is that we employ an inflow-regulated accretion model (described in Section 2.4), while the vast majority of other simulations employ some form of sink particle prescription, whereby gas residing within a fixed sink radius is removed from the simulation domain and added to the black hole. This change in approach may affect the accretion rate onto both black holes, however, we expect the impact of this difference to be small given that our model will still remain dependent on the amount of mass in the vicinity of each black hole. Hence, the thermodynamic and viscous modelling of the gas likely plays the critical role in understanding the accretion process. Our chosen cooling prescription means that heat generated via compression or shocks is not immediately radiated away. This is particularly evident at the outer edge of the cavity where the disc material puffs up (see Fig. 3). Interestingly, several previous studies of hot/thick discs have found preferential accretion onto the primary instead of the secondary (Ochi et al. 2005; Hanawa et al. 2010; Young et al. 2015), with Young et al. (2015) specifically showing that the amount of gas moving from the secondary to primary Roche lobe is dependant on gas temperature, with more gas reaching the primary in hotter systems. Additionally, Roedig et al. (2012) showed an adiabatic equation of state results in less well-defined mini-discs and streams compared to the isothermal case. They further found that  \nwhile the accretion rate onto the secondary black hole is independent of the sink radius for an isothermal gas, it scales with the sink radius in the adiabatic case, suggesting that for small enough sink radii, preferential accretion may be reversed. As well as the gas thermodynamics, the viscosity likely also plays an important role. Dittmann & Ryan (2024) recently showed that whether a constant kinematic viscosity or an 𝛼 -viscosity is used affects the level of preferential accretion on to the secondary, while Duffell et al. (2020) found that the ratio of / 𝑀 2 / / 𝑀 1 tends towards unity at low viscosity. While we reserve a detailed study to future work, in comparing these previous findings with the characteristics of our simulations, we conclude that the lack of preferential accretion on the secondary is driven by a combination of slow cooling, and hence 'hot' discs/cavity, the low intrinsic viscosity in our simulations, and to some extent the specific mass accretion prescription we have adopted. An important open question that future studies need to address is the likely range of visco-thermo-dynamical states of mini-discs representative of the entire population of black hole binaries.  \nIf accretion occurs preferentially onto the secondary black hole, it would lead to parsec-separation binary AGN being primarily powered by the secondary. A potential piece of evidence countering this situation comes from observations of large-scale radio jets that exhibit precession on timescales consistent with geodetic precession driven by a low-mass binary companion (Krause et al. 2019). Given the powerful nature of these sources, it would seem likely that the jet originates from accretion onto the primary black hole, which tentatively disfavours preferential secondary accretion. However, further observations, for example with LOFAR, SKAO, or precursor facilities, as well future GW detections with LISA are needed to definitively solve this open question.  \nFinally, LISA will allow us to build our understanding of the black hole populations spins (Amaro-Seoane et al. 2017, 2023), as such it is important to understand how black hole spins evolve during CBD phases of their evolution. The sub-grid accretion disc model of Fiacconi et al. (2018) allows us to self-consistently evolve the black hole spin. Several similar models have been developed for use in galaxy formation simulations (e.g., Dubois et al. 2014; Fiacconi et al. 2018; Bustamante & Springel 2019; Cenci et al. 2021; Sala et al. 2021; Huˇsko et al. 2022b; Massonneau et al. 2023), and have been used to determine feedback properties (e.g., Dubois et al. 2014; Beckmann et al. 2019, 2022; Talbot et al. 2021, 2022, 2024; Huˇsko et al. 2022b; Bollati et al. 2023b,a). Simulations of binary evolution in large ( ∼ 100 pc) circumnuclear discs have tracked spin evolution using similar sub-grid models to ours (Bollati et al. 2023b). However, while previous simulations of CBDs have tracked the instantaneous rotational angular momentum of the accretor (e.g., Mu˜noz et al. 2019, 2020; Dittmann & Ryan 2021), to the best of our knowledge ours is the first work to simulate parsec scale binaries and CBDs in order to self-consistently evolve sub-grid accretion discs and incorporate physical spin-alignment timescales. This allows us to not only evolve the spin magnitude but also its alignment. In general, the spins of black holes in systems in which the binary and CBD are coplanar gradually align with the global angular momentum. The picture is less clear cut, at least during early evolution, in systems where the binary and CBD are initially misaligned because the mini-disc angular momentum can be highly variable (see also, Moody et al. 2019). This can inhibit black hole spins from aligning with the binary angular momentum over long timescales, with several of our simulations undergoing non-negligible periods during which the black hole spin and binary angular momentum become further misaligned (see Fig. A1). This being said, despite some early deviations, the relative angle between black hole spins does reduce  \nduring the course of the simulations, with Δ 𝜃 bh bh differing with respect to the coplanar cases by roughly ± 30%. As such, if the system evolves unperturbed over a sufficiently long timescale such that the mini-discs are able to align with the global angular momentum of the system, then spin-alignment should be expected (Miller & Krolik 2013; Gerosa et al. 2015). However, as noted previously, the study we present here covers a restricted range of initial spin, binary and CBD orientations/properties, and the sub-grid accretion disc model does not capture all external processes driving the angular momentum evolution of the accretion disc, in particular the effect of the external tidal field of the gas and companion black hole, which could directly impact spin alignment (see e.g., Gerosa et al. 2020). Additionally, very few simulation works have studied initially misaligned systems and it would be pertinent to explore the long term spin evolution in systems covering a wide range of initial inclination angles and CBDproperties in future works, including in the regimes where tearing/breaking of the CBD (Nixon et al. 2013; Aly et al. 2015) and/or accretion discs (Nixon et al. 2011b, 2012; Gerosa et al. 2020; Nealon et al. 2022) occurs, and for cases of polar alignment of the CBD about the binary (Aly et al. 2015). Further the scenario in which the secondary black hole is embedded within the CBD could provide interesting insights into spin evolution, with recent studies of stellar mass black holes embedded in AGN accretion discs finding that turbulence can result in stochastic accretion and spin evolution (Chen & Lin 2023). As we discuss next, spin-alignment has important implications for merger remnants and GW detection, as well as feedback physics (e.g., efficiency and direction).", '4.4 Implications for black hole mergers and gravitational wave detection': 'Once at parsec scales, the evolution of the binary depends upon the environmental properties (Amaro-Seoane et al. 2023). In sufficiently gas-rich environments a CBD may form and effectively torque the binary. As discussed above, simulations have found that interactions between the binary and CBD can lead to either expansion or shrinking of the binary (Cuadra et al. 2009; Roedig et al. 2012; Roedig & Sesana 2014; Miranda et al. 2017; Tang et al. 2017; Moody et al. 2019; Mu˜noz et al. 2019, 2020; Duffell et al. 2020; Heath & Nixon 2020; Tiede et al. 2020; D\'Orazio & Duffell 2021; Franchini et al. 2021, 2022; Siwek et al. 2023b; Wang et al. 2023), which makes it non-trivial to determine the rate of GW detections. Based on our simulations we expect massive, self-gravitating CBDs to provide a plausible channel through which a binary can shrink. However, it is worth highlighting that the CBDs in reality may be more dynamic with significant scope to improve modelling of the gas thermodynamics and include the effects of star formation as well as stellar and black hole feedback. Additionally, to reach merger it will likely require numerous individual accretion events and we find that retrograde events are likely more effective than prograde events at shrinking a binary (see also, Nixon et al. 2011a; Dunhill et al. 2014; Tiede & D\'Orazio 2024).  \nAs discussed in Section 4.3, not only does the interaction of the binary with a CBD alter its orbital properties and provide a channel through which to bring the black holes closer together but also impacts the individual black hole properties that shape the eventual merger remnant properties, GW emission and the stochastic GW background. The GW strain amplitude scales with the binary chirp mass,  \nM = GLYPH<0> 𝑀 · , 1 𝑀 · , 2 GLYPH<1> 3 / 5 GLYPH<0> 𝑀 · , 1 + 𝑀 · , 2 GLYPH<1> 1 / 5 = " 𝑞 ( 1 + 𝑞 ) 2 # 3 / 5 𝑀 bin , (15)  \nwhich for a circular binary scales as M 5 / 3 (Thorne 1987). As such, for a fixed binary mass, equal mass binaries provide the highest amplitude signals. Accretion can play an important role in determining the amplitude of the GW background both in terms of simply increasing the binary mass and in setting the final mass ratio (e.g., Sesana et al. 2009b; Kelley et al. 2017; Siwek et al. 2020).  \nThe alignment (or misalignment) of black hole spins can directly impact the GW emission (Vecchio 2004; Klein et al. 2009, 2013; Lang et al. 2011; Gerosa & Moore 2016). Additionally, high magnitude misaligned spins can result in merger remnants with large recoil velocities (Campanelli et al. 2007; Gonz\'alez et al. 2007; Lousto & Zlochower 2011, 2013; Lousto & Healy 2019; Sperhake et al. 2020). This scenario could lead to black holes being ejected from their host galaxy (Redmount & Rees 1989; Volonteri et al. 2010; Gerosa & Sesana 2015), thus preventing or severely reducing possible future mergers in the galactic core and hence reducing the number of GW events (Sesana 2007; Sesana et al. 2009a; Blecha & Loeb 2008). Considering both the binary-CBD and spin alignment timescales, analytical estimates by Miller & Krolik (2013) suggest that spin alignment should occur in all but extreme mass ratio binaries. When considering coplanar systems, Gerosa et al. (2015), predict that preferential accretion onto secondary black holes inhibits spin-alignment in low mass ratio systems, although it can proceed for 𝑞 ≳ 0 . 2 systems (see also, Gerosa et al. 2020). Our simulations present somewhat different behaviour. The lack of preferential secondary accretion in our simulations means that efficient spin alignment can occur even for low 𝑞 systems and we note that longer timescales found in retrograde systems are likely a result of our inflow accretion model and a different gas cooling prescription compared to the majority of the literature. Secondly, as discussed in the previous section, additional complexity arises for systems in which there is initial misalignment between the binary and CBD. However, as long as the system evolves unperturbed such that the binary, CBD and mini-discs come into alignment, we would also expect spin alignment to occur. Additionally, if the black holes co-exist in a large-scale post-merger coherent flow, such as a nuclear gas disc, their spins may align prior to the CBD phase (e.g. Bogdanovi\'c et al. 2007; Dotti et al. 2010). Once the binary is at ≲ parsec separation and a cavity forms, the subsequent spin evolution will depend on whether the binary experiences coherent or chaotic/randomly orientated accretion. The former would likely drive and maintain high spins and their alignment. On the other hand, chaotic/randomly orientated accretion would result in misaligned accretion events (King & Pringle 2006, 2007; Nixon et al. 2013; Dunhill et al. 2014) for which we believe long-term spin-alignment requires further study. It is important to stress that high recoil velocities also require high black hole spins. While spin up could be achieved via long term coherent accretion onto the binary (which is expected to lead to spin alignment), chaotic accretion is expected to spin down black holes in the long term due to the higher specific angular momentum of retrograde events (King & Pringle 2006; King et al. 2008). This being said, for the system studied here, in which we assume the binary forms from two already highly spinning black holes, the timescale for binary spin down could be longer than the binary inspiral time, especially in the case of low mass ratio and/or retrograde systems, meaning that the black holes could still have non-negligible spins at merger. LISA will be able to put important constraints on the population of black hole spins in merging systems and therefore shed light on black hole accretion physics (Amaro-Seoane et al. 2023).', '4.5 Caveats and future work': 'There are some caveats and subsequent areas of future development that we briefly highlight here. Firstly in terms of the numerical methods, in this work, we use a hierarchical time integration (Springel et al. 2021), including hybrid force calculations that employs direction summation at small radii and a tree method at larger radii. As shown in Section 3.2.2, while this gives broadly consistent results compared to a pure direct summation calculation and maintains stable binary orbits, small differences in gravitational torques, however, would suggest that in future works direct summation between the binary components and the gas component should be used (although as in other works, to ease computational expense a tree can still be used for the gas self-gravity, e.g., Cuadra et al. 2009; Roedig et al. 2012). The black hole accretion model also differs from previous works in calculating mass inflow rates onto the sub-grid accretion disc as opposed to using a sink particle approach. A future study is important to understand the extent to which differences in these approaches affect both the binary and black hole evolution, and would allow the development of more physically realistic and robust accretion prescriptions and place limits on the required resolution to achieve this.  \nThe parameter space when modelling CBDs and binaries is vast and there are many avenues to consider when expanding on the presented simulations. Based on the above discussion, as well as extending the range of binary properties, a key area to probe is the thermal modelling of the gas. The ratio between accretion rates onto the secondary and primary black holes is very likely sensitive to the exact thermal state of the gas and hence the employed cooling prescription (see discussion in Section 4.3). Many processes govern the radiative heating and cooling of the gas around a binary, however, fully self-consistent simulations employing full radiativetransfer would be prohibitively expensive for wide parameter surveys. As such, variations in the 𝛽 -cooling, such as those studied by Wang et al. (2023) would provide a tractable and insightful avenue. Employing 𝛽 -cooling within self-gravitating CBDs can result in a system that self-regulates its temperature/thickness (as long as 𝛽 is not so small as to drive fragmentation) and as such other parameters that would be useful to investigate are the total disc mass and viscosity (see e.g., Franchini et al. 2021), including the regime in which disc breaking occurs (Nixon et al. 2013; Aly et al. 2015). Furthermore, in the currently employed 𝛽 -cooling prescription (see Equation (3)) the positions of gas cells is always calculated with respect to the binary centre of mass. However, one could alternatively treat gas within the mini-discs separately by instead using the distance to the corresponding binary component. If using a universal 𝛽 value, this would reduce cooling times of gas close to the secondary black hole in low 𝑞 systems, although we note that the choice of 𝛽 need not be the same around each black hole and in the CBD. Feedback from accreting black holes may also impact the binary evolution (del Valle &Volonteri 2018; Bollati et al. 2023b) as well as the properties of gas feeding the binary from larger scales (Dehnen & King 2013; Zubovas 2015; Bourne & Sijacki 2017; Faber & Dehnen 2018; Talbot et al. 2024). In future, we aim to include feedback prescriptions coupled to the sub-grid accretion disc model (Talbot et al. 2021, 2022, 2024, Koudmani et al., in prep) to explore this further.  \nFinally, there are some caveats related to the sub-grid modelling of the black hole-accretion disc system. Firstly we assume that the sub-grid accretion disc is a thin Shakura & Sunyaev (1973) disc. Such discs are expected to exist for high Eddington ratios, however, the low Eddington ratios seen in our simulations suggest that a thick or truncated disc may be more appropriate (Yuan & Narayan 2014;  \nKoudmani et al. 2024). In order to achieve an analytically tractable solution, the sub-grid accretion disc model of Fiacconi et al. (2018) was derived under the small-warp approximation (e.g., Martin et al. 2007; Perego et al. 2009; Dotti et al. 2013; Fiacconi et al. 2018), however, this may not be applicable when the black hole spin and accretion disc are initially highly misaligned. Additionally, depending on the disc properties, instead of attaining a steady warp the accretion disc may break and result in either reduced or enhanced rates of accretion (Nixon et al. 2011b, 2012; Gerosa et al. 2020; Nealon et al. 2022). Given that we resolve gas down to the scale of the sub-grid accretion disc, the mini-discs that we simulate should potentially also experience torquing by the spin of their host black hole. However, no such back reaction from the sub-grid model onto the resolved gas is included in our simulations and may affect the mini-disc behaviour, especially in the case of misaligned binaries. Finally, in a binary, each accretion disc is perturbed by the other black hole (Martin et al. 2009; Tremaine & Davis 2014; Miller & Krolik 2013), which can impact the spin alignment (Miller & Krolik 2013; Gerosa et al. 2020) and whether or not the accretion disc breaks (Nealon et al. 2022).', '5 CONCLUSIONS': 'Wehaveperformedhigh-resolution simulations of ∼ parsec scale separation black hole binaries surrounded by gas-rich gravito-turbulent CBDs ( 𝑀 cbd / 𝑀 bin = 0 . 1; 𝑀 bin = 2 × 10 6 M ⊙ ) including the effects of self-gravity. Employing a super-Lagrangian refinement techniques allows us to effectively resolve the streams and mini-discs that form within the cavity. By using a sub-grid model that mediates the black hole mass and spin evolution via a thin 𝛼 -disc prescription we have, for the first time, been able to self-consistently track the spin evolution of the individual black holes in such CBD simulations and to study their alignment. The full simulation suite includes twelve setups varying the initial binary mass ratio ( 𝑞 = 1, 1 / 3 or 1 / 10), eccentricity ( 𝑒 = 0 or 0 . 5) and inclination angle ( 𝑖 = 0 · , 45 · or 180 · ). Our main findings are summarised as follows:  \n- · The evolution of the CBD, cavity, streams and mini-discs depends sensitively on the binary properties (Section 3.1). Higher mass ratio and eccentric binaries result in larger cavities, while retrograde binaries result in smaller cavities. Mini-disc sizes typically scale with the black hole mass (roughly filling the Roche lobe for prograde circular binaries), although eccentric binaries have truncated mini-discs.\n- · We find that the net torque experienced by binaries drives them towards merger (Sections 3.2), with low mass ratio binaries and retrograde binaries expected to shrink more rapidly.\n- · For prograde binaries, the net gravitational torques experienced by the binary depend on the balance of net positive torques from gas at 𝑅 ≲ 𝑎 and net negative torques from gas at 𝑅 ≳ 𝑎 . For retrograde binaries, the net gravitational torque is dominated by negative torques from gas a 𝑅 ≲ 𝑎 . Accretion torques are always positive and generally sub-dominant.\n- · Unlike many previous works, accretion in our simulations does not proceed preferentially onto the secondary black hole (Section 3.3.1). This implies binaries to generally maintain a constant mass ratio up to merger, which for low mass ratio systems would result in smaller chirp masses, lower GW strain amplitudes and a weaker GW background.\n- · Weattribute the lack of preferential accretion onto the secondary black hole to a combination of using an adiabatic equation of state with slow 𝛽 -cooling, low intrinsic viscosity and possibly due to dif-  \nces between our inflow rate accretion model and the traditional sink particle model often employed in the literature.  \n- · Alignment of black hole spins and their accretion discs occurs on timescales of ≲ the binary inspiral time, with the majority of alignment times being much shorter (Section 3.3.2). The lack of preferential secondary accretion means this is true even for the primary black hole in low mass ratio systems.\n- · In coplanar systems we expect black hole spins to align with each other and the global angular momentum of the system prior to merging. Additionally, while the mini-discs, which feed the black hole accretion discs that torque the black hole spin, do not align on short timescales, especially for low mass ratio systems, black hole spins do appear to align with each other in the long term. However, we caution that wider parameter space studies, the effects of external tidal fields on the sub-grid accretion disc, and scenarios including disc breaking/tearing are necessary to fully address the question of spin alignment efficiency in CBDs. Further, in galaxies with bulges or central nuclear star clusters, stellar scattering may speed up binary shrinking on parsec scales, meaning that spin alignment is only likely to be efficient on sub-parsec scales. An understanding of the regimes in which the spins of merging black holes align with each other and the binary angular momentum (or not) is critical in determining expected black hole recoil velocities, galaxy occupation fractions and GWevent detection rates.', 'ACKNOWLEDGEMENTS': "We thank the referee for constructive comments on our original manuscript. We additionally thank Alex Dittmann, Pedro Capelo, Martin Krause, Jim Pringle and Chris Nixon for useful discussions, feedback and comments. M.A.B. acknowledges support from a UKRI Stephen Hawking Fellowship (EP/X04257X/1). M.A.B., D.F. and D.S. acknowledge support by European Research Council Starting Grant 638707 'Black holes and their host galaxies: coevolution across cosmic time'. M.A.B and D.S. additionally acknowledge support from the Science and Technology Facilities Council (STFC). The simulations were performed on the DiRAC Darwin Supercomputer hosted by the University of Cambridge High Performance Computing Service (http://www.hpc.cam.ac.uk/), provided by Dell Inc. using Strategic Research Infrastructure Funding from the Higher Education Funding Council for England and funding from the Science and Technology Facilities Council; the COSMA Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility.This equipment was funded by a BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grant ST/K00087X/1, DiRAC Operations grant ST/K003267/1 and Durham University.", 'DATA AVAILABILITY': 'The data underlying this article will be shared on reasonable request to the corresponding author.', 'REFERENCES': 'Young M. D., Baird J. T., Clarke C. J., 2015, MNRAS, 447, 2907 Yuan F., Narayan R., 2014, ARA&A, 52, 529 Zrake J., Tiede C., MacFadyen A., Haiman Z., 2021, ApJ, 909, L13 Zubovas K., 2015, MNRAS, 451, 3627 del Valle L., Volonteri M., 2018, MNRAS, 480, 439', 'APPENDIX A: SPIN ALIGNMENT TIMESCALES': "The spin alignment timescales presented in Fig. 13 and discussed in section 3.3.2 represent the median expected alignment times between the black hole spin and its accretion disc angular momentum over ∼ 400 binary orbits. In Fig. A1 we instead show the evolution of the spin alignment timescale given by Equation 12 for the primary (blue) and secondary (red) black hole for each run by the thick lines in each panel. However, the alignment timescale between the black hole spins and binary angular momentum measured directly from the simulation is potentially more relevant for considering misalignment triggered black hole recoils (Campanelli et al. 2007; Gonz'alez et al. 2007; Lousto & Zlochower 2011, 2013; Lousto & Healy 2019; Sperhake et al. 2020), and is given by 𝜏 bin bh = 𝜃 bin bh /| / 𝜃 bin bh | , where 𝜃 is the angle between the black hole spin and binary angular momentum. The evolution of 𝜏 bin bh is shown by the thin blue and red lines, for the primary and secondary black holes, respectively. Additionally, the timescale on which the black hole spins align with each other, 𝜏 bh bh = 𝜃 bh bh /| / 𝜃 bh bh | , is shown by the black line in each panel. In all cases dotted lines represent periods during which 𝜃 is increasing, i.e., the spin and binary are becoming more misaligned.  \nThere are a couple of key points to note. Firstly, for the initially coplanar systems we find that the evolution of the analytic alignment timescales given by Equation (12), the black hole spin-binary alignment timescale, 𝜏 bin bh , and the spin-spin alignment timescales, 𝜏 bh bh , follow each other very well, with only fairly minor differences. We also note that for these systems there is at most mild evolution of the alignment timescales. On the other hand, if we consider initially misaligned systems, the alignment timescales also undergo a richer and more variable evolution, with potentially significant differences between the three timescale estimates. In particular, 𝜏 bin bh exhibits more variation including periods of slow alignment and increasing misalignment, in some cases over significant periods of time, due to the relative evolution of both the black hole spin and binary angular momentum. This being said, while 𝜏 bh bh can also exhibit slow alignment and periods of increasing misalignment, this is confined to very early times and in general 𝜏 bh bh roughly follows the behaviour of 𝜏 align at later times.  \nTo consider the latter point further, in Fig. A2 we show the evolution of the change in angle between the black hole spins, Δ 𝜃 bh bh , for all prograde, initially circular binaries. Each panel shows a different binary mass ratio with black and green lines showing initially coplanar and misaligned systems respectively. We see that in general the change in the spin-spin alignment angle is fairly similar for all runs, although initially misaligned systems can undergo initial short periods of growing misalignment. By the end of the simulations, the difference in Δ 𝜃 bh bh between runs of a given mass ratio differ by at most ± 30% . Based on this behaviour, and assuming the system survives unperturbed for a sufficiently long time it is expected that all systems will eventually come into alignment (Miller & Krolik 2013; Gerosa et al. 2015). Although we point the reader to Sections 3.3.2, 4.3 and 4.4 for discussion of caveats and possible implications.  \nFigure A1. Evolution of the black hole spin alignment timescales. Each panel represents a different binary, with the grey shaded regions illustrating the first ∼ 100 binary orbits that are excluded from the mean timescale estimates used in Fig. 13. Blue and red lines represent the primary and secondary black hole in each system, respectively. Thick lines show the alignment timescale assumed by the Fiacconi et al. (2018) model as given in Equation (12), while thin lines show the alignment timescale between the black hole spin and binary angular momentum calculated directly from the simulation and given by 𝜏 bin bh = 𝜃 bin bh /| / 𝜃 bin bh | , where theta is the relevant alignment angle. The timescale on which the black hole spins align with each other 𝜏 bh bh = 𝜃 bh bh /| / 𝜃 bh bh | is additionally shown by the black line. In all cases solid and dotted lines indicate when / 𝜃 is negative or positive, i.e. when the system is aligning or misaligning, respectively. For coplanar systems, the analytic alignment timescale between the black hole spin and its accretion disc matches well to the simulated alignment timescale between the black hole spin and binary angular momentum, as one would expect. However, it is clear that differences arise between these timescales in the case of systems that are initially misaligned, indicating that global alignment could initially be inhibited in such systems.  \n<!-- image -->  \nFigure A2. Evolution of the change in angle between the spin of each black hole, Δ 𝜃 bh bh . Each panel represents a different binary mass ratio as labelled in the top right, noting that we only include prograde, initially circular binaries in order to compare initially misaligned systems (shown in green) with their coplanar (shown in black) counterparts. In general the evolution is similar for all binaries, although initially misaligned systems may initially experience periods in which Δ 𝜃 bh bh increases. For a given mass ratio, after 500 binary orbits the difference in Δ 𝜃 bh bh between initially coplanar and initial misaligned systems is at most ∼ ± 30%.  \n<!-- image -->"}