Black Hole Binaries in AGN Accretion Discs II: Gas Effects on Black Hole Satellite Scatterings

The black hole (BH) binaries in active galactic nuclei (AGN) are expected to form mainly through scattering encounters in the ambient gaseous medium. Recent simulations, including our own, have confirmed this formation pathway is highly efficient. We perform 3D smoothed particle hydrodynamics (SPH) simulations of BH scattering encounters in AGN disks. Using a range of impact parameters, we probe the necessary conditions for binary capture and how different orbital trajectories affect the dissipative effects from the gas. We identify a single range of impact parameters, typically of width $\sim0.86-1.59$ binary Hill radii depending on AGN disk density, that reliably leads to binary formation. The periapsis of the first encounter is the primary variable that determines the outcome of the initial scattering. We find an associated power-law between the energy dissipated and the periapsis depth to be $\Delta E\propto r^{-b}$ with $b=0.42\pm0.16$, where deeper encounters dissipate more energy. Excluding accretion physics does not significantly alter these results. We identify the region of parameter space in initial energy vs impact parameter where a scattering leads to binary formation. Based on our findings, we provide a ready-to-use analytic criterion that utilises these two pre-encounter parameters to determine the outcome of an encounter, with a reliability rate of>90\%. As the criterion is based directly on our simulations, it provides a reliable and highly physically motivated criterion for predicting binary scattering outcomes which can be used in population studies of BH binaries and mergers around AGN.

The formation and evolution of black hole binaries in gaseous en-vironments is not fully understood.Simulations that consider purely gas driven evolution of a single isolated black hole binary (e.g neglecting three-body scattering, GW hardening, and the additional torque of the central massive SMBH) result in different outcomes depending on the assumed initial conditions.In the 2D simulations of Tang et al. (2017), Muñoz et al. (2019), Moody et al. (2019), Tiede et al. (2020) and 3D simulations of Moody et al. (2019), which all adopted pressure conditions for a thick disc with scale height  to radius  ratio / = 0.1, a concensus that binaries should experience net positive torques emerges, leading to an outspiral.However, the opposite is true for cooler, thinner discs (Tiede et al. 2020;Heath & Nixon 2020).It was also suggested that viscosity plays a role in cases where this transition occurs between positive and negative torques with respect to / (Dittmann & Ryan 2022).The modelling of subgrid/resolution environments around the BHs was also shown to affect the simulation outcome in some cases.The binary evolution was found to depend on the sink rate (Tang et al. 2017) and softening lengths (Li et al. 2021), although more recent work found torques to have converged with respect to most numerical choices, at least for circular binaries (e.g Moody et al. 2019;Muñoz & Lithwick 2020;Duffell et al. 2020;Westernacher-Schneider et al. 2022).
In the last few years there has been a steadily growing number of hydrodynamical simulations of binaries embedded in a SMBH accretion disc.In the first such simulations of Baruteau et al. (2011) (which consider stellar binaries), gas dynamical friction was found to harden a pre-existing binary in a 2D gas disc.This occured regardless of any gap formation in the SMBH disc, where the binary is sufficiently massive to expel gas from its orbit around the SMBH (e.g.Goldreich & Tremaine 1980) more efficiently than it can be replenished through the disc's viscosity.Li et al. (2021) corroborated the binary hardening results of Baruteau et al. (2011) in the cases where the binary system orbit is retrograde with respect to its orbit around the disc/SMBH, as positive torque sources near the BHs are damped in the retrograde configuration due to disruption of the inner regions of their circum-single mini discs (CSMDs) that gravitationally deposit energy into the binary.For prograde binaries, both papers find their orbits expand over time.Kaaz et al. (2021) considered a binary in a thick disc where the binary does not open a gap in the disc using 3D wind tunnel simulations and find binary hardening in all their models.Such rapid hardening is also found in the three paper series of wind tunnel simulations of Li & Lai (2022a,b, 2023).The source of the discrepancy between the wind tunnel simulations isn't clear, but could be the choice of equation of state (EOS), where the former assume an isothermal EOS and the latter consider an adiabatic EOS.In another paper that models the global AGN disc, Li et al. (2022c) find an enhanced temperature profile in the BH minidiscs lead to more pronounced negative torques and binaries that would normaly expand under isothermal conditions instead shrink.Extending the dimensions of the problem to 3D also leads to increased negative torques (e.g Dempsey et al. 2022), although that study found that this reverses when the separation becomes 10% of their Hill radius.Other orbital elements may affect the evolution of the binary separation, such as eccentricity (Muñoz et al. 2019;D'Orazio & Duffell 2021;Siwek et al. 2023b) and mass ratio (Duffell et al. 2020;Dempsey et al. 2021;Siwek et al. 2023a).For a review of how orbital elements can affect binary evolution, see Lai & Muñoz (2023).
While the overall evolution of embedded binaries is still a largely unconsolidated problem, the formation of such binaries is an even less understood topic.Highly detailed 1D semi-analytical work by Tagawa et al. (2020a) suggested that of order ∼ 90% of binaries in AGN discs will have formed within the lifetime of the AGN disc.therefore understanding the initial orbital elements of these in situ formed binaries is essential to understanding their later evolution.These binary formation events are expected to predominantly take place through what is known as a gas assisted capture whereby upon the scattering of two BHs in the AGN disc, a complex interplay between the binary and the surrounding gas and their CSMDs leads to a net removal of the binaries' kinetic energy (e.g Goldreich et al. 2002), such that they remain bound.This was investigated very recently in semi-analytical studies (e.g.DeLaurentiis et al. 2022;Rozner et al. 2022).These studies assumed orbital energy is dissipated by an Ostriker (1999) like dynamical friction prescription between the binary BHs and the local gas.Finding higher AGN densities favour binary formation in BHs with smaller radial separations in the SMBH disc.The first full hydro simulations of this process (Rowan et al. 2023;Li et al. 2022a) corroborated the overall analytical and semianalytical expectation that scattering encounters in gas may lead to binary formation, even when their initial orbital energy is positive upon first entering each others' Hill sphere, but the details may differ from those of previous semianalitical models.In our previous paper (Rowan et al. (2023), hereafter Paper I), the subsequent evolution of the binaries was also directly simulated following their formation, showing that their evolution depends on the AGN disc mass and their orbital configuration (i.e prograde or retrograde).Inspiral was observed in some prograde cases and for retrograde binaries, eccentricities were excited to high enough values in short enough times for binary mergers to occur within timescales of only a few AGN orbits.Both Paper I and Li et al. (2022a) found inpiraling binaries that form with very high initial eccentricities, possibly explaining the discrepancy with studies that used near circular initial binaries (e.g.Baruteau et al. 2011;Kaaz et al. 2021;Li et al. 2022d,c).Li & Lai (2022a) further corroborated this with a positive correlation between eccentricity and inspiral rate.Now that the viability of the gas assisted capture mechanism has been justified, the next natural question is to quantify the efficiency of this process based on pre-encounter parameters.Our previous gasless study (Boekholt et al. 2023) identified three 'islands' in the parameter space of the impact parameter (i.e. initial radial separation relative to the central SMBH)  of black hole encounters that lead to more than one encounter, i.e. "a Jacobi capture".The dependence of the number of encounters (orbits) during the Jacobi capture on the impact parameter exhibits a fractal structure.Furthermore, the large statistics of the work allowed one to very finely sample the parameter space of the encounter to quantify the fraction of the parameter space  that lead to temporary binary formation.Here, we examine if these temporary binaries formed without gas in Boekholt et al. (2023) lead to permanent capture through gas hardening.We also studied how gas affects the fractal structure of the parameter space of initial conditions that allow for successful binary capture.In the first paper of this series Paper I and in the work of Li et al. (2022a), the gas capture process was validated using a full hydrodynamical approach for the first time.In Paper I the formation and long term evolution of each binary was modelled self consistently in one continuous simulation, so the starting conditions of the binaries were directly set by their formation and not based on any assumptions.The work demonstrated that both prograde and retrograde binaries may form in varying AGN gas densities and identified a bimodality in the eccentricity evolution where prograde binaries circularised over time and retrograde binaries became more eccentric.
In this paper we explore the regions of phase space that permit the formation of a binary from two initially isolated black holes embedded in a gaseous accretion disc of a third massive body (SMBH) and the role of gas in shaping this space using a global 3D hydrodynamical simulation based on smoothed particle hydrodynamics (SPH).
This work is the second in a series of papers investigating each element of the gas-induced binary merger mechanism in AGN.Our investigation builds on Paper I by probing the binary capture process in far greater detail using a much finer sampling (114 simulations total) in impact parameter, providing great enough statistics to probe the stochasticity of the gas-assisted formation process and develop statistical correlations between pre-encounter parameters and important physical processes during the encounter.Such correlations can then be used to generate a vastly improved analytic prescription for semi-analytic works that discuss binary formation in accretion discs like that of Tagawa et al. (2020a).We present this paper alongside a companion paper Whitehead et al. (2023) which models the same scattering problem considered here but using a local shearing box approach and employing the grid code Athena++ instead of our SPH implementation.As we will discuss later in detail, the results of both papers are largely consistent with one another.Such similarities include the relationship betweeen the initial approach trajectory of the BHs and the level of dissipation and the dependence on gas density.The numerical method and initial conditions are detailed in Sec. 2, and the results are presented in Sec. 3. We discuss our results in Sec. 4 and the caveats of the models in Sec. 5 before summarising and concluding in Sec. 6.

Hydrodynamics
We consider a three-body system consisting of a central SMBH of mass  SMBH = 4 × 10 6 M ⊙ and two stellar BHs with masses  BH =  1 =  2 = 25M ⊙ embedded in a gas disc.The BHs are inserted symmetrically about a radial encounter distance from the SMBH,  mid = 0.0075 pc.The gas disc is resolved about the BHs as an annulus of radial width  disc = 20 H where is the Hill radius of one of the satellite BHs, from an inner radius,  in , to outer radius,  out .Based on our annulus width this gives  in =  mid − 10 H and  out =  mid + 10 H .The hydrodynamics of our problem are simulated using the 3D smooth particle SPH code PHANTOM (see Price et al. 2018) utilising 25 million particles, ignoring radiative and magnetic effects which for the distances we consider from the SMBH are less relevant (Jiang et al. 2019;Davis & Tchekhovskoy 2020).Such effects may be important for modelling the CSMDs of the satellite BHs, however due to the computational expense of incorporating them at our resolution, we also ignore them here.The AGN disc, where we will simulate the encounter of two BHs, is consistent with its implementation in our previous work in Paper I, though we summarise its features here for completeness.The disc is represented by a thin Shakura-Sunyaev constant alpha disc (Shakura & Sunyaev 1973) with a viscosity parameter  SS = 0.1 and scale height to radius ratio of / in = 0.005.This value corresponds to an accretion rate at this radius of 10% Eddington (see Goodman & Tan 2004).
The enclosed disc mass as a function of radius is assumed to follow the power-law of that derived in Goodman & Tan (2004) which is defined with four additional free parameters Here   is the disc luminosity relative to the Eddington limit,  is the mean molecular mass, κ is the opacity and  is the radiative efficiency.The disc has a sound speed  s profile given by (3) where the  s,in is the sound speed at  in and q is the power-law exponent.The surface density Σ profile is similarly represented as with its own radial dependence  .The value of Σ 0 calculated through normalisation of ( 5) such that 2 That is to say the radial integration of the surface density across the annulus matches the expected mass of the annulus based on (2).This gives a 3D density distribution of

Initial Conditions
In a similar manner to our previous work of the gasless case in Boekholt et al. (2023) we examine how binary formation varies with the approach of the binary, namely the impact parameter , closest approach  min,1 and relative two-body energy in the presence of gas with disk mass  d,0 given by Eq. (2).To do this we initialise the two BH satellites on circular orbits around the SMBH with varying initial radial separations which we will denote as their impact parameter  =  2 − 1 where  1 and  2 are the respective radial positions of the inner and outer satellite in the SMBH disc.The impact parameter  is sampled using 39 evenly spaced points over the range [1.75, 4.25] H .This is an ∼ 8× increase in the number of simulations compared to Paper I which only considered 5 impact parameters per set of initial parameters.In the previous paper, this limitation prohibited any statistical analysis of the effect of varying impact parameters.Now, using far more models, we can investigate this region of the parameter space for BH binary formation.Recent literature indicates that binaries can dissipate their relative energy more efficiently in denser AGN discs (e.g.Li et al. 2022a;DeLaurentiis et al. 2022;Rowan et al. 2023).Hence, we also probe more rigorously how the capture window is tied to the density of the surrounding accretion disc by performing another sample of 45 simulations using an AGN disc with a mass of 3 d,0 .We also compare to a gasless sample with the same initial conditions for reference (see Boekholt et al. 2023, for a detailed study of the gasless case).This gives the mass of gas contained within an annulus of width  H of 0M ⊙ , 110M ⊙ and 330M ⊙ for the gasless, fiducial and 3 d,0 models respectively, indicating the high amount of gas present in our simulations.
As we are now including gas in this experimental setup, there are several additional complexities that can affect the interpretation of the results compared with the gasless case.First, BHs form circum-single mini disks (CSMDs) before their encounter.The discs enhance the gravitational attraction between the satellite BHs since their masses can be approximately added to their host BH to leading (monopole) order in a multipole expansion.Additionally, the masses of the BHs ) are the mean radius and the width of the simulated gaseous annulus,  and  set the radial dependence for the sound speed and surface density across the annulus (Eqs.4-5),  is the scaleheight ( = 0.4  H ) which also sets the pressure and temperature in the disk via Eq.( 4),  E is the Eddington ratio,  is the radiative efficiency,  is the mean molecular mass, κ is the opacity relative to the electron scattering opacity (0.4 cm 2 g −1 ).can change via accretion in the lead up to the encounter which again can alter the amount each satellite perturbs the other prior to the encounter itself.To have a more similar set of BH and CSMD masses at their encounter we scale the initial azimuthal separation Δ of each object, depending on their impact parameter , such that their approach time is approximately the same.Specifically, Δ is scaled so all simulations have the same approach time as that with  = 2.5 H and Δ = 20 • , i.e This ensures that the amount of accretion and growth of the CSMDs is approximately equal for each simulation with different .We summarise our simulation parameters and initial conditions in Table 1.

RESULTS
In this section, we present and discuss the results of our simulation suites.Starting with an overview of the shape and size of the formation parameter space in Sec.3.1, how the energy dissipation of the encounter is influenced by the closest approach is discussed in Sec.3.2.The time dependence of each mechanism that can alter the orbital energy of the binary during the encounter is described in Sec.3.3.Accretionless encounters are considered in Sec.3.4 and the effects of an increased AGN disc density are presented in Sec.3.5.

The capture cross section with gas driven dissipation
The periapsis distance of each first encounter is shown in Figure 1.These are compared directly to resimulated gasless encounters using the exact same initial conditions for the BHs.As the computational cost without gas is orders of magnitude lower, a much finer resolution of 16k (i.e. 2 14 ) models with impact parameters in the range  = [1.5, 3.5] are used.It is immediately obvious that gas leads to a much larger window in impact parameters for encounters to penetrate less than a Hill radius.Moreover, there is only one wide valley of captures in the space of , with the minimum separations varying gradually across , unlike in the gasless case where there are two separated capture windows and sharp changes in the minimum separation for very particular initial conditions (as identified in Boekholt et al. 2023).Note that the deep troughs observed in the gasless case are possibly missed by the relatively coarse sampling in impact parameters for the simulations with gas.Additionally, gas leads to binary formation at depths far shallower than in the gasless simulations, i.e at larger separation at the first closest approach,  min,1 .In order to quantify directly the size of the parameter space in  where multiple encounters may happen, we define the initial onedimensional capture cross section  as the total "length" in the 1D Minimum separation for our fiducial run (blue) compared to the gasless case (black) for the first approach as a function of the "impact parameter", i.e. the initial radial separation relative to the SMBH, denoted by .The horizontal line indicates the Hill radius of one of the objects for comparison.
The gas clearly broadens the range of impact parameters of close encounters.
In the gasless  d = 0 simulation suite, binaries with a number of encounters of   ≥ 2 and   ≥ 4 are shown in red and blue respectively.For the our hydrodynamic simulation, all binaries with two encounters also performed four and remain bound.
parameter space of the impact parameter  that allows for at least two encounters within the binary Hill sphere.Even though many gasless and gaseous encounters later decouple after this point, we examine the dissipation required for a binary to initially be bound, not whether the system is stable to the 3-body interaction with the SMBH in the long term.More generally, the cross section for encounters to have exactly  encounters within the Hill sphere is defined in Boekholt et al. (2023) as where   is the number of encounters of a simulation,  represents the simulation number,  is the Kronecker delta which vanishes unless simulations have  total encounters and Δ  = 1 2 (  i+1 −  i ) is the length of the impact parameter interval between the discrete samples around the simulation .Here we are interested in the cross section for all models which have   ≥ 2. So we replace the Kronecker delta with an indicator function   () that vanishes if case  is not part of set .For our problem, this becomes All together this gives explicitly the initial capture cross section with an associated error  of We compute this quantity based on all 2 14 gasless and all 39 of our fiducial gaseous simulations.We get a cross section for two encounters of (  ≥ 2) = (0.0407±0.0011)  H for our gasless simulations and 0.98 ± 0.13 H when gas is included.The cross section for four encounters then drops to ((  ≥ 4) = (0.00219 ± 0.00025)  H for the gasless case but remains identical for our simulations with gas.
The drop in the gasless case follows as there is no dissipation mechanism to remove energy once the binary is formed other than through its interaction with the SMBH which is entirely chaotic, thus they all inevitably decouple after usually only a short number of orbits after being disrupted by the SMBH, see Boekholt et al. (2023) for details.The fact (  ≥ 2) = (  ≥ 4) when gas is included indicates that gas can dissipate their orbital energy rapidly enough so that all binaries that undergo two encounters harden sufficiently to prevent their disruption from the SMBH.Furthermore, based on these values we conclude that for our initial conditions, there is a ∼ 25 times larger cross section to have multiple orbits when gas is included based on (  ≥ 2) with and without gas, indicating that binaries are at least temporarily captured far more efficiently in the presence of gas.Whether these binaries are permanent is examined below.Directly comparing the (  ≥ 2) and (  ≥ 4) cross sections of our gasless simulations to those of Boekholt et al. (2023) ((  ≥ 2) = 0.12 H , (  ≥ 4) = 0.0051 H ) we find slightly lower cross sections while still within the same order of magnitude for both values, though we reiterate that that paper used a far larger initial azimuthal separation so the gravitational focusing will of course be different across those simulations compared to the simulations in this investigation.The slope of  log 10 ()/  was shown in Boekholt et al. (2023) to be constant at ∼ 0.68.We calculate this slope to be ∼ 0.64 , in good agreement.

Dissipation as a function of minimum separation
Though the encounters lead to largely chaotic gas flows, we develop ready-to-use semi-analytical prescriptions based on our highresolution simulations to improve the accuracy of purely semianalytical studies such as by Tagawa et al. (2020a), DeLaurentiis et al. (2022) and Rozner et al. (2022).In this work, we define explicitly the first encounter as the period between the BHs entering within 2 H separation and either the first apoapsis or upon reaching 1 H separation after the first encounter where the choice of end point depends on whether the binary remains bound or not, respectively.If the binary becomes bound (even if temporarily) we define the second encounter as the period between first apoapsis and second apoapsis or the time of exiting the Hill sphere, and so on for higher numbers of encounters.The two-body energy of the binary is given by where  bin =  1 +  2 is the total mass of the binary,  =  1  2 / bin is the reduced mass,   ,   are the positions and velocities of satellite BH  = (1, 2) and  is the gravitational constant.
Figure 2 shows Δ bin , the change of the binary energy during the first encounter as a function of closest approach in our fiducial simulation setup expressed in the natural units of  H,0 defined as i.e. the orbital energy of a binary with a semimajor axis equal to the Hill sphere of one of the BHs, where the sign of Δ bin is denoted by different point types as indicated in the caption.Also plotted is a best-fit line for the analytic power-law profile as where  min,1 is the minimum separation of the first encounter,  and  are dimensionless fitting parameters.
Figure 2 indicates an anticorrelation trend between the depth of the first encounter and the amount of energy dissipated, i.e deeper encounters tend to dissipate more energy.1There is a large scatter where the energy dissipated in a small bin of closest separations can vary by up to two orders of magnitude for small  min,1 .This scatter may be interpreted as an error in the predicted power-law exponent  as shown.For our parameters, the binary objects intersect each others' CSMDs at Δ ∼ 0.05 H . Thus since the trend continues outside of this range, the coupling between the dissipation and periapsis changes continuously with Δ/ H , rather than akin to a step function at the intersection of the BH CSMDs as previously thought in Paper I.There is no correlation for the second and later encounters (Figure 3), as the gas morphology which determines the amount of energy dissipation becomes chaotic.The first encounter is naturally the most important for discussing energy dissipation as if the binary energy is still too high after the first encounter then no further encounters are possible.This excludes the possibility of an entirely separate second encounter of the objects as they orbit round the AGN, which is discussed in Li et al. (2022d).

Measurement of dissipation rates
To uncover the physical origin of the results in Figure 2, we consider the cumulative energy dissipation per unit mass during the encounter period for the three energy dissipation terms, defined and calculated identically to our last paper Paper I as: (i) SMBH interaction,  SMBH -energy per unit time lost or gained by the binary system through shearing tidal forces induced by the SMBH.(ii) Gas gravitational dissipation,  grav -from the gravitational interaction with the surrounding gas.(iii) Accretion,  acc -due to conservation of mass and linear momentum of accreted gas particles onto the BHs.
These dissipation mechanisms are all calculated by taking the dot product of the velocity differential  1 −  2 with the differential accelerations  1 −  2 due to the force in question, as in eq. ( 17).
Where the second term on the right-hand side is the change in specific energy per unit time due to any mass change in the binary,  bin from accretion.
From this general expression, the dissipation due to the SMBH potential is given by where the acceleration of BH  = (1, 2) induced by the BH is The dissipation by the local gas gravity follows the same form, now summing over gas particles (or more specifically the SPH kernels used in the tree algorithm) where the acceleration of BH  = (1, 2) due to the gas, from all  p SPH particles of mass  p and position   is The accretion dissipation  acc is the work done on the binary due to the linear momentum of the accreted gas.Upon accretion of a particle by a BH, the BH's mass, position, velocity and acceleration are evolved by a mass weighted average of any  acc accreted particles over the timestep Where   ,   are the accelerations and velocities of the  = (1, 2) satellite BHs before the accretion event during the current timestep.Δ  and Δv  are the changes in their acceleration and velocity due to the impulsive momentum transfer upon accretion, where   in equations ( 22) and ( 23) are the BH accelerations and velocities just prior the accretion impulse calculation (i.e after  ,SMBH and  ,gas have been incorporated into   ).The quantities  p,j and  p,j are acceleration and velocity of accreted SPH particles.Hence, from eqns.( 22) and ( 23) the dissipation due to accretion is where With the total energy dissipation as If one integrates  over some time window, this gives the cumulative or 'net total' energy exchange from the binary over the specified time window.We note that the integration matches the net change in the binary energy exactly as we deconstruct the dissipation sources and record their instantaneous value at every timestep in the simulation, so no inaccuracies arise due incomplete sampling of the timesteps.

Distinguishing orbit families
We differentiate between three different types of encounters by separating them into three families based on the form of their trajectories, determined by their impact parameters , shown in Figure 4.In the co-rotating frame of the inner satellite, the first family consists of first encounters with low impact parameters where the approaching outer BH is perturbed to the left of the inner BH into a prograde orientation, i.e. orbiting in the same direction as around the SMBH.We will refer to this family as leftsided encounters (LS for short).2 Further, encounters with intermediate impact parameters where the approaching outer BH initially passes on the inside (right) of the inner BH and is dragged into a retrograde orbit, which we label rightsided (RS) encounters. 3The third family takes an initially similar approach as the RS trajectories except they turn back on themselves and reverse to a prograde orbit, i.e their trajectories are deflected by at least 180 • in the final lead up to the encounter.These we label as turnaround (TA) encounters.These three encounter types were also identified in Boekholt et al. (2023).However, since there was no dissipation mechanism in that study, the binaries were loosely bound with low angular momentum and could periodically switch encounter family.The complexity of the three-body problem and this ability to switch encounter families gave rise to a fractal structure in the number of encounters as a function of impact parameter.In Boekholt et al. (2023), three distinctly separated islands in  allow for temporary binary formation, which we do not find here, instead we observe one large formation window as a function of impact parameter.Though binaries may flip orientations afterwards, the first encounter trajectories of the binaries in Boekholt et al. (2023) obey the same three family classifications and rough position in the space of .However the size of the successful capture region of each family is significantly smaller.More specifically, the  encounters only result in binaries near the second trough at higher  in the curve of  min,1 vs , where we find encounters lead to binaries in the whole range between the two minima.Additionally the LS and TA encounters of the gasless study have their regions of possible temporary binary formation at higher  min,1 and are much narrower in .Hence, the inclusion of gas smooths out the troughs and blends the three possible windows for formation into one larger one.The three encounter families are separated by two direct collision trajectories at the two troughs in  min,1 vs .We can conclude that the complexity is removed as there is little to no changing of orbital families (and hence prograde/retrograde orientation) even after just two orbits as the gas is efficient at altering the energy and angular momentum of the binary, shown already in our previous paper Paper I.

Time evolution of the dissipation rates
Figure 5 shows the orbit families in the parameter space of  using the same colour code as in Figure 4 and compares their cumulative energy dissipation during the first encounter.As indicated by the top panel of Figure 5, the two large troughs in the encounter depth naturally mark the transition between each family and also between prograde and retrograde binary formation, where RS encounters are retrograde and LS and TA encounters are prograde.The first trough (on the left) locates the boundary in the space of  for LS and RS encounters and the second trough marks the boundary between RS and TA encounters.At these troughs is also where theoretically one can expect infinitely close encounters or direct collisions if the impact parameters are fine tuned.Since there is a switch between prograde and retrograde orientations, there is naturally a possible set of orbits at minima of the troughs with zero angular momentum approach with an eccentricity of unity.These extremely close encounters are far better resolved in the gasless 3-body simulations.As shown in Boekholt et al. (2023) which considers the problem in a three-body framework, it is possible to fine tune the initial conditions to have arbitrary close encounters at these two impact parameters and at 3 the satellites execute a right-handed turn in this case  others for subsequent encounters provided numerical sampling is not a limitation.
The middle row of Figure 5 shows that the family of the encounter can be a good indicator of the type of energy exchange the binary will experience.RS encounters typically dissipate the most orbital energy when compared to leftsided and turnaround encounters.In all cases there is an initially positive contribution from accretion and gravitational torques as the BHs approach separations of 2 H , as observed in Paper I, which is found to be a result of a gas pileup ahead of each BH on their approach.This is initially overpowered by the shearing force on the binary from the SMBH until the final year before the closest approach, where  grav dominates, driving the spike in  tot just prior to the first closest approach ( −  enc = 0) before  acc with some addition from  SMBH quickly overpower the gas gravitational forces, removing energy from the binary rapidly.Note that ultimately  acc and  SMBH have often the same order of magnitude, see further discussion in Sec.3.4.
Deconstructing the dissipation contributions per orbital family, we find that RS encounters typically dissipate the most energy and hence are most favourable for binary formation due to both the strong accretion and the work done by the SMBH.Interestingly, the timeevolution of the RS encounters are very consistent across all dissipation mechanisms, having net positive energy transfer from the gas to the satellites, an initially positive but rapidly negative contribution from accretion, and a negative contribution from the SMBH.Both the LS and TA encounters have a range of contributions, both positive and negative, from each dissipation mechanism.Simulations with impact parameters close to the RS encounter window exhibit the same positive energy dissipation as the RS encounters.LS and TA encounters near the deep troughs in the top panel Figure 5 4 using the same colour coding.The middle panel shows the total cumulative energy change for all fiducial models.The bottom row shows the cumulative energy transfer from the local gas gravity, accretion and SMBH respectively from left to right.We artificially extend the final cumulative value from the end point of the first encounter right up to the ten year mark to more easily compare the net cumulative energy dissipation across models.Additionally, we denote unsuccessful captures with dashed lines and successful captures with solid lines.From the colour-coded results, rightsided encounters lead to the most reliably efficient binary formation which we attribute to the encounter depth dependent dissipation of Figure 2. . in this small region in the space of .As they deviate to higher or lower impact parameters from the central RS encounter region, dissipation gradually becomes less significant and flips when the impact parameter  is sufficiently high (≳ 2.8) or low (≲ 1.6) for our parameters.This is explained by the coupling of Δ bin with the close approach depth  min,1 and the behaviour of  min,1 with .Where the binaries with the closest approaches, occurring in the RS encounter region, dissipate the most energy.Binaries outside the RS encounter regions have increasingly distant close approaches and so dissipate less energy.
In the wings of the encounter outside the RS region, gas gravity becomes more negative (removing energy from the binary) while dissipation from accretion becomes more positive due to less accretion during the closest approach and dissipation from the SMBH becomes positive.The reversal of the SMBH dissipation predominantly applies to the turnaround encounters as after they initially pass each other, SMBH shear flips from slowing down the binaries' approach, to accelerating them after the first periapsis passage.This is because the shear of the SMBH acts along the direction of the radial vector from the SMBH to the binary COM.As the binary approaches this shear acts against the binary motion, removing the relative kinetic energy of the binary.After a turnaround binary executes its first pe-riapsis, the binary separation is increasing along this same vector, so the shear leads to an increase in their velocities and thus their relative energy, i.e doing positive work.

Torques during the encounter
Similarly to the dissipation, the strength and sign of torques induced on the binary depend on the impact parameter and family of the encounter.The cumulative torques over the first encounter (Figure 6) show that the net torque moves from positive to negative values as  increases.In terms of encounter families, LS encounters experience net positive torques, TA experience net negative torques and the net torque of RS encounters transitions smoothly from positive to negative.Due to the already highly eccentric nature of the binaries, we find torques highly insignificant for predicting the binary outcomes as the binary motion is almost entirely radial and so the torques ( 1 −  2 ) × ( 1 −  2 ) are minimal for the first encounter.

Accretionless encounters
In this work and in Paper I accretion appears to dissipate energy on the same order of magnitude as gravitational gas drag.Here we test 10.0 7.5 5.0 2.5 0.0 2.5 5.0 7.5 10.0 t t enc [ 6. Cumulative torque/angular momentum change in all our binaries during the first encounter from 2 H to apoapsis/exiting of Hill sphere.Results are colour-coded by encounter family; leftsided (red), rightsided (blue) and turnaround (green) encounters.The top panel shows the total cumulative angular momentum change for all fiducial models.The bottom row shows the cumulative energy transfer from the local gas gravity, accretion and SMBH respectively from left to right.We artificially extend the final cumulative value from the end point of the first encounter right up to the ten year mark to more easily compare the net cumulative torque across models.Additionally, we denote unsuccessful captures with dashed lines and successful captures with solid lines.From the colour-coded results, the torque switches sign smoothly when transitioning from leftsided to turnaround encounters.This results from the SMBH torque's dependence on the angle between the vectors from the SMBH and the binary COM and from one satellite to the other.whether captures may still occur in the absence of accretion during the close encounter.Figure 7 shows the separation vs time for our fiducial run along side 30 resimulated runs of each simulation where accretion is switched off when the BHs have a separation less than 2 H .The softening radius of 0.01 H is left unchanged.In switching off accretion at this point we preserve the initial encounter trajectory and immediate pre-encounter mass of the fiducial models so they can be directly compared.Figure 7 indicates that removing accretion leads to qualitatively identical outcomes with all successful binaries in the fiducial simu-lation suite also successfully forming in the accretionless suite and vice versa for unsuccessful captures.In addition to the scattering outcome, the depth of the initial encounter is largely unchanged, aside from the location of the two large troughs where accretionless encounters have deeper encounters for the assumed particular values of .This is likely only coincident with the sensitivity of the depth with  at these two turnaround points between prograde and retrograde encounters, where it was shown in Boekholt et al. (2023) that fine tuning  can lead to extremely deep encounters for only small changes in .Hence the small variations in the trajectory due to the exclusion of accretion can lead to very different depths at these two key points in the parameter space of .The fact the encounter depths change minimally is unsurprising as although accretion is shown to dominate overall, this is only due to a sudden reversal at the first periapsis passage,  −  enc ≈ 0. Thus the encounter trajectory is only perturbed by the SMBH and gas gravity prior to the close encounter, so the encounter depth remains largely unchanged.The capture cross section for accretionless encounters is then calculated to be  noacc = 0.92 H in the range 2.34 H <  < 3.26 H , identical to the fiducial model.
The non-negligible effect of removing accretion is in the behaviour of  grav immediately after  − enc ≈ 0, shown in Figure 8.This occurs during the time frame we would expect the extremely negative energy dissipation from accretion at the first periapsis passage,  −  enc ≈ 0. When accretion is removed, the time evolution of  acc in the fiducial simulations is reproduced in  grav , demonstrating a highly efficient period of energy removal from the binary after  −  enc ≈ 0, driving the cumulative dissipation to negative values.More formally, the time evolution of  grav without accretion mimics the rapid transition from positive to negative values observed in  acc in our accretion simulations.This is more akin to the findings of Li et al. (2022a) where gravity acts to add energy prior to periapsis then reverses afterwards, resulting in a net energy removal.This is not present when 10.0 7.5 5.0 2.5 0.0 2.5 5.0 7.5 10.0 t t enc [  Results are colour-coded by encounter family; leftsided (red), rightsided (blue) and turnaround (green) encounters.The top panel shows the total cumulative energy change for all fiducial models.The cumulative energy transfer from the local gas gravity and SMBH are shown respectively on the left and right of the second row.We artificially extend the final cumulative value from the end point of the first encounter right up to the ten year mark to more easily compare the net cumulative dissipation across models.Additionally, we denote unsuccessful captures with dashed lines and successful captures with solid lines.Turning off accretion leads to the replication of its effects in the dissipation from the local gas gravity.
accretion is enabled because accretion can alter the gas morphology close to the binary.When accretion is switched off, gas that would normally be accreted upon entering  acc can instead pass around the BH.Alternatively, in the frame of the BH, the accretion headwind felt by the BH can now pass directly through the accretion boundary and accumulate behind the BH.To visualise this, we calculate a surface density plot of the dissipation rate, averaged over the first encounter, from the local gas gravity <  grav / > and show this alongside the radially integrated cumulative dissipation over the first encounter period in Figure 9 for an accretionless and accretion enabled simulation case.
Figure 9 shows, as expected, that the morphologies are largely similar prior to the closest approach in the simulations with and without accretion.In the accretionless simulation there is a density buildup within the softening radius of the gas's gravitational effect on the satellites,  soft though this is largely axisymmetric and does not strongly affect the binary trajectory at this moment.Just before and after the periapsis passage there are dense gas clumps behind each BH in the accretionless encounters where the gas has passed around the BH and accumulated behind them, similar to the wakes caused by dynamical friction.
Figure 10 shows a zoom-in on the cumulative work done by the gas gravity on the satellites in the simulations with and without gas accretion near the first periapsis passage.The dissipation from  grav oscillates violently over the close encounter period.This was first observed in Paper I and is a result of the formation of nonaxisymetric, highly dense, inhomogeneities in the CSMDs from the BH tidal forces.As they orbit their respective BHs, the clumps tug on the BHs, resulting in a periodic oscillation in  grav .We find the oscillations have an orbital frequency which is akin to a Keplerian orbit of ∼ 0.4 min,1 .Though the strength of the oscillations in the binary energy warrants a discussion of the possibility they may have their own gravitational wave imprint on the encounter in a similar manner to other gas-induced GW signitures suggested in the literature (e.g.Kocsis et al. 2011;Yunes et al. 2011;Hayasaki et al. 2013; Barausse et al. 2014Barausse et al. , 2015;;Cole et al. 2023;Nouri & Janiuk 2023).However the characteristic frequency of the oscillations is typically of the order ∼ 10 −7 − 10 −6 Hz which is unfortunately too low to be measured even with LISA.
We conclude that the formation of these inhomogeneities is due to the tidal warping of the CSMDs being maximal at the periapsis of the encounter, this is further evidenced by the fact they form exactly during this part of the binary orbit (see Figure 10).The characteristic orbital radius of the inhomogeneities is therefore unsurprising as the BHs do not pass within 0.5 min,1 of each others CSMDs, so the gas in these regions is better retained around their original BH following periapsis.Thus gas can remain at high densities within distances smaller than the periapsis separation (i.e | −  1 | < 0.5 min,1 for the inner BH).
The oscillating dissipation phenomenon from the inner CSMDs is present in both simulations with and without accretion, however the oscillations are larger in the latter case.The cumulative effect of the rapidly varying  grav , is a net dissipation in the binary energy by the time of the first apoapsis passage.In the simulations with accretion, this is net positive in most cases, i.e. the binary gains energy and becomes less bound, in agreement with our previous results in Paper I. In contrast, in our accretionless simulations, the net dissipation is significantly negative.This is the result of the allowed permeation of what would have been the accretion radius in the previous accretion enabled simulations by the gas.Gas that would normally accrete against the BHs motion (see Figure 5) and produce the strong negative dissipation instead accumulates behind the BH.While the CSMDs perturbations are still present, this permeation leads to a net removal of binary energy via gravitational drag.
Figure 11 directly compares the net contribution of various physical processes to the binary energy dissipation during the first encounter.We find that artificially removing the accretion leads to qualitatively similar behaviour, with the RS encounters dissipating the most energy, where now it is dominated by  grav and  SMBH instead of  grav and  acc .While enabling accretion significantly alters  21) over the first encounter for a representative accretion-enabled (top) and accretionless (bottom) simulation with identical impact parameters.The map is centred on the centre of mass  CM of the binary and distances normalised to the current separation of the binary Δ.The dissipation is normalised to the max of both simulations to compare the relative strengths.The black lines show equal contours in each simulation, showing the enhanced dissipation close to the BHs in the accretionless case.Top panel: Average gravitational energy dissipation integrated from the  −  centre of mass radially outwards in the plane of the binary for the accretion-enabled case.Bottom panel: Same as the top panel but for the accretionless model.The radial and 2D maps are aligned so they can be directly compared along their axes.The binaries execute their orbit in a counter clockwise direction.
the impact parameter dependence of the integrated  grav , we find that the relative strength of the gas effects is similar across the range of impact parameters when comparing  acc +  grav to  grav in simulations with and without accretion.This suggests that although the dominant dissipation term shifts from gas gravity to accretion when accretion is enabled, it preserves the overall expected time-dependent variance and strength of the overall dissipation at least to within a factor of few.

Dependence on different accretion disc densities
We now consider BH scatterings in a 3 times higher density environment than in our fiducial model but otherwise identical initial conditions so that the disk mass is  d = 3 d,0 .This includes the modelling of accretion as in the fiducial simulations.This suite is comprised of 45 simulations, 6 more than the fiducial run to span the larger space of impact parameters for hard encounters.Figure 12 compares the encounter depth vs impact parameters.The window of captures in  is greatly enhanced with additional gas.This matches the conclusion of DeLaurentiis et al. ( 2022), Li et al. (2022a) and our sibling paper; Whitehead et al. (2023).The added gas also has the effect of shifting the 'valley' of the encounter depths to higher impact parameters, this explains why for our high disc densities in Paper I we found very different encounters compared to the lower mass simulations.There are two effects causing this shift, one is the increased mass buildup in the discs, which means the BHs are perturbed earlier along their orbit since the disc mass can effectively be added to the BH mass when they are sufficiently far away.The second is that accretion causes the Hill sphere of the objects to increase on their approach to each other.We remove the accretion dependence of the window in the bottom panel by normalising the cross sections, which are presented in units of the Hill sphere, to the Hill sphere of the BHs when they intersect twice each others Hill radii,  H 2H , instead of their initial Hill radii,  H t=0 , and denote these cross sections with a tilde: This centres the windows perfectly, showing that accretion is responsible for the shifting of the window.However, the overall window of the encounters is still larger, where the two troughs are further apart as well as the RS and turnaround encounter regions occur at lower and higher impact parameters, respectively.Additionally, captures successfully take place at shallower initial encounter depths, even when  min,1 > 0.1 unlike in the fiducial model.This evidence points to encounters in AGN discs with higher densities being more favourable for binary formation; allowing binaries more sparsely separated to have stronger encounters and affording binaries with weaker encounters a greater chance for formation.Minimum separation for our fiducial run (blue) compared to the 3 d,0 run (orange) for the first approach as a function of the initial radial separation .On the top row is the unaltered initial form of the curve.On the bottom we show the same two curves, where the minimum separations are now normalised to the size of the Hill sphere when the binaries reach 2 H in separation.The normalisation process is described by eq. ( 28).12) for successful binary formation.These values quantify formally the range in  that permits binary formation.Shown from left to right are the different simulation suites, the standard cross section , the accretion normalised cross sections λ (see eq. 28), and the ratio of the normalised cross sections to that of the fiducial simulation suite.
Figure 13 shows the contributions of the different physical processes to the satellite energy dissipation and the torques over the first encounter and per source over time.Qualitatively, the overall trend is the same initially.Dissipation is initially negative due to the SMBH, before rapidly peaking positive just prior to encounter due to gas gravity followed by a sharp drop due to strong negative dissipation from accretion.The torque behaviour is also preserved from the fiducial model with the SMBH dominating the torque in the negative direction for nearly all models, while the gas gravity and accretion induces negative torques in the turnaround encounters, gradually flipping to positive torque across the window of leftsided encounters and inducing significantly positive torques for RS encounters.While the behaviour of the dissipation mechanisms is unchanged from the fiducial model, the strength of the dissipation and torques from both accretion and gas gravity are around a factor of 3× larger (see Figure 14 for clarity).For the torques this is less impactful as the SMBH still dominates by a factor of a few.However the dissipation becomes enhanced by this factor of five as it is dominated by gravity and accretion, indicating that more energy can be dissipated in higher gas density environments.Note that this is independent of any enhanced mass gain by the BHs as our energy unit normalisation  H,c and its reduced mass  dependence is calculated for each model.The capture cross section reflects this favourability, with a value of λ3 d = 1.591.In Table 2 we summarise the normalised and unnormalised cross sections from each of our simulation suites.[yr] 0.00 0.01 0.02 acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc 10 5 0 5 10 t t enc [yr] 0.00 0.02

SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH
Figure 13.Cumulative energy and angular momentum change in all our binaries from the run with three times the AGN disc density (3 d,0 ) during the first encounter from 2 H to apoapsis/exiting of Hill sphere.Results are colour-coded by encounter family; rightsided (red), leftsided (blue) and turnaround (green) encounters.The top panel shows the total cumulative energy change for all fiducial models.The cumulative energy transfer from the local gas gravity, accretion and SMBH are shown respectively from left to right of the second row.The cumulative torque is shown on the 3rd row, with the breakdown from the three dissipation mechanisms similarly displayed in row 4. We artificially extend the final cumulative value from the end point of the first encounter right up to the ten year mark to more easily compare the net cumulative torque across models.Additionally, we denote unsuccessful captures with dashed lines and successful captures with solid lines.The enhanced gas density does not change the form of the dissipation or torques but increases their strengths.
Recall that we have shown in Figure 2 that the depth of the first encounter greatly affects the orbital energy dissipation of the binary for simulations with fiducial disc model with different impact parameters.Figure 15 shows the energy dissipation vs closest approach for all our simulations, including those with an increased disk mass and those where accretion is turned off.From the results, we find that compared to our fiducial model, the density enhancement and removala of accretion leads to a steeper dependence of the energy dissipated with the periapsis distance  ,1 .Additionally, if we compare the  coefficients we can quantify the average energy dissipation increase with density.With a value of 3.48/1.28= 2.69 ± 1.33 this corresponds to a scaling with the disc density Σ of Δ bin ∝ Σ 0.9±0.44 close to the prediction from our other paper Whitehead et al. (2023) that finds dissipation scales with density to the power of unity, albeit with a considerable error margin here.

Parameter space of captures
While analytic and semi analytic studies of BH binary population synthesis in AGN (e.g Miralda-Escudé & Gould 2000; Freitag et al. 2006;Hopman & Alexander 2006;Bartos et al. 2017b;Rasskazov & Kocsis 2019;Secunda et al. 2019;Tagawa et al. 2020a;McKernan et al. 2020b) are well suited for simulating the large BH population in AGN, they currently utilise very simplistic assumptions for the outcome of an encounter, oversimplifying the complex interactions with the gas.A revised, more physically motivated set of analytical tools for deducing the outcome of an encounter would provide far more realistic and reliable estimators for the binary formation rate which is essential for estimating the BH merger rate for GW astronomy.While we have numerically derived an expression relating dissipation and encounter depth, the encounter depth itself is affected non trivially by the gaseous effects.A more useful expression may be obtained if one ascertains a parameter space for captures using encounter parameters that are not largely effected by the complex gas morphology of the encounter itself.We consider two such parameter spaces, starting with the the binary relative energy at the start of the encounter  2H at 2 H and the impact parameter of the binary at the Hill radius  1H and .We define these quantities specifically in eq. ( 29) and eq.( 30) respectively.Where eq. ( 29) is equation eq. ( 14) evaluated when the binaries are at 2  separation.
Where as before,   and   are the positions of the satellite black holes  = (1, 2). bin =  1 +  2 is the mass of the outer binary and  =  1  2 / bin is its reduced mass.To reiterate, eq. ( 29) is evaluated when the binary separation first reduces to 2 H and eq.( 30) is evaluated when the separation reaches  H .
Before constructing our analytic tool for predicting the energy dissipated, we first have to understand how much energy must be lost.Figure 16 shows the total energy of the binary after the first encounter,  f , i.e after the work done by the three dissipation mechanisms during the encounter, Δ bin , has been added to the initial energy at 2 H ,  2H .One might consider using the energy at 2 H but the impact parameter at  H an odd choice for our parameterisation.However, as shown in Figures 5, 8 and 13 the energy exchange of the binary becomes important at scales beyond the Hill radius.Though the dissipation from 2 H is important for determining the outcome, we find no binary formation beyond just  1H ≳ 0.68 H . Thus in semi-analytical studies, one needs only consider binaries entering a single Hill radius and have access to  2H and  1H to determine the scattering outcome.
In the plot, a very clear parameter space for successful binary formation emerges.At the bottom of the parameter space there is an island of failed captures, where for  1H / H ≳ 0.68 energy is added to the binary.In the other failed capture region where  1H / H ≲ 0.68, this means the binary did not dissipate enough energy to be sufficiently hard so as not to be decoupled by the SMBH.The slope in the boundary between each island implies that before first periapsis, a binary with identical initial energy must dissipate more to remain bound if its impact parameter is higher.This is no doubt due to correlation between  1H and the encounter depth  min,1 , which we confirm in Figure 17.We define this sloped boundary as the function  f,crit (  1H ), shown below in eq. ( 31).
Next, we determine the scaling between  1H and  min,1 , which is accurately informed by our full hydrodynamical approach and therefore includes the effects and stochasticity introduced by the gas.We show our function for  min,1 in terms of  1H (see eq. 32) alongside the data in Figure 17.
Orbital energy, in units of the energy of a circular orbit at the Hill sphere  H,c , of the binaries after the first encounter (either when executing first apoapsis or upon leaving the Hill sphere if unbound),  f vs their impact parameter at one Hill radius  1H as defined in (30).The results show three distinct regions of parameter space, 1) where binaries successfully form through sufficient energy dissipation, 2) where binaries dissipate energy but not sufficiently to remain bound and 3) where binaries decouple through energy gained during the encounter.The function of the positively sloped line separating the successful from unsuccessful formations is the critical final energy,  f,crit , needed to remain bound, represented by the log-linear function in eq. ( 31).First encounter periapsis scattered against the impact parameter at one Hill radius  1H (Eq.30) for all simulations, showing that smaller impact parameters lead to closer encounters, with some scatter for lower values of  1H where the very small periapsese becomes increasingly sensitive to small changes in energy.Also shown is the log-linear line of best fit (solid red) and its error (dashed red) according to eq. ( 32).
Using our fits of the Δ bin ( min,1 ) relation to each simulation suite, we construct the parameter space of  1H −  2H that allows for binary formation in eq. ( 33) In simple terms, this equation calculates the critical final energy required to remain bound based on its impact parameter  1H and subtracts how much energy it was expected to dissipate based on the same impact parameter to calculate its expected energy at 2.Therefore, it calculates the expected energy at  2H that would allow a binary to remain bound for its impact parameter  1H .In explicit form, both the first and second term take on power laws.When all , for all our models as a function of their impact parameter as measured at  1H overlayed on our analytically derived parameter space.The red area highlights where we expect failed captures due to a lack of sufficient energy removal.The green area is the parameter space where we expect successful captures.The grey area presents a firm barrier to capture as energy is added to the binaries in this region which don't pass deep into each others' Hill sphere.The overlapping red and green area is the area of the parameter space where the outcome depends on our simulation suite.The boundary of the successful formation region for each suite is traced by the black curves and labelled in the plot, calculated via eq.( 33).
constants combined and simplified this gives: We determine the constants for this equation for all of our simulations suites collectively and individually and display them in Table 3.According to our criterion, a BH-BH scattering will result in a successfully formed binary if the following is true This new expression meets our aim of being independent of any variables that could be affected by the complex hydrodynamics of the encounter and since the energy is normalised to  H,c it is also independent of the binary mass, thus is insensitive to any mass buildup due to accretion before the encounter.We summarise the process of constructing the binary formation criterion in short steps below.We apply our criterion to the data first in Figure 18.The figure shows three distinct regions, (i) The already known region of zero binary formation at  1H ≳ 0.68.
(ii) The region of successful binary formation under the curve where the binary can remove enough energy to remain bound and unperturbed by the SMBH.(iii) The region above the curve where the binary was too energetic to dissipate enough energy to form a stable binary.
The vertical dashed line indicates the transition to where the binary must already have a negative two-body energy to remain bound for our fiducial model.The curves that separates the successful and unsuccessful captures are our criterion eq.34 calculated for each suite using the corresponding  and  values.We summarise the methodology of constructing our capture criterion in numbered steps.

Summary of Binary Formation Criterion
Step I Determine the binary's two body energy at a separation of two Hill radii,  2H , using eq.( 29).
Step II Determine the binary's impact parameter at separation of one Hill radii,  1H , using eq.( 30).
Step III Use the impact parameter  1H to determine the close approach distance  min,1 using eq.( 32).
Step V Calculate the minimum energy  f,crit required to stay bound via eq.( 31) using  1H .
Step VI Determine if the energy dissipation criteria, eq.34 is satisfied using  f,crit , Δ bin and  2H .With constants for eq.34 given by Table 3.
While this criterion is a significant improvement, it naturally still makes assumptions.These include its assumption of only 2D encounters and equal mass ratios for the two BHs.Further work is needed to determine generalisations for non-coplanar encounters misaligned with the AGN disc, and for asymmetric mass ratios.Additionally, the relationship may be depend on the mass of the BHs, since we have not verified that the mass of the minidiscs scales one to one with the BH mass.Given a large portion of the dissipation in rightsided encounters comes from colliding minidiscs, alternate ratios of minidisc to BH masses due to differing BH masses could affect the amount of energy dissipated in the same manner as changing the ratio via changing the AGN disc density.
We verify the accuracy of the criterion in Figure 19 by breaking it down by suite and also showing the errors in the curves.Counting the false positives and false negatives in each suite we find a predictive success rating of 94%, which for the errors in  and  of eq ( 16) is considerably good.Additionally, though we have shown its form to be mass dependent based on the curves in Figures 18 & 19, our relation only holds for the expected mean AGN density for our parameters and our inflated 3x density model.A more formal density scaling would allow the relation to hold for arbitrary or at least a range of AGN densities, which is the subject of our sibling paper, Whitehead et al. (2023).

Implications for mergers
Given the high efficiency with which binaries form in the simulations and that most of them are retrograde, if we combine this with the result from Paper I that retrograde binaries have a non-negligible merger rate, it would appear that overall the AGN channel is highly conducive to black hole mergers.In reality however, there is also non-zero radial and vertical velocity dispersion in the stellar mass BH population in AGN that is not encorporated in this work.The additional vertical component will naturally lead to interactions outside of the midplane where the gaseous background is less dense, leading to a lower probability of capture as indicated here, in Paper I and in Li et al. (2022a).Additionally, little is known of what gas morphology we should expect around a satellite in an AGN disc.We expect encounters with BHs that have entered from outside or at the extremities of the AGN disc will involve less gas since they cannot replenish gas lost from their Hill sphere to the SMBH as quickly or at all prior to the encounter (Kratter et al. 2010;Kocsis et al. 2011).Though this additional aspect to the problem will likely reduce the chance for binary formation and mergers, it is expected that the inclination of objects will decrease significantly over the lifetime of the AGN due to dynamical friction (e.g.Tagawa et al. 2020a;Generozov & Perets 2023) between the satellites and the SMBH disc.
Based on the density dependence of the capture criterion (eq.( 34)) favouring higher densities this would suggest a bias towards formations and merger of BH binaries closer to the SMBH according to eq. ( 5).However there is the competing factor of the velocity shear of the satellites, which is more pronounced closer to the SMBH.An increased shear increases the value of  2H which is calculated via eq.( 14) by enhancing  2 rel = ∥ 1 −  2 ∥ 2 .If we assume the radial density dependence  −0.6 of Eq. ( 5) used here, and that  2 rel ∝  −1 then whether BH binaries are more likely to form closer or further from the SMBH depends on the steepness of the exponent of eq. ( 34) with Σ (the surface gensity of the gas, eq. ( 5)).For the binary formation probability to be independent of  then it would require the dissipation in the first encounter to scale as Δ bin / H,c ∝ Σ 1.6 .Here we narrow the dependence down to Σ 0.90±0.44 .In conjunction, our wind tunnel simulations of Whitehead et al. (2023) narrow this slope down to 1.01 ± 0.04.Thus, as a rough estimate, we expect the probability, (), of a highly co-planar BH-BH strong encounter ( min,1 <  H ) forming a binary successfully to scale as () ∝  −0.6 .Hence we would expect BH binaries to form more easily closer to the SMBH.

Comparison to similar studies
Our results are a natural extension of the pure 3-body results of Boekholt et al. (2023) without gas.In Boekholt et al. (2023), threebody scatterings of BHs orbiting a SMBH were performed using a purely N-body approach.In the absense of gas, the simulations showed the formation of temporary BH binaries, where binaries with more than 2 encounters are possible for a smaller range of impact parameters, though they are short lived.Where the range of impact parameters that lead to higher encounters numbers becomes exponentially smaller for each additional encounter.Whether they have a deep enough encounter for GWs to induce mergers is dependent on the chaotic evolution of the 3-body system, which can only rarely lead to these chance extremely close encounters.Aside from these extreme cases the lack of a dissipation mechanism leads to their eventual disruption by the SMBH, with probability of surviving multiple  encounters dropping off rapidly with  (e.g Li et al. 2022d;Boekholt et al. 2023).As shown in this work, the inclusion of gas leads to the formation of reliably hardened binaries in a single large window in impact parameter space, destroying the fractal structure observed in the paramater space of the impact parameter and number of encounters in Boekholt et al. (2023).This results from the inability of the binaries to switch orientations reliably after even one encounter  18 but for all simulations including those with an increased disk mass and those where accretion is turned off as a function of their impact parameter as measured at  1H overlayed on our analytically derived parameter space.The red area highlights where we expect failed captures due to a lack of sufficient energy removal.The green area is the parameter space where we expect successful captures.The grey area presents a firm barrier to capture as energy is added to the binaries in this region which don't pass deep into each others' Hill sphere.The overlapping red and green area is the area of the parameter space where the outcome depends on our simulation suite.The boundary of the successful formation region for each suite is traced by the solid colour-coded curves according to the simulation suite.The dashed curves span the error range of the boundary region.Out of 115 simulations, we find only 5 outliers.We highlight the point where accretionless encounters that already have a negative two body energy still need to dissipate energy further with the horizontal and vertical dashed lines.
due to the efficiency of the gas-induced energy dissipation enabling the binary to rapidly harden.DeLaurentiis et al. ( 2022) verify that at exceptionally low ambient gas densities, the fractal structure in the number of encounters vs impact parameter reappears in their 3-body scatterings that include semi-analytic dynamical friction.The recovery of the fractal structure at sufficiently low densities is also hinted in Whitehead et al. (2023).Where lower density encounters have a greater tendency to flip orientation, most commonly at the boundary between families.Though, the sampling is too course to make any reliable comment on the recovery of the fractal structure.
A closely related work to ours is the gas scattering experiments of Li et al. (2022a).In their 2D hydro study of embedded BH scatterings, which uses multiple density values, they show a roughly linearly increasing window in their impact window size as a function of initial ambient density.Instead of varying the radial separation they instead vary the initial azimuthal separation of their satellite BHs, which is possibly why they find nearly all binaries to be formed in a retrograde configuration.In this work we varied the initial radial separation in the AGN disc and scaled the azimuthal separation so the approach time for each simulation was identical.Using this alternate setup, we find both prograde and retrograde binaries, though likewise find the cross section to be larger for retrograde encounters, which make up the RS encounter region, see Sec. 3.3.2.Comparing the outcomes, there is good agreement in the dissipation mechanisms being positively correlated with the local gas density.There is also agreement that for lower AGN densities, the maximum allowed encounter energy is lower, as highlighted in the final figure of our paper, Figure 19.
Comparing to the recent semi-analytical studies, DeLaurentiis et al. ( 2022) consider three-body scatterings with an Ostriker dynamical friction prescription to account for gaseous effects.DeLau-rentiis et al. (2022) find similar "island" regions of  to our purely gravitational study of Boekholt et al. (2023), where binaries may successfully form.Whereas with our inclusion of gas we find only a single central valley of captures.Given they switch off their dynamical friction at the closest approach they neglect the subsequent dissipation between periapsis and their approach back towards the Hill sphere, which can be significant, as shown in Figure 5 for example.This extra dissipation could harden the extremely low angular momentum encounters at the troughs that have the largest apoapses and naturally are more prone to being disrupted by the SMBH.Additionally they find increasing  d leads to the formation window moving to lower  values which is the opposite to our findings.Similarly to Li et al. (2022a), this discrepancy in the approach trajectories could be due to their azimuthal separation being far smaller than in this study.
In a sibling paper to this (Whitehead et al. 2023), we investigated the same scenario posed in this paper using a shearing box model using the Eulerian grid code Athena++.The use of the new code allows us to compare the two different numerical approaches to the hydrodynamics and by using a shearing box, maintain a higher resolution inside the Hill sphere for less computational cost.We find excellent agreement on several aspects of the BH-BH encounters.This includes the order of magnitude of dissipation expected for gas assisted encounters, which lies around 10 − 100 times the energy of the binary upon reaching a separation of one Hill sphere.Also corroborated is the identification of the three encounter families that lie at low, moderate and high impact parameters, which have the orientations prograde, retrograde and prograde respectively (relative to the orbit about the AGN).The slope of Δ/ H,c as a function of  min,1 / H is also corroborated to within 10% percent of the values shown here.A larger sampling of  d is performed in (Whitehead et al. 2023) which finds that sufficiently low densities the single formation window splits down the centre into two separate islands of capture.Though the shape of the capture window in (Whitehead et al. 2023) is identical to those here (i.e a single window with three types of first encounter), the cross section and location of the window tends to be at lower values of  for the same AGN disc densities.However, comparing the size of the formation window of the global simulations here to the straight wind tunnel of (Whitehead et al. 2023) directly is likely not reliable since the linear geometry of the shearing box model will lead to different gravitational focusing during the initial approach.

CAVEATS
• In this paper we restricted attention to strictly co-planar encounters in the midplane of a 3D gaseous disk, but in reality there is also a 3D distribution of BHs around AGN with some vertical velocity dispersion.Inclination has been shown by Li et al. (2022d) to reduce the likelihood of encounter, thus inclination may alter the size of the capture windows shown here.Though they do not include gas in their study.In a follow up to our 2D study in preparation, we extend the previous work in Boekholt et al. (2023) to 3D.
• We assume gas heating/thermal shocks do not significantly increase the gas temperature anywhere in the entire simulation domain.While this may be a good assumption for the main body of the annulus distant from the encounter, the high speed intersections of the accretion discs of the BHs will in practice lead to significant heating of their material due to their high density and relative velocity.Tagawa et al. (2022) show that even single black holes can provide a heating mechanism to the local gas through accretion driven jet outflows.Therefore we encourage future studies to include such effects if feasible.
• Despite the number of simulations shown here, we are still dealing with far fewer statistics than 3-body studies and recommend further simulations to further probe additional parameters that affect the dissipation processes.For example, can we constrain analytically the slope of the disspation -encounter depth relation and improve the errors on the gradient overall.
• In addition to non inclined orbits, our BHs are initialised with only circular orbits around the SMBH.A more accurate nuclear BH population would include also include a radial velocity dispersion, i.e orbits around the SMBH would have their own eccentricity as well, which would affect the energy of a binary going into an encounter.
• We have shown gravitational dissipation and torques from the CSMDs during the encounter can be reproduced by accretion, despite the accretion radius being far larger than its true value by many orders of magnitude.However there is still a factor of a few difference in the strength of these effects between some models.Additionally, as in Paper I, we find the accretion to be super-Eddington by several orders of magnitude.It is possible that a more accurate handling of the accretion leads to results different to both the accretionless and accretion enabled results shown in this paper.However, accretion at such small scales necessitates extreme resolutions and the inclusion of radiative effects which becomes extremely computationally expensive.

SUMMARY AND CONCLUSIONS
Using a hydrodynamic numerical simulations, we performed a quantitative analysis of the binary formation process in AGN discs, with the aim of constraining the likelihood of binary formation depending on encounter parameters and providing the community with physically motivated analytical tools for predicting binary formation.We performed a total of 114 hydrodynamic simulations comprised of three suites, i) a fiducial suite with the expected AGN density for a thin Shakura-Sunyaev disc ii) a suite with ambient density three times higher and iii) a suite identical to the fiducial where accretion is turned off.Captures occur in each of our suites with varying success rates.We summarise our findings into key points below.
(i) The addition of the gaseous AGN discs leads to captures from a wider range of radial impact parameters  in the disc and leads to consistently close (Δ < 0.1 H ) encounters within one large window, unlike the gasless case which has only two very narrow, separate regions.This results from increased gravitational focusing in the lead up to the encounter from the mass of the CSMDs and strong dissipation from gas effects in the immediate lead up to the first periapsis.(ii) Binary scatterings in higher density environments form more easily in an even larger range of impact parameters.Additionally, they also allow easier formation of binaries that have larger close separations, owed to increased focusing from more massive CSMDs and stronger gaseous energy dissipation.(iii) In our simulations with accretion, dissipation during the encounter is initially dominated by the SMBH whose shear removes energy from the binary.Once the binaries get closer (Δ ≲ 0.5 H ) energy is then deposited into the binary from primarily the gas gravity and a little accretion as the CSMDs are perturbed by the opposing BH.
During the periapsis passage, when the CSMDs collide, accretion dissipation quickly switches sign and removes energy from the binary very efficiently, leaving it bound.(iv) When accretion is switched off, during the close approach gas can flow past the BHs where they would normally be accreted, leading to a net accumulation of gas behind the BH and inducing a drag.When accretion is enabled, this drag is instead reflected in the accretion dissipation where there is a direct 'headwind' on the BHs.We find little change in the trajectory of the BHs when accretion is disabled and similar energy dissipation, indicating accretion mimics the gas drag well on the small scales around the BHs.(v) A detailed understanding of the chronology and characteristics of the encounter was developed.Three characteristic orbital trajectories, determined by their impact parameter  were identified (see Figure 4) which demonstrate specific behaviours in their dissipation processes.
In the frame of the inner BH orbiting counterclockwise around the SMBH the trajectory of the three families are as follows: (a) Leftsided encounters -encounters with low  that pass to the left of the inner BH which form prograde binaries.(b) Rightsided encounters -encounters with moderate  values that pass to the right of the inner BH and form retrograde binaries (c) Turnaround encounters -encounters with high  that pass initially to the right before looping back on themselves and going round the left of the inner BH, forming prograde binaries.
We find that rightsided encounters are most favourable for successful binary formation as they are located in the centre of the space of  for successful binary formation.leftsided and Turnaround encounters impact parameters that place them further from the rightsided undergo encounters with increased minimum separations and are less likely to remain bound due to less accretion and gravitationally induced energy dissipation from the local gas.(vi) Strong oscillations in the dissipation from gas gravity, first identified in our last paper Paper I, are again observed in all our simulations where the binaries have deeper encounters due to non-cynlindrically symmetric perturbations to the CSMDs arising from tides induced by the opposing binary BH.These perturbations dominate the gas dissipation during the close approach in both simulations with and without accretion and accretion, though the net dissipation is found to be positive in the former and negative in the latter case.The strength of the oscillations in the binary energy could leave a detectable GW signature, however the orbital frequency of these inhomogeneities is of the order ∼ 10 −7  which puts them out of the low frequency sensitivity range of even LISA.(vii) We find the dissipation during the first encounter in all simulation suites follows a power-law with the minimum separation.This slope is found to be steeper when the AGN disk density is increased.(viii) We construct an analytic criterion for binary formation, eq. ( 34).
The relation is derived directly from the fully hydrodynamic simulations in this study yet only depends on pre-encounter properties.Thus, it is physically well informed but does not require the user to implement any such physics in order for it to operate, making it well suited to improving population studies such as (e.g.Tagawa et al. 2020a).Using the two-body energy at two Hill spheres separation and the impact parameter measured at one Hill sphere separation, the criterion can predict with a success rate >90% whether the binary will remain bound, for the AGN densities shown in our paper.The accuracy of the predictor is affirmed in Figure 19.

Figure 2 .
Figure 2. The change in the binary energy during the first encounter, Δ bin normalised to  H,c , as a function of the first periapsis passage  min,1 .The sign of Δ bin is indicated by the symbol type: circular green points indicate where energy is removed from the binary while the black triangular points represent models where energy is added.Only binaries that pass within  H are shown since all encounters outside the Hill sphere result only in flyby encounters.The red solid and dashed lines show respectively the best-fit power-law relation to the dotted data points and the 1 errors on the slope.

Figure 3 .
Figure 3. Similar to Figure 2 but showing the change in the binary energy during the second encounter, normalised to  H,c as a function of the closest approach during the second encounter  min,2 .Note that there are now naturally fewer data points as many models decouple after only their first encounter.

Figure 4 .
Figure 4. Classification of the trajectories of the first encounters between two sattelites moving on initially Keplerian circular orbits around a central SMBH.The trajectory of the outer satellite is shown in the frame of the inner satellite.A leftsided encounter is shown in red, rightsided in blue and turnaround in green.The dashed lines indicated the two failed encounters adjacent to each side of the capture window, to visualise how encounters far from the centre (far centre) of the capture window proceed.Their specific position on the impact parameter space is highlighted in Figure5.

Figure 5 .
Figure5.Cumulative energy change in all our binaries from the fiducial run during the first encounter from 2 H to apoapsis/exiting of Hill sphere.Results are colour-coded by encounter family; leftsided (red), rightsided (blue) and turnaround (green) encounters.The top panel re-illustrates which encounters belong to which family in the simulation suite.The dashed vertical lines indicate the two failed encounters adjacent to each side of the capture window and the bold lines highlight the three other selected models from each family used in Figure4using the same colour coding.The middle panel shows the total cumulative energy change for all fiducial models.The bottom row shows the cumulative energy transfer from the local gas gravity, accretion and SMBH respectively from left to right.We artificially extend the final cumulative value from the end point of the first encounter right up to the ten year mark to more easily compare the net cumulative energy dissipation across models.Additionally, we denote unsuccessful captures with dashed lines and successful captures with solid lines.From the colour-coded results, rightsided encounters lead to the most reliably efficient binary formation which we attribute to the encounter depth dependent dissipation of Figure2. .

Figure 7 .
Figure 7. Minimum separation for the first approach for our fiducial run (blue) and the same run where accretion is turned off when the BH satellites get within 2 Hill radii of each other (grey) as a function of the initial radial separation .Horizontal line indicates the Hill radius size of one of the objects.The outcome of all the simulations remains unchanged after switching off accretion.

Figure 8 .
Figure 8. Cumulative energy change change in all our binaries from the no-accretion run during the first encounter from 2 H to apoapsis/exiting of Hill sphere.Results are colour-coded by encounter family; leftsided (red), rightsided (blue) and turnaround (green) encounters.The top panel shows the total cumulative energy change for all fiducial models.The cumulative energy transfer from the local gas gravity and SMBH are shown respectively on the left and right of the second row.We artificially extend the final cumulative value from the end point of the first encounter right up to the ten year mark to more easily compare the net cumulative dissipation across models.Additionally, we denote unsuccessful captures with dashed lines and successful captures with solid lines.Turning off accretion leads to the replication of its effects in the dissipation from the local gas gravity.

Figure 9 .
Figure9.Central two panels: 2D average dissipation per area from gravitational interaction with the local gas (21) over the first encounter for a representative accretion-enabled (top) and accretionless (bottom) simulation with identical impact parameters.The map is centred on the centre of mass  CM of the binary and distances normalised to the current separation of the binary Δ.The dissipation is normalised to the max of both simulations to compare the relative strengths.The black lines show equal contours in each simulation, showing the enhanced dissipation close to the BHs in the accretionless case.Top panel: Average gravitational energy dissipation integrated from the  −  centre of mass radially outwards in the plane of the binary for the accretion-enabled case.Bottom panel: Same as the top panel but for the accretionless model.The radial and 2D maps are aligned so they can be directly compared along their axes.The binaries execute their orbit in a counter clockwise direction.

Figure 10 .
Figure10.A zoom-in on the cumulative energy change of the simulations with gas accretion being enabled and turned off, respectively, in Figure9around the first periapsis passage.Artificially turning accretion off leads to far stronger oscillations in the work done on the binary, due to the stronger gravitational tug of the gas minidiscs.

Figure 11 .
Figure11.Net energy change of the binaries over their first encounters, comparing the fiducial simulations where gas accretion is possible (left) to the simulations where accretion is turned off (right).The net contributions from different physical processes are shown with different colours.In the accretion-enabled panel, we also show the strength of  grav +  acc to compare with  grav in the accretionless simulations.When accretion is turned off, the sign of the net gravitational dissipation switches from positive to negative.Though accretion is the most efficient remover of energy, its omission only changes the dissipation by roughly a factor two.This indicates that it fairly accurately models the net dissipation from  acc contributions near the accreting boundary.In otherwords one may neglect accretion and retain fairly similar results qualitatively.

Figure 14 .
Figure 14.Net energy change of the satellites over their first encounters, comparing the fiducial simulations (left) to the 3 d,0 enhanced disc density simulations (right).Also shown are the net contributions from each dissipation mechanism.The enhancement of the dissipation from the increased ambient density is clear in the rightsided encounters (see Figure 5 for definition) in the central flat plateau of the encounter window, showing on average approximately three times greater energy dissipation.

Figure 15 .
Figure15.Energy change in the binaries, normalised to  H,c , during the first encounter as a function of the first periapsis  min,1 as in Figure2but for all simulations including those with an increased disk mass and those where accretion is turned off.Dotted points indicate where energy is removed to the binary while triangular points represent models where energy is added.Only binaries that pass within  H are shown.Also shown in the same colour as the raw data are the power-law fits for the fiducial  d,0 (blue), 3 d,0 , accretionless  d,0 and whole sample (black).
Figure17.First encounter periapsis scattered against the impact parameter at one Hill radius  1H (Eq.30) for all simulations, showing that smaller impact parameters lead to closer encounters, with some scatter for lower values of  1H where the very small periapsese becomes increasingly sensitive to small changes in energy.Also shown is the log-linear line of best fit (solid red) and its error (dashed red) according to eq. (32).

Table 2 .
1D cross sections  in the space of the impact parameter  as calculated by eq. ( tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc acc SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH SMBH H dt tot H dt tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot tot totH dt grav grav grav grav grav grav grav grav grav grav grav grav grav grav grav grav grav grav grav grav grav grav grav grav grav grav grav grav grav grav grav grav grav grav grav grav grav grav grav grav grav grav Encounter energy at two Hill radii,  2H

Table 3 .
Fit parameters for eq.34 for our  d =  d,0 ,  d = 3 d,0 and accretionless simulations, also shown is the average fit using all our simulations.
Encounter energy at two Hill radii,  2H , as in Figure