The Evolution of Inclined Binary Black Holes in the Disks of Active Galactic Nuclei

The accretion disks that fuel active galactic nuclei (AGNs) may house numerous stars and compact objects, formed in situ or captured from nearby star clusters. Embedded neutron stars and black holes may form binaries and eventually merge, emitting gravitational waves detectable by LIGO/VIRGO. AGN disks are a particularly promising environment for the production of high-mass gravitational-wave events involving black holes in the pair-instability mass gap, and may facilitate electromagnetic counterparts to black hole binary mergers. However, many orders of magnitude separate the typical length scales of binary formation and those on which gravitational waves can drive binary inspirals, making binary mergers inside the disk uncertain. Previous hydrodynamical simulations of binaries have either been restricted to two dimensions entirely, or focused on binaries aligned with the midplane of the disk. Herein we present the first three-dimensional, high-resolution, local-shearing-box, inviscid hydrodynamical simulations of disk-embedded binaries over a range of orbital inclinations. We find that retrograde binaries can shrink up to 4 times as quickly as prograde binaries, and that all binaries not perfectly aligned (or anti-aligned) with the AGN disk are driven into alignment. An important consequence of this is that initially retrograde binaries will traverse the inclinations where von Zeipel–Lidov–Kozai oscillations can drive binary eccentricities to large values, potentially facilitating mergers. We also find that interactions with the AGN disk may excite eccentricities in retrograde binaries and cause the orbits of embedded binaries to precess.


INTRODUCTION
The accretion disks around supermassive black holes (SMBH), which power active galactic nuclei (AGN) (e.g.Lynden-Bell 1969;Ho 2008), are a potential host site for merging stellar-mass black holes.Stellar mass black holes embedded in these disks will grow via gas accretion, and potentially through repeated mergers facilitated by the deep potential well of the SMBH (e.g.McKernan et al. 2012;Gerosa & Fishbach 2021).Thus, AGN disks are readily able to produce black holes within the pair instability mass gap, such as those inferred to have caused GW190521 (85 +21 −14 M ⊙ and 66 +17 −18 M ⊙ at 90% confidence, Abbott et al. (2020); see however Nitz & Capano (2021)).Additionally, the presence of gas introduces the possibility of electromagnetic counterparts to stellar-mass black hole binary mergers (e.g.Stone et al. 2017;McKernan et al. 2019;Tagawa et al. 2023).Flaring AGN have been associated with the localization error volumes of not only GW190521 (Graham et al. 2020), but also eight other gravitational wave events (Graham et al. 2023) during the third observing run of Advanced LIGO and Advanced VIRGO (The LIGO Scientific Collaboration et al. 2021), although these associations may not be statistically robust (e.g.Nitz & Capano 2021;Palmese et al. 2021).
Numerous investigations have estimated the rates of black hole binary mergers within AGN disks using semianalytical methods or one-dimensional simulations (e.g.Stone et al. 2017;Bartos et al. 2017;Tagawa et al. 2020;McKernan et al. 2020b,a), and it has been estimated that up to ∼ 80% of the black hole binary mergers de-tected by LIGO/VIRGO may occur within AGN disks (Ford & McKernan 2022).
The progenitors of these GW events -stars and compact objects -may be captured into the disk from a nuclear star cluster through, for example, torques from the gaseous disk or dynamical friction (e.g.Syer et al. 1991;Artymowicz et al. 1993;Rauch 1995;MacLeod & Lin 2020).Additionally, gravitational instability in the outer regions of AGN disks may lead to star formation (e.g.Kolykhalov & Syunyaev 1980;Goodman 2003).Furthermore, stellar evolution simulations suggest that most stars captured into or formed within AGN disks will leave behind black holes and neutron stars at the end of their lives (Cantiello et al. 2021;Jermyn et al. 2021;Dittmann et al. 2021Dittmann et al. , 2023;;Ali-Dib & Lin 2023), and neutron stars within the disk will likely accrete from the disk to the point of collapse (e.g.Perna et al. 2021;Pan & Yang 2021).Therefore, it is plausible that AGN disks host a substantial population of stellar-mass black holes. 1 growing body of work has studied binary formation within AGN disks using multidimensional N-body simulations (e.g.Secunda et al. 2019Secunda et al. , 2020;;Li et al. 2022a) as well as hydrodynamical simulations (e.g.Li et al. 2023;Rowan et al. 2023a;Whitehead et al. 2023;Rowan et al. 2023b).An assumption that pervades the literature on black hole binary mergers within AGN disks is that gas-binary interactions help drive binaries towards inspiral.This assumption stems from early analytical (e.g.Pringle 1991;Artymowicz et al. 1991) and numerical (e.g.MacFadyen & Milosavljević 2008) work on isolated binary-disk interaction.However, subsequent higherresolution and longer-duration simulations of circumbinary accretion have demonstrated that binary inspirals are far from universal (e.g.Moody et al. 2019;Muñoz et al. 2019Muñoz et al. , 2020;;Tiede et al. 2020;Dittmann & Ryan 2021, 2022;Wang et al. 2023, see also the recent review Lai & Muñoz 2023).
Studies of viscous circumbinary disks inclined relative to their binaries have found that the binary and disk are typically brought into alignment in the circular case (e.g.Lubow et al. 2015;Moody et al. 2019), but also that disks may equilibriate to a polar configuration about eccentric binaries (e.g.Martin & Lubow 2017).
Compared to the study of isolated binaries, the study of embedded binaries is nascent.Early studies pursued global two-dimensional simulations, identifying the importance of resolution (c.f.Baruteau et al. 2011;Li et al. 2021), binary orientation (Li et al. 2021), and disk thermodynamics (Li et al. 2022b).More recently, twodimensional shearing-sheet simulations have surveyed binary eccentricity and mass ratio (Li & Lai 2022a), typically finding that binaries shrink and eccentricities are damped.Recent two-dimensional studies have also surveyed the gas equation of state and binary separation (Li & Lai 2022b).
Heretofore the only three-dimensional simulations of binary stellar-mass black holes embedded in AGN disks have been those presented in Dempsey et al. (2022).Using isothermal shearing-box simulations that directly evolve the three-body orbital evolution of the system, that work illustrated the dependence of binary orbital evolution on the binary separation; balance between gravitational forces and gas pressure; and dimensionality of the problem, finding that three dimensional flows dramatically differ in structure compared to those found in two-dimensional simulations.
All previous simulations of binaries embedded in AGN disks have focused on binary orbits aligned or antialigned with the midplane of the SMBH disk, while in reality binaries may form with a wide range of inclinations and eccentricities.In particular, the twodimensional simulations of gas-assisted binary formation by Li et al. (2023) suggest that newly-formed binaries should form with high eccentricities, predominantly in retrograde configurations.Thus motivated, we focus in the present study on the orbital evolution of inclined embedded binaries, and will present the evolution of eccentric binaries in a companion paper, in both cases investigating both prograde and retrograde binaries.
We review AGN disk models and the orbits of embedded binaries in Section 2, and describe our simulation methodology in Section 3. We present the results of our simulations in Section 4, including that retrograde binaries can inspiral more than four times faster than prograde binaries and binary inclinations are gradually driven to zero.We discuss some shortcomings of the present study, as well as its implications for mergers in AGN disks and various dynamical phenomena, in Section 5. We conclude in Section 6, and provide additional technical details in Appendix A. We present some convergence tests in Appendix B and tabulate some of our results in Appendix C.

Objects embedded in AGN disks
The accretion disks that fuel AGN are typically modeled by finding solutions to the equations of viscous hydrodynamics, often under the assumptions of azimuthal symmetry, a steady state, and the limit of a thin disk such that the disk pressure scale height H at a given radial distance r • is H ≪ r • (e.g.Shakura & Sunyaev 1973;Novikov & Thorne 1973;Frank et al. 2002;Sirko & Goodman 2003;Thompson et al. 2005), although models which relax some of these assumptions have also been developed (e.g.Abramowicz et al. 1988;Sądowski 2009;Gilbaum & Stone 2022).A common feature of many disk models is that at large distances from the center the free-fall timescale of the gas disk becomes shorter than the dynamical timescale leading to gravitational instability (e.g.Toomre 1964) and the formation of stars (e.g.Kolykhalov & Syunyaev 1980;Derdzinski & Mayer 2023).Some disk models have accounted for this by imposing marginal stability arbitrarily (e.g.Sirko & Goodman 2003) or by invoking feedback from stars or black holes embedded in the disk (e.g.Thompson et al. 2005;Dittmann & Miller 2020;Gilbaum & Stone 2022).
We provide examples of several AGN disk models from the literature in Figure 1 (see also Figure 1 of Dempsey et al. 2022), specifically their aspect ratios H/r • , surface densities (Σ), the ratio of local disk mass to object mass, ratio of the Hill radius to disk scale height, and characteristic timescales for stellar mass objects to change their orbits.
Restricting ourselves to thin disks (H ≪ r • ) which are vertically isothermal, the density profile as a function of height above the disk midplane (z) at a given radius is where ρ 0 is the density in the disk midplane and H is the pressure scale height of the disk.Furthermore, H = c s /Ω, where c s is the sound speed within the disk at a given radius, Ω = GM • /r 3 • is the angular velocity of the disk, and M • is the SMBH mass.Considering objects 2 of mass m ≪ M • embedded within the AGN disk, the characteristic length scale within which gas is bound to the embedded object rather than the SMBH is the Hill radius (Hill 1878), The local disk mass (M d ) within the Hill sphere of an object is given, assuming only vertical density variation, 2 Within this subsection the difference between individual objects and tightly-bound binaries is immaterial so we use 'object' in place of either.by The ratio R H /H is a measure of the strength of gravitational forces relative to gas pressure, as well as the degree to which the disk is vertically stratified within the vicinity of the embedded object.In the limit R H ≪ H, the gas surrounding an embedded object has a near-constant density and the local disk mass is M d ≈ 4πρ 0 R 3 H /3. In the limit R H ≫ H, the gas about an embedded object is strongly stratified and the lo- As illustrated in Figure 1, objects embedded in AGN disks can occupy both extremes.Previous studies have demonstrated that R H /H significantly influences the orbital evolution of binaries embedded within accretion disks, with moreembedded binaries inspiraling more quickly at a given orbital separation (e.g.Dempsey et al. 2022).
The ratio M d /m gauges whether the self-gravity of the gas disk is dynamically relevant to the orbital evolution of embedded objects.Throughout most of the AGN disk models shown in Figure 1 the self-gravity of the AGN disk should have a negligible effect on the dynamics of embedded objects, although it may be marginally important in the outer regions of the disks modeled in Dittmann & Miller (2020).Therefore, we are somewhat justified neglecting in our simulations the disk self-gravity and the back-reaction of gas upon the black holes.
The orbit of an embedded object about the SMBH can be characterized by a semi-major axis a, an eccentricity e, and an inclination angle relative to the disk midplane i.As embedded objects orbit within the disk, they excite waves that can transfer orbital energy and angular momentum between the disk and object.Because the mass of an embedded object is typically much smaller than the mass of the corresponding SMBH, we approximate the timescales to change the orbital elements τ a for the semi-major axis and τ w for the eccentricity and inclination) using the classical Type I values (Goldreich & Tremaine 1980;Tanaka et al. 2002;Tanaka & Ward 2004) where q ≡ m/M • and τ orb ≡ Ω −1 is the local orbital timescale.Because astrophysical disks generally have H/r • < 1, τ w will generally be shorter than τ a .
The ratio between τ a and τ orb measures the rate of change of the semi-major axis of an embedded object.If τ a /τ orb ≫ 1, then the physical properties of the AGN disk are roughly constant on the orbital timescale of the embedded object.In this case the local AGN disk properties in the vicinity of an embedded objects are approximately constant, whereas if τ a ∼ τ orb the disk properties change relatively rapidly.In our simulations of embedded binaries, we assume that the properties of the AGN disk are constant in time, which is justified when τ a /τ orb ≫ 1. Considering the disk models shown in Figure 1, the steady-state assumption in our simulations is reasonable in all AGN disk regions save for the outer regions of the disk models calculated in Dittmann & Miller (2020).
Objects formed within the disk will typically have inclinations i ≲ H/r • , while objects captured into the disk may have larger initial inclinations (i > H/r • ) that are eventually damped to i ≲ H/r • .Additionally, gravitational scattering events between embedded objects may excite orbital eccentricities and inclinations over time (e.g.Stewart & Ida 2000).For example, Stone et al. (2017) estimated that over 10 8 years, typical inclinations exceeded H/r by factors of order unity far from the SMBH and more than an order of magnitude closer in.Assuming that binaries form within the disk with orbital separations a b ∼ R H (e.g., Li et al. 2023), in regions of the disk where R H /H ≲ 1 and τ w /τ orb ≳ 1 binaries may form with virtually any orientation relative to the disk midplane.On the other hand, where τ w /τ orb ≲ 1, binaries will likely form in the disk midplane, either aligned or anti-aligned with the orbit of the binary about the SMBH, although the linear analysis performed in Tanaka & Ward (2004) breaks down in the τ w /τ orb ≲ 1 limit.In regions with τ w /τ orb > 1, binary-single interactions may excite inclination even in binaries which are initially coplanar with the disk midplane (e.g.Tagawa et al. 2020).

Binary Orbits in AGN disks
In this section, we introduce basic definitions related to black hole binaries (BHBs) in AGN disks, drawing from Chapter 6 of Danby (1988).Subsequently, we briefly review relevant dynamical phenomena.

Definitions
The BHB is characterized by both its specific energy and its specific angular momentum vector, where m b = m 1 + m 2 is the binary mass, r = r 2 − r 1 is the relative position vector and v = v 2 − v 1 is the relative velocity vector.From E and h ≡ √ h • h, we define the orbital semi-major axis, a b , and eccentricity, e b , using the well-known relations In addition to the scalar eccentricity, we also define the eccentricity vector which, in the orbital plane, points toward the periapse, x y A schematic of a binary orbit oriented in threedimensional space.The reference plane is shown in gray, and the portion of the orbit below the plane is plotted using a dashed line.The ascending node, the point at which r passes from below to above the reference plane, is marked as k.
Orienting an inclined and eccentric BHB in space requires three additional angles (i b , ω b , Ω b ) which are respectively the orbit inclination, argument of periapse, and longitude of the ascending node.For reference, we illustrate this connection in Figure 2. From the schematic, the orientation angles are related to h and e via, where because ω b and Ω b are each confined to a plane they can be oriented without ambiguity, i b simply takes values between 0 and π,3 and we have defined the magnitude of the in-midplane component of the binary angular momentum vector to be h ⊥ ≡ h 2 x + h 2 y .We note that the role of Ω b in the expressions for ω b is simply a rotation in the x − y plane, and that Ω b is undefined when h ⊥ = 0, in which case only Ω b + ω b is a meaningful quantity and we set Ω b = 0.
Gravitational interaction with the surrounding gas changes the binary's total energy, E, and angular momentum, J.These are related to the specific values through the reduced mass The power delivered to the binary, Ė, and the total torque on the binary, J, together with the mass accretion rates ṁ1,2 , induce changes in a b and e b at the rates, where J = J • J/J.In addition to ȧb and ėb , the orientation angles change with the rates In Appendix A we provide further details on how we measure changes in the binary energy, angular momentum, and other quantities; how we derive average rates of change of binary orbital elements using these measurements; and list additional mathematical relations omitted here for brevity.In Section 4, we will make use of a coordinate system aligned with the BHB made up of the triad h, r, and h × r.

Dynamics
We briefly draw attention to a few flavors of dynamical interactions between binaries, such as those described above, and a distant, perturbing third body.First, binaries with large separations ( a b ≳ 0.4 R H ) are known to evolve chaotically; tidal forces from the SMBH are likely to disrupt such loosely bound binaries (e.g.Eggleton & Kiseleva 1995;Mardling & Aarseth 2001).Another process that may operate, notable for functioning even in purely coplanar triple systems, is the evection resonance (Touma & Wisdom 1994).This process relies on the orbital angular frequency of the outer black hole binary about the SMBH (n • ) being commensurable to the rate of change of the longitude of periapsis n • ≈ πb , where ϖ b = ω b + Ω b .The resonance may be able to drive black hole binaries migrating inward through AGN disks towards high eccentricities (and thus faster gravitational wave-driven mergers, e.g.Peters 1964), although this mechanism is thought to preferentially apply to binaries in AGN disks composed of an intermediate-mass black hole and a stellar-mass black hole (Muñoz et al. 2022;Bhaskar et al. 2022).If binaries acquire appreciable eccentricities through the evection resonance or other processes, their inclination may also be excited through the eviction resonance (Touma & Wisdom 1994) at the expense of the eccentricity of the binary.

NUMERICAL METHODS
We simulate binaries embedded in AGN disks using the Athena++ code (Stone et al. 2020), generally following the procedure described in Dempsey et al. (2022).We assume that the gas is isothermal and inviscid, and use the shearing-box approximation (Hawley et al. 1995;Stone et al. 1996), which expands the equations of hydrodynamics about a reference point in the disk, in our case the BHB center of mass.We illustrate this methodology schematically in Figure 3, which highlights both the quasi-global features captured by our simulations, such as the horseshoe orbits followed by certain fluid elements and the spiral arms excited in the AGN disk, and the high fidelity with which the flow of gas around the binary is captured in three dimensions.
The equations of compressible, isothermal hydrodynamics in the shearing box frame, along with source terms describing the interaction between the fluid and the BHB, are where Φ is the gravitational potential of the BHB, S Σ is a mass sink term, and S p is a momentum sink term.Each sink term S Σ = i S Σ,i , S p = i S p,i is localized to the vicinity of the BHB point masses.As illustrated in Figure 3, we place the shearing-box at some orbital radius, R 0 , in the AGN disk.The coordinate basis vectors in this frame rotate in time and are x, which points from the SMBH to the center of the shearing-box; ẑ, the direction normal to AGN disk midplane; and ŷ which is defined as ŷ = ẑ × x to complete the orthonormal triad.The rotation rate of the box is set to the Keplerian value, Ω 0 = GM • /R 3 0 , and the binary center of mass is initially placed at the shearing box center, but we note that the binary center of mass is not fixed to this location.In the remainder of the paper we use Ω 0 and n • interchangeably since n • does not change ap-preciably over the course of a simulation.We close the equations of motion with an isothermal equation of state such that P = c 2 s ρ, where the sound-speed is c s = H 0 Ω 0 , and H 0 is a scale height which we take to be constant throughout the domain.The Ω 2 0 zẑ term represents the vertical gravitational acceleration in the thin disk due to the SMBH.As discussed in Section 2, we neglect both the self-gravity of the disk and the force of the disk on the black holes.
We model the gravitational potential of each black hole using a spline function which is exactly Keplerian outside of a softening radius r s , within which the potential is softened to avoid singularities (Springel et al. 2001).We set the softening length to r s = 0.04a b . 4We apply torque-free sink terms S Σ and S p (Dempsey et al. 2020;Dittmann & Ryan 2021) for gas within r s from either point mass as described in Appendix A of Dempsey et al. (2022).We gradually introduce the black hole binary to the system over a timescale of 0.5n −1 • , only allowing accretion after the conclusion of that period.We orient the binary as described in Section 2.2.1, using the simulation x-axis as a reference direction and x−y plane as a reference plane for defining Ω b , i b , and ω b .
Our simulation domain extends from −24H to 24H in x and y, and from −4H to 4H in z, resolving the scale height with five cells in each direction.Additionally, we utilize seven levels of static mesh refinement, allowing Athena++ to automatically place additional refinement regions in order to maintain that neighboring cells differ by at most one refinement level.Accordingly, we resolve the gas around the binary with more than 1200 cells per a b .Our simulations used 2nd-order spatial reconstruction (van Leer 1974), the HLLE approximate Riemann solver (Harten et al. 1983;Einfeldt 1988), and a second-order strong stability preserving Runge-Kutta time integrator (RK2 Gottlieb & Shu 1998).We found that the Runge-Kutta integrator leads to a much more stable time step than the 2nd-order van Leer integrator (Stone & Gardiner 2009)  Our simulations were initialized according to the background equilibrium solution, with v y = −3Ω 0 x/2, v x = v z = 0, and ρ = ρ 0 exp [−(z/H) 2 /2].Because we neglect the self-gravity of the disk and disk-induced binary motion, our simulations are scale-free and the value of ρ 0 is set to unity without loss of generality.
We used shear-periodic boundary conditions in the x−direction and periodic boundary conditions in the y−direction.In the z−direction we applied standard outflow boundary conditions to the velocities, but extrapolated the density so that vertical hydrostatic equilibrium could be maintained.

N-body evolution
Previous studies of embedded BHBs either use fixed binary orbits with ad hoc precession terms (e.g.Li & Lai 2022a) or solved the 3-body equations of motion under simplifying assumptions (Dempsey et al. 2022;Whitehead et al. 2023).We found that the approximate N-body evolution scheme employed in Dempsey et al. (2022), while sufficient for circular i b = 0 binaries, did not properly conserve angular momentum for inclined binaries and resulted in ZKL cycles of erroneous magnitude, motivating the present treatment.The difficulty lies in that fact that even for mildly hierarchical triples (e.g., a b ∼ R H /4) such as the ones presented in Section 4, the orbital evolution of the binary is non-trival.We illustrate this by the evolution of the orbital elements of an i b = 165 • binary in Figure 4.
To ensure that we evolve the BHB accurately, we have developed an N-body extension to Athena++ that fully couples the evolution and interaction of N gravitating and accreting bodies to the hydrodynamics solver.Integration of the gravitational system is done with a full REBOUND (Rein & Liu 2012) simulation embedded in Athena++ and uses the IAS15 integrator, an implementation of a 15th-order Radau quadrature (Radau 1880;Rein & Spiegel 2015).During one Athena++ time step, REBOUND evolves the positions and velocities of all objects in the simulation to the appropriate end time of each stage of the integrator.For example, for the RK2 integrator used here, the two-stage update equations for the schematic system dU/dt = F (U ) are such that the final updated values are, and where U (0) , U (1) , and , (2) are the values of U at the beginning of the step, at the end of the first stage, and at the end of the second stage, respectively.For our simulations U is the vector of conserved variables, as well as the position and velocities of all gravitating bodies in the REBOUND simulation.
During stage 1, we first measure the force between each body and the gas for the positions and velocities at the beginning of the step (i.e., evaluating F (U (0) ) and then evolve the N-body system to the end of the step.During stage 2, we measure the force between each BH and the gas for the positions and velocities at from the end of stage 1 (i.e, evaluating F (U (1) )).For each stage we add the force on each BH to a running time average with the appropriate weights.With this method we are able to achieve 2nd order in time coupling between rebound and Athena++ for the RK2 integrator.We have also tested our method on the RK1, VL2, and RK3 time integrators available in Athena++.The REBOUND simulation is always performed in the "global" frame of the AGN disk, but the Athena++ simulation is performed in the "local" shearing-box frame.We thus need to rotate the REBOUND results into the local frame when computing the gas forces.
With REBOUND taking care of the numerical integration of the N-body system between each hydrodynamical time step we are able to accurately model many systems of interest for embedded Black holes in AGN disks that were previously too numerically difficult to evolve.Examples include highly eccentric and inclined binaries, small clusters of objects, arbitrarily complicated inner and outer binary configurations, and binary formation.Our sub-stage coupling method also allows for accurate book-keeping of the forces on the BHs allowing for highfidelity time averages of orbital evolution quantities.1952), over the course of our simulations.Retrograde (i b > 90 • ) binaries accrete at higher rates than prograde (i b < 90 • ) binaries.We use shades of red and orange to plot the accretion rate for prograde binaries, while we plot the accretion rate for retrograde binaries using shades of blue, in each case using lighter shades for binaries closer to 90 • .

RESULTS
We present six simulations over a range of binary inclinations, • , or ∼ 3 orbits of the binary around the SMBH, which was sufficient for the accretion rates onto the binary to reach a quasi-steady state, as illustrated by Figure 5.We perform all of our measurements of orbital evolution at times later than 14n −1 • , well after initial transients have died out.We focus first on somewhat large scales (∼ R H ) in Section 4.1, at which the influence of the binary begins causing deviations in the flow compared to a single black hole of the same mass.We investigate the flow of gas through the accretion torus (minidisk or circum-single disk, 'CSD') surrounding each black hole in Section 4.2, and binary orbital evolution in Section 4.3.

Large-Scale Structures
On scales much larger than the binary semi-major axis, the presence of the binary has a largely negligible affect on the flow compared to a single black hole of the same mass.For example, the fluid density slices in Figures 3 and 6 demonstrate the formation of a spiral arm in the larger disk.The streamlines in the aforementioned figures illustrate the horseshoe orbits followed by fluid elements that are deflected away from the binary, and the separatrices in the midplane velocity field where fluid is captured by the binary.
However, because the binary separation is only a quarter the size of the Hill radius, the binary has a nonnegligible influence on the gas within the Hill sphere, as illustrated for binaries at four different inclinations in Figure 6.Turning first to the density distribution in the disk midplane, it is clear that the binary creates m = 2 spirals at larger scales -while for prograde binaries these extend down and into the CSDs, these spirals are not as well defined within r ∼ a b about retrograde binaries.Contrasting the simulations of aligned (i b = 0 • ) and anti-aligned (i b = 180 • ) binaries, while the gas just within the Hill sphere is prograde in both cases, the gas in the minidisks of the aligned binary orbits in a prograde sense while the minidisks of the anti-aligned binary orbit in a retrograde sense.The minidisks around retrograde binaries are significantly smaller than those of their prograde counterparts, largely due to the intense ram pressure experienced by gas orbiting along with the binary against gas falling in from larger scales in a prograde sense, which we investigate further in Section 4.2.
Turning to the slices through the x−z and y−z planes displayed in the second and third columns of Figure 6, we observe large-scale circulation, especially along the axis connecting the binary and SMBH.Although the precise dynamics depend on the binary inclination, matter more than ∼ 2a b from the binary tends to be pushed away from the binary and subsequently lifted upwards to |z| ≳ a b , joining the meridionally circulating flow which accretes onto the binary from higher latitudes.Much of the accretion onto the binary occurs along the vertical direction, necessitating three-dimensional study of these accretion flows.The circulation pattern observed here is qualitatively similar to that observed in some threedimensional simulations of planets embedded in circumstellar disks (e.g.Szulágyi et al. 2014;Fung & Chiang 2016), as well as the previous simulations of aligned embedded binaries (Dempsey et al. 2022).
One peculiar feature of retrograde binaries is their higher accretion rates relative to their prograde counterparts, as illustrated in Figure 5: prograde binaries typically accrete at about 7.2% of the Bondi rate, whereas retrograde binaries typically accrete at about 8.6% of the Bondi rate.A similar trend was noted in (Li & Lai 2022a), which attributed the increased accretion rate onto retrograde binaries to those binaries lacking CSDs, suggesting that the presence of minidisks around prograde binaries decreases the accretion rate.Unlike (Li & Lai 2022a), we observe CSDs around retrograde bi-naries,6 so their existence or lack thereof cannot dictate these trends in accretion rate.Moreover, our i b = 150 binary has both the smallest minidisks and lowest accretion rate out of the retrograde binaries studied here.
We have found that throughout the Hill sphere of each binary, to distances ≳ R H = 4a b , the average fluid density is lower and average radial velocity towards the binary barycenter is higher around retrograde binaries.Because gas typically enters the Hill sphere orbiting in a prograde sense, the prograde binaries will have, on average, longer-duration gravitational interactions with gas, potentially doing more work and hindering the binding of gas to the binary.Regardless, the higher accretion rates onto retrograde binaries begin at large distances from the binary and do not depend on the size of the minidisks.

Minidisks
The accretion disks that form around each black hole in our simulations display a number of notable characteristics: these mindisks exhibit meridional circulation, gas outflowing along the midplane and inflowing at higher latitudes, similar the flows on the much larger scales of the Hill sphere shown in in Figure 6.The accretion disks that form around the black holes in prograde binaries orbit in a prograde sense, whereas the accretion disks that form around retrograde binaries form in a retrograde sense (that is, with the same handedness of the binary: l b • l gas > 0); and the minidisks themselves are typically aligned with the AGN disk midplane -that is to say the angular momentum vector of gas in a minidisk around a prograde binary is nearly parallel to z, and the angular momentum vector of gas comprising the CSD of a retrograde binary is nearly parallel to −z.Particularly for inclined binaries, these CSD alignments could appreciably alter the effective spin distribution of AGN-channel black hole binary mergers (e.g McKernan & Ford 2023).As an example, a volume rendering of the minidisks of a i b = 30 • binary are displayed in Fig- ure 7. Notably, while the gas most proximate to each black hole is aligned with the AGN disk midplane, the lower-density gas at larger distances, corresponding to the spirals and outflows seen in Figure 6, is more aligned with the orbital plane of the binary.We caution that this result does not hold at all binary inclinations: preliminary simulations of more misaligned binaries, which we plan to present in a subsequent paper, suggest that the minidisks of those binaries will align more with the • .The darkest colors highlight gas at densities of ρ ∼ ρ0/10, and the brightest colors highlight gas at densities of ρ ∼ 10 4 ρ0, with intermediate colors spaced log-uniformly in density.Even though the orbit of the binary is inclined relative to the AGN disk midplane, the minidisks remain oriented along with the AGN disk rather than the plane of the binary orbit.Notably, the lower-density gas flowing away from the binary is in closer alignment with the binary orbital plane.
binary orbital plane than do the minidisks of binaries in the present investigation.
It is worth comparing the alignment of these minidisks to the disk-binary alignment observed in earlier simulations of viscous, isolated circumbinary disks (e.g.Moody et al. 2019): in such simulations the binary accretes from a circumbinary disk much thinner than the orbit of the binary, and (in a frame oriented with the circumbinary disk) out-of-plane velocities in the circumbinary disk are negligible compared to the orbital motion of the binary.However, embedded binaries accrete gas vertically, and over the course of a binary orbital period just as much matter falls on from above as below.Gas flows onto isolated binaries from their circumbinary disks, whereas gas falls onto embedded binaries ballistically (Dempsey et al. 2022).Additionally, Avara et al. (2023) recently identified minidisk tilt oscillations on the order of the disk aspect ratio.
The gas constituting the minidisks of embedded binaries is also not particularly long-lived.While the typical viscous inflow timescale of material in the minidisks of isolated binaries can be tens to hundreds of binary orbital periods, 7 we find that the mass-weighted average inflow timescales within the Roche lobe of each black 7 A simplistic estimate follows from treating each CSD as a circular Keplerian disk with a radius matching the size of the Roche lobe of each black hole, in which case r/vr ≈ 91n −1 • ≈ 15t orb for an unrealistically high viscosity of ν = 10 −3 a 2 b Ω b and equal-mass binary.
hole are on the order of ∼ 0.2 − 4 binary orbital periods t orb in the minidisks themselves and ≲ 0.1t orb in the polar regions around each black hole.We note, however, that for gas to form a CSD around a black hole, it must necessarily be moving with a similar bulk velocity to that black hole, lest it be left behind or captured by the other black hole; thus the gas forming each CSD has significantly torqued the black hole around which it orbits, even if the orbital plane of that gas does not align with the binary orbital plane.
The flow of gas onto and throughout the minidisks over a range of inclinations is illustrated in Figure 8, which plots the azimuthally averaged gas density (in the frame of one of the black holes) along with velocity streamlines in the same frame.Although quasiballistic accretion occurs onto prograde binaries most prominently at higher latitudes, such accretion occurs over a wider range of angles in retrograde binaries until reaching the minidisk.Intriguingly, the CSDs surrounding the black holes comprising the i b = 150 • binary are visibly smaller than those around the i b = 180 • binary, which are in turn smaller than those of the black holes in prograde binaries.In fact, the CSD of the i b = 150 • binary is only about twice as large as the size of our sinks.However, the CSDs themselves display meridional circulation, with high-latitude inflows and midplane outflows, regardless of their size.
Another view of the flow of gas throughout and around the minidisks is provided in Figure 9, which plots the time-and azimuth-averaged profiles of the gas specific angular momentum and radial mass flux.Patterns of meridional circulation and midplane outflows are again evident in the distribution of radial mass flux, which also highlights accretion through the polar regions.Even though the CSDs around the black holes in the i b = 150 • binary are very small (≲ a b /10 in extent), they still display the same circulation patterns.Although the specific angular momentum profile is not particularly informative for prograde binaries, around the black hole in retrograde binaries it very clearly demarcates which gas is flowing in a retrograde sense (and thus with the same handedness as the orbit of the binary) compared to gas orbiting in a prograde sense.
A slightly more condensed view of the velocity profile throughout the Roche lobe of each black hole is shown in Figure 10.The top panel therein illustrates the specific angular momentum profile in comparison to a reference Keplerian curve.Prograde binaries have sub-Keplerian disks at r ≳ a b /10, and the time-averaged flow does not depend substantially on the binary inclination.The bottom panel illustrates the azimuthal and radial Mach numbers (M ϕ and M r respectively) of Figure 10.One-dimensional profiles, averaged over spherical shells centered on each black hole, of the fluid specific angular momentum (top panel) and both radial and azimuthal fluid Mach numbers (bottom panel).Warm colors display results for prograde binaries and cool colors display results for retrograde binaries.In the top panel, the dashed black like displays the specific angular momentum profile of a Keplerian disk.In the bottom panel the radial Mach number is plotted using dashed lines and the azimuthal Mach number is plotted using solid lines.We have calculated the characteristic Mach number at a given radius by integrating the radial and azimuthal momentum within each shell, dividing by mass in that shell and the sound speed, and taking the absolute value of the result.The vertical dotted line indicates the radius within which we apply sink terms and gravitational softening.
the gas around each black hole.As characteristic of a disk, the azimuthal flow is supersonic throughout, varying as M ϕ ∝ r 1/2 , and in the prograde disks the radial Mach number is everywhere both below unity and orders of magnitude below the azimuthal Mach number.The Mach number profiles of CSDs around prograde binaries show very little variation with inclination.
The flow of gas around the black holes in retrograde binaries, however, is markedly different.The top panel of Figure 10 illustrates how the size of the CSD around each black hole varies with binary inclination, smallest for the i b = 150 • binary and largest for the i b = 180 • binary, and how the disks are still Keplerian at small radii despite their smaller size.The bottom panel shows how fluid motion is generally azimuthally subsonic prior to forming a disk, but plunges in a radially supersonic manner up until that point.The gas initially falls in ballistically, with M r ∝ r −1/2 , before forming an accretion disk.
The flow of gas onto, off of, and between the minidisks in each binary is displayed in Figure 11, which plots a projection of the fluid density along the binary angular • .; in each case the angular momentum vector of the binary is pointing out of the page, and in this frame the binary is moving in a counterclockwise sense.Because the CSDs are aligned with the AGN disk midplane rather than the orbital plane of the binary, the minidisks of inclined binaries appear slightly prolate in these projections.We can clearly see the smaller minidisks around the retrograde binaries, as well as the streams of gas being stripped off of the retrograde binaries.Major ticks are placed 0.5a b apart, with minor ticks every 0.25a b .momentum vector.The CSDs around prograde binaries appear very similar to those observed in simulations of standard circumbinary disks (e.g.Ryan & MacFadyen 2017;Muñoz et al. 2019;Dittmann & Ryan 2021): spiral waves are present throughout the disks, some matter flows onto the binary from larger distances, and streams of gas flow between the two black holes.Although some very faint spirals can be seen in the CSDs of retrograde binaries, they are far less prominent (cf.Li et al. 2021).As presaged by Figures 9 and 10, ample gas appears to be stripped away from the retrograde binaries, con-tributing to their smaller sizes; this trailing gas appears to be more substantial for the binaries that orbit outside of the AGN disk midplane, in line with the observed trends in CSD size.

Orbital evolution
Based on our examinations of the flow of gas onto the binary and around the black holes constituting it, we can draw a number of qualitative inferences.As a first approximation, we can assume that gas accreting onto each black hole does so on average with the same velocity as that black hole.Thus, the accretion of gas onto the binary should add to its specific energy largely through changing the mass of the binary and thus its specific gravitational binding energy, and for circular binaries with |r| ≈ a b , Ė/E ≈ 2 ṁb /m b .However, gas must journey deep into the potential well of an individual black hole before accreting, potentially loosing energy to the binary in the process and helping unbind it.For more tightly bound binaries, more work must be done on gas before it can accrete, and for sufficiently bound binaries this may cause binaries to outspiral as found in Dempsey et al. (2022).
The above reasoning is directly applicable to prograde binaries, but must be adapted slightly for retrograde binaries.For example, if CSDs do not form at all (e.g.Li & Lai 2022a), then the assumption that gas accretes with the same velocity as a given black hole breaks down, and the kinetic contribution to the time derivative of the binary specific energy will become nonzero.On the other hand, if CSDs form around each black hole, as in our simulations, then the velocity of the gas must be significantly altered before it can accrete, as shown in Section 4.2, and thus the binary will do a significant amount of work on the gas.Similarly, as shown in Figure 9, the angular momentum of gas must change dramatically for it to become part of the minidisk, in which case the binary will experience a strong torque, greatly reducing the magnitude of its angular momentum, although it is not clear a priori if this will predominantly affect the binary eccentricity or semi-major axis.And similarly, for gas in the AGN disk which is on average aligned with the AGN disk midplane to accrete onto either black hole, that gas must gain significant out-of-plane angular momentum, and thus damp the inclination of the binary.
With the aforementioned notions in mind, we turn to our measurements of the rate of change of the specific energy and angular momentum of each binary in our simulations in Figure 12.Turning first to the energy, our prior expectations are confirmed: in all cases the contribution to Ė/E due to accretion is ∼ 2, and while the contribution due to gravity is small in the prograde case, it is quite large for retrograde binaries, suggesting that the binaries should inspiral rapidly.Also in line with our qualitative expectations, accretion does not change the specific angular momentum of any binary.8As argued above, retrograde binaries experience significant torques, especially contributing to changes in the z−component of their angular momenta.
The orbital evolution of each binary follows from the above discussions.First, in the top row of Figure 13, we see that although all binaries contract (at this separation, see Dempsey et al. (2022)), retrograde binaries do so at a much higher rate.Similarly to previous studies we find that eccentricity is excited in retrograde binaries (e. Figure 13.The rates of change (due to gravity, accretion, and in total) of the binary semi-major axis, squared eccentricity, inclination, argument of the ascending node, and argument of periapsis.All quantities were averaged in time between 14 n −1 • and 20 n −1 • .When reporting values of ib , ωb , and Ωb , we use radians rather than degrees.We only report the contributions to orbital elements due to accretion and gravitational interactions with the surrounding gas, which are independent of choices like the reference plane for defining Ω b .
grade binaries (Muñoz et al. 2019;Zrake et al. 2021).Although ib is not well-defined for i b = 0 and i b = 180, we see as expected that binary inclinations are damped.Otherwise put, this result suggests that any binary at all misaligned with the midplane of the AGN disk should be gradually aligned with the disk, in a prograde sense, over time.Thus, it appears that i b = 0 is a stable equilibrium while i b = 180 • is unstable.If the trends we observe here continue to higher-inclination systems, initially near-retrograde binaries could realign with the AGN disk within a few mass-doubling timescales.
The time derivative of ωb is not particularly welldefined for near-circular binaries (due to the terms in the denominator of Equation ( 18) that go to zero as e b → 0), so we have normalized those results by the time-averaged values of e b ; for these effectively circu-lar binaries, we simply emphasize that ωb is generally nonzero.The time derivative of Ω b is more well defined, at least for binaries with h ⊥ ̸ = 0, and we find that it is negative for prograde binaries and positive for retrograde binaries.These signs relate to the handedness of the binaries passage through the disk midplane and the background shear of the AGN disk, and as we discuss in Section 5.2, the important conclusion is that this diskinduced precession rate is nonzero.

DISCUSSION
In this section, we briefly discuss the implications of our results on the manner in which accretion onto embedded black holes proceeds, the dynamics of jets which may be formed during the accretion process, and the orbital evolution of embedded BHBs.We caution that the present work has eschewed magnetic fields, the effects of radiation, gas self-gravity, and general deviations from isothermality, and thus some caution must be taken when applying our results to real AGN disks.

Accretion
Although the gravitational radii of embedded black holes were under-resolved by orders of magnitude in our simulations, our study achieved sufficient scale separation to observe that accretion onto each black hole does not occur in a spherically symmetric manner, but rather through an accretion disk.Additionally, accretion rates onto the embedded binaries in our simulations are about ∼ 8% of the Bondi accretion rate ṀB = 4πρ 0 G 2 m 2 b c −3 s , which is reasonable given that our simulations are vertically stratified rather than uniform in density, and that the binary Hill radius is comparable to the Bondi radius R B = Gm b /c 2 s such that tidal forces may limit accretion onto the binary.This is largely consistent with the results of previous three-dimensional simulations of accretion onto embedded binaries (Dempsey et al. 2022), but as noted in Section 4.1, retrograde binaries accrete at a slightly higher rate than prograde ones due to their disparate gravitational influence on gas entering the Hill sphere.
However, for typical gas densities in AGN disks such as those modeled in Figure 1, ρ 0 ∼ 10 −18 −10 −9 g cm −3 , the Bondi accretion rate may be substantially larger than the Eddington accretion rate.Fixing R H /H, as in our simulations, the Bondi accretion rate can be expressed as We can compare this to the Eddington-limited accretion rate onto a binary of the same mass, (26) where ϵ is the fraction of rest mass radiated during accretion onto the black holes, c is the speed of light, and κ is the opacity of the ambient gas.For brevity, we assume that κ = 0.4 cm 2 g −1 , the value for fully ionized hydrogen gas.Then, the ratio of the Bondi accretion rate to the Eddington accretion rate is Thus, for a wide range of parameters, the accretion rates onto the binaries in our simulations, at ∼ 8% of the Bondi rate, are highly super-Eddington.If the accretion rate onto each black hole is radially constant, the accretion process onto these embedded black holes may lead to significant feedback to larger scales in the form of jets and winds.However, such feedback would be driven on scales many orders of magnitude smaller than the smallest resolved in our simulations, and will be the subject of future work.Furthermore, thermodynamic effects neglected in our isothermal simulations may suppress the accretion rate (possibly as Ṁ ∝ r 1/2 , see Guo et al. 2023), resulting in a much smaller accretion rate near the event horizon.
If either embedded black hole has an appreciable spin, jets will likely be launched into the AGN disk.Because jet efficiency strongly depends on the angular momentum of the gas accreting onto a given black hole (e.g.Kwan et al. 2023), it seems probable that prograde binaries could support jets; on the other hand, because the CSD size of retrograde binaries appears to depend strongly on softening, if they exist at all (Li & Lai 2022a), retrograde binaries may not be able to launch strong jets.Although on small scales jets are typically aligned with the spin of the black hole by which they are launched, on larger scales jets can be aligned with the accretion disk surrounding them (McKinney et al. 2013;Liska et al. 2018).It remains to-be-determined whether jets launched by embedded black holes would reach scales comparable to the size of their minidisks, although if they do they may be re-collimated and oriented along the normal to the AGN disk midplane.

Orbital Evolution
The two most robust results, in terms of orbital evolution, of our study of stellar-mass black hole binaries embedded in AGN disks are that retrograde binaries inspiral more quickly and that binaries are gradually aligned AGN disk midplane.Specifically, anti-aligned binaries inspiral more than four times faster than aligned binaries.Furthermore, the rate at which binary inclinations decreases depends on how far out of the AGN disk midplane their orbital planes lie.Thus, binaries which form in a high-inclination retrograde configuration will both have their inclinations reduced and have their semimajor axes rapidly decreased through interactions with gas.Recent studies (e.g.Li et al. 2023) have suggested that a large fraction of binaries should form in retrograde configurations.
Additionally, we have observed that both ω b and Ω b precess due to interactions with the AGN disk.Although the former is not well-defined due to the low eccentricities in our simulations, and the latter is generally small, depending on the ambient gas density, the precession of these angles may affect various resonances which could operate on embedded binaries.It is wellknown that precession of the argument of periapsis can both prevent ZLK cycles (Holman et al. 1997) and decrease the maximum achievable eccentricities in cases where ZLK cycles still occur (Miller & Hamilton 2002).However, because of the small binary eccentricities in our simulations, the degree to which ZLK cycles will be affected is uncertain, although our results suggest that ZLK cycles should be less common for embedded binaries than those in vacuum.Evection resonances are more complicated still, due to their dependence on n • as well as πb .Interactions with the gas disk may either help or hinder evection resonances, although our low-i b , lowe b simulations are not sufficient to predict how evection resonances should be affected in a general sense.
The comparative rates of change of the semi-major axis and inclination aside, we expect that all binaries with inclinations i b < 180 • will be gradually driven towards i b = 0 • .Thus, our results suggest that many retrograde binaries will pass through the range of inclinations where ZLK cycles occur as their inclinations are gradually damped.Although the rate of orbital evolution depends strongly on disk properties, if it is slow compared to the timescale of ZLK cycles, the binary will undergo ZLK oscillations in inclination and eccentricity as long as they are not disrupted by disk-induced precession.Moreover, binaries near i b ∼ 90 • may approach e b ∼ 1 via ZLK oscillations, in which case the dissipation of energy and angular momentum due to gravitational waves could become extreme near pericentre leading to a large population of merging black holes with polar alignment relative to the AGN disk.

CONCLUSIONS
We have conducted a series of three-dimensional hydrodynamical simulations of stellar-mass binaries embedded in AGN disks.
Whereas previous threedimensional studies of embedded binaries focused on those with prograde orbits aligned with the midplane of the AGN disk, we have explored a range of binary inclinations including both prograde and retrograde configurations.Our work has uncovered qualitatively new aspects to both the accretion flows around and orbital evolution of inclined embedded binaries.However, we caution that our simulations have used a simplified isothermal equation of state, neglected magnetic fields, and focused on nearly-circular binaries misaligned by not more than 30 • from the AGN disk midplane.
Binaries with higher inclinations tend to accrete at higher rates, in all cases accreting at less than the Bondi rate.Retrograde binaries have substantially smaller minidisks, although the black holes constituting the i b = 180 • binary possessed larger minidisks than those constituting the i b = 150 • or i b = 165 • binaries, due at least in part to ram pressure stripping.We find that the minidisks around each black hole are typically aligned with the AGN disk midplane.On the larger scales of the minidisks and AGN disk there are patterns of midplane outflows, meridional circulation, and accretion along the surface of gaseous tori.Some gas also flows directly onto each black hole along the polar directions.
Retrograde binaries inspiral up to four times as quickly as their prograde counterparts due to gravitational interactions with the disk.Additionally, we find that binaries with orbital planes outside of the AGN disk midplane experience inclination damping, which is small for binaries only slightly misaligned with the disk midplane, but grows larger for more misaligned binaries.We have also found that eccentricities can be excited in nearcircular inclined binaries, although we cannot comment on the maximum eccentricities achievable through diskbinary interactions due to the limited range of orbital parameters surveyed herein.Similarly, we have found evidence for disk-driven precession of both the binary argument of periapsis and longitude of the ascending node, which may affect various the susceptibility of embedded binaries to various resonances.Because a significant fraction of binaries in AGN disks may form with retrograde configurations, our results paint an optimistic, albeit incomplete, picture of gas interactions assisting binary inspirals down to to scales where gravitational waves can drive their merger, resulting in gravitational wave signals such as those observed by LIGO/VIRGO.SOFTWARE matplotlib (Hunter 2007), cmocean (Thyng et al. 2016) • ., in analogy to Figure 11.However, each point mass in the simulations pictured here has a gravitational softening length (and sink length) of 0.12a b , three times larger than the simulations presented in the main text.Many features on larger scales are completely consistent between these simulations, such as the streams of gas being stripped off of the retrograde binaries and gas flowing between the members of the prograde binaries.However, in these simulations the minidisks around the black holes in retrograde binaries are appreciably larger than in simulations using smaller softening lengths.Major ticks are placed 0.5a b apart, with minor ticks every 0.25a b .
To illustrate how CSD size depends on the softening length, we plot in Figure 14 projections of the gas density along the binary angular momentum vector analogous to Figure 11, except for simulations using sink and softening lengths of 0.12a b , the same settings as Dempsey et al. (2022).The minidisks are virtually unchanged for prograde binaries.However, the CSDs of retrograde binaries are visibly larger in Figure 14  Although the sizes of the minidisks appear to differ based on the gravitational softening length, we have found that this results in only very minor changes in binary orbital evolution.To illustrate this, we have plotted the time-averaged rates of change of the binary orbital elements from a suite of simulations using a softening length of 0.08a b , twice that used in the simulations comprising the main body of our study.The rates of change of the binary semi-major axis, eccentricity, and inclination are virtually identical to the results from smaller-softening simulations.However, we find that the time-averaged values of Ωb and ωb differ by factors of order unity between the 0.04a b -and 0.08a b -softening simu-   et al. (2022).However, Dempsey et al. (2022) used not only a softening length of 0.12a b in comparison to the default value of 0.04a b in this work, but also a different method of solving the equations of motion of the binary.Thus, we have repeated our i b = 0 simulation using softening lengths of 0.12a b , 0.08a b , 0.04a b , and 0.02a b , using the same value for the sink length in each case.All of these simulations used our standard resolution settings except for the smallest-softening simulations, for which we added an additional level of refinement around the binary to achieve a resolution of ∼ 2400 cells per a b .Because of the great cost associated with the highest-resolution simulation, it was only carried out until t = 5.5n −1 • .Although not enough to reach a converged state or achieve robust statistics, this was sufficient for the binary to reach a quasi-steady accretion rate.Because the accretion contribution to the evolution of the semi-major axis is independent of softening length, we focus on the gravitational contribution to ⟨ Ė⟩/⟨E⟩.For the aforementioned softening lengths, we found values of ⟨ Ė⟩/⟨E⟩, averaged from t = 4.5n −1 • to t = 5.5n −1 • , of ∼ 2.5, 0.78, 0.57, 0.76.Thus, it appears that our fiducial softening choice is sufficiently converged, at least for prograde binaries, although our results can be assumed uncertain at the level of a few percent.

C. QUANTITATIVE RECORDS OF ORBITAL EVOLUTION
We list in Table 1 a number of the quantities plotted in Figure 13, as well as the accretion rate onto each binary, in all cases averaged from t = 14n −1 • to t = 20n −1 • .

Figure 1 .
Figure1.Four example AGN disk models: the left and center panels display AGN disk models (those ofSirko & Goodman (2003) andThompson et al. (2005) respectively) around M• = 10 8 M⊙ SMBHs computed inBellovary et al. (2016); and the right column displays models for disks around 4 × 10 6 M⊙ SMBHs computed inDittmann & Miller (2020), based onThompson et al. (2005) but making different assumptions about feedback and gas opacity.The first two rows display the disk aspect ratio and surface density.The bottom four rows display the range of M d /m, RH /H, τa/τ orb , and τw/τ orb for objects with masses ranging from 5 M⊙ to 300 M⊙.Horizontal dashed gray lines plot the line y = 1 where applicable.

Figure 3 .
Figure3.A schematic diagram of our simulation methodology.We do not simulate the global AGN disk (left), but instead consider a shearing box which orbits the SMBH along with the binary.The upper right panel plots a density slice through the midplane during a simulation of an inclination i = 30 • binary, illustrating the quasi-global features captured by our simulation such as the excitation of spiral arms in the AGN disk and the horseshoe orbits of fluid elements which remain unbound from the binary.The bottom panel depicts opaque density iso-contours in the immediate vicinity of the binary in the same simulation, demonstrating that the flow of gas in three dimensions is resolved well by our methodology.In both plots of simulation outputs, brighter colors indicate higher densities, which are colored on a logarithmic scale.

Figure 4 .
Figure 4. Evolution of the semi-major axis (blue dashed), inclination (orange solid), and eccentricity (green dotted) of our i b = 165 • binary over the course of one orbit about the SMBH.

Figure 5 .
Figure 5.The accretion rate onto each binary, normalized to the Bondi accretion rate ṀB = 4πρ0G 2 m 2b /c 3 s(Bondi 1952), over the course of our simulations.Retrograde (i b > 90 • ) binaries accrete at higher rates than prograde (i b < 90 • ) binaries.We use shades of red and orange to plot the accretion rate for prograde binaries, while we plot the accretion rate for retrograde binaries using shades of blue, in each case using lighter shades for binaries closer to 90 • .
We focus on a single binary configuration, with mass ratio m b /M • = 1.536 × 10 −6 , separation a b = R H /4 and a disk scale height of H = 0.01R • such that R H /H = 0.8 (c.f.Dempsey et al. 2022).In each case, the mean motion of the binary was n b ∼ 13.9n • .Each binary was initialized with e b = 0, ω b = 0, and Ω b = 0.The simulations were run for at least 20n −1

Figure 6 .
Figure 6.Slices of density along each coordinate plane, along with velocity streamlines, at t = 14.5n −1 • in our i b = {0 • , 30 • , 150 • , 180 • } simulations.The first, second, and third columns plot data in the x − y, x − z, and y − z planes respectively.Curved light blue arrows illustrate streamlines of the in-plane velocity field in each slice.In all cases, it is evident that accretion along the vertical direction plays a significant role.Additionally, there is evidence for an m = 2 spiral in each case, although these spirals are weaker when i b > 0 • , and the deflection of fluid elements on approximate horseshoe trajectories is clearly visible through the streamlines in the x − y plane.The axis ticks are spaced 2a b = RH /2 apart, and each panel roughly encompasses a slice through the Hill sphere of each binary.

Figure 7 .
Figure 7.A volume rendering of the gas around the black holes in our i b = 30 • simulation at t = 14.5n −1• .The darkest colors highlight gas at densities of ρ ∼ ρ0/10, and the brightest colors highlight gas at densities of ρ ∼ 10 4 ρ0, with intermediate colors spaced log-uniformly in density.Even though the orbit of the binary is inclined relative to the AGN disk midplane, the minidisks remain oriented along with the AGN disk rather than the plane of the binary orbit.Notably, the lower-density gas flowing away from the binary is in closer alignment with the binary orbital plane.

Figure 8 .
Figure 8. Profiles, averaged in time and azimuth, of the gas density and velocity in a cylindrical polar coordinate system in the frame of one of the black holes from each of our i b ∈ {0 • , 30 • , 150 • , 180 • } simulations.Streamlines illustrate meridional circulation in the minidisks and vertical accretion onto the binary through the direction of the azimuthally-averaged velocity field, and color indicates the azimuthally-averaged fluid density.Time-averaging used 27 individual snapshots over the course of three binary orbits.The dashed white circle denotes the region where our sink particle operates and the gravitational potential is softened.

Figure 9 .
Figure 9. Azimuthally averaged profiles of radial mass flux (top row) and the specific angular momentum of the gas (bottom row) in a cylindrical polar coordinate system in the frame of one of the black holes from each of our i b ∈ {0 • , 30 • , 150 • , 180 • } simulations.Time-averaging used 27 individual snapshots over the course of three binary orbits.The dashed white circle denotes the region where our sink particle operates and the gravitational potential is softened.

Figure 11 .
Figure 11.Projections of the fluid density along the orbital angular momentum vector of each binary at t = 18.4n −1• .; in each case the angular momentum vector of the binary is pointing out of the page, and in this frame the binary is moving in a counterclockwise sense.Because the CSDs are aligned with the AGN disk midplane rather than the orbital plane of the binary, the minidisks of inclined binaries appear slightly prolate in these projections.We can clearly see the smaller minidisks around the retrograde binaries, as well as the streams of gas being stripped off of the retrograde binaries.Major ticks are placed 0.5a b apart, with minor ticks every 0.25a b .

Figure 12 .
Figure 12.The top row plots the rate of change of the binary specific energy normalized by the average binary specific energy and accretion timescale m b / ṁb for each of our simulations (the total rate shown by blue stars, the contribution due to accretion by orange crosses, and the contribution due to gravity by green circles).The second row plots the average rate of change of the magnitude of the binary specific angular momentum, the third the rate of change of the z-component, and the bottom row shows the rate of change of the perpendicular component, all normalized by the binary specific angular momentum and accretion timescale.In the first row, the gray dashed line indicates the critical value Ė/E = ṁb /m b above which equal-mass binaries are driven to inspiral.All quantities were averaged in time between 14 n −1 • g. Schnittman & Krolik 2015; Tiede & D'Orazio 2023), while eccentricity is damped in near-circular pro-

Figure 14 .
Figure 14.Projections of the fluid density along the orbital angular momentum vector of each binary at t = 18.4n −1• ., in analogy to Figure11.However, each point mass in the simulations pictured here has a gravitational softening length (and sink length) of 0.12a b , three times larger than the simulations presented in the main text.Many features on larger scales are completely consistent between these simulations, such as the streams of gas being stripped off of the retrograde binaries and gas flowing between the members of the prograde binaries.However, in these simulations the minidisks around the black holes in retrograde binaries are appreciably larger than in simulations using smaller softening lengths.Major ticks are placed 0.5a b apart, with minor ticks every 0.25a b .

Figure 15 .
Figure15.The same quantities shown in Figure13for our simulations with the same resolution but gravitational softening parameters and sink rates twice as large, in this case 0.08a b compared to the value of 0.04a b used in the body of the paper.
lations, though in all cases these fluid-induced precession rates could hinder ZLK cycles.Notably, the i b = 0 simulations suggest an appreciably slower rate of inspiral for binaries with a b /R h = 1/4 binaries than the equivalent simulations in Dempsey et al. (2022): ⟨ ȧb ⟩/⟨a b ⟩ ≈ −1.76⟨ ṁb ⟩/m b in simulations compared to ⟨ ȧb ⟩/⟨a b ⟩ ≈ −3.1⟨ ṁb ⟩/m b in Dempsey