ABSTRACT

There are indications that stellar-origin black holes (BHs) are efficiently paired up in binary black holes (BBHs) in active galactic nucleus (AGN) disc environments, which can undergo interactions with single BHs in the disc. Such binary–single interactions can potentially lead to an exceptionally high fraction of gravitational-wave mergers with measurable eccentricity in LIGO/Virgo/KAGRA. We here take the next important step in this line of studies by performing post-Newtonian N-body simulations between migrating BBHs and single BHs set in an AGN disc-like configuration, with a consistent inclusion of the central supermassive black hole (SMBH) in the equations of motion. With this set-up, we study how the fraction of eccentric mergers varies in terms of the initial size of the BBH semimajor axis relative to the Hill sphere, as well as how it depends on the angle between the BBH and the incoming single BH. We find that the fraction of eccentric mergers is still relatively large, even when the interactions are notably influenced by the gravitational field of the nearby SMBH. However, the fraction as a function of the BBH semimajor axis does not follow a smooth functional shape, but instead shows strongly varying features that originate from the underlying phase-space structure. The phase-space further reveals that many of the eccentric mergers are formed through prompt scatterings. Finally, we present the first analytical solution to how the presence of an SMBH in terms of its Hill sphere affects the probability for forming eccentric BBH mergers through chaotic three-body interactions.

1 INTRODUCTION

The observations of gravitational-wave (GW) sources with LIGO/Virgo/KAGRA (LVK) are ongoing and have to date revealed sources of both binary black holes (BBHs; Abbott et al. 2016a,b,c, 2017a,b; Venumadhav et al. 2020), binary neutron stars (BNSs; Abbott et al. 2017c), as well as possible black hole (BH) neutron star mergers (Hamers et al. 2021; Vynatheya & Hamers 2022; Abbott et al. 2023). However, how and where these sources form in our universe remain major unsolved problems. Many formation channels have been proposed, including field binaries (Dominik et al. 2012, 2013; Kinugawa et al. 2014; Dominik et al. 2015; Belczynski et al. 2016a,b; Murguia-Berthier et al. 2017; Silsbee & Tremaine 2017; Rodriguez & Antonini 2018; Schrøder, Batta & Ramirez-Ruiz 2018), dense stellar clusters (Portegies Zwart & McMillan 2000; Banerjee, Baumgardt & Kroupa 2010; Lee, Ramirez-Ruiz & van de Ven 2010; Tanikawa 2013; Bae, Kim & Lee 2014; Rodriguez et al. 2015; Rodriguez, Chatterjee & Rasio 2016a; Rodriguez et al. 2016b,b; Askar et al. 2017; Park et al. 2017; Samsing 2018; D’Orazio & Samsing 2018; Samsing & D’Orazio 2018; Samsing et al. 2020; Rozner & Perets 2022), active galactic nucleus (AGN) accretion discs (Bartos et al. 2017; Stone, Metzger & Haiman 2017; McKernan et al. 2018; Tagawa, Haiman & Kocsis 2020; Rozner, Generozov & Perets 2023), galactic nuclei (O’Leary, Kocsis & Loeb 2009; Hong & Lee 2015; Antonini & Rasio 2016; Stephan et al. 2016; VanLandingham et al. 2016; Hamers et al. 2018; Hoang et al. 2018), very massive stellar mergers (Loeb 2016; Woosley 2016; Janiuk et al. 2017; D’Orazio & Loeb 2018), and single-single GW captures of primordial BHs (Bird et al. 2016; Carr, Kühnel & Sandstad 2016; Cholis et al. 2016; Sasaki et al. 2016). To disentangle these different scenarios using GWs, we have to understand how the observable parameters differ between the different channels. For example, studies indicate that one can distinguish dynamically induced mergers from isolated binaries, by analysing the relative spin orientation of the merging BHs (Rodriguez et al. 2016c), the orbital eccentricity at some reference GW frequency (Gültekin, Miller & Hamilton 2006; Samsing, MacLeod & Ramirez-Ruiz 2014; Samsing & Ramirez-Ruiz 2017; Rodriguez et al. 2018; Samsing 2018; Samsing & D’Orazio 2018; Samsing & Ilan 2018; Samsing, MacLeod & Ramirez-Ruiz 2018a; Samsing, Askar & Giersz 2018b; Zevin et al. 2019; Samsing, Hamers & Tyles 2019b; Samsing et al. 2020), as well as the mass spectrum (Zevin et al. 2017). The environment in which the BBH was formed and merged can also be imprinted in the GW form, showing up as, for example, a GW phase shift (e.g. Inayoshi et al. 2017; D’Orazio & Loeb 2020; Hendriks, Zwick & Samsing 2024; Samsing et al. 2024). Other probes of formation include, for example, stellar tidal disruptions (e.g. Kremer et al. 2019b; Lopez Martin et al. 2019; Samsing et al. 2019a). From these studies it generally follows that dynamically formed mergers tend to have mass ratios near one (e.g. Rodriguez et al. 2018), random spin orientations (e.g. Rodriguez et al. 2016c), as well as a non-negligible fraction of sources with measurable eccentricity in LISA (D’Orazio & Samsing 2018; Samsing & D’Orazio 2018; Kremer et al. 2019a), DECIGO/Tian-Qin (e.g. Chen & Amaro-Seoane 2017; Samsing et al. 2020), and LVK (Samsing 2018). This is in contrast to isolated binary mergers, which are more likely to have a non-random spin distribution (e.g. Kalogera 2000; Hotokezaka & Piran 2017; Zaldarriaga, Kushnir & Kollmeier 2018; Piran & Piran 2020), larger mass ratios, and eccentricity |$\approx 0$| near LISA, DECIGO/Tian-Qin, and LVK.

However, some of the observed GW mergers start to challenge these classical formation channels. For example, GW190521 (Abbott et al. 2020b,c) seems to have masses that are either above the so-called BH mass gap (e.g. Fishbach & Holz 2020), or one above and one below (Fishbach, Farr & Holz 2020; Nitz & Capano 2021), which certainly was not generally expected from any models prior to this observation. In addition, GW190521 suggests better consistency with an eccentric waveform with |$e\ge 0.1$| at 10 Hz, rather than a quasi-circular one (e.g. Romero-Shaw et al. 2020; Gayathri et al. 2022; Romero-Shaw, Lasky & Thrane 2022), and a possible corresponding electromagnetic (EM) counterpart (e.g. Graham et al. 2020). These features have led to the proposal that GW190521 could have formed in an AGN disc, where BHs can grow to high masses through gas accretion (e.g. McKernan et al. 2012; Gilbaum & Stone 2022) or repeated mergers (Tagawa et al. 2021). They can also encounter each other in disc-like configurations that have been shown to produce up to two orders of magnitude more eccentric mergers compared to isotropic environments (Samsing et al. 2022, from hereon referred to as S22), as well as interact with the surrounding gas to create possible EM counterparts (e.g. Graham et al. 2023).

This observation, the underlying rich physics that brings together stellar-origin BHs, SMBHs, gas, dynamics and disc-accretion physics, have opened up a wealth of studies on how BBHs might form and merge in such AGN–disc environments (Tagawa, Haiman & Kocsis 2020). The general picture for this formation channel is that BHs can be formed through different mechanisms. They can be captured by the disc (e.g. Bartos et al. 2017; Fabj et al. 2020; MacLeod & Lin 2020; Generozov & Perets 2023; Nasim et al. 2023; Wang, Zhu & Lin 2024), form in the disc as a result of star formation in the outskirts (Stone et al. 2017; Artale et al. 2019), or they can be the result of core collapse of stars that undergo rapid mass accretion in high-density regions of the accretion disc (Cantiello, Jermyn & Lin 2021; Dittmann, Cantiello & Jermyn 2021; Jermyn et al. 2021, 2022; Dittmann, Jermyn & Cantiello 2023).

We envision this population of BHs to then migrate through the disc during which they may pair up to form binaries. These binaries are either brought to merge in the disc through a combination of gas drag and GW radiation (Tagawa et al. 2020; Dempsey et al. 2022; Li & Lai 2022, 2023, 2024; Calcino et al. 2024; Dittmann, Dempsey & Li 2024), or surviving inside the disc for them to later encounter either other BBHs or single BHs (Tagawa et al. 2020; Samsing et al. 2022). In S22, it was argued that three-body interactions, i.e. interactions between assembled BBHs and single BHs in the disc dominate the formation of BBH mergers, as each interaction hardens the BBH, possibly increase its eccentricity, and thereby bring the BBHs to merge on a short time-scale. In S22, the first AGN-disc like three-body simulations were performed with the inclusion of post-Newtonian (PN) radiation terms (e.g. Blanchet 2006, 2014), which showed that such interactions lead to an exceptionally high fraction of eccentric mergers.

Other studies focusing on the formation process of the BBH itself under the influence of an SMBH, through processes known as Jacobi captures, have also found that these BHs tend to pair up and possibly merge with very high eccentricity (Boekholt, Rowan & Kocsis 2023). Recent work (Trani, Quaini & Colpi 2024) has performed complementary studies on how BHs interact in stellar-remnant disc environments that do not necessarily move on Keplerian orbits as expected in AGN discs, and find that eccentric mergers are formed in large fraction. Recent studies have started to explore the effect of a gaseous disc in the problem of BH pairing with the inclusion of gas (DeLaurentiis, Epstein-Martin & Haiman 2023; Li et al. 2023; Rowan et al. 2023, 2024; Rozner et al. 2023; Whitehead et al. 2024). While these studies are encouraging, the problem of gas friction exerted by the accretion disc when modelling the PN three-body scattering problem in AGNs is a rather challenging task. Therefore, in this work we neglect the role of gas and focus on the other prominent environmental effect in AGNs: the tidal field exerted by the SMBH.

In this paper, we take the next step in building up our understanding of how BBHs might be brought to merger in disc-like environments, by performing PN simulations of single BHs interacting with BBHs in disc-like configurations with the inclusion of the central massive BH in the equation of motion. This is a necessary extension to the work by S22, who did not include the effect from the tides due to the SMBH in the scatterings. In addition, we adopt a set of initial conditions (ICs) that are closer to the ones that are expected in disc-like environments, where the single and binary slowly approach each other until they interact through their common Hill sphere. The question that naturally arises is whether the tidal field and the more constrained ICs drastically decrease the high number of eccentric mergers found in the aforementioned study, or whether instead eccentric mergers are a robust indicator of BBHs assembled in disc-like environments. In this paper, we address the previous questions, present results that are directly comparable to the ones from S22, and introduce new important elements to acquire a better understanding of this problem.

With these motivations, we start in Section 2.1 by reviewing the theory of assembling eccentric mergers through binary–single interactions taking place in disc-like environments. We then in Section 2.2 extend this theory to include the effect from a tidal boundary caused by the presence of the nearby SMBH. After this, we move on to performing PN interactions between a BBH and single BH all orbiting an SMBH (Section 4.1), for which we present outcome distributions (Sections 4.2 and 4.3) and phase-space diagrams in Section 4.4. We especially focus on the outcome of probabilities for eccentric mergers as a function of BBH semimajor axis (SMA) relative to the size of the Hill sphere, as well as how the results depend on the inclination angle between the BBH and the incoming single BH (Section 4.5). Finally, we conclude and highlight future directions in this problem.

2 THREE-BODY INTERACTIONS IN A TIDAL FIELD

The three-body scattering problem has been studied in great detail with different physical motivations, from the context of classical cluster evolution (Heggie 1975) to including PN correction terms to follow the formation of GW driven mergers of compact objects (Samsing et al. 2014; Samsing 2018). While several theoretical studies have been put forward for describing especially the fraction of eccentric merges forming in chaotic three-body interactions (e.g. Samsing et al. 2014; Rodriguez et al. 2016a; Samsing 2018), only a few have extended this theory to include how a nearby perturber, such as an SMBH, affects the range, fraction, and nature of outcomes (Leigh et al. 2018; Trani et al. 2019, 2024; Ginat & Perets 2021; Rom, Sari & Lai 2024 ). In S22, it was shown analytically and numerically that co-planar scatterings, such as those that might take place in AGN–disc environments, naturally give rise to an exceptionally high fraction of eccentric GW mergers due to the controlled geometric set-up, however, neglecting the effect of the SMBH. We start by reviewing the relevant three-body scattering theory, after which we derive how the presence of a nearby object to leading order has an effect on the fraction of eccentric mergers forming during the chaotic state of the three-body system. This not only adds an important component to this theoretical framework, but also sets the stage and provides motivation for PN scatterings with an SMBH that we present later in the paper.

2.1 Formation of eccentric mergers

We start by considering Fig. 1, which shows the evolution of an isolated binary–single scattering, i.e. without a nearby perturber, that goes through several resonances, or intermediate binary–single states (IMS), before one of them undergoes a GW inspiral merger while being bound to the remaining single. This outcome we refer to as a three-body merger (Samsing 2018). The condition for this to happen is that the inspiral time of the IMS binary, |$t_{b}$|⁠, has to be less than the time for the remaining single to return, |$T_{bs}$|⁠. In the high eccentricity limit, the merger time of the IMS binary can be approximated by (Peters 1964)

(1)

where |$a_b$|⁠, m, and |$r_m$|⁠, denote the initial IMS binary SMA, individual BH mass (throughout the paper we assume the triple system to be equal mass), and pericentre distance, respectively. The time-scale of the binary–single orbital time |$T_{bs}$| scales with the orbital time of the initial binary such that

(2)

up to a constant factor that does not play a role here for our purpose. By setting |$t_m$| equal to |$T_{bs}$|⁠, one can isolate for the critical pericentre distance, |$r_m$|⁠, that the two IMS binary objects need to have upon formation for their merger time to be short enough for them to merge while still being bound to the single. This leads to

(3)

where |$\mathscr {R}$| here denotes the Schwarzschild radius of one of the BHs with mass m. To get a sense of the scales, one finds that |$r_m/\mathscr {R} \approx 100$|⁠, for |$m = 20{\rm M}_{\odot }$| and |$a_b = 1\ \rm au$|⁠, i.e. if the IMS binary is formed with a pericentre distance |$\lt 100 \times \mathscr {R}$| then it is likely to undergo a GW capture merger for this set-up (Samsing 2018).

Illustration of binary–single interaction propagating from left to right resulting in a merger between two of the three objects while the remaining single is still bound (see Section 2.1). The system goes through several intermediate states (IMS) before GW capture of two of the three objects. For illustrative purposes, we include horizontal dashed lines to indicate the size of a Hill sphere ($R_{\rm H}$) in case the interaction takes place in the presence of a tidal field. If the interaction extends beyond $R_H$, the system becomes unbound and undergoes tidal breakup. The tidal boundary therefore limits the number of interactions $\mathcal {N}$ and it increases the probability of breaking up the system before natural completion. This consequently leads to a reduction in the probability for three-body and eccentric mergers. The derived correction factor is shown in equation (22).
Figure 1.

Illustration of binary–single interaction propagating from left to right resulting in a merger between two of the three objects while the remaining single is still bound (see Section 2.1). The system goes through several intermediate states (IMS) before GW capture of two of the three objects. For illustrative purposes, we include horizontal dashed lines to indicate the size of a Hill sphere (⁠|$R_{\rm H}$|⁠) in case the interaction takes place in the presence of a tidal field. If the interaction extends beyond |$R_H$|⁠, the system becomes unbound and undergoes tidal breakup. The tidal boundary therefore limits the number of interactions |$\mathcal {N}$| and it increases the probability of breaking up the system before natural completion. This consequently leads to a reduction in the probability for three-body and eccentric mergers. The derived correction factor is shown in equation (22).

To estimate the actual probability that a given IMS binary undergoes such a merger, we have to consider the distribution of SMA and eccentricity for a given IMS, as |$r_m$| is directly linked to these quantities through |$r_m = a_b(1-e_b)$|⁠. In the equal mass case, the binary SMA, |$a_{b}$|⁠, is of order the initial SMA, |$a_0$|⁠, and we (for this part) therefore assume that

(4)

The IMS binary eccentricity, |$e_b$|⁠, on the other hand can vary greatly (Monaghan 1976; Valtonen & Karttunen 2006; Stone & Leigh 2019) and in the co-planar chaotic limit follows the distribution,

(5)

This distribution implies that the probability that an IMS binary in the co-planar limit forms with a pericentre distance |$\lt r$| is given by

(6)

which in the case of an IMS binary merger set by the condition |$r \lt r_m$| translates to the relation,

(7)

where we have used the notation from equation (3). Now, this is the probability that a given IMS undergoes a merger in the co-planar limit. To get the full probability that the entire binary–single interaction results in a three-body merger, we have to take into account that the interaction goes through several IMS, a number we denote |$\mathcal {N}$| (Samsing et al. 2014). In the limit where the probability for a single IMS to undergo a merger is |$\ll 1$|⁠, the total probability for a three-body merger is simply the probability for one IMS binary to undergo a merger, i.e. |$P(r \lt r_m)$|⁠, times the number of intermediate states, |$\mathcal {N}$|⁠,

(8)

where |$\mathcal {N} \approx 20$| as discussed further in the sections below.

For a three-body merger that is always highly eccentric at formation, the corresponding GW frequency where most of the power is emitted, also referred as the GW peak frequency (D’Orazio & Samsing 2018), is related to the binary pericentre |$r_f$| distance as

(9)

which implies that if an IMS binary forms with a pericentre less than

(10)

it consequently shows up eccentric at GW frequency that is higher than |$f_{\rm gw}$|⁠. It might be that the merger forms with a GW peak frequency that is below the observable band of LVK (Samsing et al. 2014). In this case, one has to evolve the binary until it reaches the band through the emission of GWs; however, at this point the eccentricity might be too low to be measured. Therefore, three-body mergers are not necessarily representative of the fraction of observable highly eccentric mergers in e.g. LVK as discussed in Samsing (2018). For an IMS merger to clearly show up as an eccentric merger, its GW peak frequency at formation has to be near the observable band.

With the critical distance from equation (10), one can easily estimate the probability for the outcome in a similar way to what we did in equation (8), which is |$\approx \mathcal {N} \times \sqrt{2r_f/a_0}$|⁠. In the more general case for the IMS binary to appear at GW frequency |$f_{\rm gw}$| with an eccentricity |$e_f$|⁠, the probability should instead be based on the critical pericentre distance,

(11)

as discussed in S22. Here, the first term is simply |$r_f$| from equation (10), which implies that the ratio |$r_{e,f}/r_f$| only depends on the eccentricity threshold at |$f_{\rm gw}$|⁠. For example, for |$e_f = 0.1$|⁠, |$r_{e,f}/r_f \approx 3$|⁠, meaning that if the IMS binary is assembled with a pericentre distance that is |$\approx 3$| times larger than |$r_f$|⁠, it appears with an eccentricity 0.1 at GW peak frequency |$f_{\rm gw}$|⁠. Generally, one should note that for an eccentric source to also merge the pericentre has to be less than |$r_{m}$|⁠, which indeed is the case for most astrophysical systems near the LVK band.

With relevant astrophysical scales, we can summarize these analytical estimations as

(12)

and

(13)

as further explained in S22. The question is how these general results, which also have shown excellent agreement with numerical simulations (S22), change when the chaotic three-body interactions take place in a tidal field from a nearby perturber, such as an SMBH. We estimate how this addition affects the outcomes.

2.2 Effects from tidal fields

The effect from a nearby object shows up (radially) to leading order in the three-body dynamics as a tidal field across the spatial scattering domain. Although in this analysis we perform actual three-body PN simulations with a nearby SMBH, we first explore analytically how a tidal field might impact the results from the above section. For this we imagine the interaction to be taking place in a tidal field that we assume can be represented by the size of the Hill radius,

(14)

where |$R_{\rm CM}$| is the distance of the binary centre of mass (CM) from the central BH (as later described in Section 3.1), m is the stellar-origin BH mass, and |$M_{\rm BH}$| is the SMBH mass (as we consider equal mass stellar-origin BHs). By definition, outside the sphere the tidal forces from the perturber separate the interacting objects, which then lead to a termination of the scattering. An example of such a sphere, or distance, is also shown in Fig. 1 (see figure caption). Assuming the interaction inside the Hill sphere is unaffected by the perturber, the effect from a tidal field on the three-body scattering is therefore to terminate the interaction before it would reach its natural end-state (defined to be the outcome in the isolated case). If we translate this into the quantities we use in the theoretical model from Section 2.1, this means that a tidal field reduces the number of IMS binaries |$\mathcal {N}$|⁠, which would therefore lead to a reduction in the probability of both three-body mergers and eccentric GW mergers. In other words, there are less chances, i.e. |$\mathcal {N}$|⁠, to undergo a three-body merger when tides can terminate the binary–single interaction. In the following, we estimate how tides entering in this Hill sphere picture affects the number |$\mathcal {N}$|⁠, and thereby the resultant merger probabilities.

The dynamics giving rise to the three-body system temporarily splitting into a successive series of |$\mathcal {N}$| IMS binaries and bound singles, relate to the energy exchange between the binary and the incoming single. In the initial configuration, the incoming single is unbound with respect to the binary, but after the first close interaction it most likely happens that the single delivers some of its dynamical energy to the binary, which then expands slightly as a response. If the energy exchange is larger than the initial energy between the single and the binary, the single gets temporally bound to the binary. This bound state is what we call an IMS. The IMS can essentially be considered as a regular ‘binary’, with one object being the ‘IMS binary’ (with mass |$= 2m$|⁠) and the other the ‘bound single’ (with mass |$= 1m$|⁠), with an SMA |$a_{bs}$| that is given by (Samsing et al. 2014)

(15)

where |$|E_0|$| and |$|E_b|$| here denote the initial energy of the three-body system and the energy of the IMS binary, respectively. Relating this to the Hill sphere boundary, we now note that the interaction terminates if |$a_{bs} \gt R_{H}$|⁠. To calculate the correction from the Hill sphere, we therefore have to calculate the probability that the three-body system splits into an IMS with |$a_{bs} \gt R_{H}$|⁠. For this, we first use that the total initial energy of the three-body system in the considered hard binary limit is about the binding energy of the initial binary, i.e. |$|E_0| \approx Gm^2/(2a_0)$|⁠. In this limit, the energy of a given IMS binary relates to the binary + single SMA |$a_{bs}$| by

(16)

where we have used equation (15). By now setting |$a_{bs} = R_{H}$|⁠, we find

(17)

which defines the IMS binary energy at which the system splits into a binary + single state with |$a_{bs} = R_{\rm H}$|⁠. To relate this characteristic energy threshold to a probability, we first have to consider the distribution of |$|E_b|$| for the assembly of each IMS binary. For this we make use of the theory presented in Valtonen & Karttunen (2006) and Stone & Leigh (2019), which through arguments related to statistical mechanics derives that in the co-planar case the binary energy |$|E_b|$| approximately follows the distribution,

(18)

To get a better understanding of how these energies relate to different three-body outcomes, we now consider Fig. 2, which shows |$P(|E_b|)$|⁠, with three characteristic limits. Generally, the more compact the binary in the IMS is, i.e. the larger |$|E_b|$| is, the larger the orbit of the single with respect to the binary is. As shown in the illustrations at the bottom of the figure, this exchange of energy gives rise to three distinct states and outcomes: (1) ‘Bound State’: Here, the energy of the binary |$|E_b|$| is so low that the single is on an orbit that has a comparable SMA, |$a_{bs}$|⁠, to the SMA of the binary, |$a_{b}$|⁠. The limit |$|E_{21}|$| therefore denotes the minimum energy the binary can have for the system to be described as an IMS. (2) ‘Tidal Breakup’: At the critical energy |$|E_H|$|⁠, the SMA of the bound single with respect to the binary, |$a_{bs}$|⁠, is similar to the Hill sphere distance, |$R_H$|⁠. If |$|E_b|$| is greater than this value, the system experiences a tidal breakup. Therefore, if |$|E_{21}|\lt |E_b|\lt |E_H|$|⁠, then the system continues interacting as an IMS, whereas if |$|E_b|\gt |E_H|$| the interaction terminates. (3) ‘Isolated Unbound’: When there is no tidal field, the system naturally breaks up when the energy of the single with respect to the binary is positive, i.e. when it is being sent out on an unbound orbit by the binary. In the isolated case this happens when |$|E_b|\gt |E_0|$| (in the hard binary limit). The ratio between the different areas enclosed by the energy limits (⁠|$|E_{21}|,|E_{H}|,|E_{0}|,|E_{\infty }|$|⁠), relates to the outcome probabilities for a given formed IMS binary. We use this picture to calculate how the effect from a tidal field of size |$R_{\rm H}$|⁠, as a tidal breakup limits the possibility for forming three-body and eccentric mergers by breaking the resonance chain of IMS binaries.

Illustration of how the distribution of binary energies, $P(|E_b|)$, as a function of $|E_b|$ relates to different interaction and outcome states of the three-body system when the interactions take place in a tidal field characterized by a Hill sphere of size, $R_{\rm H}$. As further described in Section 2.2, $|E_{21}|$ denotes the minimum energy the binary can have for the three-body system to still be defined as a binary with a bound single, i.e. as a ‘2 + 1’ state. $|E_{H}|$ is the IMS binary energy at which the SMA of the binary + single, $a_{bs}$, equals the Hills sphere, $R_{\rm H}$. $|E_{0}|$ is the initial energy of the IMS binary, and defines the limit of escape in the isolated binary–single interaction problem. The ratio between the different areas enclosed by the energy limits is related, for a given formed IMS binary, to the outcome probabilities. This is illustrated in equations (19) and (20).
Figure 2.

Illustration of how the distribution of binary energies, |$P(|E_b|)$|⁠, as a function of |$|E_b|$| relates to different interaction and outcome states of the three-body system when the interactions take place in a tidal field characterized by a Hill sphere of size, |$R_{\rm H}$|⁠. As further described in Section 2.2, |$|E_{21}|$| denotes the minimum energy the binary can have for the three-body system to still be defined as a binary with a bound single, i.e. as a ‘2 + 1’ state. |$|E_{H}|$| is the IMS binary energy at which the SMA of the binary + single, |$a_{bs}$|⁠, equals the Hills sphere, |$R_{\rm H}$|⁠. |$|E_{0}|$| is the initial energy of the IMS binary, and defines the limit of escape in the isolated binary–single interaction problem. The ratio between the different areas enclosed by the energy limits is related, for a given formed IMS binary, to the outcome probabilities. This is illustrated in equations (19) and (20).

With this formalism we can now estimate the number of IMS, |$\mathcal {N}$|⁠, with and without a tidal field. In the absence of the tidal field, |$\mathcal {N}$| is given by the ratio between the bound-state area, and the escaper-area,

(19)

where we have used the distribution from equation (18). In the case of system constrained by a Hill sphere |$R_{\rm H}$|⁠, the corresponding number of IMS, |$\mathcal {N}_H$|⁠, is instead given by the ratio:

(20)

Since the probability for the merger types we consider is |$\propto \mathcal {N}$|⁠, as seen in, for example, equation (8), we can now deduce that when a system is tidally limited, then the probability for merger during the interaction is reduced by the factor,

(21)

where we have used the relation from equation (17). From this we can conclude that when the interaction is limited by the Hill sphere, |$R_H$|⁠, the merger probabilities are corrected in the following way:

(22)

If we ask at what |$R_{\rm H}/a_0$| the probability has decreased by a factor of 2, one finds that happens when |$R_{\rm H}/a_0 \sim 7$|⁠. Therefore, if the binary is within a few times the Hill sphere to start out with, then the outcomes from isolated interactions without the tidal perturber are expected to return probabilities that are accurate within a factor of order of unity. With a solid theoretical understanding and models to test, in the section below we treat this problem with PN N-body simulations.

3 MODEL AND METHODS

Having described and considered various theoretical aspects of how tides might impact the chaotic interaction of binary–single scatterings and their outcomes, we now turn to exploring this problem in more detail using PN N-body simulations. As the key interest is related to how eccentric mergers are produced at high rates in disc-like environments, where the population of co-planar objects is significant, we start here by focusing on the co-planar set-up. After this we move on to out-of-plane configurations in Section 4.5. In the following sections, we present our adopted model, numerical methods, and results from N-body scatterings including a nearby SMBH.

3.1 Initial configuration of binary–single system orbiting the SMBH

In this work, we directly expand on the analysis of S22, most notably on the result that in general, co-planar scatterings result in order-of-magnitude more eccentric BBH mergers compared to the normal isotropic case taking place in, for example, stellar clusters (Samsing 2018). However, the results in S22 were based on scatterings without a nearby SMBH, and did as well include outcomes from a range of impact parameters that do not necessarily reflect the conditions expected in an AGN disc. Motivated by this, we here extend their analysis to now also include a nearby SMBH and ICs that likely are closer to the ones taking place in AGN accretion-disc-like environments. More specifically, we here perform scattering experiments between a BBH and a single BH while they orbit an SMBH with mass |$M_{\rm BH}$| on a circular orbit with their CM placed at a radius |$R_{\rm CM}$|⁠, as further discussed and shown in Fig. 3. Since we expect embedded orbiters to move subsonically, we set the initial eccentricity of the binary system |$e_{b}=0$|⁠, as we imagine that the binary in the subsonic regime has been circularized through (previous) interactions with gas (e.g. O’Neill et al. 2024). The BBH and the single BH are placed at the beginning of each simulation on their own circular orbits around the SMBH. The single BH is put at the gravitational sphere of influence of the triple with respect to the SMBH, which is set by the Hill radius shown by blue circle in Fig. 3 and its size is set by equation (14). This estimate for the tidal influence of the central BH on the three-body system is only used for setting up the ICs, as the evolution of the few-body systems naturally includes effects from the SMBH through the N-body solver.

Overall (left) and close-up (right) schematic representation of the scatterings configuration. The SMBH is at the centre, and the BBH CM is placed at a distance $R_{\rm CM}$ from the SMBH (red line on the left side). The triple BH system (BBH formed by $m_1$ and $m_2$, and incoming single BH $m_3$) is here on a co-planar circular orbit with respect to the SMBH (2D configuration). The circle represents the Hill sphere with radius defined in equation (14). We test two different configuration set-ups. Setup I (blue): single BH is approaching the binary from the right; Setup II (pink): the single BH is placed closer to the SMBH and approaches therefore the binary from the left. We test both BBH CR and CT with respect to the orbit around the SMBH. For a fixed SMA, we vary the initial binary phase angle indicated by $\phi _0$.
Figure 3.

Overall (left) and close-up (right) schematic representation of the scatterings configuration. The SMBH is at the centre, and the BBH CM is placed at a distance |$R_{\rm CM}$| from the SMBH (red line on the left side). The triple BH system (BBH formed by |$m_1$| and |$m_2$|⁠, and incoming single BH |$m_3$|⁠) is here on a co-planar circular orbit with respect to the SMBH (2D configuration). The circle represents the Hill sphere with radius defined in equation (14). We test two different configuration set-ups. Setup I (blue): single BH is approaching the binary from the right; Setup II (pink): the single BH is placed closer to the SMBH and approaches therefore the binary from the left. We test both BBH CR and CT with respect to the orbit around the SMBH. For a fixed SMA, we vary the initial binary phase angle indicated by |$\phi _0$|⁠.

The illustration to the right in Fig. 3 shows a close-up of the initial configuration. For each scattering experiment, we vary the initial phase of the binary angle |$\phi _0$|⁠. In addition, we explore co-rotating (CR, green arrow) and counter-rotating (CT, orange arrow) configurations, both of which are believed to take place in AGN disc environments (Secunda et al. 2021). The definition of CR and CT is naturally set by the direction of motion of the triple system with respect to the SMBH. With this, we explore two different set-ups as shown in the right panel of Fig. 3:

  • Set-up I (blue): The single BH is initially placed at the outer boundary of the Hill sphere w.r.t. the SMBH, so that its distance from the SMBH is |$R_{\rm CM} + R_{\rm H}$|⁠.

  • Set-up II (pink): The single BH is initially placed at the inner boundary of the Hill sphere w.r.t. the SMBH, so that its distance from the SMBH instead is |$R_{\rm CM}-R_{\rm H}$|⁠.

The most likely initial configuration depends on the astrophysical settings, e.g. if the binary and the single meet due to different migration speeds inside the disc, Set-up II is the most likely as the heavier BBH ‘sweeps up’ the lighter single BH as they both migrate through the disc. On the other hand, if the interactions are taking place in a migration trap, Set-up I is the most likely, as the single BH approaches the BBH from the outer parts (Secunda et al. 2019). If one allows for significantly different masses, both Set-up I and Set-up II can take place throughout the disc. It is important to note that as we vary the position of the incoming single BH, the relative orbital velocity between the BBH and the BH is varied as well (indicated by the green arrows on the left side of Fig. 3). With these more physical realistic ICs, we now explore the possible differences between the results from S22 and the results from the scatterings.

3.2 Merger classification

The end-state of few-body scatterings can generally be classified and divided into a range of distinct outcomes depending on what environments the scatterings are taking place in, on the probed observables, and on how finite sizes as well as energy dissipative terms (GW radiation, tides, gas, etc.) are included (e.g. Samsing et al. 2017). In this work, we only report results for BBH mergers that have a notable eccentricity, e, at frequency |$f_{\rm gw}$| as these carry unique information about the nature of the formation channel. In addition, when focusing on eccentric BBH mergers in the LVK bands, the merger time is so short that we do not have to complicate the study by introducing and working with several outcome types. For example, if we considered eccentric LISA sources instead, the relevant population would be the BBHs that do not undergo a three-body merger but instead are left unbound after the interaction with a notable eccentricity (D’Orazio & Samsing 2018; Samsing & D’Orazio 2018). However, in the AGN disc, tracking this population is very difficult, as the evolution of such binaries after the interaction is likely to either be greatly affected by the gas, or by other objects embedded in the disc. In stellar clusters this is in contrast relatively easy as was demonstrated in Samsing (2018). Focusing on eccentric LVK sources in the adopted AGN set-up is more robust, as they form and merge promptly. For these reasons, in this paper we mainly focus on the prompt BBH mergers that form in our simulations and appear with notable eccentricity (>0.1) in the LVK bands (⁠|$\sim 10\ \rm Hz$|⁠).

3.3 Numerical methods and PN corrections

We perform simulations with the same extensively tested N-body code as in Samsing et al. (2014, 2017, 2022) and Samsing (2018), which adopts a solver based on lsode with a relative and absolute error set to |$10^{-12}$| and an adaptive time-step with a scale based on the relative position, velocity, and acceleration of the interacting bodies following standard procedures. The code includes GR effects by the use of the PN formalism (Blanchet 2006, 2014). In this formalism, the effects from GR are added to the Newtonian acceleration as an expansion in |$v/c$|⁠,

(23)

where |$\boldsymbol{a}_0$| denotes the newtonian acceleration (0PN order), |$c^{-2} \boldsymbol{a}_2 + c^{-4} \boldsymbol{a}_4$| (1PN + 2PN order) are energy conserving and lead to orbital precession (Blanchet 2006), and finally |$\boldsymbol{a}_5$| is the dissipative term that in this formalism is the leading order term describing the energy and momentum loss due to GW radiation. When doing an orbit average, the effect from this term is essentially equivalent to the well-known Peters (1964) formulae. In this work, we do not include the conservative 1PN, 2PN corrections, as these terms are essential for evolving hierarchical systems where the interactions are taking place over many orbits, but have been shown not to play a significant role in the type of chaotic scatterings we are exploring here. We focus instead on the leading effect from GW energy dissipation for assembling BBH mergers, i.e. the terms we include in the N-body code are |$\boldsymbol{a} = \boldsymbol{a}_0 + c^{-5} \boldsymbol{a}_5$|⁠. This PN-acceleration is applied pairwise between any two objects ‘1’ and ‘2’ such that the acceleration terms can be written as

(24)

and

(25)

where bold symbols indicate vectors. Using the same notation as in Blanchet (2006), the separation vector is defined here as |$\boldsymbol{r}_{12}= \boldsymbol{r}_{1}-\boldsymbol{r}_{2}$| with corresponding relative velocity vector |$\boldsymbol{v}_{12}$|⁠, and unit vector |$\hat{\boldsymbol{r}}_{12}$| given by |$\boldsymbol{r}_{12}/r_{12}$|⁠. We note that PN terms have been derived up to higher orders; however, it is not clear how the series converges, and for stability we restrict ourselves not to include higher order terms.

4 RESULTS

We here present the main results on the formation of eccentric BBH mergers forming as a consequence of binary–single interactions taking place near an SMBH in a disc-like set-up. We start by showing a few illustrative cases where the tidal field from the SMBH clearly plays a role in the interaction and outcome. We then show the probability for forming eccentric mergers as a function of the BBH SMA, relative to the Hill sphere of the SMBH, for both CR and CT interactions. These results are compared to the scatterings by S22.1 The new features we find are then explored by considering the phase-space, or topology, of the scatterings. Finally, we show results for out-of-plane scatterings, as the binary–single orbital inclination is believed to be critical in the production of eccentric mergers as shown in S22.

4.1 Escape versus capture binary–single interactions

Fig. 4 shows two different examples of three-body scatterings for Set-up I, where the single BH is approaching the BBH from the outside, as earlier described in Fig. 3. Both scatterings are plotted with respect to the binary + single CM. The top panel displays an interaction that is terminated before its natural outcome as a result of the tidal field from the SMBH, which breaks up the binary as the interaction here leads to an IMS that extents outside of the Hill sphere (blue circle). For this scattering, the initial binary SMA is 1 au and the distance from the SMBH is chosen such that the size of the Hill radius is 10 au. The initial binary phase angle |$\phi _0$| is 0|$^{\circ }$|⁠. For the experiment in the bottom panel, the SMA and Hill radius are kept the same as in the left, but the binary angle is instead 135|$^{\circ }$|⁠. As shown, this variation in the binary phase leads to a completely different outcome, namely the merger between two of the three objects (pink, dark blue) while the single BH is still bound, i.e. the outcome we denote a ‘three-body merger’. This outcome is not possible in a purely Newtonian code, as the inspiral of the two objects undergoing the merger is entirely driven by the radiation arising from the 2.5PN term. Comparing the two shown scatterings, we clearly see that some interactions almost promptly undergo a merger (bottom panel), whereas some enter the resonating state as we saw in Fig. 1. As described in Section 2.2, when the system enters these chaotic states the tidal field from the SMBH is likely to split it apart before, for example, a three-body merger is formed. If gas is present one can imagine the drag forces acting continuously on the bottom panel interaction, which would gradually transfer energy out of the entire system. Such processes could essentially protect the system from undergoing a tidal breakup as it would harden over time. In this simple picture, the gas in the accretion disc could therefore help increase the number of three-body mergers that otherwise could have terminated earlier by the SMBH tidal field.

Examples of equal-mass binary–single interaction under the influence of the SMBH resulting in an escape from the Hill sphere (top) and a three-body GW merger (bottom) for Set-up I. The dots indicate the initial positions. The encounter is plotted with respect to the three-body CM frame, where the single BH is placed at the Hill sphere indicated by the circle with a radius set to be 10 au. The initial binary SMA is 1 au for both experiments. The initial binary phase in the top panel is set to $0^{\circ }$ and $135^{\circ }$ for the bottom panel, indicating that with fixed initial SMA and $R_{\rm CM}$, varying the BBH phase can easily change the scattering outcome.
Figure 4.

Examples of equal-mass binary–single interaction under the influence of the SMBH resulting in an escape from the Hill sphere (top) and a three-body GW merger (bottom) for Set-up I. The dots indicate the initial positions. The encounter is plotted with respect to the three-body CM frame, where the single BH is placed at the Hill sphere indicated by the circle with a radius set to be 10 au. The initial binary SMA is 1 au for both experiments. The initial binary phase in the top panel is set to |$0^{\circ }$| and |$135^{\circ }$| for the bottom panel, indicating that with fixed initial SMA and |$R_{\rm CM}$|⁠, varying the BBH phase can easily change the scattering outcome.

4.2 Probability of eccentric black hole mergers

Fig. 5 shows the probability for the considered scatterings to result in an eccentric merger (⁠|$e\gt 0.1$|⁠) with a peak frequency |$f_{\rm gw} \ge 10$| Hz (LVK band), as a function of the SMA of the target BBH. For this figure we have assumed |$m = 100 \, {\rm M}_{\odot }$|⁠, |$M_{\rm BH} = 10^{8} \, {\rm M}_{\odot }$|⁠, and that the single encounters the binary from a distance equal to a Hills sphere |$R_{\rm H}$| = 10 au, which corresponds to |$R_{\rm CM} = 10^3$| au. For each point in the figure we vary the phase angle |$\phi _0$| from 0 to 2|$\pi$| in N steps, where |$N = 2 \times 10^5$| for the results presented here. The probability is estimated by taking the total number of outcomes resulting in an eccentric inspiral (⁠|$e\gt 0.1$| at |$f_{\rm gw} \ge 10$| Hz) to the total number of scatterings per SMA (⁠|$N = 2 \times 10^5$|⁠).

Probability for forming a GW merger with $e\gt 0.1$ at peak frequency $f_{\rm gw} \ge$ 10 Hz, for CR (red) and CT (blue) BBHs in set-up I, for $m = 100 {\rm M}_{\odot }$ triple BH system. The large black dots are the replicated results from S22. The grey black dots represent the correction factor due to the tidal field according to equation (22). At lower SMA (up to $\approx$0.3 au) the trend of the probability follows approximately a power law, that then breaks after which the curve becomes more fluctuating as the size of the SMA approaches the size of $R_{\rm H}$. Further analysis in Section 4.3.
Figure 5.

Probability for forming a GW merger with |$e\gt 0.1$| at peak frequency |$f_{\rm gw} \ge$| 10 Hz, for CR (red) and CT (blue) BBHs in set-up I, for |$m = 100 {\rm M}_{\odot }$| triple BH system. The large black dots are the replicated results from S22. The grey black dots represent the correction factor due to the tidal field according to equation (22). At lower SMA (up to |$\approx$|0.3 au) the trend of the probability follows approximately a power law, that then breaks after which the curve becomes more fluctuating as the size of the SMA approaches the size of |$R_{\rm H}$|⁠. Further analysis in Section 4.3.

Generally, we do not evolve each binary until merger, therefore to classify which BBHs are falling into the category of having an |$e\gt 0.1$| at |$f_{\rm gw} \ge 10$| Hz, we take the binary output from the N-body code, which essentially corresponds to the end-state binary SMA and eccentricity, and evolve it using Peters (1964) equations. Specifically, first we estimate what pericentre distance, |$r_f$|⁠, corresponds to a given GW peak-frequency, |$f_{\rm gw }$|⁠, using equation (10). After this, we take the binary SMA and eccentricity outputted from the code, here denoted by |$a_{\rm out}$| and |$e_{\rm out}$|⁠, respectively, and derive the constant |$c_{\rm out}$|⁠,

(26)

We then solve for the eccentricity at which the binary reaches a pericentre distance of |$r_f$|⁠, using the relation,

(27)

For the sources that form with a GW peak frequency above the imposed threshold, e.g. above 10 Hz, we label these as eccentric mergers as well (S22). Theoretically, these mergers belong to the category of being eccentric in the observable band, whereas from an evolutionary point of view, they could have very different kinds of GW waveforms. For example, mergers with extremely high eccentricity show up as burst sources with a GW signal similar to, for example, GW190521 (Abbott et al. 2020a), whereas less eccentric mergers gradually inspiral throughout the observable band. For the results in Fig. 5, we map both CR (red) and CT (blue) binaries from Set-up I. As one of the main goals with our analysis is to explore if the findings from S22, stating eccentric mergers are produced at exceptionally high probability in AGN environments, also holds for more realistic ICs and with the inclusion of an SMBH, we here also show data from S22 in black. From this we can now conclude that in terms of the overall scaling and dependence on SMA, the results with an SMBH and with the single approaching from the Hill sphere is similar to the results found in S22, which were based on isolated three-body scatterings (and no Hill’s boundary). This especially implies that the co-planar geometry that likely is facilitated by the AGN disc, is still highly effective in creating eccentric GW sources, despite the influence from the nearby SMBH. Looking closer at the scattering results, we do see some unique features in the data when the SMA of the BBH approaches the Hill sphere |$R_{\rm H}$|⁠, i.e. in this case when |$a_0$| approaches 10 au. In general, the more we increase the SMA, i.e. the binary size becomes comparable to the Hill radius, the more the probability trend becomes less linear and much more fluctuating. This is not due to numerical noise, but a real feature of the scatterings. There are even regions in which the probability increases with increasing binary size (e.g. after 1 au, where |$\rm a_0/\rm R_{\rm H}=0.1$|⁠), to then start decreasing again (after 3 au). Although a direct comparison cannot be made, Boekholt et al. (2023) also found clear (fractal) structures in their two-body Jacobi captures, which also results from keeping track and varying the ICs in a systematic way as we are also doing, but here in the three-body problem.

Finally, Fig. 6 shows the actual distribution of GW peak frequency returned from the scatterings at the time we stop the simulation, where the stopping criterion corresponds to either a GW capture inside the Hill sphere, or to the interruption of the simulation if the outcome is an escaper, where all three objects have left the Hill sphere (same as the two outcomes shown in Fig. 4). The chosen value for the SMA is 0.3 au, for which |$\phi _0$| is varied |$10^4$| times. The blue distribution indicates the binaries that underwent a GW capture inside the Hill sphere, while the pink distribution represents the population of binaries that survived the interaction. Both eccentric LISA and LVK mergers are forming, which would play a major role when multiband GW observations become possible with LISA, LVK, and possible 3G instruments operating simultaneously. Starting with the LVK sources, we see that a significant fraction of mergers are assembled everywhere from 10 to |$10^3$| Hz, implying that many of the mergers we are finding are closer to burst sources or direct plunge sources than to classical eccentric mergers that almost adiabatically evolve from high eccentricity to circular inspiral. Burst and plunge searches are therefore equally important to perform as eccentric searches for probing the contribution of disc-like environments. We note that such searches have been done (Klimenko et al. 2016; Abbott et al. 2020a) with no clear detections yet, except maybe for GW190521 (Abbott et al. 2020c).

Distribution of GW peak frequency of binaries formed through binary–single interactions near an SMBH evaluated at the time of assembly (see Section 4.2). The initial SMA is $a_0 = 0.3\,\mathrm{ au}$. The blue distribution represents the binaries that merged inside the Hill sphere (three-body mergers), while the pink distribution indicates binaries that survived the interaction. Overplotted are the LISA and LVK peak frequency bands (dashed vertical black lines), showing that a large fraction of sources enter both LISA and LVK bands at formation.
Figure 6.

Distribution of GW peak frequency of binaries formed through binary–single interactions near an SMBH evaluated at the time of assembly (see Section 4.2). The initial SMA is |$a_0 = 0.3\,\mathrm{ au}$|⁠. The blue distribution represents the binaries that merged inside the Hill sphere (three-body mergers), while the pink distribution indicates binaries that survived the interaction. Overplotted are the LISA and LVK peak frequency bands (dashed vertical black lines), showing that a large fraction of sources enter both LISA and LVK bands at formation.

4.3 Effects of periodic orbital encounters

To explore the nature of the fluctuating scattering results seen in Fig. 5, we first consider Fig. 7, which further expands on the results by showing the probability for an eccentric merger for two different BH masses (15 |${\rm M}_{\odot }$|⁠, top panels and 100 |${\rm M}_{\odot }$|⁠, bottom panels) and for Set-up I (left panels) and Set-up II (right panels). All four configurations show similar behaviour to the one described in Fig. 5, i.e. at lower SMA, the probability mostly follows a linear trend, while it becomes increasingly more fluctuating as we widen the BBH SMA. The general reason is that our set-up, where the single approaches the BBH from the Hill sphere, is so constrained that we are seeing real features in the underlying phase-space or topology of the three-body system (Samsing et al. 2018a). Before studying the topology of the scatterings, we first explore analytically if we can estimate when the scatterings go from being a near smooth power law, to greatly varying. For this, we note that the scatterings most likely have some underlying periodic features, as we keep the initial distance between the binary and single fixed (⁠|$R_{\rm H}$| does not change with the BBH SMA) while changing the SMA of the target BBH. This set-up gives rise to periodic interactions, where the single falls from the Hill sphere to interact with the BBH after the BBH has undergone n half-revolutions as the SMA decreases. In other words, we might expect periodic or repeating features when the following condition is met for integer n:

(28)

where |$T_{\rm int}$| is the time it takes the single to reach the BBH (interaction time), and |$T_{\rm BBH}$| is the orbital time of the BBH. Starting with |$T_{\rm int}$|⁠, the natural time-scale associated with this interaction time, scales with the parameters of the problem as

(29)

where for this analysis we assume the BHs to have equal mass, m. Correspondingly, the time-scale associated with the BBH is simply

(30)

By using equation (28), we can now isolate for the BBH SMA that corresponds to integer values of n,

(31)

The values of |$a_n$| for integer values of n is plotted over the probability distributions shown in Fig. 7 with black vertical solid lines. The vertical lines line up well across the different combinations of ICs, with the fluctuating and repeating features in the probability curves. Furthermore, because of the scaling from equation (31), the BBH takes more and more turns, and the spacing between the lines becomes smaller and smaller, with decreasing SMA, which explains why in this limit the probability curves have less features and more closely follows a simple power law. The more the BBH turns and the higher line density for lower binary SMA imply that the system increasingly loses information about the exact ICs, in which case it approaches the more chaotic limit that, for example, was shown in S22 to follow a smoother curve. Correspondingly, when the SMA increases towards |$R_{\rm H}$|⁠, the impact from the exact ICs is important and the underlying phase-space features start to show up as features in the probability curves. These observations and the theoretical model from Section 2.2 are all complementary, and clearly indicate that, at least for this case, the largest visible effect in the probability when introducing the SMBH is not ‘tidal breakup’, but instead fluctuations from the underlying phase-space. However, it might be that the real astrophysical environment gives rise to more chaotic ICs, e.g. if the single BH does not slowly approach the BBH from the edge of the Hill sphere, in which case the effect from tidal breakup could be a more dominant effect. Continuing to consider results from the scatterings shown in Fig. 7, we find that the critical n for which the system can be considered chaotic is for |$n \sim 10$|⁠.

Probability of mergers at $e\gt 0.1$ and $f_{\rm gw} \gt 10 \, \rm Hz$ as a function of initial inner binary SMA (similar to Fig. 5). The left column displays results for Set-up I for a 15 ${\rm M}_{\odot }$ (top) and 100 ${\rm M}_{\odot }$ (bottom) equal mass triple system, while the right column corresponds to results with Set-up II for the same BH masses. As in Fig. 5, the binary SMA ranges from values of 0.1 to 6 au both CR (red) and CT (blue). The black vertical lines represent solutions to equation (31), where n equals an integer number of half-orbits. As the size of the SMA approaches $R_{\rm H}$, the fluctuating regions become more spread out and the trend of the probability becomes less linear and more fluctuating, as described in Section 4.3.
Figure 7.

Probability of mergers at |$e\gt 0.1$| and |$f_{\rm gw} \gt 10 \, \rm Hz$| as a function of initial inner binary SMA (similar to Fig. 5). The left column displays results for Set-up I for a 15 |${\rm M}_{\odot }$| (top) and 100 |${\rm M}_{\odot }$| (bottom) equal mass triple system, while the right column corresponds to results with Set-up II for the same BH masses. As in Fig. 5, the binary SMA ranges from values of 0.1 to 6 au both CR (red) and CT (blue). The black vertical lines represent solutions to equation (31), where n equals an integer number of half-orbits. As the size of the SMA approaches |$R_{\rm H}$|⁠, the fluctuating regions become more spread out and the trend of the probability becomes less linear and more fluctuating, as described in Section 4.3.

We continue by exploring the scattering phase-space, or topology, in greater detail to further understand the probability features and where and how the eccentric BBH mergers are forming.

4.4 Topology of end-states

As we find very similar features for the phase space distribution of scatterings for Set-up I and Set-up II, we focus on the comparison of the CR and CT cases in relation to the probability features, rather than on the comparison of the two configurations. We therefore in this section, only show examples from Set-up I. Fig. 8 shows the phase space distribution of outcomes (CR for the left panel, and CT for the right panel), in terms of the initial binary SMA |$a_0$| and phase angle |$\phi _0$|⁠, as defined in Fig. 3. Similarly to Samsing & Ilan (2018), the figure is colour-coded in terms of final end-state:

  • BS[ij]: Outcome where the BBH and single BH interact and both escape from the Hill sphere without merging before. The possible combinations can be 12 (yellow), 13 (blue), and 23 (pink), where as shown in Fig. 3 the initial binary components are labelled by ‘1’ and ‘2’, and the incoming single by ‘3’.

  • GW (black): The outcome of the scattering is a three-body merger, where two of three objects undergo a GW inspiral that leads to a GW merger while all three objects are still within the Hill sphere.

Phase-space distribution of a selected region for Set-up I in the CR (left) and CT (right) case. We fix the impact parameter of the incoming single at the Hill radius and vary the inner BBH SMA and phase angle. This particular section of the phase space is between SMAs of 1 and 5 au. The plot is colour-coded in terms of final outcome as explained in Section 4.4: Escaper from Hill sphere 12 (yellow), 13 (blue), and 23 (pink). The black region corresponds to a three-body GW merger where two of the three objects merge inside the Hill sphere.
Figure 8.

Phase-space distribution of a selected region for Set-up I in the CR (left) and CT (right) case. We fix the impact parameter of the incoming single at the Hill radius and vary the inner BBH SMA and phase angle. This particular section of the phase space is between SMAs of 1 and 5 au. The plot is colour-coded in terms of final outcome as explained in Section 4.4: Escaper from Hill sphere 12 (yellow), 13 (blue), and 23 (pink). The black region corresponds to a three-body GW merger where two of the three objects merge inside the Hill sphere.

We only show a portion of the overall analysed parameter space to look more closely at what happens when the SMA approaches the size of |$R_{\rm H}$|⁠. From the phase space distribution, a few elements stand out. In terms of |$\phi _0$|⁠, for both the CR and CT cases the phase-space structure has a trivial periodicity, where the pattern repeats but in opposite colours after |$180^{\circ }$|⁠. In terms of the structure along the SMA, we see that as we approach the |$R_{\rm H}$| limit, the regions widen and for some angular intervals, e.g. in the CR case (left), there is really only one possible outcome. For example, at |$200^{\circ }$| in the CR case, the system can only result in outcome ‘BS13’ as the binary has no time to undergo significant orbital evolution before it encounters the single, which in turn turns into this restricted outcome space. At lower SMA, the binary can for just small changes in SMA undergo several revolutions, before the single arrives and initiates the interaction (see equation 31), which give rise to the seen band-like structure in both the CR and the CT case. Another interesting behaviour that we further analyse in Fig. 8, is that in the CR case (left) there appears to be a significant higher fraction of non-chaotic regions than chaotic ones compared to the CT case (right). Intuitively, this is to be expected from the set-up, as the single, which arrives from the Hill sphere, generally encounters the binary from ‘below’ (since the binary in this set-up is close to the SMBH than the single), as seen in Fig. 3, and therefore initially orbits the binary in the same direction as the CT set-up. The objects therefore encounter each other with a lower relative velocity, which generally leads to a temporary dynamical capture of the single. This effect starts the resonating interaction, which corresponds to the seed to enter a chaotic state. This is in contrast to the CR case, where the single encounters one of the binary objects with a relatively high velocity, which is more likely to result in a prompt exchange, which might bring the objects out of the Hill sphere after the first encounter. What is not expected is that the resonating interaction and corresponding chaotic regions are so rare in the CR case.

Furthermore, the phase space figures provide important hints as to how and where the eccentric BBH mergers form, which is also illustrated in Fig. 9. Here, we highlight mergers formed in the chaotic (left side) and in the prompt regions (right side). The interactions taking place in the resonating and chaotic regions are much more prone to additional effects such as gas, if present, and are therefore likely to be greatly affected, as they have much more time to lose energy through drag over the duration of the resonating state. The opposite can be stated regarding prompt mergers, although they could turn into resonating interactions if enough drag is acting on them before they leave the Hill sphere.

Close-up view of parameter space (Set-up I) with values of initial SMA ranging between 1 and 2 au, along with four examples of possible interaction scenarios depending on the location in the parameter space. The phase-space appears as divided into bands, alternating between chaotic and prompt regions. The number of interactions for the three-body encounters in the chaotic region is significantly higher than in the prompt regions.
Figure 9.

Close-up view of parameter space (Set-up I) with values of initial SMA ranging between 1 and 2 au, along with four examples of possible interaction scenarios depending on the location in the parameter space. The phase-space appears as divided into bands, alternating between chaotic and prompt regions. The number of interactions for the three-body encounters in the chaotic region is significantly higher than in the prompt regions.

Finally, a more direct comparison between the outcome of probabilities (top) and phase-space (bottom) is shown in Fig. 10 for a smaller portion of the SMA range (1–5 au). The vertical dashed lines indicate the corresponding region of the probability in the phase space. Fig. 10 highlights how the band-like configuration observed in the phase space is connected to the merger probability. As highlighted in Section 4.3, for larger values of |$a_0/R_{\rm H}$|⁠, the encounters become more controlled as the single BH is starting off closer to the binary. As a result, the regions of orbital resonance shown in Fig. 7 become more spread-out causing large fluctuations in between these regions. Consequently, as the system becomes more sensitive to the ICs, the outcome of the scattering (prompt merger or short interaction and escape) strongly depends on the combination of |$\phi _0$| and |$a_0$|⁠.

Comparison of probability for an eccentric inspiral (top) and the phase space distribution (bottom) for a 15 ${\rm M}_{\odot }$ equal mass triple system in Set-up I. The purple vertical dashed lines in both figures are plotted as an indicator of which region of the parameter space the probability corresponds to, showing for which SMAs the probability hits a low or high number of inspirals and how the fluctuation becomes more apparent as the value of $a_0/R_{\rm H}$ increases.
Figure 10.

Comparison of probability for an eccentric inspiral (top) and the phase space distribution (bottom) for a 15 |${\rm M}_{\odot }$| equal mass triple system in Set-up I. The purple vertical dashed lines in both figures are plotted as an indicator of which region of the parameter space the probability corresponds to, showing for which SMAs the probability hits a low or high number of inspirals and how the fluctuation becomes more apparent as the value of |$a_0/R_{\rm H}$| increases.

4.5 Out-of-plane interactions

Until now we have focused on properties and outcomes of the co-planar set-up, where all four objects (SMBH, binary + single) are interacting in the same plane, as this configuration gives rise to the most distinct outcomes relative to the well studied isotropic case (Samsing et al. 2014; Samsing 2018). However, there are different factors that can result in a non-fully 2D distribution of BHs embedded in the disc. When BHs in the nuclear star cluster are captured by the disc due to drag force processes they do not necessarily get aligned in fully co-planar orbits (e.g. Fabj et al. 2020; Nasim et al. 2023; Wang et al. 2024), dynamical heating will occur if a larger population of the disc is present as discussed in Stone et al. (2017), and the orbital evolution of objects inside the disc is believed to be subject to, for example, gas turbulence (Whitehead et al. 2024). All these processes ultimately result in mutual interactions that are not co-planar, which we refer to as out-of-plane interactions.

We quantify the out-of-plane interaction by the inclination angle, i, which refers to the angle between the incoming single relative to the orbital BBH and disc plane. As shown in the top part of Fig. 11, we place the single incoming BH at the Hill sphere at angle i above the disc plane, with a velocity vector that is equivalent to the co-planar set-up, i.e. we consider the situation where the single BH initially has zero velocity in the direction perpendicular to the disc. A similar set-up was explored in S22, but in their case the single BH was assumed to be coming from infinity in a plane with inclination i relative to the BBH orbital plane, in which case the single BH passes through the Hill sphere approximately with the same angle i, but with a slightly different angular momentum. In S22, it was shown that the outcome probability of eccentric mergers strongly depends on i, with nearly one order-of-magnitude decrease when the angle has opened to |$1^{\circ }$| relative to the |$0^{\circ }$| co-planar case. This strong dependence on the orbital plane inclination angle naturally questions how the results actually depend on the assumption that the single BH approached from infinity as normally assumed in scattering systems. It could be that if the object is initiated at the Hill sphere, and not propagated to the Hill sphere from far away, the dependence on i is much weaker, as what really matters is how much angular momentum the single BH brings. With this motivation, we now consider the results shown in Fig. 11, which shows the probability for eccentric mergers (⁠|$e\gt 0.1$| with a peak frequency |$f_{\rm gw} \ge 10$| Hz), resulting from our out-of-plane Set-up I where the single starts at the Hill sphere at an angle i above the plane. In Fig. 11, the BH masses are equal and set to |$15 \, {\rm M}_{\odot }$|⁠. We show in total four different setups characterized by two different SMAs (1.7 au, cross symbols and 1 au, circles) and for CR (purple) and CT (orange) binaries. As seen, the general trend is consistent with S22, i.e. at angle |$10^{-2}$| the probabilities converge to the value found in the co-planar case shown in Fig. 7, with a rapid decrease by almost 2–3 orders of magnitude once i passes |$1^{\circ }-10^{\circ }$|⁠. Even in this set-up, the probability for eccentric mergers therefore still seems to be very sensitive to the inclination angle, which adds further motivation for studying how the gas has an impact on the alignment of the objects in the disc. This is particularly important when considering the binaries surviving the interaction, as they are scattered both out of the disc plane and have their orbital momentum tilted (S22). Many of these binaries are likely to undergo subsequent interactions inside the disc, but how aligned they are in the disc at the time of interaction depends on the alignment time-scales relative to the interaction time-scale, as well as the properties of the binary after interaction. The subsequent question is naturally what the expected out-of-plane angle could be in AGN disc like environments. Assuming the binary orbital angular moment is perpendicular to the disc angular momentum, a simple estimate can be made in the case the Hill sphere is larger than the disc height, in which case the maximum opening angle is around the AGN disc height divided by the Hill radius, i.e. |$\mathrm{ max}(i) \approx h/R_H \approx (h/R) \times \left(M/m\right)^{1/3}$|⁠. The maximum opening angle can span a range of different values, depending on both aspect ratio and mass ratio. For example, for a system described by |$m=10^{1} {\rm M}_{\odot }$|⁠, |$M=10^{7} {\rm M}_{\odot }$|⁠, then |$\mathrm{ max}(i) \approx h/R_{\rm H} \times 10^{2}$|⁠, i.e. for thin-disc models with |$h/R \sim 10^{-3}$| like Shakura & Sunyaev (1973), the characteristic angle could very well be |$\lt 1 ^{\circ }$|⁠. Reading from the bottom panel in Fig. 11, this would still mean that eccentric mergers are likely to form in overabundance in disc environments compared to isotropic environments. However, for more realistic disc models, this simple estimate could vary throughout the disc, reaching values of |$h/R \sim 0.1$| (e.g. Sirko & Goodman 2003; Thompson, Quataert & Murray 2005). This structure would make some regions more prone to forming eccentric mergers than others.

Top figure: Illustration of the set-up we use for the out-of-plane interactions we explore in Section 4.5, which is similar to the set-up shown in Fig. 3, but now also with inclination angle i. An example of an out-of-plane scattering is also shown in the large blue circle indicating the Hill sphere. We here initiate the single object at the Hill sphere at an angle i above the plane with the same velocity vector that we used in the co-planar case, i.e. we assume that the velocity perpendicular to the plane is equal to zero when the simulation starts. This set-up is used to explore and quantify how the vertical oscillations and relative orbital fluctuations there are in the disc impact the formation of eccentric mergers. Bottom figure: Probability of eccentric mergers as a function of i for the out-of-plane scatterings in both CR (purple) and CT (orange) configurations for two SMAs (1.7 and 1 au). The dashed lines represent the 2D (disc) and 3D (isotropic) limits for eccentric mergers.
Figure 11.

Top figure: Illustration of the set-up we use for the out-of-plane interactions we explore in Section 4.5, which is similar to the set-up shown in Fig. 3, but now also with inclination angle i. An example of an out-of-plane scattering is also shown in the large blue circle indicating the Hill sphere. We here initiate the single object at the Hill sphere at an angle i above the plane with the same velocity vector that we used in the co-planar case, i.e. we assume that the velocity perpendicular to the plane is equal to zero when the simulation starts. This set-up is used to explore and quantify how the vertical oscillations and relative orbital fluctuations there are in the disc impact the formation of eccentric mergers. Bottom figure: Probability of eccentric mergers as a function of i for the out-of-plane scatterings in both CR (purple) and CT (orange) configurations for two SMAs (1.7 and 1 au). The dashed lines represent the 2D (disc) and 3D (isotropic) limits for eccentric mergers.

5 CONCLUSIONS AND CAVEATS

In this paper, we have explored the outcomes from binary–single interactions between migrating BBHs and single BHs in AGN disc-like environments in the presence of an SMBH. In particular, we have studied how the fraction of eccentric mergers depends on the BBH SMA relative to the Hill sphere created by the tidal field of the SMBH, the BBH direction of rotation relative to the disc, as well as the relative inclination angle between the BBH and the single BH as they approach each other through their Hill sphere. We have further presented the first analytical solution to how a tidal field in general terms has an impact on the evolution of resonating binary–single scatterings, which are one of the main pathways for producing eccentric mergers. In the following list we summarize our main conclusions and discuss the possible implications of adding the effect of a gaseous disc.

  • Theory of three-body interactions in a tidal field. When binary–single objects interact, they often undergo a large number of resonant states that are characterized by long excursions of the single relative to the binary. Each of these IMS states has been shown to be one of the main pathways to form eccentric mergers in systems where such interactions can take place. In Section 2.2, we have extended the theory of how such eccentric three-body mergers form in resonant interactions, by including the effect from a tidal field characterized by a Hill sphere boundary. By the use of statistical theory for how the energy and angular momentum is distributed in chaotic encounters, based on works by Valtonen & Karttunen (2006) and Stone & Leigh (2019), we showed that the probability for three-body mergers, i.e. mergers forming during the interactions, are generally reduced by a factor of |$(1-2 a_0 /R_{\rm H})^2$|⁠, where |$a_0$| is the initial BBH SMA and |$R_{\rm H}$| is the Hill sphere (see equation 14). While true chaotic interactions are difficult to achieve in environments with a strong nearby perturber, we identify this as the first indication of how a nearby SMBH might impact the fraction of eccentric mergers. Fig. 5 shows this correction factor applied on the results from S22.

  • PN simulations of eccentric mergers under the influence of the SMBH. We performed a systematic study of PN scatterings between a BBH and a single incoming BH taking place near the SMBH. In this set-up, we especially explored the formation probability of eccentric BBH mergers (e > 0.1 at |$f_{\rm gw}$| > 10 Hz), and how this depends on the SMA of the initial BBH |$(a_0)$| relative to the effective Hill sphere of the SMBH (⁠|$R_{\rm H}$|⁠), its orbital angular momentum relative to its direction around the SMBH (co- and counter-rotating), the location of the single relative to the BBH (Set-up I and Set-up II, see Fig. 3), as well as the initial inclination between the BBH and single BH. In all scatterings, we initiated the single BH at the Hill sphere, to emulate the likely process that takes place when a BBH and a single BH meet in an AGN-disc environment through disc migration. S22 showed that three-body scatterings in disc-like configurations with no SMBH included, produce a high fraction of eccentric mergers. With more realistic set-ups we here find that eccentric mergers are indeed also found to be produced in large numbers even when the interactions take place close to the SMBH. The reason is partly that many of the near co-planar interactions produce mergers that are formed through prompt interactions (see Fig. 9), that are not greatly influenced by the SMBH upon formation. Finally, as shown in Fig. 6, many of the three-body mergers do actually start with a GW peak frequency that can be far above 10 Hz, hinting that a significant fraction might show up as bursts or prompt mergers, rather than classical eccentric inspirals, where the system evolves from eccentric to circular over hundreds of cycles. This result is in agreement with the analysis of Trani et al. (2024) and Rom et al. (2023), where for the latter they explore the effects of tidal forces coming from the SMBH on binary capture in AGNs. These findings further motivate searches for burst-like GW signals.

  • Merger types and phase-space distributions. We studied the distribution of outcomes as a function of SMA, and orbital phase angle |$\phi _0$| for different set-ups, and found a clear structure that propagates to our main results of the probability for eccentric mergers as shown in Fig. 10. Generally, we found that when the SMA of the BBH is much smaller than the Hill sphere, the scatterings nearly lose information about the ICs, and the results therefore approach the one found by S22 that by construction explored outcomes from a much less constrained set-up (the single coming in isotropically in a plane, where we e.g. start always at the same point at the Hill sphere). On the other hand, when the SMA increases and starts to approach the size of |$R_{\rm H}$|⁠, the outcomes strongly depend on the exact angular phase relative to the SMA, with an approximate periodicity arising from considering integer values of the interaction time relative to the BBH orbital time (see equation 31 and Fig. 7). Despite these strong dependencies when approaching the Hill sphere, we find that the probability for eccentric mergers in this region as well approximately matches the one found by S22, but with significant fluctuations coming from the underlying phase-space structure. In addition, both CT and CR for Set-up I and Set-up II all showed similar outcome distribution.

  • Importance of orbital inclinations. No interactions are expected to be perfectly co-planar, studying how our results depend on the initial inclination angle between the BBH and single when the single enters the Hill sphere (angle i as seen in Fig. 11), is therefore of major importance. In Section 4.5, we explored the fraction of eccentric BBH mergers as a function of the initial inclination angle i, and found that across the different constrained set-ups, the probability does falls off quickly with i as was also pointed out in S22. All AGN-disc environments are different from each other, but we argue that for thin disc models, in regions where the maximum opening angle is |$i \sim 1^{\circ }$|⁠, a relatively high fraction of eccentric mergers are likely to be produced.

  • Potential impact of gas in the accretion disc. As discussed, we do not include gaseous effects in our analysis but acknowledge their potential impact on the presented results and on the overall fraction of eccentric mergers. A relevant question to be explored is if the additional gas drag force term keeps the objects contained in the Hill sphere, ultimately leading to more chaotic interactions in the co-rotating and counter-rotating case as a contrast to the high fraction of prompt interaction regions we observe from, for example, Figs 8 and 9. In addition, a full dynamical model of migrating BHs inside AGN-disc-like environments is necessary to state what relevant values of |$R_{\rm H}$| and SMA are to be expected, and therefore how many interactions can be considered chaotic or not. As far as the nature and outlook of mergers is concerned, the scattering seen in the bottom panel of Fig. 4 might in comparison not be much affected by gas, as the binary undergoing merger is formed promptly. Therefore, we speculate that gas might not have a great impact on the formation of such prompt mergers, which could hint that our results might also hold when gaseous effects are included. There are, however, several reasons to believe that our considerations and assumptions from above need to be adjusted, as the relevant setup that we consider involves binaries that encounter singles as a result of migration through an AGN-disc (e.g. Secunda et al. 2019, 2020; Tagawa et al. 2020). Lastly, whether gas friction is expected to damp or pump the eccentricity of the binary is a non-trivial answer, as its evolution depends on a number of factors of the parameter space, such as mass ratio or whether binary is in prograde or retrograde motion with respect to the gas (e.g. Li & Lai 2022, 2023, 2024; Dittmann, Dempsey & Li 2024); Valli et al. 2024).

ACKNOWLEDGEMENTS

The authors thank the anonymous referee for the useful suggestions and comments. GF thanks the GW-Astro group at NBIA, in particular Martin Pessah and Daniel J. D’Orazio, for insightful discussions and constant support, as well as the TIDY-NYC group for interesting discussions. The Tycho supercomputer hosted at the SCIENCE HPC Center at the University of Copenhagen was used for performing the simulations presented in the paper. This work was funded by the Villum Fonden grant no. 29466 and European Research Council starting grant no. 101043143 — BlackHoleMergs.

DATA AVAILABILITY

The data underlying this article will be shared on request to the corresponding author.

Footnotes

1

We note that in Trani et al. (2024), a large fraction of eccentric mergers were also found in their set-up with an SMBH, but the ICs and underlying physical set-up (inclusion of stellar-disc velocity dispersion) are different from ours, which makes it difficult to compare side-by-side.

REFERENCES

Abbott
B. P.
et al. ,
2016a
,
Phys. Rev. X
,
6
,
041015

Abbott
B. P.
et al. ,
2016b
,
Phys. Rev. Lett.
,
116
,
061102

Abbott
B. P.
et al. ,
2016c
,
Phys. Rev. Lett.
,
116
,
241103

Abbott
B. P.
et al. ,
2017a
,
Phys. Rev. Lett.
,
118
,
221101

Abbott
B. P.
et al. ,
2017b
,
Phys. Rev. Lett.
,
119
,
141101

Abbott
B. P.
et al. ,
2017c
,
Phys. Rev. Lett.
,
119
,
161101

Abbott
B. P.
et al. ,
2020a
,
Class. Quantum Gravity
,
37
,
055002

Abbott
R.
et al. ,
2020b
,
Phys. Rev. Lett.
,
125
,
101102

Abbott
R.
et al. ,
2020c
,
ApJ
,
900
,
L13

Abbott
R.
et al. ,
2023
,
Phys. Rev. X
,
13
,
041039

Antonini
F.
,
Rasio
F. A.
,
2016
,
ApJ
,
831
,
187

Artale
M. C.
,
Mapelli
M.
,
Giacobbo
N.
,
Sabha
N. B.
,
Spera
M.
,
Santoliquido
F.
,
Bressan
A.
,
2019
,
MNRAS
,
487
,
1675

Askar
A.
,
Szkudlarek
M.
,
Gondek-Rosińska
D.
,
Giersz
M.
,
Bulik
T.
,
2017
,
MNRAS
,
464
,
L36

Bae
Y.-B.
,
Kim
C.
,
Lee
H. M.
,
2014
,
MNRAS
,
440
,
2714

Banerjee
S.
,
Baumgardt
H.
,
Kroupa
P.
,
2010
,
MNRAS
,
402
,
371

Bartos
I.
,
Kocsis
B.
,
Haiman
Z.
,
Márka
S.
,
2017
,
ApJ
,
835
,
165

Belczynski
K.
,
Holz
D. E.
,
Bulik
T.
,
O’Shaughnessy
R.
,
2016a
,
Nature
,
534
,
512

Belczynski
K.
,
Repetto
S.
,
Holz
D. E.
,
O’Shaughnessy
R.
,
Bulik
T.
,
Berti
E.
,
Fryer
C.
,
Dominik
M.
,
2016b
,
ApJ
,
819
,
108

Bird
S.
,
Cholis
I.
,
Muñoz
J. B.
,
Ali-Haïmoud
Y.
,
Kamionkowski
M.
,
Kovetz
E. D.
,
Raccanelli
A.
,
Riess
A. G.
,
2016
,
Phys. Rev. Lett.
,
116
,
201301

Blanchet
L.
,
2006
,
Living Rev. Relativ.
,
9
,
4

Blanchet
L.
,
2014
,
Living Rev. Relativ.
,
17
,
2

Boekholt
T. C. N.
,
Rowan
C.
,
Kocsis
B.
,
2023
,
MNRAS
,
518
,
5653

Calcino
J.
,
Dempsey
A. M.
,
Dittmann
A. J.
,
Li
H.
,
2024
,
ApJ
,
970
,
107

Cantiello
M.
,
Jermyn
A. S.
,
Lin
D. N. C.
,
2021
,
ApJ
,
910
,
94

Carr
B.
,
Kühnel
F.
,
Sandstad
M.
,
2016
,
Phys. Rev. D
,
94
,
083504

Chen
X.
,
Amaro-Seoane
P.
,
2017
,
ApJ
,
842
,
L2

Cholis
I.
,
Kovetz
E. D.
,
Ali-Haïmoud
Y.
,
Bird
S.
,
Kamionkowski
M.
,
Muñoz
J. B.
,
Raccanelli
A.
,
2016
,
Phys. Rev. D
,
94
,
084013

D’Orazio
D. J.
,
Loeb
A.
,
2018
,
Phys. Rev. D
,
97
,
083008

D’Orazio
D. J.
,
Loeb
A.
,
2020
,
Phys. Rev. D
,
101
,
083031

D’Orazio
D. J.
,
Samsing
J.
,
2018
,
MNRAS
,
481
,
4775

DeLaurentiis
S.
,
Epstein-Martin
M.
,
Haiman
Z.
,
2023
,
MNRAS
,
523
,
1126

Dempsey
A. M.
,
Li
H.
,
Mishra
B.
,
Li
S.
,
2022
,
ApJ
,
940
,
155

Dittmann
A. J.
,
Cantiello
M.
,
Jermyn
A. S.
,
2021
,
ApJ
,
916
,
48

Dittmann
A. J.
,
Jermyn
A. S.
,
Cantiello
M.
,
2023
,
ApJ
,
946
,
56

Dittmann
A. J.
,
Dempsey
A. M.
,
Li
H.
,
2024
,
ApJ
,
964
,
61

Dominik
M.
,
Belczynski
K.
,
Fryer
C.
,
Holz
D. E.
,
Berti
E.
,
Bulik
T.
,
Mandel
I.
,
O’Shaughnessy
R.
,
2012
,
ApJ
,
759
,
52

Dominik
M.
,
Belczynski
K.
,
Fryer
C.
,
Holz
D. E.
,
Berti
E.
,
Bulik
T.
,
Mandel
I.
,
O’Shaughnessy
R.
,
2013
,
ApJ
,
779
,
72

Dominik
M.
et al. ,
2015
,
ApJ
,
806
,
263

Fabj
G.
,
Nasim
S. S.
,
Caban
F.
,
Ford
K. E. S.
,
McKernan
B.
,
Bellovary
J. M.
,
2020
,
MNRAS
,
499
,
2608

Fishbach
M.
,
Holz
D. E.
,
2020
,
ApJ
,
904
,
L26

Fishbach
M.
,
Farr
W. M.
,
Holz
D. E.
,
2020
,
ApJ
,
891
,
L31

Gayathri
V.
et al. ,
2022
,
Nat. Astron.
,
6
,
344

Generozov
A.
,
Perets
H. B.
,
2023
,
MNRAS
,
522
,
1763

Gilbaum
S.
,
Stone
N. C.
,
2022
,
ApJ
,
928
,
191

Ginat
Y. B.
,
Perets
H. B.
,
2021
,
MNRAS
,
508
,
190

Graham
M. J.
et al. ,
2020
,
Phys. Rev. Lett.
,
124
,
251102

Graham
M. J.
et al. ,
2023
,
ApJ
,
942
,
99

Gültekin
K.
,
Miller
M. C.
,
Hamilton
D. P.
,
2006
,
ApJ
,
640
,
156

Hamers
A. S.
,
Bar-Or
B.
,
Petrovich
C.
,
Antonini
F.
,
2018
,
ApJ
,
865
,
2

Hamers
A. S.
,
Fragione
G.
,
Neunteufel
P.
,
Kocsis
B.
,
2021
,
MNRAS
,
506
,
5345

Heggie
D. C.
,
1975
,
MNRAS
,
173
,
729

Hendriks
K.
,
Zwick
L.
,
Samsing
J.
,
2024
,
preprint
()

Hoang
B.-M.
,
Naoz
S.
,
Kocsis
B.
,
Rasio
F. A.
,
Dosopoulou
F.
,
2018
,
ApJ
,
856
,
140

Hong
J.
,
Lee
H. M.
,
2015
,
MNRAS
,
448
,
754

Hotokezaka
K.
,
Piran
T.
,
2017
,
ApJ
,
842
,
111

Inayoshi
K.
,
Tamanini
N.
,
Caprini
C.
,
Haiman
Z.
,
2017
,
Phys. Rev. D
,
96
,
063014

Janiuk
A.
,
Bejger
M.
,
Charzyński
S.
,
Sukova
P.
,
2017
,
New Astron.
,
51
,
7

Jermyn
A. S.
,
Dittmann
A. J.
,
Cantiello
M.
,
Perna
R.
,
2021
,
ApJ
,
914
,
105

Jermyn
A. S.
,
Dittmann
A. J.
,
McKernan
B.
,
Ford
K. E. S.
,
Cantiello
M.
,
2022
,
ApJ
,
929
,
133

Kalogera
V.
,
2000
,
ApJ
,
541
,
319

Kinugawa
T.
,
Inayoshi
K.
,
Hotokezaka
K.
,
Nakauchi
D.
,
Nakamura
T.
,
2014
,
MNRAS
,
442
,
2963

Klimenko
S.
et al. ,
2016
,
Phys. Rev. D
,
93
,
042004

Kremer
K.
et al. ,
2019a
,
Phys. Rev. D
,
99
,
063003

Kremer
K.
,
Lu
W.
,
Rodriguez
C. L.
,
Lachat
M.
,
Rasio
F. A.
,
2019b
,
ApJ
,
881
,
75

Lee
W. H.
,
Ramirez-Ruiz
E.
,
van de Ven
G.
,
2010
,
ApJ
,
720
,
953

Leigh
N. W. C.
et al. ,
2018
,
MNRAS
,
474
,
5672

Li
R.
,
Lai
D.
,
2022
,
MNRAS
,
517
,
1602

Li
R.
,
Lai
D.
,
2023
,
MNRAS
,
522
,
1881

Li
R.
,
Lai
D.
,
2024
,
MNRAS
,
529
,
1

Li
J.
,
Dempsey
A. M.
,
Li
H.
,
Lai
D.
,
Li
S.
,
2023
,
ApJ
,
944
,
L42

Loeb
A.
,
2016
,
ApJ
,
819
,
L21

Lopez
Martin J.
,
Batta
A.
,
Ramirez-Ruiz
E.
,
Martinez
I.
,
Samsing
J.
,
2019
,
ApJ
,
877
,
56

MacLeod
M.
,
Lin
D. N. C.
,
2020
,
ApJ
,
889
,
94

McKernan
B.
,
Ford
K. E. S.
,
Lyra
W.
,
Perets
H. B.
,
2012
,
MNRAS
,
425
,
460

McKernan
B.
et al. ,
2018
,
ApJ
,
866
,
66

Monaghan
J. J.
,
1976
,
MNRAS
,
177
,
583

Murguia-Berthier
A.
,
MacLeod
M.
,
Ramirez-Ruiz
E.
,
Antoni
A.
,
Macias
P.
,
2017
,
ApJ
,
845
,
173

Nasim
S. S.
et al. ,
2023
,
MNRAS
,
522
,
5393

Nitz
A. H.
,
Capano
C. D.
,
2021
,
ApJ
,
907
,
L9

O’Leary
R. M.
,
Kocsis
B.
,
Loeb
A.
,
2009
,
MNRAS
,
395
,
2127

O’Neill
D.
,
D’Orazio
D. J.
,
Samsing
J.
,
Pessah
M. E.
,
2024
,
ApJ
,
974
,
216

Park
D.
,
Kim
C.
,
Lee
H. M.
,
Bae
Y.-B.
,
Belczynski
K.
,
2017
,
MNRAS
,
469
,
4665

Peters
P. C.
,
1964
,
Phys. Rev.
,
136
,
1224

Piran
Z.
,
Piran
T.
,
2020
,
ApJ
,
892
,
64

Portegies Zwart
S. F.
,
McMillan
S. L. W.
,
2000
,
ApJ
,
528
,
L17

Rodriguez
C. L.
,
Antonini
F.
,
2018
,
ApJ
,
863
,
7

Rodriguez
C. L.
,
Morscher
M.
,
Pattabiraman
B.
,
Chatterjee
S.
,
Haster
C.-J.
,
Rasio
F. A.
,
2015
,
Phys. Rev. Lett.
,
115
,
051101

Rodriguez
C. L.
,
Chatterjee
S.
,
Rasio
F. A.
,
2016a
,
Phys. Rev. D
,
93
,
084029

Rodriguez
C. L.
,
Haster
C.-J.
,
Chatterjee
S.
,
Kalogera
V.
,
Rasio
F. A.
,
2016b
,
ApJ
,
824
,
L8

Rodriguez
C. L.
,
Zevin
M.
,
Pankow
C.
,
Kalogera
V.
,
Rasio
F. A.
,
2016c
,
ApJ
,
832
,
L2

Rodriguez
C. L.
,
Amaro-Seoane
P.
,
Chatterjee
S.
,
Kremer
K.
,
Rasio
F. A.
,
Samsing
J.
,
Ye
C. S.
,
Zevin
M.
,
2018
,
Phys. Rev. D
,
98
,
123005

Rom
B.
,
Sari
R.
,
Lai
D.
,
2024
,
ApJ
,
964
,
43

Romero-Shaw
I.
,
Lasky
P. D.
,
Thrane
E.
,
Calderón Bustillo
J.
,
2020
,
ApJ
,
903
,
L5

Romero-Shaw
I.
,
Lasky
P. D.
,
Thrane
E.
,
2022
,
ApJ
,
940
,
171

Rowan
C.
,
Boekholt
T.
,
Kocsis
B.
,
Haiman
Z.
,
2023
,
MNRAS
,
524
,
2770

Rowan
C.
,
Whitehead
H.
,
Boekholt
T.
,
Kocsis
B.
,
Haiman
Z.
,
2024
,
MNRAS
,
527
,
10448

Rozner
M.
,
Perets
H. B.
,
2022
,
ApJ
,
931
,
149

Rozner
M.
,
Generozov
A.
,
Perets
H. B.
,
2023
,
MNRAS
,
521
,
866

Samsing
J.
,
2018
,
Phys. Rev. D
,
97
,
103014

Samsing
J.
,
D’Orazio
D. J.
,
2018
,
MNRAS
,
481
,
5445

Samsing
J.
,
Ilan
T.
,
2018
,
MNRAS
,
476
,
1548

Samsing
J.
,
Ramirez-Ruiz
E.
,
2017
,
ApJ
,
840
,
L14

Samsing
J.
,
MacLeod
M.
,
Ramirez-Ruiz
E.
,
2014
,
ApJ
,
784
,
71

Samsing
J.
,
MacLeod
M.
,
Ramirez-Ruiz
E.
,
2017
,
ApJ
,
846
,
36

Samsing
J.
,
MacLeod
M.
,
Ramirez-Ruiz
E.
,
2018a
,
ApJ
,
853
,
140

Samsing
J.
,
Askar
A.
,
Giersz
M.
,
2018b
,
ApJ
,
855
,
124

Samsing
J.
,
Venumadhav
T.
,
Dai
L.
,
Martinez
I.
,
Batta
A.
,
Lopez
M.
,
Ramirez-Ruiz
E.
,
Kremer
K.
,
2019a
,
Phys. Rev. D
,
100
,
043009

Samsing
J.
,
Hamers
A. S.
,
Tyles
J. G.
,
2019b
,
Phys. Rev. D
,
100
,
043010

Samsing
J.
,
D’Orazio
D. J.
,
Kremer
K.
,
Rodriguez
C. L.
,
Askar
A.
,
2020
,
Phys. Rev. D
,
101
,
123010

Samsing
J.
et al. ,
2022
,
Nature
,
603
,
237

Samsing
J.
,
Hendriks
K.
,
Zwick
L.
,
D’Orazio
D. J.
,
Liu
B.
,
2024
,
preprint
()

Sasaki
M.
,
Suyama
T.
,
Tanaka
T.
,
Yokoyama
S.
,
2016
,
Phys. Rev. Lett.
,
117
,
061101

Schrøder
S. L.
,
Batta
A.
,
Ramirez-Ruiz
E.
,
2018
,
ApJ
,
862
,
L3

Secunda
A.
,
Bellovary
J.
,
Mac Low
M.-M.
,
Ford
K. E. S.
,
McKernan
B.
,
Leigh
N. W. C.
,
Lyra
W.
,
Sándor
Z.
,
2019
,
ApJ
,
878
,
85

Secunda
A.
et al. ,
2020
,
ApJ
,
903
,
133

Secunda
A.
,
Hernandez
B.
,
Goodman
J.
,
Leigh
N. W. C.
,
McKernan
B.
,
Ford
K. E. S.
,
Adorno
J. I.
,
2021
,
ApJ
,
908
,
L27

Shakura
N. I.
,
Sunyaev
R. A.
,
1973
,
A&A
,
24
,
337

Silsbee
K.
,
Tremaine
S.
,
2017
,
ApJ
,
836
,
39

Sirko
E.
,
Goodman
J.
,
2003
,
MNRAS
,
341
,
501

Stephan
A. P.
,
Naoz
S.
,
Ghez
A. M.
,
Witzel
G.
,
Sitarski
B. N.
,
Do
T.
,
Kocsis
B.
,
2016
,
MNRAS
,
460
,
3494

Stone
N. C.
,
Leigh
N. W. C.
,
2019
,
Nature
,
576
,
406

Stone
N. C.
,
Metzger
B. D.
,
Haiman
Z.
,
2017
,
MNRAS
,
464
,
946

Tagawa
H.
,
Haiman
Z.
,
Kocsis
B.
,
2020
,
ApJ
,
898
,
25

Tagawa
H.
,
Kocsis
B.
,
Haiman
Z.
,
Bartos
I.
,
Omukai
K.
,
Samsing
J.
,
2021
,
ApJ
,
908
,
194

Tanikawa
A.
,
2013
,
MNRAS
,
435
,
1358

Thompson
T. A.
,
Quataert
E.
,
Murray
N.
,
2005
,
ApJ
,
630
,
167

Trani
A. A.
,
Spera
M.
,
Leigh
N. W. C.
,
Fujii
M. S.
,
2019
,
ApJ
,
885
,
135

Trani
A. A.
,
Quaini
S.
,
Colpi
M.
,
2024
,
A&A
,
683
,
A135

Valli
R.
et al. ,
2024
,
A&A
,
688
,
A128

Valtonen
M.
,
Karttunen
H.
,
2006
,
The Three-Body Problem
,
Cambridge University Press
,
Cambridge, UK

VanLandingham
J. H.
,
Miller
M. C.
,
Hamilton
D. P.
,
Richardson
D. C.
,
2016
,
ApJ
,
828
,
77

Venumadhav
T.
,
Zackay
B.
,
Roulet
J.
,
Dai
L.
,
Zaldarriaga
M.
,
2020
,
Phys. Rev. D
,
101
,
083030

Vynatheya
P.
,
Hamers
A. S.
,
2022
,
ApJ
,
926
,
195

Wang
Y.
,
Zhu
Z.
,
Lin
D. N. C.
,
2024
,
MNRAS
,
528
,
4958

Whitehead
H.
,
Rowan
C.
,
Boekholt
T.
,
Kocsis
B.
,
2024
,
MNRAS
,
533
,
1766

Woosley
S. E.
,
2016
,
ApJ
,
824
,
L10

Zaldarriaga
M.
,
Kushnir
D.
,
Kollmeier
J. A.
,
2018
,
MNRAS
,
473
,
4174

Zevin
M.
,
Pankow
C.
,
Rodriguez
C. L.
,
Sampson
L.
,
Chase
E.
,
Kalogera
V.
,
Rasio
F. A.
,
2017
,
ApJ
,
846
,
82

Zevin
M.
,
Samsing
J.
,
Rodriguez
C.
,
Haster
C.-J.
,
Ramirez-Ruiz
E.
,
2019
,
ApJ
,
871
,
91

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.