Prompt gravitational-wave mergers aided by gas in active galactic nuclei: the hydrodynamics of binary-single black hole scatterings

ABSTRACT

Black hole binary systems embedded in active galactic nucleus (AGN) discs have been proposed as a source of the observed gravitational waves (GWs) from LIGO-Virgo-KAGRA. Studies have indicated binary-single encounters could be common place within this population, yet we lack a comprehensive understanding of how the ambient gas affects the dynamics of these three-body encounters. We present the first hydrodynamical simulations of black hole binary-single encounters in an AGN disc. We find gas is a non-negligible component of binary-single interactions, leading to unique dynamics, including the formation of quasi-stable hierarchical triples. The gas efficiently and reliably dissipates the energy of the three-body system, hardening the triple provided it remains bound after the initial encounter. The hardening time-scale is shorter for higher ambient gas densities. Formed triples can be hardened reliably by |$2-3$| orders of magnitude relative to the initial binary semimajor axis within less than a few AGN orbits, limited only by our resolution. The gas hardening of the triple enhances the probability for a merger by a minimum factor of |$3.5-8$| depending on our assumptions. In several cases, two of the black holes can execute periapses of less than 10 Schwarzschild radii, where the dynamics were fully resolved for previous close approaches. Our results suggest that current time-scale estimates (without gas drag) for binary-single induced mergers are an upper bound. The shrinkage of the triple by gas has the prospect of increasing the chance for unique GW phenomena such as residual eccentricity, dephasing from a third object and double GW mergers.

gravitational waves, hydrodynamics, binaries: general, galaxies: nuclei, black hole mergers

1 INTRODUCTION

The astrophysical mechanisms driving compact object mergers and observed gravitational wave (GW) sources (Abbott et al. 2016, 2019, 2020a, b, c, d, 2022, 2023, 2023a, b; Venumadhav et al. 2020) remain an open question. Possible environments that would allow black holes (BHs) to reach small enough separations to merge within a Hubble time include isolated stellar binary evolution (e.g. Lipunov, Postnov & Prokhorov 1997; Belczynski et al. 2010b, 2016; Dominik et al. 2012, 2013, 2015; Giacobbo & Mapelli 2018; Tagawa & Umemura 2018), dynamical three-body or even four-body interactions in globular clusters (e.g. Portegies Zwart & McMillan 2000; Miller & Hamilton 2002; Mouri & Taniguchi 2002; Downing et al. 2010; Rodriguez et al. 2015, 2018; Rodriguez, Chatterjee & Rasio 2016; Samsing & D’Orazio 2018; Fragione & Kocsis 2019; Zevin et al. 2019; Di Carlo et al. 2020; Hamers & Samsing 2020; Arca Sedda, Li & Kocsis 2021; Liu & Lai 2021; Antonini et al. 2023; Samsing et al. 2024; Hendriks et al. 2024b) and in active galactic nuclei (AGNs; e.g. O’Leary, Kocsis & Loeb 2009; Bartos et al. 2017; Stone, Metzger & Haiman 2017; Hoang et al. 2018; Tagawa & Umemura 2018; Yang et al. 2019; Gröbner et al. 2020; McKernan et al. 2020; Rowan, Whitehead & Kocsis 2024a; Tagawa, Haiman & Kocsis 2020a; Tagawa et al. 2020b, 2021a, b; Ford & McKernan 2022; Delfavero et al. 2024; Rom, Sari & Lai 2024; Vaccaro et al. 2024; Wang, Zhu & Lin 2024).

AGNs present a unique environment due to the presence of a highly dense and thin accretion disc, in addition to a large spherical stellar component, i.e the nuclear stellar cluster (NSC). These two components play host to many different dynamical processes. In the NSC, we have two-body, binary-single, and binary-binary scatterings (e.g. Antonini & Perets 2012; Antonini & Rasio 2016; Trani, Quaini & Colpi 2024). Objects within the NSC that cross the AGN disc can become embedded through a combination of accretion/gas dynamical drag (e.g. Bartos et al. 2017; Panamarev et al. 2018; Yang et al. 2019; Fabj et al. 2020; Wang et al. 2024) and vector resonant relaxation (Rauch & Tremaine 1996; Kocsis & Tremaine 2015; Szölgyén & Kocsis 2018; Gruzinov, Levin & Zhu 2020; Magnan et al. 2022; Máthé, Szölgyén & Kocsis 2023). The disc-like geometry of this population, compared to a spherical component, may increase the frequency of the aforementioned BH interactions due to the enhanced number density of objects and through disc migration (e.g. Secunda et al. 2019, 2021). In addition to the same kind of encounters within the NSC, the high gas density of the disc facilitates the efficient formation of black hole binaries (BBHs) through two-body scatterings by dissipating the relative two-body energy of the BHs through gas drag, known as ‘gas-captured’ binaries (Tagawa et al. 2020a). This BBH hydrodynamical formation mechanism has been studied through hydrodynamical simulations of the global disc (e.g. Rowan et al. 2023, 2024b), shearing box simulations (e.g. Li et al. 2023; Whitehead et al. 2024a, b) and N-body simulations using analytical gas drag treatments (e.g. DeLaurentiis, Epstein-Martin & Haiman 2023; Rozner, Generozov & Perets 2023; Dodici & Tremaine 2024).

In this paper, we examine the influence of gas on binary-single interactions which are expected to be common during the evolution of disc-embedded binaries if the typical gas-driven inspiral time-scale is greater than the collisional time-scale for three-body interactions (as suggested by Tagawa et al. 2020a). Binary-single interactions play a dominant role in shaping the merging mass distribution and the spin–orbit angle distribution due to exchange interactions and by reorientating the binaries’ orbital plane. The time-scale for the inspiral of embedded BBHs is still highly uncertain, as the spatial scales involved to accurately model the hydrodynamics from the BBH’s formation to merger can be over seven orders of magnitude, with estimates between |$10^{2}\,\mathrm{ and}\,10^{6}$| yr (e.g. Baruteau, Cuadra & Lin 2011; Dittmann & Cantiello 2025; Ishibashi & Gröbner 2024; Vaccaro et al. 2024). There is a growing amount of evidence that retrograde binaries (orbiting in the opposite direction with respect to the orbit around the AGNs) experience eccentricity pumping, which may accelerate the merger time-scale (e.g. Dittmann & Ryan 2022; Rowan et al. 2023; Calcino et al. 2024). The expected merger and eccentricity pumping time-scale scales inversely with the density of the AGN disc (e.g. Baruteau et al. 2011; Ishibashi & Gröbner 2020; Li, Lai & Rodet 2022b; Rowan et al. 2023). Therefore, given the immense range in expected AGN disc densities across the various initial parameters in AGN disc models, this can lead to highly varying time-scales regardless of the method used to estimate the binary torques.

The dynamics of three body scatterings are inherently chaotic (e.g. Aarseth & Zare 1974; Heggie 1974; Hut 1983; Samsing & Ilan 2018), yet three (or more) body systems are abundant throughout the Universe. Examples include asteroids (Kozai 1962; Lidov 1962), planetary binaries (Nagasawa, Ida & Bessho 2008; Perets & Naoz 2009), stellar triples (Eggleton & Kiseleva-Eggleton 2001; Duchêne & Kraus 2013; Toonen, Hamers & Portegies Zwart 2016), stellar compact object triples (Thompson 2011), black hole binaries in AGNs (e.g. Leigh et al. 2018; Trani et al. 2019, 2024; Ginat & Perets 2021; Samsing et al. 2022; Fabj & Samsing 2024; Rom et al. 2024) and supermassive black hole (SMBH) systems (Blaes, Lee & Socrates 2002).

The stability of a three body system depends on the relative separations between the objects. In a hierarchical triple, one object exists at a large separation from a more tightly bound binary system. Provided the separation is sufficiently large, the outer object may only weakly perturb the inner binary orbit and vice versa over secular time-scales (e.g. Ford, Kozinsky & Rasio 2000; Toonen et al. 2016). If there is sufficient inclination between the inner and outer binary orbit, then secular oscillations in the eccentricity and inclination driven by the Kozai–Lidov (Kozai 1962; Lidov 1962) oscillations could alter the periapsis of the inner binary, speeding up the time-scale for merger in the case of BHs (e.g. Naoz et al. 2013; Antonini, Murray & Mikkola 2014; Naoz 2016; Silsbee & Tremaine 2017). If the separations are comparable, the evolution is non-secular and the objects continually exchange energy and angular momentum, until one object is ejected from the system (in the absence of any dissipation mechanism). Fabj & Samsing (2024) frame this case as a sequence of binary-single like states, where the eccentricities and semimajor axes may be sampled from probability distributions each time the less bound third object executes its periapsis. The continual ‘re-shuffling’ of the binary eccentricity during the evolution of non-hierarchical BH systems has profound implications for the expected eccentricity distributions of GW signals, where we expect binaries may be able to merge with residual eccentricity through dynamical BH triple interactions in AGNs (e.g. Ishibashi & Gröbner 2020) and globular clusters (e.g. Samsing, MacLeod & Ramirez-Ruiz 2014; Dall’Amico et al. 2024). Constraining the rates of these highly eccentric mergers from both theory and GW observations will shed light on the proportion of BH mergers coming from dynamical channels versus isolated binary evolution.

Although analytical and numerical N-body studies present exciting prospects for unique observational GW features and possibly shorter BH merger time-scales in AGNs, they currently neglect the AGN disc component. As indicated by the efficiency of the gas-capture mechanism during single-single interactions, the gas may also significantly affect the dynamics during a binary-single encounter in the AGN disc. In this work, we perform the first simulations of binary-single encounters in a gaseous accretion disc that include the effects of the gas, using a hydrodynamical approach. We examine how gaseous drag alters the dynamics of the initial encounter as well as the evolution of formed black hole triple systems. In later sections, we calculate the GW merger probabilities and their enhancement when gas is included to quantify the significance of gaseous effects. We describe the simulations and their initial conditions in Section 2, before presenting our results in Section 3. The implications of our results are discussed in Section 4. We summarize and conclude our results in Section 5.

2 METHODS

This work extends our previous studies of single-single black hole encounters (Rowan et al. 2023, 2024b; Whitehead et al. 2024a, b) in AGNs to the binary-single case. We simulate the encounters and the hydrodynamics using the Eulerian GRMHD code Athena++ (Stone et al. 2020). We utilize a second-order accurate van Leer predictor–corrector integrator with a piecewise linear method (PLM) spatial reconstruction and Roe’s linearised Riemann solver. See Stone et al. (2020) for a detailed discussion of the integration scheme used in Athena++. We simulate a 2D portion of the AGN disc in a rectangular shearing box (e.g. Goldreich & Tremaine 1978; Hawley, Gammie & Balbus 1994, 1995). The important length-scales for the system are the single, binary, and three-body Hill spheres: |$r_{\rm H,i}, r_{\rm H,b}$|⁠, and |$r_{\rm H,t}$|⁠, respectively:

$$\begin{eqnarray} r_{\rm H,i} =\, R_0\bigg (\frac{M_{\rm i}}{3M_\bullet }\bigg)^{1/3} , \end{eqnarray}$$

(1)

$$\begin{eqnarray} r_{\rm H,b}=\, R_0\bigg (\frac{M_\mathrm{bin}}{3M_\bullet }\bigg)^{1/3} , \end{eqnarray}$$

(2)

$$\begin{eqnarray} r_{\rm H,t}=\, R_0\bigg (\frac{M_\text{trip}}{3M_\bullet }\bigg)^{1/3} , \end{eqnarray}$$

(3)

where |$M_{\rm bin} = M_1+M_2$| is the mass of the binary, |$M_{\rm trip} = M_\text{bin}+M_3$| is the mass of the triple, |$M_\bullet$| is the supermassive black hole (SMBH) mass, and |$R_0$| is the radial position of the shearing box from the SMBH. We define |$i=(1,2)$| as the BHs present in the initial (pre-encounter) binary and |$i=3$| as the incoming single BH.

2.1 Nomenclature

Concerning the language of triple encounters, use of the term close triple encounter will specifically refer to an encounter where the separations of the three objects are comparable and within |$r_{\rm H,t}$|⁠, during the three-body interaction (which may occur multiple times). During a triple encounter, the terms single or BH pertain to the currently least bound object to the three-body centre of mass (COM). The terms binary or BBH refer to the currently most energetically bound pair of BHs. The semimajor axes, eccentricities, and periapses for the binary and single components during the three-body scattering are denoted |$a_\mathrm{bin}$|⁠, |$e_\mathrm{bin}$|⁠, and |$r_\mathrm{p,bin}$| for the binary and |$a_\mathrm{sin}$|⁠, |$e_\mathrm{sin}$|⁠, and |$r_\mathrm{p,sin}$|⁠, respectively. As utilized in Sections 2.5 and 3.7, the binary quantities are calculated assuming purely 2-body dynamics and the single quantities are similarly calculated assuming the binary is a point mass with mass |$M_i+M_j$| and |$i,j\in \lbrace 1,2,3\rbrace$|⁠.

2.2 The shearing box

We utilize a comparable shearing box setup to Whitehead et al. (2024a, b), where the encounter takes place in a non-inertial reference frame that co-rotates with the AGN disc at a fixed radius |$R_0$| and angular frequency |$\Omega _0=\sqrt{GM_\bullet /R_0^{3}}$|⁠. The Cartesian coordinate system of the shearing box |$\lbrace x,y\rbrace$| can be translated to a position in the global AGN disc in cylindrical coordinates |$\lbrace R,\phi \rbrace$| as

$$\begin{eqnarray} \boldsymbol{r}=\bigg ({\begin{array}{*{10}c}R \\ \phi \end{array}}\bigg) =\bigg ({\begin{array}{*{10}c}R_0 + x \\ \Omega _0 t+y/R_0 \end{array}}\bigg). \end{eqnarray}$$

(4)

The gas and BHs are subject to accelerations from the Coriolis and centripetal forces. These can be added together as

$$\begin{eqnarray} \boldsymbol{a}_{\text{SMBH}}= 2 \boldsymbol{u} \times \Omega _0 \boldsymbol{\hat{z}} + 2q\Omega _0^2 (x-x_\mathrm{C})\hat{x}, \end{eqnarray}$$

(5)

where |$\boldsymbol{u}$| represents here the gas or BH velocity in the corotating frame, |$x_\mathrm{C}$| is the x position of co-rotation and |$q=-\frac{\mathrm{d}\ln \Omega }{\mathrm{d}\ln R}=\frac{3}{2}$| is the shear parameter for a Keplerian disc. Thus, the equilibrium trajectories where equation (5) is zero are

$$\begin{eqnarray} \boldsymbol{u}_\mathrm{eq}=\bigg ({\begin{array}{*{10}c}u_{\rm x} \\ u_{\rm y} \end{array}}\bigg) =\bigg ({\begin{array}{*{10}c}0 \\ -q\Omega _0 (x-x_\mathrm{C}) \end{array}}\bigg)\, . \end{eqnarray}$$

(6)

These trajectories represent Keplerian circular orbits of varying radii in the global picture, but straight lines in the tangential |$\boldsymbol{\hat{y}}$| direction in the shearing frame.

The shearing box has a radial width of |$0.24R_{0}$| and azimuthal extent |$0.96R_{0}$|¹, corresponding to |$\sim 18r_{\rm H,t}$| and |$\sim 72r_{\rm H,t}$|⁠, respectively (see Section 2.5.2). Defining |$x_{\rm C}$| as the x position of the shear in the box about the midpoint |$(x=0)$|⁠, the boundary conditions at the upper (⁠|$y=y_\mathrm{max}$|⁠) and lower edge (⁠|$y=y_\mathrm{min}$|⁠) of the box are:

$$\begin{eqnarray} (y=y_\mathrm{max}):\, \, \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\text{outflow} & x< x_\mathrm{C},\\ \text{refill} & x>x_\mathrm{C},\\ \end{array}\right. \end{eqnarray}$$

(7)

$$\begin{eqnarray} (y=y_\mathrm{min}):\, \, \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\text{refill} & x< x_\mathrm{C},\\ \text{outflow} & x>x_\mathrm{C}.\\ \end{array}\right. \end{eqnarray}$$

(8)

The boundary conditions at the x boundaries are set to outflow, where the flow in the ghost cells match the flow properties at the boundary (i.e. a zero-gradient boundary). The refill regions assume the gas velocities and densities at the start of ambient flow, i.e |$\lbrace u_\mathrm{x},u_\mathrm{y}\rbrace =\lbrace 0,-q\Omega _0 (x-x_\mathrm{C})\rbrace$| and |$\Sigma =\Sigma _{0}$| (see Section 2.5.1), such that there is a continual inflow of gas from the refilling ghost cells.

2.3 Gas dynamics

Athena++ solves the fluid equations in Eulerian form through the extended Navier–Stokes equations, in 2D:

$$\begin{eqnarray} \frac{\partial \Sigma }{\partial t} + \nabla \cdot \left(\Sigma \boldsymbol{u}\right) = 0, \end{eqnarray}$$

(9)

$$\begin{eqnarray} \frac{\partial \left(\Sigma \boldsymbol{u}\right)}{\partial t} + \nabla \cdot \left(\Sigma \boldsymbol{u} \boldsymbol{u} + P \boldsymbol{I} + \boldsymbol{\Pi }\right) = \Sigma \left(\boldsymbol{a}_{\text{SMBH}}-\nabla \phi _\text{BH}\right)\, . \end{eqnarray}$$

(10)

Here, |$\Sigma$|⁠, |$\boldsymbol{u}$|⁠, P, |$\Pi$|⁠, are the cell gas surface density, velocity, pressure, and viscous stress tensor

$$\begin{eqnarray} \Pi _{ij} = \Sigma \nu \left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}-\frac{2}{3}\delta _{ij}\nabla \cdot \boldsymbol{u}\right), \end{eqnarray}$$

(11)

for a viscosity |$\nu$|⁠. The viscosity is modelled using the standard |$\alpha$|-disc approach, where |$\nu =\alpha c_{\rm s}H$| and |$H=c_{\rm s}/\Omega$| is the disc scale height. In all simulations here, |$\alpha =0.1$|⁠. We assume an isothermal equation of state where |$P=\Sigma c_{\rm s}^{2}$| and the sound speed |$c_{\rm s}$| is fixed. We ignore any magnetic effects and the self-gravity of the gas.

The remaining |$\nabla \phi _{\rm BH}$| term is the acceleration from the potential of all (⁠|$n_\mathrm{BH}=3$|⁠) stellar mass BHs,

$$\begin{eqnarray} -\nabla \phi _\text{BH}(\boldsymbol{r}) = \sum ^{n_\text{BH}}_{i=1} M_i \, g\left(\frac{\boldsymbol{r}-\boldsymbol{r}_n}{h}\right), \end{eqnarray}$$

(12)

where |$g(s)$| is the gas gravitational softening kernel, see Price & Monaghan (2007)

$$\begin{eqnarray} g(\boldsymbol{s}) = -\frac{G}{h^2}\hat{\boldsymbol{s}} \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\frac{32}{3}s-\frac{192}{5}s^3 + 32s^4 & 0 < s \le \frac{1}{2} \\ -\frac{1}{15s^2} + \frac{64}{3}s-48s^2 + \frac{192}{5}s^3-\frac{32}{3}s^4 & \frac{1}{2} < s \le 1 \\ \frac{1}{s^2} & s >1\, . \end{array}\right.\\ \end{eqnarray}$$

(13)

As three-body encounters can often lead to chaotic trajectories, we cannot rule out very close BH–BH separations. Therefore, we set the softening length h to be smaller than our previous work at |$h_{i}=0.005r_{\rm H,i}$|⁠. This allows us to trust the dynamics of these close encounters as the enclosed mass within |$r< h$| (where gas gravitation is not accurately described) is very low and unlikely to alter the trajectories of the BHs significantly should they execute trajectories within h. Note that only the gas–BH interactions are softened, but the BH–BH interactions are not (see Section 2.5.2).

2.4 Mesh refinement

We apply an adaptive mesh refinement (AMR) procedure, allowing us to resolve the area around the BHs to a high degree, whilst minimizing compute time and resources. The mesh closer to the location of each BH becomes more refined, down to a minimum refinement level. In all simulations shown here, we maintain a base resolution of 256×1024, with 8 mesh refinement levels. For our box size, this gives a maximal and minimal cell size of |$\delta _\mathrm{max}\simeq 0.1r_{\rm H,t}$| and |$\delta _\mathrm{min}\simeq 0.00039r_{\rm H,t}$|⁠, The softening length is then resolved by |$h/\delta _\mathrm{min}\gtrsim 10$| cells across one dimension. The AMR scheme is centred on the positions of the BHs and moves with them throughout the simulation, ensuring the necessary resolution around the BHs at all times. We set the spatial scales for minimal/maximal refinement to 0.2 and 0.5 |$r_\mathrm{H,i}$|⁠, respectively, i.e the domain is maximally refined within 0.2|$r_\mathrm{H,i}$|⁠, see Whitehead et al. (2024a) for more detail.

2.5 Initial conditions

2.5.1 The AGN disc

To model an accurate AGN disc, we determine the ambient surface density |$\Sigma _{0}$| and sound speed |$c_{\rm s,0}$| of the gas in the shearing box from AGN disc profiles generated using pAGN (Gangardt et al. 2024). In this paper, we consider a fiducial setup assuming an AGN disc with an Eddington fraction |$L_\mathrm{\epsilon }=0.1$|⁠, radiative efficiency |$\epsilon =0.1$|⁠, and hydrogen/helium fractions X/Y = 0.7/0.3. We consider an AGN with a SMBH mass at the peak of the anticipated merger rate distribution for gas-captured binaries predicted in Rowan et al. (2024a), where |$M_\bullet =10^{7}\mathrm{M}_{\odot }$|⁠. The shearing box radius is set to |$R_0=0.01$|pc. For these parameters, this gives |$\Sigma _0\simeq 1.6\times 10^{5}$| kg m|$^{-2}$|⁠, |$c_{\rm s,0}\simeq 11.9$| km s|$^{-1}$| and disc thickness ratio |$H/R\simeq 0.0057$|⁠. We investigate how the AGN disc density affects the triple dynamics by considering two additional values such that we have three disc densities |$\Sigma =\lbrace 0.1,0.25,1\rbrace \Sigma _0$|⁠, where we take |$\Sigma =\Sigma _0$| to be our fiducial value.

2.5.2 The black holes

For simplicity, we consider three representative equal mass BHs (⁠|$M_1=M_2=M_3=25\mathrm{M}_{\odot }$|⁠). This gives the mass scales for a single BH and triple system: |$q_1=M_i/M_\bullet =2.5\times 10^{-6}$| and |$q_3 = \sum _\mathrm{i=1}^{3}M_i/M_\bullet =7.5\times 10^{-6}$|⁠, respectively. The BHs are represented by non-accreting point-like particles. The motion of these BHs is determined purely by gravity through direct summation of the forces between the BH sinks, background SMBH frame forces (equation (5)), and the gas in each cell. Like the gas, the gravitational accelerations due to cells within h of a BH is softened according to equation (13). The BH–BH interactions remain unsoftened. The gravitation of the BHs is purely Newtonian and hence dissipation from GWs is not considered explicitly, though we discuss its implication in Section 3.7.

2.5.3 The encounter

The binary and single system are initialized with zero inclination and marginal eccentricity of |$e_0=0.02$|⁠, as a result all scatterings will also be co-planar. In this 2D geometry the scattering is parametrized by the radial impact parameter |$p=x_\mathrm{s}-x_\mathrm{b}$|⁠, where |$x_\mathrm{s}$| and |$x_\mathrm{b}$| are the x positions of the single and the binary COM, respectively, at the moment the single object is introduced. The initial binary semimajor axis |$a_0$| is set to |$0.1r_\mathrm{H,b}$|⁠. We fix the semimajor axis of the binary at the encounter and therefore turn off the gravitation of the BHs due to the gas until they reach a separation of |$2r_{\rm H,t}$|⁠. Not doing so would introduce variable hardening rates on the binary for the different initial gas densities considered here, leading to a correlation between the assumed ambient density and the initial conditions (binary semimajor axis) of the encounter, complicating the analysis.

In our previous shearing box simulations of single-single BH scatterings (Whitehead et al. 2024a, b), the time required for the gas to form a steady morphology around each object was short enough to permit their construction naturally within a reasonable azimuthal extent of the shearing box. As we are now dealing with a binary, the time required for the gas around the binary to form a steady state is longer. Therefore, we initialize only the binary in the simulation for |$\sim 50$| binary orbits before inserting the single BH, to allow the formation of a well defined and steady circumbinary discs. After this settling time, we initialize the single with impact parameter p and an azimuthal displacement that is normalized such that the approach time for each encounter is the same as that for |$p=1.5r_{\rm H,t}$| and an initial azimuthal distance |$\Delta y=30r_{\rm H,t}$|⁠, i.e |$\Delta y/r_{\rm H,t}=30 p/1.5r_{\rm H,t}$|⁠. Normalizing the azimuthal distance to maintain the same approach time ensures the binary and single have the same amount of time to accumulate mass in their Hill spheres and develop their wakes. These initial conditions allow the binary |$\sim 120$| orbits to fill its Hill sphere and reach a steady morphology before the encounter with the single.

We perform simulations with 24 different impact parameters in the range |$p=[1.2,1.9]r_{\rm H,t}$| for each value of the density |$\Sigma$| for a total of 72 simulations. An initial linear sweep over the full range was performed with 15 simulations. We then doubled the sampling uniformly in the parameter range, which leads to chaotic triple encounters (⁠|$p=[1.375,1.775]r_{\rm H,t}$|⁠, see Section 3.6) with a further 9 simulations.

2.5.4 Shifting the shearing frame

At the beginning of each simulation, we anticipate the location of the scattering using the three-body COM. We initialize the initial binary’s position such that the three-body COM lies in the origin when the single is inserted. In this work, we consider the scattering of a binary on an inner orbit compared to the single, therefore the binary is located in a position towards the lower left of the simulation domain. To allow the gas around the binary to form a steady morphology before inserting the single, we fix the co-rotating position |$x_\mathrm{C}$| initially on the binary (i.e |$x_\mathrm{C,0}=-2p/3$|⁠) so it does not have any motion in the shearing box. After five AGN orbital time-scales (⁠|$5\Omega _0^{-1}$|⁠), we insert the single and shift |$x_\mathrm{C}$| to the origin and boost the velocity of all sinks and gas cells to match the new co-rotation frame with the velocity boost given by |$\Delta v_\mathrm{y}=q\Omega _0 x_\mathrm{C,0}=-\frac{2}{3}q\Omega _0 p$|⁠. The single’s mass is initially zero and linearly increased to its final value (25|$M_\odot$|⁠) over a period of 0.1|$\Omega _0^{-1}$| for stability.

2.5.5 Reaching steady state

The most influential parameter for the dissipation during embedded object encounters is the mass within each object’s Hill sphere |$m_\mathrm{H}$|⁠. If the encounter occurs before the mass has time to build-up and reach a steady value, the energy/angular momentum exchange between the objects and the gas will typically be lower. In Fig. 1, we show the gas mass |$m_\mathrm{H}$| and angular momentum |$L_\mathrm{H}$| enclosed within the Hill sphere of both the binary and the single prior to the encounter. The binary Hill mass and angular momentum reaches a steady state after approximately 50 inner binary orbits (⁠|$P_\mathrm{bin}\approx 0.12\Omega _0^{-1}$|⁠), remaining approximately flat until the encounter at |$120P_\mathrm{bin}$|⁠. The gas within the Hill sphere of the single reaches a steady state after roughly |$40P_\mathrm{bin}$|⁠. As the single nears the binary (⁠|$\sim 100P_\mathrm{bin}$|⁠), the gas flow into their respective Hill spheres changes due to interaction with the undersense regions associated with the gas morphology around each perturber, leading to a reduction in |$m_\mathrm{H}$| prior to encounter. The approach time for the encounter was set to the minimum that allows |$m_\mathrm{H}$| in the single to reach a steady state before this reduction takes place.

$Time evolution of gas properties in the vicinity of the binary and single prior to encounter for the simulation with $p=1.7r_{\rm H,3}$ and $\Sigma =0.1\Sigma _0$. Top: the gas mass $m_\mathrm{H}$ contained within the Hill sphere of the binary ($r< r_{\rm H,12}$) and single ($r< r_{\rm H,3}$), normalized to the mass initially enclosed within their respective Hill radii $m_\mathrm{H,0}$. Bottom: the angular momentum of the enclosed gas $L_\mathrm{H}$ (as measured from the single and the binary COM) normalized to the angular momentum of the binary $L_\mathrm{bin}$. Both quantities are shown as a function of time in inner binary orbits $P_\mathrm{bin}$. The binary Hill mass reaches a steady state after $50 P_\mathrm{bin}$ and the single after $40 P_\mathrm{bin}$, respectively.$

Figure 1.

Time evolution of gas properties in the vicinity of the binary and single prior to encounter for the simulation with |$p=1.7r_{\rm H,3}$| and |$\Sigma =0.1\Sigma _0$|⁠. Top: the gas mass |$m_\mathrm{H}$| contained within the Hill sphere of the binary (⁠|$r< r_{\rm H,12}$|⁠) and single (⁠|$r< r_{\rm H,3}$|⁠), normalized to the mass initially enclosed within their respective Hill radii |$m_\mathrm{H,0}$|⁠. Bottom: the angular momentum of the enclosed gas |$L_\mathrm{H}$| (as measured from the single and the binary COM) normalized to the angular momentum of the binary |$L_\mathrm{bin}$|⁠. Both quantities are shown as a function of time in inner binary orbits |$P_\mathrm{bin}$|⁠. The binary Hill mass reaches a steady state after |$50 P_\mathrm{bin}$| and the single after |$40 P_\mathrm{bin}$|⁠, respectively.

Open in new tab Download slide

A visual timeline of one of our simulations is shown in Fig. 2 through visualizations of the surface density |$\Sigma$| in the full extent of the shearing box. All density visualizations in this work are normalized to the initial ambient value |$\Sigma _{\infty }$|⁠. In the figure, we show the initial fixed position of the binary as it begins to form a circum-binary disc, followed by the insertion of the single, the approach of the binary and single, the encounter and a close-up of the encounter. We also compare the resolution to the scale of the binary, demonstrating the binary is well resolved by |$\sim 250$| cells over its semimajor axis. We note that the streamers of the binary only settle their final configuration once the circum-binary disc mass/angular momentum plateaus. The remnants of the unsettled streamers are present in panel III of Fig. 2 and quickly exit the domain well before the encounter.

$The timeline of a binary-single scattering simulation ($\Sigma =0.1\Sigma _0$, $p=1.4r_{\rm H,t}$), visualized through the surface density $\Sigma$ normalized to the initial ambient value $\Sigma _\mathrm{\infty }$. Panel I: the binary is initialized in its offset position with the shear centred on its COM. Panel II: the single is injected (green dot) and the radial location where the shear vanishes $x_\mathrm{C}$ is shifted to the origin, which is the COM of the three-body system. Panel III: the single forms a circum-single disc as both the single and binary approach the encounter location at the origin. Panel IV: With both components now in steady state, they enter the underdense gas region just prior to encounter (see Section 2.5.5). Panel V: the system undergoes a close encounter. V*: A zoom in of the close encounter as the colliding gas forms a shock. The small grey line indicates the length-scale of 100 times the smallest cell size $\delta _\mathrm{min}$.$

Figure 2.

The timeline of a binary-single scattering simulation (⁠|$\Sigma =0.1\Sigma _0$|⁠, |$p=1.4r_{\rm H,t}$|⁠), visualized through the surface density |$\Sigma$| normalized to the initial ambient value |$\Sigma _\mathrm{\infty }$|⁠. Panel I: the binary is initialized in its offset position with the shear centred on its COM. Panel II: the single is injected (green dot) and the radial location where the shear vanishes |$x_\mathrm{C}$| is shifted to the origin, which is the COM of the three-body system. Panel III: the single forms a circum-single disc as both the single and binary approach the encounter location at the origin. Panel IV: With both components now in steady state, they enter the underdense gas region just prior to encounter (see Section 2.5.5). Panel V: the system undergoes a close encounter. V*: A zoom in of the close encounter as the colliding gas forms a shock. The small grey line indicates the length-scale of 100 times the smallest cell size |$\delta _\mathrm{min}$|⁠.

Open in new tab Download slide

3 RESULTS

3.1 Fiducial example

As a qualitative example, we show the trajectories and separations of the three BHs in addition to the gas morphology at four intervals from one of our simulations (⁠|$p=1.7$|⁠, |$\Sigma =0.1\Sigma _0$|⁠) in Fig. 3. The first panel shows the initial shock formed when the gas in the binary and single Hill spheres come into contact. Rapidly after this point, both components lose their initially well-defined discs and the morphology becomes highly complex and variable. In this example, the interaction leads to a binary swap (⁠|$t\sim 15.1\Omega _0^{-1}$|⁠), visible in the separation between each of the three BHs. The bottom left panel shows the formation of the temporary binary system formed when the single swaps place with one of the other BHs. This is then reversed around |$t=15.6\Omega _0^{-1}$|⁠, where another close triple encounter leads to the re-pairing of the original binary BHs and the third is ejected, leaving the binary hardened compared to both the original and the temporary binary formed after the first encounter. While the three-body dynamics play out, the gas leads to a net hardening of the temporary binary and also continues to harden the final binary after the third is ejected. This fiducial example demonstrates simultaneous hardening of a BH through a three-body exchange and gas drag.

$An example encounter with $\Sigma =0.1\Sigma _0$ and $p=1.7r_{\rm H,t}$. Top left: the gas morphology at the moment the minidiscs of the binary and single collide, the colours represent the surface density relative to the ambient value of the simulation $\Sigma _{\infty }$. Top centre: the morphology of the gas as the single begins plunging into the binary. Bottom left: the gas morphology following the first major encounter of the third object. Here, one binary object has been replaced by the single via an exchange interaction. Bottom centre: the gas morphology after one object has been removed from the system, hardening the remaining binary. Top right: the trajectories of the three BHs in the three-body COM. Bottom right: the separation between each BH $\Delta r_{ij}$ as a function of time. The vertical dashed lines represent the timestamps of the $\Sigma$ plots.$

Figure 3.

An example encounter with |$\Sigma =0.1\Sigma _0$| and |$p=1.7r_{\rm H,t}$|⁠. Top left: the gas morphology at the moment the minidiscs of the binary and single collide, the colours represent the surface density relative to the ambient value of the simulation |$\Sigma _{\infty }$|⁠. Top centre: the morphology of the gas as the single begins plunging into the binary. Bottom left: the gas morphology following the first major encounter of the third object. Here, one binary object has been replaced by the single via an exchange interaction. Bottom centre: the gas morphology after one object has been removed from the system, hardening the remaining binary. Top right: the trajectories of the three BHs in the three-body COM. Bottom right: the separation between each BH |$\Delta r_{ij}$| as a function of time. The vertical dashed lines represent the timestamps of the |$\Sigma$| plots.

Open in new tab Download slide

3.2 Dissipation

The closest analogue to the binary-single encounters in gas is the single-single gas-capture mechanism. During such an encounter, the two-body energy of the BHs is dissipated to the gas through a combination of accretion (direct linear momentum transfer) or gas dynamical drag (dynamical linear momentum transfer). The amount of energy that can be transferred is approximately proportional to the gas contained within the Hill sphere of each approaching object (e.g. Whitehead et al. 2024a; Rowan et al. 2024b) as this effectively serves as the available reservoir to deposit the energy of the BHs. For a three body system, the energy of BHs is given by the three body energy

$$\begin{eqnarray} E_{\rm trip} = \sum _{i}^{3}\frac{1}{2}M_i v_i^{2}-\sum _{j}^{3}\sum _{i>j}^{3}\frac{G M_i M_j}{r_{i,j}}, \end{eqnarray}$$

(14)

where |$v_i$| is the velocity of BH i relative to the three-body COM and |$r_{i,j}=\Vert \boldsymbol{r}_i-\boldsymbol{r}_j\Vert$| is the separation of BHs i and j. In the two_body case with BHs i, j the energy is given by

$$\begin{eqnarray} E_{\rm bin} = \frac{1}{2}\mu \Vert v_i-v_j\Vert ^{2}-\frac{G M_i M_j}{r_{i,j}}, \end{eqnarray}$$

(15)

where |$\mu =M_1 M_2/M_\text{bin}$| is the reduced mass of the binary.

We show how |$E_\mathrm{trip}$| and |$E_\text{bin}$| for the most bound pair of BHs changes throughout the fiducial model in Fig. 4. We define a natural energy scale for the system |$E_\mathrm{H,3}$| as the absolute energy of a two-body system with masses |$M_3$| and |$M_1+M_2$|⁠, and a separation of |$r_\text{H,t}$|⁠, i.e.

$The evolution of the three-body energy $E_{\rm trip}$ in the fiducial simulation. Top: the three body energy of the system (equation (14)) and energy of the most bound BH pair (equation (15)) as a function of time in units of $\Omega _{0}^{-1}$. Bottom: the separation between each of the three objects $\Delta r_{ij}$. The vertical grey lines border the time period where all three objects are within $r_\mathrm{H,t}$.$

Figure 4.

The evolution of the three-body energy |$E_{\rm trip}$| in the fiducial simulation. Top: the three body energy of the system (equation (14)) and energy of the most bound BH pair (equation (15)) as a function of time in units of |$\Omega _{0}^{-1}$|⁠. Bottom: the separation between each of the three objects |$\Delta r_{ij}$|⁠. The vertical grey lines border the time period where all three objects are within |$r_\mathrm{H,t}$|⁠.

Open in new tab Download slide

$$\begin{eqnarray} E_\mathrm{H,3}= \frac{G(M_1+M_2)M_3}{2r_\mathrm{H,t}} . \end{eqnarray}$$

(16)

At the first encounter, a large amount of energy is dissipated from the system by the gas. This initial large energy loss is consistent with single-single scatterings, owed to the large gas mass contained in the Hill sphere of the binary and single that they accumulated before the encounter. Following the first encounter and the binary object exchange, the newly formed triple system is hardened slightly (visible in the more gradual decline in |$E_\mathrm{trip}$|⁠) before another close three-body interaction occurs and one BH is ejected. At the second close encounter, there is another spike in the dissipation. We attribute this to the single BH re-accumulating gas while executing its large apoapsis before the second close encounter. The remaining binary is hardened (observe the second gradual decline in |$E_\mathrm{bin}$| and |$E_\mathrm{trip}$| following the second encounter). When the ejected BH leaves |$r_\mathrm{H,t}$| the SMBH potential cannot be ignored and |$E_\mathrm{trip}$| becomes unreliable as a metric for the system’s hardness.

3.3 Alternate encounter outcomes

The dynamics of the encounters and their outcomes can be divided into characteristic families.

Glancing encounters – where the minimum separation of the third body remains of the order of or larger than |$r_{\rm H,t}$| and only a weak interaction takes place.
Hierarchical encounters – when a single approaches and is successfully captured through gas into a hierarchical system where the single orbits the original binary and evolves in a secular fashion.
Temporary chaotic encounters – where the periapsis of the single passes close to the binary and gas dissipation removes an insufficient amount of energy, leading to a chaotic three-body interaction which ends with one object being ejected.
Hardened chaotic encounters – where the periapsis of the single passes close to the binary and gas dissipation is highly efficient, leading to a prolonged three-body encounter as gas continually removes energy from the system.

We show examples of the trajectories and object separations over time for each encounter type in Fig. 5. In the glancing encounter case, the BBH and BH do not pass deeply into each other’s accretion discs and have typically high approach velocities, therefore dissipation is not efficient and the single does not become bound. As the periapsis of the encounter is large, the BBH is only mildly perturbed. In Fig. 5, one can observe the encounter induces an eccentricity change from |$e\simeq 0.05$| to |$e\simeq 0.14$| in the BBH. In the hierarchical encounter, the periapsis of the encounter is still large compared to the BBH semimajor axis, but the single is successfully captured through gas drag. Provided the single has sufficient angular momentum that subsequent encounters have periapses greater than the BBH semimajor axis, the single BH remains stable in its orbit about the BBH, forming a quasi-stable hierarchical triple system (see panel 2 in Fig. 5).

$The trajectories and object separations for the four characteristic types of triple encounters in gas. Left to right: Glancing encounter, hierarchical encounter, temporary chaotic encounter, hardened chaotic encounter. Top row: the separation between each pair of BHs $\Delta r_{ij}$ as a function of time. Top row: the trajectories in x and y of each BH during the encounter in the three-body COM frame. In the hierarchical panel, we turn off gas effects at $t=16\Omega ^{-1}$ (vertical dashed line), observing the stagnation of the single BHs semimajor axis.$

Figure 5.

The trajectories and object separations for the four characteristic types of triple encounters in gas. Left to right: Glancing encounter, hierarchical encounter, temporary chaotic encounter, hardened chaotic encounter. Top row: the separation between each pair of BHs |$\Delta r_{ij}$| as a function of time. Top row: the trajectories in x and y of each BH during the encounter in the three-body COM frame. In the hierarchical panel, we turn off gas effects at |$t=16\Omega ^{-1}$| (vertical dashed line), observing the stagnation of the single BHs semimajor axis.

Open in new tab Download slide

If the periapsis of the single is approximately equal to the initial BBH semimajor axis, the three-body dynamics become chaotic. If the gas dissipation is less efficient, then the typical apoapses of the objects (with respect to the three-body COM) will be larger in accordance with the energy of the system (equation (14)). This gives a larger probability that one of the BHs will be launched with an apoapsis large enough to escape |$r_\mathrm{H,t}$| and become ionized by the SMBH. Upon ionization, the remaining binary is hardened relative to the penultimate temporary BBH in all our simulations. This is expected as it follows that the energy gained to have |$r_\mathrm{p,sin}>r_\mathrm{H,t}$| must have come from an energy exchange with the BBH if the previous orbit of the single had |$r_\mathrm{p,sin}< r_{\rm H,t}$| and |$a_\mathrm{sin}< r_{\rm H,t}/2$|⁠.

The final scenario, hardened chaotic encounters, are akin to the temporary chaotic encounters but with more efficient gas drag during the initial interaction and in the following evolution. In these simulations, energy is removed rapidly from the three-body system, quickly shrinking the mean separation of each of the three BHs, creating a far more compact three-body system that remains chaotic. As the three-body system loses energy, the required energy gain of one of the BHs to remove it from the system becomes increasingly large, requiring extremely close encounters (and correspondingly massive binary hardening) to eject a BH. Since this mechanism requires efficient drag to shrink the three-body system quickly (so that an object is not ejected in the early stages of the three-body interaction where |$a_\mathrm{sin}$| is typically larger) this scenario becomes more common for higher ambient gas densities, as we will see in Section 3.6.

As a simple test, we switch off gas in our hierarchical example at |$t=16\Omega _0^{-1}$| (Fig. 5, top row, second column). When gas is switched off, we observe the evolution of |$a_\mathrm{sin}$| and |$a_\mathrm{bin}$| immediately stagnates. Modulations still exist in the separation of the binary + single separation due to accelerations from the SMBH.

3.4 Hierarchical triples, in more detail

For encounters that have a typically larger periapsis, the binary may still be treated as a single object during the encounter. In such cases, the BH and gas dynamics are akin to the single-single case. After the initial encounter, the binary and single orbit each other, whilst shedding gas through spiral outflows. The gas dissipation then allows the entire three-body system to remain bound as a hierarchical triple. We show an example of this scenario in Fig. 6. In the example, we have the first approach (top panel) of the single where it is captured into a counter clockwise orbit (prograde with the initial circum-binary disc) around the binary via gaseous dissipation. The resulting orbit of the single is relatively circular (⁠|$e_\mathrm{sin}\le 0.1$|⁠). The circular orbit and the prograde rotation of the single lead to the generation of large spiral wakes (middle panel). At this point, there is a steady and efficient removal of energy from the system, hardening both |$a_\mathrm{sin}$| and |$a_\mathrm{bin}$|⁠, with |$e_\mathrm{sin}$| remaining small. In all hierarchical encounters identified, |$a_\mathrm{sin}$| shrinks faster than |$a_\mathrm{bin}$|⁠. If we consider the BBH and the BBH-BH systems in isolation, this is reminiscent of the |$\mathrm{d}a/\mathrm{d}t\propto a$| behaviour expected in binary systems (e.g. Ishibashi & Gröbner 2020), since |$a_\mathrm{sin}>a_\mathrm{bin}$|⁠. This behaviour leads to all hierarchical systems gradually moving back towards chaotic type encounters in their late evolution.

$The gas morphology during the formation of and evolution of an embedded hierarchical triple system ($\Sigma =\Sigma _0$, $p=1.8r_{\rm H,t}$). Top: the intersection of the circum-binary and circum-single discs just prior to the first close encounter. Middle: the production of strong spiral outflows that persist for the first $\sim 15-20$ orbits of the single around the binary. Bottom: Late-stage evolution of the system, where the gas has refiled the three-body Hill sphere to form a circum-triple disc.$

Figure 6.

The gas morphology during the formation of and evolution of an embedded hierarchical triple system (⁠|$\Sigma =\Sigma _0$|⁠, |$p=1.8r_{\rm H,t}$|⁠). Top: the intersection of the circum-binary and circum-single discs just prior to the first close encounter. Middle: the production of strong spiral outflows that persist for the first |$\sim 15-20$| orbits of the single around the binary. Bottom: Late-stage evolution of the system, where the gas has refiled the three-body Hill sphere to form a circum-triple disc.

Open in new tab Download slide

We note that the rotation of the gas around the binary component remains prograde with the binary. The maintenance of a low eccentricity in the effective binary system of the hierarchical binary-single system is consistent with the eccentricity damping in prograde embedded binary simulations (e.g. Dittmann & Ryan 2022; Li et al. 2022b; Rowan et al. 2023). After |$\sim 15-20$| orbits of the single, the gas relaxes and replenishes the three-body Hill sphere, forming a circum-triple disc. We still identify the presence of a circum-single disc around the single and a circum-binary disc around the binary, with each component of the binary now also having their own miniscule circum-single discs!

We find hierarchical systems are only formed at the beginning of binary-single encounters and no such systems form following an initially chaotic scattering. We suspect this is likely due to a lower local gas mass following the first close encounter and the messy and chaotic gas flows not having time to form the morphology of the hierarchical system shown in Fig. 6 (which follows directly from the interaction of the BBH and BH accretion discs). Therefore the gas is unable to sufficiently circularize the orbit of the single and it undergoes another close three-body scattering.

3.5 Hardened chaotic triples

Here, we examine the aforementioned hardened chaotic encounters. When the initial periapsis of the single is of the same order as the binary semimajor axis, gas dissipation is highly efficient, forming a tightly bound three-body system that evolves in a chaotic manner. These systems are arguably the closest hydrodynamical analogue to the gasless three-body simulations of e.g. Samsing et al. (2022), Fabj & Samsing (2024), and Trani et al. (2024), where the separations between each BH are constantly changing. We show how the three-body energy and BH separations evolve as an example in Fig. 7. In this example, there is a small initial amount of dissipation following a more shallow first encounter. This softens the initial binary such that the system is no longer a hierarchical system. This is followed by a strong three-body encounter that efficiently removes energy from the system, reducing the mean separation between all BHs. Energy is then efficiently and continuously removed from the system. We observe fluctuations where energy is injected rather than removed from the triple or temporary softenings of the most bound BH pair. However, in all hardened encounters there is a net removal of energy. As another test, we switch off the gas just prior to the first strong encounter (bottom panel of Fig. 7). When gas is switched off, the system retains a large mean separation between the three objects, i.e. |$a_\mathrm{bin}$| and |$a_\mathrm{sin}$| remain large during each binary-single state. The absence of gas also leads little change in |$a_\mathrm{bin}$| as the single executes its larger orbit, where we would before expect gradual inspiral from gas drag.

$Example of a hardened chaotic encounter ($p=1.7r_{\rm H,t}$ and $\Sigma =\Sigma _0$). Top: the gas morphology after the second triple encounter, where the triple becomes strongly bound. Black markers have been added to better indicate the positions of the BHs. Middle: the energy of the triple system $E_\mathrm{trip}$ (equation (14)) and most bound BH pair $E_\mathrm{bin}$ (equation (15)) and the BH separations as a function of time. The vertical dash indicates the time the hydrodynamical snapshot is taken. Bottom: the BH separations where gas is turned is off at $t=15.1\Omega _0^{-1}$ (dashed vertical line). The faded lines indicate the triple evolution when gas gravity is switched off. h and $\delta _\text{min}$ are represented by the horizontal lines.$

Figure 7.

Example of a hardened chaotic encounter (⁠|$p=1.7r_{\rm H,t}$| and |$\Sigma =\Sigma _0$|⁠). Top: the gas morphology after the second triple encounter, where the triple becomes strongly bound. Black markers have been added to better indicate the positions of the BHs. Middle: the energy of the triple system |$E_\mathrm{trip}$| (equation (14)) and most bound BH pair |$E_\mathrm{bin}$| (equation (15)) and the BH separations as a function of time. The vertical dash indicates the time the hydrodynamical snapshot is taken. Bottom: the BH separations where gas is turned is off at |$t=15.1\Omega _0^{-1}$| (dashed vertical line). The faded lines indicate the triple evolution when gas gravity is switched off. h and |$\delta _\text{min}$| are represented by the horizontal lines.

Open in new tab Download slide

What separates the temporary from hardened encounters is how quickly the system can contract relative to the number orbits of the single. At each close three-body encounter, the re-shuffling of |$a_\mathrm{sin}$| presents another opportunity for the single BH to escape. If the system contracts quickly, then the region of the parameter space for the single (⁠|$a_\mathrm{sin}$| and |$e_\mathrm{sin}$|⁠) that will facilitate its escape quickly becomes smaller. Since the dissipation scales with the density |$\Sigma$|⁠, we find hardened chaotic encounters are more common in our higher density simulation suites.

3.6 The parameter space of binary-single encounters

To better understand the dependence of the encounter type on the initial encounter parameters, we show the minimum separation of the single and the binary COM as a function of the single’s initial impact parameter p and colour-code the data points according to the encounter type in Fig. 8. The ‘w’ shape of the curve is reminiscent of the periapsis-impact parameter distribution of single-single BH scatterings (e.g. Boekholt, Rowan & Kocsis 2023; Rowan et al. 2024b; Whitehead et al. 2024a). The periapsis of the first encounter, |$r_\mathrm{p,sin,1}$|⁠, in the centre of our impact parameter p domain is more stochastic in the binary single case. This is due to the BBH not being another point mass, so the separation between the BBH COM and the single becomes hard to predict when the single reaches a separation of |$a_\mathrm{bin}$| from the BBH. As the original binary hardens slightly less during the approach of the single for lower |$\Sigma$|⁠, this stochasticity is larger.

$The closest approach of the single $r_\mathrm{p,sin,1}$ as measured from the binary COM as a function of impact parameter p for each simulation suite in density $\Sigma$. The data points are colour coded by the type of encounter: No close interaction (F), glancing encounter (G), hierarchical encounter (H), temporary chaotic encounter (TC), and hardened chaotic encounter (HC). The results indicate increased AGN disc densities produce more HC encounters. The crossed markers represent triple interactions where there is no second encounter. The triangular markers indicate runs were the separation of two BHs is small enough to reasonably expect a merger while the system was resolved (see Section 3.7.1).$

Figure 8.

The closest approach of the single |$r_\mathrm{p,sin,1}$| as measured from the binary COM as a function of impact parameter p for each simulation suite in density |$\Sigma$|⁠. The data points are colour coded by the type of encounter: No close interaction (F), glancing encounter (G), hierarchical encounter (H), temporary chaotic encounter (TC), and hardened chaotic encounter (HC). The results indicate increased AGN disc densities produce more HC encounters. The crossed markers represent triple interactions where there is no second encounter. The triangular markers indicate runs were the separation of two BHs is small enough to reasonably expect a merger while the system was resolved (see Section 3.7.1).

Open in new tab Download slide

Similar to the single-single case, the rotation of the single about the three-body COM is prograde (counter clockwise) for low impact parameters left of the first trough and higher parameters beyond the second trough. For mergers in-between, the single executes a retrograde trajectory, where the transition between prograde and retrograde encounters occurs on the troughs, representing the cases where the angular momentum of the single as it enters |$r_{\rm H,t}$| tends towards zero at periapsis.

As mentioned, hierarchical and glancing encounters occur in prograde encounters with larger initial periapses. For systems in the central retrograde portion of the |$p-r_\mathrm{p,sin,1}$| domain, we observe either temporary or hardened chaotic encounters. Comparing the |$\Sigma =\Sigma _0$| to the |$\Sigma =0.1\Sigma _0$| runs, we find an increase larger than 200 per cent in the number of hardened chaotic encounters (13 versus 5) despite a comparable number of at least temporary formed triples. (20 versus 17), owed to typically more efficient energy dissipation. We summarize the scattering outcomes for each value of the density |$\Sigma$| suite in Table 1.

Table 1.

Open in new tab

The number of encounter types per simulation suite with density |$\Sigma$|⁠. The encounter types are labelled as: failed encounters (F), i.e no close encounter, glancing encounters (G), hierarchical (H), temporary chaotic encounters (TC), and hardened chaotic encounters (HC).

\|$\Sigma /\Sigma _0$\|	F	G	H	TC	HC
1	4	0	4	3	13
0.25	5	3	0	9	6
0.1	5	2	1	11	5

Table 1.

Open in new tab

\|$\Sigma /\Sigma _0$\|	F	G	H	TC	HC
1	4	0	4	3	13
0.25	5	3	0	9	6
0.1	5	2	1	11	5

3.7 Gravitational wave prospects

3.7.1 Constraining the merging time-scale at the resolution limit

Following our findings that gas reliably hardens the three-body system, we consider the prospects for gravitational wave sources originating from embedded binary-single scatterings. Given the hardened chaotic encounters manage to reach the resolution limit of our simulations, we can infer the fate of the system by anticipating the least favourable outcome for merger. Here, we make pessimistic assumptions for the evolution of the system to give a lower bound on the merger probability.

If we assume that the three-body system will continue to efficiently harden at scales smaller than the resolution |$\delta _\mathrm{min}$|⁠, then the ejection of the third object is the only means to increase the time-scale for a BH merger. This follows as the evolution of embedded binaries with separations on the scale of |$\delta _\mathrm{min}$| has not yet been simulated self-consistently from the larger length-scales shown here, therefore we remain agnostic to their hydrodynamical evolution. The ejection of an object also removes the chance of the chaotic three-body system to result in a very close flyby of two BHs that may induce a merger (see Section 3.7.2). We estimate the hardening of the remaining binary upon the ejection of a BH by considering the energy exchange required to soften |$a_\mathrm{sin}=\delta _\mathrm{min}$| to |$a_\mathrm{sin}=r_{\rm H,t}$|⁠. The required gain in energy to escape to a distance |$r_2$| from some initial value |$r_1$| is given by

$$\begin{eqnarray} E_\mathrm{esc}=\frac{G M_3 (M_1+M_2)}{2}\bigg (\frac{1}{r_1}-\frac{1}{r_2}\bigg). \end{eqnarray}$$

(17)

Taking the energy of the binary as |$E_\mathrm{bin}=-G M_1 M_2/2a_\mathrm{bin}$|⁠, the hardening of the binary is then related to the energy exchange by

$$\begin{eqnarray} \frac{\Delta a_\mathrm{bin}}{a_\mathrm{bin}}=\frac{-E_\mathrm{esc}}{E_\mathrm{bin}}=\frac{M_3(M_1+M_2)}{M_1M_2}\bigg (\frac{a_\mathrm{bin}}{r_2}-\frac{a_\mathrm{bin}}{r_1}\bigg)\, . \end{eqnarray}$$

Setting |$r_1=\delta _\mathrm{min}$|⁠, |$r_2=r_{\rm H,t}$| and noting here |$M_1=M_2=M_3$|⁠, this simplifies to

$$\begin{eqnarray} \frac{\Delta a_\mathrm{sin}}{a_\mathrm{sin}}=2\bigg (\frac{a_\mathrm{bin}}{r_{\rm H,t}}-\frac{a_\mathrm{bin}}{\delta _\mathrm{min}}\bigg)\, . \end{eqnarray}$$

(18)

Taking |$a_\mathrm{bin}=a_\mathrm{sin}/2=\delta _\mathrm{min}/2$|⁠, (i.e. the binary is still very soft with respect to the triple system) this gives

$$\begin{eqnarray} \frac{\Delta a_\mathrm{sin}}{a_\mathrm{sin}}=2\bigg (\frac{\delta _\mathrm{min}}{2r_{\rm H,t}}-\frac{1}{2}\bigg)\approx -1\, \, \, \, (\delta _\mathrm{min}<< r_{\rm H,t})\, . \end{eqnarray}$$

(19)

This means the resulting binary will have a final semimajor axis of |$\delta _\mathrm{min}/4\approx 0.0019$| au (0.4|$R_\odot$|⁠) under these pessimistic assumptions. We refer to Fabj & Samsing (2024) for a more detailed explanation of this analytical approximation. Assuming all binaries have this maximal semimajor axis, we compute the GW inspiral time as a function of the eccentricity |$e_\mathrm{bin}$| by numerically integrating the orbital evolution equations of Peters (1964)

$$\begin{eqnarray} \bigg \langle \frac{\mathrm{d}a}{\mathrm{d}t}\bigg \rangle =-\frac{64}{5}\frac{G^3 M_1 M_2 (M_1 + M_2)}{c^5 a^3(1-e^2)^{7/2}}\bigg (1+\frac{73}{24}e^2+\frac{37}{96}e^4\bigg), \end{eqnarray}$$

(20)

$$\begin{eqnarray} \bigg \langle \frac{\mathrm{d}e}{\mathrm{d}t}\bigg \rangle =-\frac{304}{15}\frac{G^3 M_1 M_2 (M_1 + M_2)}{c^5 a^4(1-e^2)^{5/2}}\bigg (e+\frac{121}{304}e^3\bigg)\, . \end{eqnarray}$$

(21)

We compare the inspiral time-scale |$t_\mathrm{GW}$| with the viscous time-scale for gas at the outer edge of a circum-binary disc of size.² |$R_\mathrm{disc} = r_{\rm H,12}$|/2, |$\tau _\mathrm{visc}=R_\mathrm{disc}^2/\nu$| as a measure of time to re-establish a typical circum-binary gas morphology. The relative time-scales are shown in Fig. 9 for final binaries with |$a_\mathrm{bin}=\lbrace 1,1/2,1/4,1/8,1/16\rbrace \delta _\mathrm{min}$|⁠. The results show that the length-scales of our simulations come close to those where GWs could begin to dominate the evolution of the binary. We can be more confident in the fate of more eccentric binaries, owed to the stronger GW emission at periapsis. From the |$1/a^3$| dependence of equation (20), just a factor 4 decrease in separation of |$a_\mathrm{bin}=\delta _\mathrm{min}/16$| from the fiducial |$a_\mathrm{bin}=\delta _\mathrm{min}/4$| value leads to reliable mergers independently of e. Assuming the eccentricity of the final binary follows the theoretical co-planar derived distribution of |$P(e_\mathrm{bin})=e_\mathrm{bin}/\sqrt{1-e_\mathrm{bin}^2}$| (e.g. Monaghan 1976), the percentage of binaries that meet the criterion |$\tau _\mathrm{GW}<\tau _\mathrm{visc}$| is {28 per cent, 38 per cent, 61 per cent, 92 per cent, 100 per cent} for |$a_\mathrm{bin}=\lbrace 1, 1/2,1/4,1/8,1/16\rbrace \delta _\mathrm{min}$|⁠, implying a non-negligible amount of systems may merge. We stress that this conclusion contains several pessimistic assumptions:

Most notably, for this calculation we assume gas does not further harden the triple system before a BH is ejected, which may still be significant below the resolvable scales in our simulations.
We assume gas does not further harden the resulting binary upon the ejection of the single. If the binary experiences net energy removal at these separations, this would imply all such systems would merge on a time-scale shorter than predicted on the right-hand axis of Fig. 9 regardless of whether |$\tau _\mathrm{GW}<\tau _\mathrm{visc}$| or not.
We assume a BH is ejected with the minimum amount of energy to reach |$r_{\rm H,t}$|⁠. In reality the ejected BH can carry more energy away from the system, leaving behind a tighter binary. Where the remaining energy of the binary approximately follows the distribution P(⁠|$|E_\mathrm{bin}|)\propto |E_\mathrm{bin}|^{-3}$|⁠, see Valtonen & Karttunen (2006) and Stone & Leigh (2019).
As we do not modify the viscosity |$\nu$| beyond the assumptions of a steady AGN disc, the viscous time will be underestimated owed to the scale height |$H=c_\text{s,0}/\Omega =c_\text{s,0}\sqrt{\frac{R^{3}}{GM_\text{bin}}}$| in the circum-binary disc being smaller than that of the AGN disc (maintaining the isothermal assumption).
We ignore the prospect for GW dissipation to halt the removal of the single during its final periapsis, which would similarly harden the triple system and potentially move the system into the regime where any resulting binary falls within the regime required for reliable merger.

Note that the hydrodynamical evolution of the triple will become inaccurate before reaching the resolution |$\delta _\mathrm{min}$|⁠. If we take the softening length of our sinks to be the limiting scale of the hydrodynamics³ for a three-body interaction, then this corresponds to |$\sim 10\delta _\mathrm{min}$|⁠. Therefore we must also trust that the system will harden beyond this point. We note that we find no triple systems are softened by the gas over the range |$a=a_0$| to |$a=10\delta _\mathrm{min}$|⁠, therefore we expect this to be a fair assumption.

$The maximal time-scale until a BBH merger, $t_\mathrm{GW}$ for varying initial binary separations around the resolution limit of our simulations $a_\mathrm{bin}=\lbrace 1,1/2,1/4,1/8,1/16\rbrace \delta _\mathrm{min}$. These values assume a three-body (or binary) system is not further hardened by gas, which may significantly reduce these values. The time-scale is shown in the units of the viscous time-scale for the circum-binary disc $\tau _\mathrm{visc}=R_\mathrm{disc}^{2}/\nu$. Whether the binary will merge within the viscous time of the disc depends sensitively on $e_\mathrm{bin}$ and $a_\mathrm{bin}$.$

Figure 9.

The maximal time-scale until a BBH merger, |$t_\mathrm{GW}$| for varying initial binary separations around the resolution limit of our simulations |$a_\mathrm{bin}=\lbrace 1,1/2,1/4,1/8,1/16\rbrace \delta _\mathrm{min}$|⁠. These values assume a three-body (or binary) system is not further hardened by gas, which may significantly reduce these values. The time-scale is shown in the units of the viscous time-scale for the circum-binary disc |$\tau _\mathrm{visc}=R_\mathrm{disc}^{2}/\nu$|⁠. Whether the binary will merge within the viscous time of the disc depends sensitively on |$e_\mathrm{bin}$| and |$a_\mathrm{bin}$|⁠.

Open in new tab Download slide

3.7.2 Highly eccentric dynamical mergers

In several simulations (9 of 65 initially bound triple systems formed), we find BH separations of |$\Delta r_{ij}< 10r_\bullet$| where |$r_\bullet =GM_\mathrm{BH}/c^{2}$| is the Schwarzschild radius of one of the BHs. These cases are the result of the inherently unstable three-body configuration, where each encounter randomizes the energy and angular momentum of the objects. While we are limited in our simulations to a resolution limit of |$\delta _\mathrm{min}\simeq 3\times 10^{4}r_\bullet$|⁠, these cases often occur while the separations of the BHs were thus far fully resolved. As we don’t expect the trajectories of the BHs to significantly change in the very small time-frame of the approach from |$\Delta r=\delta _\mathrm{min}$| to |$\Delta r< 10r_\bullet$|⁠, we expect the close periapsis to result in strong GW emission and possible prompt merger of the BHs involved. Following a similar methodology to Fabj & Samsing (2024), we estimate the probability enhancement for a binary merger with observable residual eccentricity due to random ultra-close encounters during hardened chaotic interactions.

The inspiral time of a BBH in the high eccentricity limit (assuming the presence of the more distant third body can be ignored) is given by

$$\begin{eqnarray} t_\mathrm{GW,e}=\frac{5\sqrt{2}c^{5}}{64G^{3}}\frac{a_\mathrm{bin}^{1/2}r_\mathrm{p,bin}^{7/2}}{M_1 M_2 (M_1+M_2)}\, . \end{eqnarray}$$

(22)

Assuming the current three-body dynamics are un-hierarchical and the triple system has a characteristic scale |$r_\mathrm{t}$| (i.e. |$a_\mathrm{bin}\sim a_\mathrm{sin}\sim r_\mathrm{t}$|⁠), then the orbital period of the single around the binary is of the order of

$$\begin{eqnarray} T_\mathrm{sin}\approx 2\pi \sqrt{\frac{r_\mathrm{t}^{3}}{G(M_1+M_2+M_3)}}\, . \end{eqnarray}$$

(23)

For the binary to merge before the single as the chance to change the binary dynamics, we must have |$t_\mathrm{GW,e}< T_\mathrm{sin}$|⁠. Again assuming only the equal mass case |$M=M_1=M_2=M_3$| this gives a maximal periapsis for the binary to reliably merge of

$$\begin{eqnarray} r_\mathrm{mrg}(r_\mathrm{t}) \lesssim \bigg (\frac{256\pi G^{5/2}}{5\sqrt{6}c^{5}}\bigg)^{2/7}r_\mathrm{t}^{2/7}M^{5/7}, \end{eqnarray}$$

(24)

where we note a small factor |$3^{-2/7}$| correction to equation (3) in Fabj & Samsing (2024). Expressing this as an eccentricity |$e_\mathrm{mrg} = 1-r_\mathrm{mrg}/a_\mathrm{bin}$|⁠, we have

$$\begin{eqnarray} e_\mathrm{mrg}(r_\mathrm{t}) \gtrsim 1-\bigg (\frac{256\pi G^{5/2}}{5\sqrt{6}c^{5}}\bigg)^{2/7}r_\mathrm{t}^{-5/7}M^{5/7}\, . \end{eqnarray}$$

(25)

In the absence of gas, the scale of the triple system is set by the initial binary separation |$a_0$|⁠, as the system cannot contract since there is no dissipation.⁴ Therefore the minimal⁵ enhancement |$\eta$| to the probability of a merger due to gas is the ratio of the probability that |$e>e_\mathrm{mrg}(r_\mathrm{t})$|⁠, i.e.

$$\begin{eqnarray} \eta =\frac{P(e>e_\mathrm{mrg}(r_\mathrm{t}))}{P(e>e_\mathrm{mrg}(a_0))}\, . \end{eqnarray}$$

(26)

Recalling |$P(e)\approx e/\sqrt{1-e^2}$| in the co-planar limit, we compute these probabilities as a function of |$r_\mathrm{t}$| in Fig. 10. We also show the probability that a three-body system will undergo a merger as a function of the number of encounters and |$r_\mathrm{t}$|⁠, assuming the probability scales as |$P_\mathrm{mrg}(r_\mathrm{t})=1-[1-P(e>e_\mathrm{mrg}(r_\mathrm{t}))]^{N_\mathrm{enc}}$|⁠, where |$N_\mathrm{enc}$| is the number of encounters. For triples hardened to |$r_\mathrm{t}=h$| and |$r_\mathrm{t}=\delta _\mathrm{min}$| through gaseous drag, we find the probability per encounter to have |$e>e_\mathrm{crit}$| is enhanced by a factor |$\eta \approx 3.5$| and 8, respectively, compared to |$r_\mathrm{t}=r_{\rm H,t}$|⁠. As we have remained pessimistic on potential triple hardening beyond our resolution limits, these are lower bounds on the value of |$\eta$|⁠. However, we are limited in its convergence due to the immense length/resolution requirements of these simulations. Adopting the mean number of encounters of |$N_\mathrm{enc}\approx 20$| as suggested by Fabj & Samsing (2024), we find the probability for systems that harden to |$r_\mathrm{t}=h$| and |$\delta _\mathrm{min}$| to merge is |$\sim 0.4$| and |$\sim 0.65$|⁠. We note that the addition of gas often leads to longer three-body exchanges in addition to prolonged periods where the identification of a binary-single like state is not possible, although whether the latter result enhances or reduces the likelihood of a highly eccentric encounter between two BHs is uncertain.

$The enhancement of the merger probability due to gas hardening. Top: the minimum critical eccentricity $e_\mathrm{mrg}$ (equation (25)) required for two BHs to merge during a triple exchange as a function of the compactness of the triple $r_\mathrm{t}$ in Schwarzschild radii. Also shown is the probability enhancement per encounter $\eta$ (equation (26)) that this eccentricity is achieved relative to the initial binary compactness. Bottom: the probability of a merger $P_\mathrm{mrg}$ as a function of the number of binary-single like states created during a three-body exchange $N_\mathrm{enc}$ as a function of the compactness of the triple. We mark the contours of $P_\mathrm{mrg}=\lbrace 0.5,0.9\rbrace$. In both plots, we indicate the length-scales of the resolution $\delta _\mathrm{min}$, the softening h and initial binary semimajor axis $a_0$.$

Figure 10.

The enhancement of the merger probability due to gas hardening. Top: the minimum critical eccentricity |$e_\mathrm{mrg}$| (equation (25)) required for two BHs to merge during a triple exchange as a function of the compactness of the triple |$r_\mathrm{t}$| in Schwarzschild radii. Also shown is the probability enhancement per encounter |$\eta$| (equation (26)) that this eccentricity is achieved relative to the initial binary compactness. Bottom: the probability of a merger |$P_\mathrm{mrg}$| as a function of the number of binary-single like states created during a three-body exchange |$N_\mathrm{enc}$| as a function of the compactness of the triple. We mark the contours of |$P_\mathrm{mrg}=\lbrace 0.5,0.9\rbrace$|⁠. In both plots, we indicate the length-scales of the resolution |$\delta _\mathrm{min}$|⁠, the softening h and initial binary semimajor axis |$a_0$|⁠.

Open in new tab Download slide

4 DISCUSSION AND IMPLICATIONS

Our results suggest that the combination of gas drag and the complex three-body dynamics of binary-single scatterings in AGN discs could lead to the rapid merger of previously secularly evolving embedded BBHs. This would indicate that predictions for the number of required binary-single encounters to harden a system to merger from N-body or analytical studies that do not include gas are an upper bound. For example, Leigh et al. (2018) predict a mean of |$\sim$|10 such encounters to dynamically harden a binary enough to inspiral via GWs. The degree of this overestimation is sensitive to the density of the AGN disc, therefore we stress the need for constraints on this quantity, which can vary wildly depending on which AGN disc model and/or viscosity is assumed e.g. Thompson (Thompson, Quataert & Murray 2005) or Sirko-Goodman discs (Sirko & Goodman 2003; Goodman & Tan 2004) and their physical parameters (e.g. SMBH accretion rate).

Though we only consider co-planar binary-single encounters with no velocity dispersions, the three-body simulations of Trani et al. (2024) around an SMBH find that the distribution of final eccentricities and semimajor axes are largely unaffected provided the dispersion |$\sigma$| is |$\sigma < 10^{-4}v_\mathrm{K}$|⁠, where |$v_\mathrm{K}$| is the Keplerian orbital velocity about the SMBH. Therefore we do not expect the three-body dynamics to differ significantly within this threshold (for our fixed parameters, e.g. object mass ratios and equation of state). At higher dispersions, they predict the remaining binary eccentricities will follow a more thermal rather than super-thermal distribution, reducing the merger likelihood. The scale height |$H/R \approx 0.0056\sim 10^{-2}$| provides the length-scale where the gas density is approximately unchanged and where we expect only small changes in the magnitude of gaseous triple hardening. However, we suspect there will be a stronger inclination dependence on the level of dissipation for non-negligible inclinations (⁠|$i\gtrsim H/R$|⁠) during the first encounter (where the mini-discs are still preserved) depending on the relative inclination of the encounter to the circum-single/binary discs, which warrants further study.

The reliable shrinkage of the triple system that we report has several implications not considered in the main body of this work. The reduction in the mean value of |$a_\mathrm{bin}$| and more importantly |$a_\mathrm{sin}$| increases the likelihood that a third object may be in close proximity to a merging binary during the exchange. This implies that triples in AGN could produce unique dephasing features in the GW waveform of the inspiralling binary (e.g. Samsing et al. 2024; Hendriks, Zwick & Samsing 2024a; Hendriks et al. 2024b). By a similar argument, a tighter three-body system increases the chance that a tight binary with a small enough separation and large enough eccentricity to have a measurable residual eccentricity. Such cases could potentially be detected as more sensitive third-generation detectors come online, with Cosmic Explorer (e.g. Reitze et al. 2019; Evans et al. 2023; Borhanian & Sathyaprakash 2024) and the Einstein Telescope (e.g. Hild, Chelkowski & Freise 2008; Franciolini et al. 2023) design sensitivities able to detect |$e\gtrsim 5\times 10^{-3}$| and |$e\gtrsim 10^{-3}$|⁠, respectively, at a GW frequency of 10 Hz (e.g. Saini 2024). The hardening of embedded triples also increases the binding energy between a merging binary and the tertiary. This increases the prospect for the now 2-body system to be retained after the remnant receives a post-merger kick (e.g. Campanelli et al. 2007; González et al. 2007; Varma et al. 2022), allowing for possible back to back or ‘double’ mergers (e.g. Samsing & Ilan 2019).

In principle, the gas drag during a chaotic encounter of a binary and a single BH directly relates to the scenario of binary-binary scatterings (e.g. Mikkola 1984; Bacon, Sigurdsson & Davies 1996; Heggie 2000; Arca Sedda et al. 2021), which serve to reduce the binary fraction on average (e.g. Sweatman 2007; Sollima 2008; Ryu, Leigh & Perna 2017). There is no reason to suggest that gas drag will not dissipate energy during the initial encounter of two binaries and during the chaotic exchange that follows, hardening the system. We therefore posit that possible reductions in the merger rate due to binary-binary (and binary-single) encounters in an AGN disc will be mitigated to some degree when gas is included, since the merger probability will be enhanced. Similarly to the hierarchical triple scenario, it is plausible that gas could aid in the formation of highly exotic hierarchical (2 + 2) 4-body systems, provided the impact parameter is relatively large compared to the binaries’ semimajor axes.

5 SUMMARY AND CONCLUSIONS

We performed 72 hydrodynamical simulations of binary-single black hole interactions in an AGN disc covering 3 disc densities and 24 impact parameters, examining the role of gas on the fate of the triple system. We summarize our key findings below:

Similar to single-single BH encounters, binary-single interactions feature significant energy dissipation during the first encounter, allowing for the formation of an energetically bound triple system.
We identify four characteristic interactions for close binary-single encounters in gas (Section 3.3): Glancing encounters, where the single does not become bound but modifies the eccentricity and/or semimajor axis of the original binary. Hierarchical encounters, where the single is bound through gas dissipation into a low eccentricity orbit around the original binary, forming a quasi-stable hierarchical triple system. Temporary chaotic encounters, where the periapsis of the single is on par with the initial binary separation, leading to a chaotic encounter that ends with one object being thrown out. Hardened chaotic encounters, where the dissipation during a chaotic encounter is efficient enough to quickly reduce the mean separation between the three BHs, leading to prolonged three-body encounters.
For the chaotic encounters (Section 3.5), the gas morphology is highly dynamic and volatile, leading to fluctuations in the energy of the three-body system. Although, in all cases, net energy is removed from the triple.
Assuming the system undergoes a close encounter with separations less than the three-body Hill sphere, chaotic encounters are the most common encounter type for our range of AGN disc densities |$\Sigma$| (Section 3.6).
The type of encounter is primarily set by first periapsis distance of the single object |$r_\mathrm{p,sin,1}$| and the gas density |$\Sigma$|⁠. Glancing and hierarchical encounters occur at larger |$r_\mathrm{p,sin,1}$| and |$\Sigma$|⁠, temporary chaotic encounters are more likely at small |$r_\mathrm{p,sin,1}$| and |$\Sigma$|⁠. Chaotic hardened encounters occur at small |$r_\mathrm{p,sin,1}$| and large |$\Sigma$| (Fig. 8).
Provided the triple remains bound, gas continually removes energy from the three-body system, causing it to contract. For hardened chaotic encounters, this shrinks the typical semimajor axes of the BBH and BBH–BH components during the establishment of binary-single like states during the three-body interaction. For hierarchical encounters, we find the semimajor axis of the BBH–BH system shrinks faster than the BBH, gradually bringing the system to a more chaotic state. The efficiency of the dissipation in all cases scales with the density of the AGN disc.
Dissipation is most efficient in hardened chaotic encounters, where the three-body system shrinks to our resolution limit |$\delta _\mathrm{min}$| in a few AGN orbits, a length-scale two orders of magnitude smaller than the initial binary separation. We evolved the orbital elements of binaries at the resolution limit, assuming a BH is immediately ejected (Section 3.7.2). Sampling over the eccentricity distribution, we find 28 per cent of our chaotic systems are expected to merge through secular GW emission.
Depending on the compactness |$r_\mathrm{t}$| of the triple, the critical minimum eccentricity that two BHs must have in order to merge during a binary + single state is given in equation (25). The relation predicts an |$\sim$|3.5 and 8 times increase in the probability (during a single encounter) for the three body system to merge after being hardened to the softening h and resolution |$\delta _\mathrm{min}$| lengths, respectively (Section 3.7.1). For a fiducial value of 20 binary + single states, the probability of merger is |$\sim 0.4$| and |$\sim 0.65$| for |$r_\mathrm{t}=h$| and |$r_\mathrm{t}=\delta$|⁠, respectively (Section 3.7.2). This suggests that while they are potentially rare, binary-single encounters in AGN discs could reliably induce BH mergers.
The ability for gas to shrink the size of the three-body system will increase the likelihood of unique BH–BH merging scenarios. This includes potential observable eccentricity in GWs, dephasing from the tertiary object and repeated mergers. As we are limited by our resolution limit we leave the quantitative analysis of these sub-grid phenomena to future research.

We have made several assumptions and simplifications in this work. First, we have fixed the semimajor axis of the binary until the encounter in order to better understand the effect of the ambient disc density on the evolution of the triple. If the embedded binary inspiral time is non-negligible compared to the triple encounter time-scale then there will be a bias towards smaller initial semimajor axes in the binary during the triple interaction. We posit that this will lead to more hierarchical encounters in denser mediums (see Section 3.4).

We have restricted ourselves to co-planar binary-single encounters. In practice there will be some vertical and radial velocity dispersion, which may alter the outcomes of the initial scattering and subsequent evolution (e.g. Trani et al. 2024). We neglect the accretion of mass by the stellar BHs which may be significant (Li & Lai 2023; Rowan et al. 2023), although we previously found that softening lengths of comparable radii to the accretion boundary lead to similar dissipation (e.g. Rowan et al. 2024b). We only consider equal mass BHs here, whereas the true range may lie in the region of |$\sim 5\!-\! \, 40\, \mathrm{M}_{\odot }$| (depending on the assumed metalicity of the nuclear stellar cluster, e.g. Belczynski et al. 2010a). For computational reasons, we assumed an isothermal equation of state. Our previous simulations of single-single scatterings utilizing a non-isothermal equation of state (Whitehead et al. 2024b, see also the work of Li & Lai 2023) reported non-negligible mass reduction in the Hill sphere of each BH before and during the first encounter. Since the dissipation scales with the local density, we expect that evolving the fluid energy equation will produce results more akin to our lower density runs (i.e. more temporary than hardened chaotic encounters compared to the isothermal case).

AGN discs host a large range of dynamical and hydrodynamical problems. As in single-single scatterings, we conclude that gas drag is an important component of binary-single encounters and the evolution of embedded three-body systems, noticeably increasing the likelihood of a merger during their evolution. Current approaches to incorporate the effects of gas dynamical drag in few-body systems are usually based on employing gas drag formulas which are valid only for rectilinear (e.g. Ostriker 1999) or circular (e.g. Kim & Kim 2007; Kim, Kim & Sánchez-Salcedo 2008) trajectories of the massive bodies even when more complex motions are involved. Recent studies hold the promise to consider the effects of gas under more general conditions (e.g. O’Neill et al. 2024; Suzuguchi et al. 2024). Our findings here, combined with our previous work on single-single scatterings (Rowan et al. 2023, 2024b; Whitehead et al. 2024a, b), provide motivation and guidance for the future development of more accurate analytical, or semi-analytical, models to capture gas effects on single-single and binary-single black hole evolution and mergers in AGNs. These effective models could then be used in population studies to produce more accurate estimations for the GW observables arising from black hole mergers in AGN discs accounting explicitly for the effects of gas.

ACKNOWLEDGEMENTS

(i) This work was supported by the Science and Technology Facilities Council Grant Number ST/W000903/1 attributed to Bence Kocsis.

(ii) The research leading to this work was supported by the Independent Research Fund Denmark via grant ID 10.46540/3103-00205B attributed to Martin Pessah.

(iii) This was work was supported by the Villum Fonden grant No. 29466, and by the ERC Starting Grant no. 101043143 – BlackHoleMergs attributed to Johan Samsing.

(iv) This work was performed using resources provided by the Cambridge Service for Data Driven Discovery (CSD3) operated by the University of Cambridge Research Computing Service (www.csd3.cam.ac.uk), provided by Dell EMC and Intel using Tier-2 funding from the Engineering and Physical Sciences Research Council (capital grant EP/T022159/1), and DiRAC funding from the Science and Technology Facilities Council (www.dirac.ac.uk).

DATA AVAILABILITY

The data used to create this work will be shared upon reasonable request.

Footnotes

Note that the global dynamics of the disc may not be as accurately captured here compared to a global disc setup as the azimuthal extent of the disc is large. However, we must maintain this azimuthal extent in order to reach steady state in the minidiscs prior to encounter. We expect that any changes in the dynamics will arise from deviations in the gas mass within |$r_{\rm H,t}$| during the triple encounter. As such, this is unlikely to significantly affect our findings as this will only affect the dissipation time-scale, see later in Section 3.6.

Once the single is ejected, the binary Hill sphere becomes the appropriate length-scale for the formation of a disc.

Given the encounters are typically highly eccentric, the majority of the orbit will be spent at higher separations.

In principle chance encounters that are very close could dissipate energy via GWs but not induce merger, but we assume this to be overall far less efficient than gas on these larger scales.

The value is a minimum as it assumes the number of encounters is unaffected by gas, whereas in reality gas will induce a drag on objects that are on initially escape trajectories, increasing |$N_\mathrm{enc}$|⁠. We do not have enough statistics here to comment on this directly.

REFERENCES

Aarseth

S. J.

Zare

1974

Celest. Mech.

185

Month:	Total Views:
April 2025	41
May 2025	10

Article Contents

Prompt gravitational-wave mergers aided by gas in active galactic nuclei: the hydrodynamics of binary-single black hole scatterings

ABSTRACT

1 INTRODUCTION

2 METHODS

2.1 Nomenclature

2.2 The shearing box

2.3 Gas dynamics

2.4 Mesh refinement

2.5 Initial conditions

2.5.1 The AGN disc

2.5.2 The black holes

2.5.3 The encounter

2.5.4 Shifting the shearing frame

2.5.5 Reaching steady state

3 RESULTS

3.1 Fiducial example

3.2 Dissipation

3.3 Alternate encounter outcomes

3.4 Hierarchical triples, in more detail

3.5 Hardened chaotic triples

3.6 The parameter space of binary-single encounters

3.7 Gravitational wave prospects

3.7.1 Constraining the merging time-scale at the resolution limit

3.7.2 Highly eccentric dynamical mergers

4 DISCUSSION AND IMPLICATIONS

5 SUMMARY AND CONCLUSIONS

ACKNOWLEDGEMENTS

DATA AVAILABILITY

Footnotes

REFERENCES

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only