Topology of Black Hole Binary-Single Interactions

TL;DR

This study explores how the outcomes of binary-single black hole interactions depend on initial conditions, revealing the potential for multiple gravitational wave mergers, including scenarios observable by LIGO, especially in co-planar configurations.

Contribution

The paper introduces a comprehensive analysis of the topology of black hole interactions using extensive N-body simulations that include gravitational wave emission, highlighting the possibility of sequential GW mergers.

Findings

Co-planar interactions often lead to rapid sequential GW mergers.

Inclusion of GW emission in simulations reveals new merger pathways.

Some interactions produce events observable by LIGO within short timescales.

Abstract

We present a study on how the outcomes of co-planar binary-single black hole (BH) interactions distribute as a function of the orbital initial conditions. We refer to this distribution as the topology. Using an $N$-body code that includes BH finite sizes and gravitational wave (GW) emission in the equation-of-motion (EOM), we perform more than a million three-body scatterings to explore the topology of both the classical and the relativistic limit. We further describe how the inclusion of GW emission in the EOM naturally leads to scenarios where the binary-single system undergoes two successive GW mergers. The time span from the first to the second GW merger is shortest in co-planar interactions, and are for some configurations even short enough for LIGO to observe both events. We note that co-planar interactions could be frequent in environments such as active galactic nuclei discs.

Full text

11affiliationtext: Department of Astrophysical Sciences, Princeton University, Peyton Hall, 4 Ivy Lane, Princeton, NJ 08544, USA22affiliationtext: Einstein Fellow

1. Introduction

Dynamical interactions of black holes (BHs) in dense stellar environments, are likely to play an important role for the assembly of binary black hole (BBH) mergers observable by gravitational wave (GW) interferometers such as LIGO (e.g. Abbott et al., 2016a; Rodriguez et al., 2016; Antonini & Rasio, 2016). Major effort is especially being devoted to understand the formation and evolution of BHs in globular cluster (GC) systems (e.g. Samsing et al., 2014; Rodriguez et al., 2016; Wang et al., 2016; Askar et al., 2017; Samsing & Ramirez-Ruiz, 2017; Samsing et al., 2017), where the derived BBH merger rate interestingly seems to be comparable with that inferred from the current sample of LIGO detections (e.g. Abbott et al., 2016c, b, d, 2017). Related few-body studies are therefore often focused on deriving cross sections and rates assuming the orbital and angular momentum initial conditions (ICs) follow an isotropic distribution at infinity, which would be the case for most dense stellar systems (e.g. Hut & Bahcall, 1983; Gültekin et al., 2006; Samsing et al., 2014; Antognini & Thompson, 2016; Samsing et al., 2016; Samsing & Ramirez-Ruiz, 2017).

However, not all few-body interactions take place under isotropic conditions. Dynamical interactions of BHs in active galactic nuclei (AGN) discs provide one example (e.g. Bellovary et al., 2016; Stone et al., 2017; Bartos et al., 2017; McKernan et al., 2017), where gas torques force objects to move and interact in approximately the same orbital plane as the disc itself (e.g. Ostriker, 1983; Syer et al., 1991; Goodman & Tan, 2004; Levin, 2007; Kocsis et al., 2011; McKernan et al., 2011; Just et al., 2012). Although the number of BHs in AGN discs is currently unknown, it has interestingly been pointed out that if BHs do exist in such environments, then they are likely to pile up and undergo frequent interactions in migration-traps similar to those found in planetary discs (e.g. Horn et al., 2012; McKernan et al., 2012; Bellovary et al., 2016; McKernan et al., 2017). If the BHs have different masses, then the corresponding difference in migration time can also lead to interactions throughout the disc (e.g. McKernan et al., 2017). While this picture is still incomplete and poorly understood, it does provide at least one example of an environment where few-body interactions do not follow an isotropic distribution.

Motivated by these recent ideas, we explore in this paper the dynamics and outcomes of co-planar binary-single interactions. For this, we perform a detailed study of the binary-single phase-space, a representation we broadly refer to as the binary-single topology in analogy with Hut (1983b). Our work should partly be considered as a continuation of that by Hut (1983b), where the topology of co-planar binary-single interactions in the Newtonian point-particle limit was studied. In our work we expand on that study by including General Relativity (GR) corrections in the N-body equation-of-motion (EOM) as well as BH finite sizes. These corrections break the scale-free nature of the Newtonian point-particle binary-single problem (e.g. Hut, 1983b; Hut & Bahcall, 1983), leading to several new effects including the formation of GW captures and BH collisions (Samsing et al., 2014). Parts of our presented study describe highly idealized situations, however, we see this as a necessary first step, as the three-body problem with GR effects is largely unexplored for anisotropic distributions; a proper understanding of the BBH outcome distribution and how it relates to the binary-single ICs is crucial for connecting observations of GW events to their astrophysical origin.

We further introduce a new unexplored dynamical pathway for generating two successive GW merger events originating from the same binary-single interaction. Here the first GW merger happens during the interaction while the three BHs are still in a bound state (e.g. Samsing et al., 2014), where the second GW merger happens between the BH formed in the first merger and the remaining bound BH. Interestingly, we find the time span from the first to the second GW merger in co-planar interactions to occasionally be short enough for LIGO to observe both GW mergers. This suggests that AGN discs might be a likely birthplace for such double GW merger events.

The paper is organized as follows. In Section 2 we introduce the binary-single interaction channel, together with the ICs and outcomes that are relevant for this study. After this follows in Section 3 a detailed study of the binary-single topology for both the classical and the relativistic limit. In Section 4 we explore in greater detail how the inclusion of GW emission in the N-body EOM affects the resultant topology and range of outcomes. Results from zoom-in simulations performed to study micro-topological structures are given in Section 5. Finally, in Section 6 we present our proposed channel for generating double GW merger events through binary-single interactions. Conclusions are given in Section 7

2. Black Hole Binary-Single Interactions

In this work we study interactions between a circular binary and an incoming single, and how the corresponding outcomes distribute as a function of the ICs. Throughout the paper we assume the three objects to be identical BHs each with mass $m=10M_{\odot}$, and corresponding Schwarzschild radius $R_{\rm s}$. All interactions are performed with our $N$-body code presented in Samsing et al. (2016), which includes GR corrections using the post-Newtonian (PN) formalism (e.g. Blanchet, 2014). For the results presented in this paper, we include the 2.5PN term for describing the dynamical loss of energy and angular momentum through GW emission. The lower order energy preserving precession terms are not included for simplicity.

In this section we first summarize a few dynamical properties of binary-single interactions with GW emission and finite sizes included in the EOM. We then list the ICs relevant for our considered binary-single study together with the range of possible outcomes. The last part contains a description of what we refer to as the topology. We refer the reader to Hut (1983b) for further details.

2.1. Three-Body Energetics and Interactions

The range of possible outcomes is broadly determined by the orbital velocity of the initial binary, $v_{\rm orb}$, relative to the velocity of the incoming single at infinity, $v_{\infty}$. In the equal mass case, one finds that the initial orbital binary-single energy is negative for $v<1$ and positive for $v>1$, where $v\equiv v_{\infty}/(v_{\rm orb}\sqrt{3})$, as further described in Hut & Bahcall (1983). These two limits are usually referred to as the hard binary (HB) limit and the soft binary (SB) limit, respectively (e.g. Heggie, 1975). We describe a few characteristics of these two limits below. We will argue that only the HB limit is relevant for our study.

2.1.1 Soft-Binary Limit

In the SB limit the binary-single interaction is always prompt (Hut, 1983a), any dynamical corrections to the $N$-body EOM will therefore rarely have a chance to affect the standard Newtonian distribution of outcomes. As a result, the inclusion of GW emission, as well as collisions and tidal interactions, does not play a significant role in this limit (e.g. Samsing et al., 2014). The topology of the SB limit was studied by Hut (1983b), and will therefore not be discussed further in this paper.

2.1.2 Hard-Binary Limit

In the HB limit the binary-single system either undergoes a direct interaction (DI) or a resonant interaction (RI), depending on the exact ICs (e.g. Samsing et al., 2014). For equal mass interactions, these two interaction types appear statistically in near equal numbers. However, the RI channel has a much higher probability than the DI channel of entering a state where GR effects and finite sizes can affect the three-body dynamics, as described in the following.

A DI is characterized by a duration time that is similar to the orbital time of the initial binary, which makes it somewhat similar to a SB limit interaction. The kinematics of a DI is often such that the incoming single promptly pairs up with one of the binary members through a sling-shot interaction that ejects the remaining member. The short nature of this interaction type implies that GR effects and collisions rarely have a chance to affect the outcome.

A RI lasts instead from a few to several thousand of initial orbital times (e.g. Hut & Bahcall, 1983). The inclusion of GR effects and finite sizes in the EOM therefore have a relative high chance of affecting the dynamics during this type of interaction. How this in particular leads to an increase in the number of BH mergers can be understood as follows. As envisioned by Samsing et al. (2014), a RI can be described as a series of intermediate states (IMS), each of which is characterized by an IMS binary with a bound single. Due to pair-wise semi-chaotic exchanges of energy and angular momentum during the resonating evolution, each IMS binary has a finite probability for being formed with a high eccentricity, even when the initial binary is circular. If the eccentricity is high enough, GR will affect the EOM of the IMS binary (e.g. Peters, 1964), which most notably leads to an inspiral driven by GW emission that is followed by a prompt merger (Gültekin et al., 2006; Samsing et al., 2014). For even higher eccentricities the IMS binary members will instead undergo a prompt collision (e.g. Samsing et al., 2016, 2017).

2.2. Initial Conditions

The IC parameters that are relevant for our presented binary-single study, are listed and described in the following. An accompanying illustration is shown in Figure 1.

Semi-major axis

($a_{0}$). The semi-major axis of the initial target binary is denoted by $a_{0}$. As described later, both the classical and the relativistic limit can be studied by simply varying the value of $a_{0}$.

Velocity

($v$). We denote the rescaled relative velocity between the initial binary COM and the incoming single at infinity by $v$, as described in Section 2.1. In this work we set $v=0.01$ for all our scatterings, however, none of our conclusions depend on the exact value of $v$, as long as $v\ll 1$.

Binary phase

($f$). The orbital phase of the binary when the single is exactly $20\times a_{0}$ away from the binary COM is denoted by $f$. The value of $f$ is defined relative to an axis that is parallel to the velocity vector of the incoming single at infinity.

Impact parameter

($b$). We define the rescaled impact parameter at infinity by $b\equiv b_{\infty}/b_{c}$, where $b_{\infty}$ is the physical value and $b_{c}=a_{0}/v$. The value of $b$ can be positive or negative, where a positive (negative) $b$ here corresponds to that the angular momentum vectors of the binary and the incoming single, respectively, have the same (opposite) sign. The initial three-body angular momentum in the co-planar case will therefore equal zero for $b=-\sqrt{3/4}\approx-0.87$ (e.g. Hut, 1983b). We further note that the pericenter distance of the single w.r.t. the binary COM, $r_{\rm p}$, relates to $b$ as $r_{\rm p}/a_{0}=3b^{2}/8$.

Orbital plane angle

($\gamma$). The angle between the orbital planes of the initial binary and the incoming single is denoted by $\gamma$. In this notation, a co-planar interaction will have $\gamma=0$.

2.3. Outcomes and Kinematical Properties

The set of binary-single outcomes and kinematical properties we keep track of in this work, is listed and described below. Further details on defining and identifying outcomes in $N$-body simulations can be found in, e.g., Fregeau et al. (2004); Samsing et al. (2014, 2016).

2.3.1 Outcomes

Binary-single

(BS[ij]). The most common endstate is a binary with an unbound single. We denote this endstate by BS[ij], where ‘BS’ is short for ‘Binary-Single’, and [ij] refers to the object pair in the final binary. Throughout the paper we denote the two initial binary members by ‘1’ and ‘2’, and the incoming single by ‘3’.

GW inspiral

(GW). A close passage between any two of the three interacting objects, occasionally leads to a GW capture, where the two binary members undergo an inspiral through the emission of GWs (e.g. Gültekin et al., 2006; Samsing et al., 2014). We generally refer to this endstate by a ‘GW inspiral’, but will in our figures simply use ‘GW’ to keep the notation short.

Collision

(Coll). If two objects pass each other at a distance smaller than the sum of their radii without undergoing a GW inspiral first, then the outcome is denoted a ‘collision’, or in short by ‘Coll’. In our case, the collisional distance is simply $2R_{\rm s}$.

Very long interaction

(VLint). A finite fraction of the interactions can in theory last infinitely long time. This happens generally when the single in an IMS is sent out on a nearly unbound orbit. In this work we denote an interaction that did not finish within $2500\times T_{\rm orb}$ by ‘VLint’ , where $T_{\rm orb}$ is the orbital time of the initial target binary – with the only exception of the results shown in Figure 2, where the limit is set to $1000\times T_{\rm orb}$.

2.3.2 Kinematical Properties

Number of intermediate states

($N_{\rm IMS}$). For each interaction we count the number of times the binary-single system splits into an IMS. The total number from initial interaction to final outcome is what we denote $N_{\rm IMS}$.

Minimum distance

($R_{\rm min}$). For all interactions ending as a BS[ij] endstate, we save the minimum distance any two of the three objects undergo during the interaction. This distance is denoted by $R_{\rm min}$.

Binary-single energy

($E_{\rm BS}$). At the end of each interaction, we calculate the orbital energy of the single and the COM of the remaining pair. We denote this energy by $E_{\rm BS}$. In this notation, a BS[ij] endstate will have $E_{\rm BS}>0$, where collisions and GW inspirals generally have $E_{\rm BS}<0$, as they preferentially form during the interaction.

2.4. Topological Mappings

The graphical representation of the distribution of binary-single outcomes as a function of the ICs, is what we loosely refer to as the topology, in analogy with Hut (1983b). As an example, Figure 2 shows the distribution of binary-single endstates as a function of $f,\ b$, for $a_{0}=10^{-2}$ AU and $\gamma=0$. In our terminology, we refer to this distribution as the endstate topology, or the endstate-map. Building on the work of Hut (1983b), we focus in this study on topological maps derived from varying $f$ and $b$, while keeping $v$ and $\gamma$ fixed.

3. Topology of Binary-Single Interactions

In this section we explore the binary-single topology on large-scales. For this, we first study the classical limit, which represents the idealized case where the objects are all point-like particles interacting through Newtonian gravity only. After this follows a study of the relativistic limit, in which GR corrections and finite size effects start to affect the classical topology. In the last part we briefly explore how the topological maps change for varying $\gamma$.

For all scatterings performed in this section we keep finite sizes and GW emission in the EOM for consistency; the classical and the relativistic limit are therefore studied by simply varying $a_{0}$. We note that $f$ conveniently happens to approximately be the phase of the binary as the incoming single enters the binary. For example, at $f\approx 0,\ b\approx 0$, we find that all three objects are approximately on a line as the single enters. In our following descriptions we will therefore simply refer to $f$ instead of the exact binary phase as the single enters the binary. Furthermore, we often refer to just a single value of $f$, however, due to the rotational symmetry of the binary, results are similar for any phase $f+n\pi$, where $n$ is an integer.

3.1. Classical Limit

If the three interacting objects are point-like and only Newtonian gravity is included in the EOM, the system is said to be scale-free, as the topology and the relative number of each endstate do not depend on the SMA $a_{0}$ and mass $m$, as long as $v$ remains constant (Hut, 1983b; Hut & Bahcall, 1983). Although such an idealized system does not exist in nature, one approaches this classical scale-free limit as $a_{0}$ increases relative to the radius of the interacting objects. In this section we study the topology for $a_{0}=1$ AU, which provides a reasonable representation of the classical limit. Results are shown in Figure 3, and discussed below.

3.1.1 Endstate Topology

The endstate-map (Figure 3, upper left) shows a clear large-scale topological pattern. We find this to emerge from the distribution of two different small-scale topological pattern types, that we in the following refer to as the homogeneous-type and the random-type. The homogeneous-type is characterized by that neighboring points share the same endstate, where the random-type is characterized by that neighboring points have a semi-random distribution of endstates. As seen on the figure, the random-type forms one fully connected large-scale structure, which divides the homogeneous-type into a few disconnected islands. This is a feature of the ultra HB limit for which $v\ll 1$, as the random-type structure starts to separate into isolated parts as $v$ approaches $1$ (e.g. Hut, 1983b).

Comparing the $N_{\rm IMS}$-map (Figure 3, lower left) with the endstate-map, shows that the random-type generally originates from RIs and the homogeneous-type from DIs. This is a consequence of the chaotic nature of RIs, which makes the outcomes very sensitive to the ICs. Neighboring points on the $f,\ b$ grid in RI regions are therefore likely to have a semi-random distribution of endstates. In contrast DIs are always prompt. Small variations to the ICs will therefore not have a chance to undergo any significant growth during this type of interaction. Neighboring points on the $f,\ b$ grid in DI regions therefore often have similar endstates.

The $E_{\rm BS}$-map (Figure 3, upper right) relates closely to the other maps, as it is the orbital energy of the single w.r.t. the COM of the first binary that forms (endstate binary or IMS binary), denoted here by $E_{\rm 1BS}$, that determines if the interaction proceeds as a DI or a RI. By definition, $E_{\rm 1BS}>0$ leads to DI, in which case $E_{\rm BS}=E_{\rm 1BS}$, where $E_{\rm 1BS}<0$ leads to a RI. As the distribution of $E_{\rm 1BS}$ seems to vary smoothly across the $f,b$ space, we conclude that the boarders separating DIs from RIs – or equivalently the coastline of the DI islands – are simply the path for which $E_{\rm 1BS}=0$. As $E_{\rm 1BS}$ approaches zero along these boarders, the associated binary-single interaction time correspondingly approaches infinity, which explains why the ‘VLint’ endstate generally appears between DI and RI regions.

As seen, there is a clear asymmetry in the large-scale topology along the $b$-axis, despite corresponding symmetry in the positional configuration of the three objects. This asymmetry results from the way the binary rotates relative to the incoming single, where $b>0$ and $b<0$ correspond to prograde and retrograde motion, respectively. For $v\ll 1$, the incoming single almost has the same velocity as the binary members as it passes the binary. In the case of prograde orbits this leads to situations where the three objects almost appear stationary w.r.t. each other. As a result, for certain combinations of $f,b$ the incoming single does not directly encounter any of the binary members, which give rise to the large wave-like BS[12] structures seen for $b>0$. In contrast, negative values of $b$ always result in encounters between the single and a binary member due to the retrograde motion. A comparison with the high velocity results shown in Hut (1983b) confirms that the relative orbital motion indeed is the origin of the $b$-axis asymmetry. For example, the results in Hut (1983b) clearly show how the topology along the $b$-axis becomes increasingly symmetric as $v$ increases, or in other words, as the velocity of the binary members become less important.

3.1.2 Kinematical Examples

Having described a few topological relations in the sections above, we now turn to a more detailed kinematical analysis. However, we are here only able to provide quantitive descriptions, as we find that both RIs and DIs generally undergo relative complex orbital evolutions. We note that an analytical solution for describing the dynamics of the SB limit is possible (e.g. Hut, 1983a), where a statistical approach has been shown to provide valuable insight in the HB limit (e.g. Samsing et al., 2014). We divide our analysis below into the study of the topology around two orbital configurations: Configuration A ($f\approx 0$) and Configuration B ($f\approx\pi/2$). All results presented below are found by visually inspecting a representative set of three-body interactions.

Configuration A:

We first consider the change in topology for varying $b$ at $f\approx 0$. For this configuration, we refer to the binary member that is closest to the incoming single as it enters by ‘1’, the remaining member by ‘2’, and thereby the incoming single by ‘3’ (this notation matches the endstate labels on our figures for $f\approx 0+\pi$). This configurative labeling is also used in Figure 4.

Starting at $b\approx 0$, we see that the topology changes rather abruptly from a BS[23] DI region to a BS[13] DI region. The kinematics leading to these two neighboring regions is as follows. For slightly negative values of $b$ object 3 first encounters 1 from below which leads to an upwards deflection of 3 and a downwards deflection of 1. Along its new orbit 1 then interacts with 2, which sends 1 up towards 3 while 2 is being ejected downwards. As a result, [13] escapes as a binary. This interaction is illustrated to the left (labeled ‘E1’) in Figure 4. For slightly positive $b$ object 3 instead encounters 1 from above, which give rise to an interaction where the two objects almost perfectly replace each other, which leads to 1 flying to the left and 3 moving upwards. Object 1 then interacts with 2 from above, which sends 1 downwards and 2 upwards where it pairs up with 3. As a result, [23] escapes as a binary. This interaction is illustrated to the right (labeled ‘E2’) in Figure 4.

At $b\approx-1.5$ we see a region that is somewhat similar to the one at $b\approx 0$, but with the two BS[23] and BS[13] regions switched. For the region at $b\lesssim-1.5$ we find that 3 first interacts with 2 from below after which the binary has rotated $\approx\pi/2$. This interaction sends 2 downwards and 3 up towards 1. 1 and 3 then undergo a near parabolic hook-interaction, which sends 3 back towards 2 by ejecting 1 further upwards. As a result, [23] escapes as a binary. For the region at $b\gtrsim-1.5$ the first interaction between 3 and 2 is simply so strong that 2 is immediately kicked out by 3 which then takes its place in the binary. As a result, [13] ends up as a binary.

Configuration B:

We now consider the change in topology for varying $b$ near $f\approx\pi/2$. For this configuration, we refer to the binary member that moves towards the incoming single by ‘1’, the remaining member by ‘2’, and thereby the single by ‘3’. This configuration is similar to that shown in Figure 4, but for the binary rotated clockwise by $\pi/2$.

We first notice that for any value of $b\approx 0\pm 0.5$, the endstate will always be of the same type, in contrast to Configuration A considered above. The kinematics near $b\approx 0$ are as follows. As object 3 enters the binary, it first accelerates 1 in the upwards direction along its orbital tangent. Object 3 then interacts with 2 from below, which results in 2 being ejected downwards and 3 being sent up towards 1. As a result, [13] escapes as a binary. This basically happens for all $b$ for which 3 enters the binary on an orbit that goes in between 1 and 2.

Across larger variations of $b$, significant differences in the topology appear. For example, interactions near $b\approx 1.5$ all evolve as a RI, whereas interactions near $b\approx-1.5$ instead evolve as a DI. We find the kinematical differences to be as follows. At $b\approx-1.5$ object 3 no longer moves in between 1 and 2, but instead passes the binary from below. As it passes, 3 and 1 first interact, pulling 1 out of its circular motion along an orbit pointing towards the initial motion of 3. Object 3 then interacts with 2 which sends 3 back towards 1 by ejecting 2 in the opposite direction. As a result, [13] escapes as a binary. At $b\approx 1.5$, the kinematics is similar to that of $b\approx 0$, however, instead or 3 ejecting 2 into an unbound orbit, 2 remains bound after which the interaction then proceeds as a RI.

3.1.3 Close Encounter Topology

The $R_{\rm min}$-map (Figure 3, lower right), shows a new structural pattern, that is characterized by connected wavy bands. These bands map the closest distance any two of the three objects undergo before the system splits into its first binary-single state (BS[ij] or IMS). For example, the band that intersects $f\approx 0,\ b\approx 0$ and rises until $f\approx\pi/2,\ b\approx 1.0$, originates from the first close encounter between object 1 and 3, where ‘1’ here refers to the binary member that is closest to the incoming single ‘3’ at $f\approx 0,\ b\approx 0$. Since the trajectory of 3 is not perturbed significantly before it first encounters 1, the value of $b$ for which 1 and 3 undergo a close encounter is therefore approximately that of the shadow image of 1 along an axis parallel to $b$. The first part of the considered band therefore takes the shape of a sinusoidal curve, which explains why it first rises and then peaks at $f\approx\pi/2$. However, when it goes back towards $b\approx 0$ near $f\approx\pi$ the sinusoidal curve is strongly perturbed, because object 2 comes in between 1 and 3, which leads to deflections of 3 before it encounters 1. For velocities $v>1$ object 3 is hardly deflected, and the bands will therefore trace out a near perfect sinusoidal form.

The $R_{\rm min}$-map also shows a population of low $R_{\rm min}$ values inside the RI regions. This population originates from IMS binaries that are formed with high eccentricity and corresponding small pericenter distance.

3.2. Relativistic Limit

We now move to a study of the topology in the relativistic limit where finite sizes and GR corrections start to play a role. For this, we perform binary-single scatterings with the initial SMA set to $a_{0}=10^{-4}$ AU, which naturally leads to a significant number of very close pairwise BH encounters. This unrealistic low value of $a_{0}$ has only been chosen for illustrative purposes, however, the following results do apply for any $a_{0}$. Results are shown in Figure 5, and discussed below.

3.2.1 Black Hole Collisions

Looking at the endstate-map (Figure 5, top plot) and the $N_{\rm IMS}$-map (Figure 5, bottom plot), we see that BH collisions form through both DIs and RIs. These two collision populations are described in the following.

The DI-collisions trace the wavy $R_{\rm min}$ bands shown in Figure 3, with only small differences due to the modified dynamics resulting from including GW emission in the EOM. The bands broaden as the value of $b$ increases. This is due to the corresponding decrease in the relative velocity between the binary member in question and the incoming single, which leads to an increase in the gravitational focusing, and thereby to the seen broadening in $b$. As the value of $b$ decreases, the relative velocity increases, which makes the gravitational focusing less efficient. The orbit of the incoming single therefore has to be relatively fine-tuned to collide with a binary member, which explains the thinner bands at negative $b$.

The RI-collisions originate from the population of IMS binaries that are formed with a pericenter distance that is smaller than the collisional distance $2R_{\rm s}$. Despite the apparent randomness of the RI regions, the collision endstates do not seem to be completely randomly distributed in these regions. For example, the relative number of collisions clearly decreases for $b\gtrsim 2$. This partially results from the corresponding increase in the total angular momentum, which naturally leads to fewer low angular momentum orbits.

Further insight into the collision endstates distribution is provided in Figure 6, which shows what pair that collides. Besides bands originating from collisions with the incoming single (Coll[13], Coll[23]), we interestingly see a band passing through $f\approx 0.7,\ b\approx-0.5$ that is formed by collisions of the binary members themselves (Coll[12]). We find that the corresponding dynamics is such that the incoming single, during its first passage, simply perturbs the two binary members into a prompt collisional orbit. Although this seems like a fine-tuned problem, the corresponding bands appear relative broad, which indicates a non-negligible formation probability.

3.2.2 GW Inspirals

A clear population of GW inspirals is seen near the boarders of the wavy collision bands. The origin of this inspiral population is not surprising, as the collision boarders are characterized by near collisional orbits. When GW emission is included in the EOM, these orbits therefore naturally result in strong GW emission, which then leads to the observed populations of GW inspirals. Just as the population of collisions in the wavy bands forms through DIs, so do the GW inspirals forming along the bands, as seen in the bottom plot in Figure 5.

We note that there are a few collisional bands located in the DI regions with no GW inspirals forming along their edges. We find here that the dynamics leading to the prompt collision at the same time ejects the third BH into an unbound orbit. The endstates populating these bands are therefore characterized by a collision with an unbound companion. This explains why no GW inspirals form along their edges, as the binary-single system splits up into a BS[i,j] state before the GW inspiral is complete. The two binary objects will of course still undergo the GW inspiral, but it will happen without a bound companion.

Moving away from the wave-like bands, we see how GW inspirals clearly form throughout the RI regions. We note here that collisions and GW inspirals do not have to follow the same distribution, as the formation of a GW inspiral is not related to a particular value of the IMS binary pericenter distance unlike a collision. Instead, the formation of a GW inspiral relates to how much energy it can lose through GW emission during the orbital time it is isolated from the bound single (Samsing et al., 2014). As a result, GW inspirals can form from IMS binaries with both small and large pericenter distances, depending on the corresponding isolation time, and can therefore appear in regions where collisions are relative rare. For example, in the upper plot in Figure 5 there are many GW inspirals forming for $b>2.5$ but only relative few collisions. A more detailed picture is provided in Figure 6, which also shows what pair that undergoes the GW inspiral.

3.3. Out-of-Plane Interactions

For completeness, we here study the change in topology for varying orbital-plane angle $\gamma$. Endstate-maps for $\gamma=\pi/4$ and $\gamma=\pi/2$ are shown in Figure 7, and briefly discussed below.

The most notable difference between the co-planar case ($\gamma=0$) and our presented out-of-plane cases ($\gamma=\pi/4,\ \pi/2$), is a breakup of the connected wave-like collision and GW inspiral bands seen for $\gamma=0$. To understand this change, we note that the incoming single can only undergo a prompt collision or a prompt GW inspiral near regions where its orbital plane intersects with the binary orbit. When $\gamma=0$, the incoming single intersects for sinusoidal variations in $f,b$, however, when $\gamma>0$, the orbital plane of the single is only able to intersect the binary orbit at maximum two points, namely at $\{-\pi/2,+\pi/2\}$. This correspondingly result in isolated regions in the topological map instead of connected bands.

As with the co-planar case, the DI collision regions for $\gamma>0$ are surrounded by GW inspirals that form through a prompt capture. However, a large population of both DI and RI collisions and GW inspirals also form beyond these DI regions. In fact, we find that the total number of collisions and GW inspirals is not strongly dependent on $\gamma$, although the distribution of orbital parameters at merger might depend on $\gamma$ due to the corresponding change of angular momentum. We save such a study for future work.

4. The Role of GW emission

We now perform a slightly more controlled study on how the inclusion of GW emission in the EOM affects the topology on different scales. For this, we explore two different topological regions, with and without GW emission included in the EOM. For all interactions we set $a_{0}=10^{-4}$ AU and $\gamma=0$ to have GW emission effects stand out clearly. Results are discussed below.

4.1. Large-Scale Effects

We start by analyzing how the inclusion of GW emission affects the large-scale topology. Results are shown Figure 8.

When GW emission is not included in the EOM (top plot), we see that the collision regions simply appear superimposed on the classical endstate-map, with a distribution that follows the $R_{\rm min}$ distribution as expected. Other than occupying part of the map, the inclusion of collisions does not affect the large-scale topological structure, as the EOM is unchanged.

When GW emission is included (bottom plot), we see that this not only leads to the formation of GW inspirals, but also a slight modification of the topology. For example, at $f\approx 0.75$, $b\approx 0.75$, where the collision band intersects the RI region in the classical map, the inclusion of GW emission clearly smoothes out the transition region, by turning a part of the BS[1,3] endstates formed in the DI region into a semi-random mix of outcomes. Generally, one finds that the largest topological changes indeed are found near regions where collision endstates, DIs, and RIs intersect. This results from two factors: the GW emission is strong (the objects are close to colliding) and the the orbital energy $E_{\rm 1BS}$ is near zero (the boarder between a DI and a RI is characterized by $E_{\rm 1BS}=0$). This combination makes it possible for the GW emission to change a classical DI into a RI, by emitting the orbital energy out of the system that otherwise would have ejected the single into an unbound orbit. In fact, in Samsing et al. (2016) it was briefly discussed if this way of ‘capturing’ the single into a RI through GW emission would be able to change the relative frequency of endstates. We have here shown and argued that this indeed happens, but only in a very small part of the available phase space. The effect is therefore negligible from a statistical perspective.

4.2. Small-Scale Effects

We now study how the topology changes on small scales as GW emission is included in the EOM. Results are shown in Figure 9.

In the top plot, where GW emission is not included, we first notice that the considered RI region has a surprisingly rich and complex micro-topology. This will be studied in greater detail in Section 5. As in the large-scale example, we see here how the collision regions simple appear as sharp shadows superimposed on top of the classical map.

We now turn to the bottom plot where GW emission is included. As seen, the inclusion drastically changes the topology close to the classical collision bands. Not only is a clear population of GW inspirals forming along the bands, but the GW emission further results in band-like structures that alter between collisions and GW inspirals. To understand the kinematics leading to this alternation we visually inspected a set of binary-single interactions for $f\approx\pi/2$ and varying $b$. For this configuration we here denote the binary member moving away from the incoming single by ‘2’, the remaining member by ‘1’, and thereby the incoming single by ‘3’. We find the kinematics to be as follows. First, the large collisional band at the bottom ($b\approx 1.15$) is formed by 2 directly colliding with 3 as it enters the binary, as also described in Section 3.1.3. Above this is an inspiral band ($b\approx 1.3$), which we find consists of interactions where 3 instead of directly colliding with 2, slightly grazes it, which then leads to a prompt GW inspiral between the two. The subsequent collision band ($b\approx 1.4$) is formed by 3 first being captured into a relative wide low angular momentum orbit by 2, after which 3 and 2 collide at their first pericenter passage along their new orbit. The following narrow BS[23] band ($b\approx 1.42$) arises from a rather fine-tuned dynamical configuration where 1 becomes kinematically unbound just before 2 and 3 undergo a collision or a GW inspiral. This band is therefore sensitive to the tidal threshold used for labeling a BS[ij] endstate, and is therefore not a fundamental feature. The next few bands form in a similar way to the previous collision band, only with 2 and 3 changing between undergoing a collision or a GW inspiral. To conclude, these topological changes all relate to the capture of object 3 into a bound orbit partly through GW emission. As the capture happens at low angular momentum, and angular momentum further is radiated away during the capture, the system evolves in a state that is likely to result in either a prompt collision or a GW inspiral.

5. Topological Microstructures

We now explore the topology of the microstructures that surprisingly seem to appear within RI regions. For this, we perform scatterings for $a_{0}=1$ AU, and $\gamma=0$, with collisions and GW emission included in the EOM. However, we do note that our considered small scale topology is likely to change slightly if additional PN terms are added to the EOM, we therefore limit ourself to general conclusions. Results from two selected zoom in regions, referred to as ‘Zoom Region 1’ and ‘Zoom Region 2’, respectively, are shown in Figure 10, and discussed below.

5.1. Zoom Region 1

The top plot shows a zoom in on one of the larger RI regions. As seen, the region does not only consist of semi-randomly distributed endstates, as it appears to on large scales, but in fact hosts clear micro-topological structures of various kinds.

To gain insight into the dynamics of this region, we visually investigated a representative set of binary-single interactions. From this exercise we found that most of the interactions near the center ($f\approx 1.3,\ b\approx 1.2$) start out in a similar way, where the objects first interact non-chaotically with their orbits tracing out a well defined pattern with a clear degree of symmetry. A schematic illustration of this primary orbital interaction assuming perfect symmetry is shown in Figure 11. Although the interaction in this example promptly ends as a BS[23] endstate, the outcomes do vary depending on the exact ICs. In a few cases we even observed that the system after this primary interaction entered a quasi-stable state which it stayed in for 5-10 cycles. We note here that the more stable a three-body state is, the less likely it is to form (e.g. Hut, 1983b), which implies that fully stable three-body states are not assembled through our binary-single experiments (Newtonian three-body states can not be stable at plus infinity if they are not stable at minus infinity). Of course the inclusion of GR in the EOM immediately changes the situation, as the implied GW emission destabilizes the system within a finite time. Further studies on the possible assembly of quasi-stable three-body states with and without GR is beyond this paper, but still of great theoretical interest (e.g. Šuvakov & Dmitrašinović, 2013).

Although the state illustrated in Figure 11 is not stable, it certainly classifies as a special interaction due to its symmetry and interesting orbital topology. However, the classification of such interactions is not entirely clear and to our knowledge still unexplored. A systematic search through the $f,\ b$ space could reveal if this is a member of a larger family of peculiar interactions, or just a single interesting incident.

5.2. Zoom Region 2

The bottom plot shows a zoom in on a transition zone between RIs (lower left) and DIs (upper right). As seen, the only pattern that is clearly visible consists of nested rings with a spacing that decreases as one moves closer to the DI boundary from the RI side. Focusing on the BS[13] endstate bands, we find that the corresponding IMS binary configuration is such that when the single returns from its first excursion, the resultant IMS binary-single interaction leads to a prompt BS[13] endstate. This happens in each band, with the only difference being how many times the IMS binary rotates before the single returns. As a result, there is a band for each integer $n$, where $n$ is the number of times the IMS rotates before the single returns (by visual expectation we find that the band in the lower left corner corresponds to $n=10$). As $E_{\rm 1BS}$ approaches zero as one moves towards the DI region, the number $n$ will correspondingly approach infinity. As a result, the spacing between the blue bands becomes smaller and smaller, approaching zero near the boundary. This explains why substructures with infinitely small spacings exist and where they are located. The precess of such structures implies that the true richness of micro topological structures cannot be studied with a finite grid approach as the one we have presented here.

6. Double Gravitational Wave Mergers

All GW inspirals form during the binary-single interaction while the three objects are still bound to each other in an IMS (Samsing et al., 2014). Ignoring dynamical kicks caused by asymmetric GW emission (e.g. Baker et al., 2006), this interestingly implies that the BH formed from the GW inspiral merger and the remaining single are still bound to each other after the three-body interaction. A compelling question is how long time it takes for this binary to merge through GW emission, or equivalently, what the time span is from the first to the second GW merger. If the time span is short enough, it might be possible for LIGO to observe both GW mergers, a scenario that undoubtful would be extremely exiting.

In this section we explore for which combinations of $f$ and $b$ the time span from first to second GW merger, denoted by ${\Delta}t_{12}$, takes its shortest value. A more detailed study of this double GW merger scenario is reserved for a separate paper, and the following discussions are therefore kept short. An illustration of our proposed double GW merger channel is shown in Figure 12, and results presented below.

6.1. Time from First to Second GW Merger

The distribution of $\log({\Delta}t_{12}/T_{0})$, where $T_{\rm 0}$ is the orbital time of the initial target binary, as a function of $f,\ b$ for $\gamma=0$, is shown in Figure 13. For this we have included all GW inspirals and BH collisions that form with a bound companion, where ${\Delta}t_{12}$ is estimated using the equations from Peters (1964). To keep this primary study simple and clear, we have not incorporated effects related to relativistic mass losses or asymmetric GW emission; corrections that in general would lead to an increase in ${\Delta}t_{12}$.

As seen on the figure, the minimum value for ${\Delta}t_{12}$ is highly sensitive to $b$, but not $f$. The reason is that only $b$ relates to the total angular momentum, which in turn directly relates to the corresponding GW inspiral time (Peters, 1964). This further explains why the value for $b$ that minimizes ${\Delta}t_{12}$ is simply the value for which the total three-body angular momentum equals zero. As derived in Section 2.2, the angular momentum of the single and the binary cancels out for $b\approx-0.87$, which clearly is where ${\Delta}t_{12}$ takes its minimum as confirmed by the simulation results shown in the figure. Close to this value of $b$ the three objects can, as a result of the vanishing angular momentum, promptly merge without any significant use of GW emission. This leads to values of ${\Delta}t_{12}$ that are of order the orbital time of the initial binary, instead of its initial GW inspiral time.

The angular momentum cancelation is less efficient in out-of-plane interactions, which explains why the co-planar interaction case is the relevant channel to consider for generating double GW mergers. This strongly motivates further studies with GR effects included in the $N$-body EOM of co-planar few-body interactions taking place in AGN environments.

7. Conclusions

In this paper, we explored how the outcomes of binary-single BH interactions distribute as a function of the ICs. For this we studied the corresponding graphical representation of the distribution (Figure 2) – a representation we loosely refer to as the topology in analogy with Hut (1983b). Our study shows for the first time how this binary-single topology changes when GR corrections and BH finite sizes are consistently included in the $N$-body EOM. Such corrections give rise to a non-negligible distribution of BH collisions and GW inspirals that is directly observable by LIGO (e.g. Samsing et al., 2014; Samsing & Ramirez-Ruiz, 2017).

Using zoom-in simulations, we resolved for the first time the emergence of micro-topological structures residing inside RI regions (Figure 10). By visually inspecting orbits near these structures, we found that many of the interactions (in the RI region centered around $f\approx 1.3,\ b\approx 1.25$) start out in a similar dynamical way characterized by a particular ordered and symmetrical orbital evolution (Figure 11). To our knowledge, the classification and assembly of such symmetrical orbital evolutions is still unexplored, and a systematic search for these throughout the $f,\ b$ space could therefore be of potential interest.

Finally, we pointed out that the inclusion of GW emission in the $N$-body EOM naturally leads to interactions resulting in two GW mergers. In this regard, we described how the co-planar case can lead to scenarios where the time span from the first to the second GW merger is short enough for LIGO to observe both events. This double GW merger channel could have profound astrophysical consequences, and will therefore be explored in greater detail in an upcoming paper.

From an astrophysical perspective, our most relevant result is that the distribution of outcomes of binary-single interactions is highly dependent on the ICs. While this certainly applies to small variations to the ICs, due to the chaotic nature of the three-body problem, we surprisingly find that also large scale variations to the ICs give rise to corresponding large scale differences in the outcome distribution. For example, in the co-planar case ($\gamma=0$) the distribution is clearly different between prograde ($b>0$) and retrogade ($b<0$) interactions (Figure 2). This also applies to how the relative number of BH mergers distributes with the ICs (Figure 6). Clear differences are also seen when the orbital-plane angle $\gamma$ is varied (Figure 7). This means that if properties of these distributions, such as the eccentricity distribution at merger, could be observed by LIGO, then one would be able to roughly deduce the underlying distribution of ICs, which could give important hints about the astrophysical origin. As an example, in AGN environments (e.g. McKernan et al., 2012; Bellovary et al., 2016; McKernan et al., 2017), the interactions will not only be close to co-planer, but they will also likely be unequally distributed between prograde or retrograde interactions, depending on the formation history of the interacting BHs and their interaction with the disc. Our idealized study indicates that the outcome distribution is clearly different for these cases. A natural extension to our presented work could therefore be to derive observables related to BH collisions, GW inspirals, and high eccentricity BH mergers, for different scenarious including retrograde interactions, prograde interactions, in-plane interactions, out-of-plane interactions, and for varying mass hierarchy.

Although our paper presents details to a level that is not directly observable, we do note that similar studies on how the distribution of BH mergers and related observables change with the ICs, are likely to play a key role in distinguishing different astrophysical GW merger channels from each other. An exercise that will become highly relevant when the number of observed GW mergers starts to grow.

Support for this work was provided by NASA through Einstein Postdoctoral Fellowship grant number PF4-150127 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060.