Dissipative Evolution of Unequal Mass Binary-Single Interactions and its Relevance to Gravitational Wave Detections

TL;DR

This paper develops an analytical framework to estimate the formation rates of compact binaries through energy-loss mechanisms like gravitational wave emission during three-body interactions, with implications for gravitational wave detection.

Contribution

It introduces a new formalism to calculate cross sections of two-body captures in binary-single interactions considering energy dissipation, enhancing understanding of dynamical binary formation.

Findings

Two-body gravitational wave captures can form during three-body interactions.

The estimated rates are relevant for Advanced LIGO detections.

Formation of eccentric binaries has observable consequences.

Abstract

We present a study on binary-single interactions with energy loss terms such as tidal dissipation and gravitational wave emission added to the equation-of-motion. The inclusion of such terms leads to the formation of compact binaries that form during the three-body interaction through two-body captures. These binaries predominantly merge relative promptly at high eccentricity, with several observable and dynamical consequences to follow. Despite their possibility for being observed in both present and upcoming transient surveys, their outcomes are not firmly constrained. In this paper we present an analytical framework that allows to estimate the cross section of such two-body captures, which permits us to study how the corresponding rates depends on the initial orbital parameters, the mass hierarchy, the type of interacting objects, and the energy dissipation mechanism. This formalism is applied here to study the formation of two-body gravitational wave captures, for which we estimate absolute and relative rates relevant to Advanced LIGO detections. It is shown that two-body gravitational wave captures should have compelling observational implications if a sizable fraction of detected compact binaries are formed via dynamical interactions.

Equations

$${\sigma }_{X_{\mathrm{ij}}}=\frac{N_{X_{\mathrm{ij}}}}{N_{\mathrm{tot}}}\times \pi b_{\max }^2,$$

$$b_{\max }=\sqrt{\frac{2G{m}_{\mathrm{bs}}{r}_{\max }}{v_{\infty }^2}},$$

$${\sigma }_{X_{\mathrm{ij}}}=\frac{N_{X_{\mathrm{ij}}}}{N_{\mathrm{tot}}}\times \frac{2\pi G{m}_{\mathrm{bs}}{r}_{\max }}{v_{\infty }^2}.$$

$$v_{\mathrm{c}}=\sqrt{\frac{G{m}_1{m}_2}{m_3(m_1+m_2)}\frac{m_{\mathrm{bs}}}{a_0}}.$$

$${\sigma }_{\mathrm{CI}}=\mathcal{C}\frac{2\pi G{m}_{\mathrm{bs}}{a}_0}{v_{\infty }^2}.$$

$${\sigma }_{X_{\mathrm{ij}}}\approx P(X_{\mathrm{ij}}|\mathrm{CI})\times {\sigma }_{\mathrm{CI}},$$

$${\Gamma }_{X_{\mathrm{ij}}}^{(V)}\approx \frac{N_{\mathrm{bin}}{N}_{\sin }{v}_{\infty }}{V}\int {\sigma }_{X_{\mathrm{ij}}}{P}_{\mathrm{bin}}(a_0)d{a}_0.$$

$$\frac{d{E}_{\mathrm{orb}}}{dt}\approx \frac{\Delta E_{\mathrm{p}}}{T_{\mathrm{orb}}}=\frac{\Delta E_{\mathrm{p}}}{G{m}_{\mathrm{ij}}}\frac{\sqrt{2}}{\pi }{\mu }_{\mathrm{ij}}^{-3/2}{E}_{\mathrm{orb}}^{3/2},$$

$$t_{\mathrm{insp}}=2\pi \sqrt{G{m}_{\mathrm{ij}}}{\mu }_{\mathrm{ij}}\frac{\sqrt{a}}{\Delta E_{\mathrm{p}}}.$$

$$\Delta E_{\mathrm{p}}=\mathcal{E}\frac{G{M}^2}{\mathcal{R}}{\left(\frac{\mathcal{R}}{r_{\mathrm{p}}}\right)}^{\beta },$$

$$t_{\mathrm{insp}}=2\pi r_{\mathrm{p}}^{\beta }\sqrt{a}\left[\frac{m_{\mathrm{ij}}{\mu }_{\mathrm{ij}}{\mathcal{R}}^{1-\beta }}{\mathcal{E}{M}^2\sqrt{G{m}_{\mathrm{ij}}}}\right].$$

$$t_{\mathrm{iso}}=2\pi \sqrt{\frac{a_{\mathrm{bs}}^3}{G{m}_{\mathrm{bs}}}},$$

$$\frac{m_1{m}_2}{2a_0}=\frac{m_{\mathrm{i}}{m}_{\mathrm{j}}}{2a}+\frac{m_{\mathrm{ij}}{m}_{\mathrm{k}}}{2a_{\mathrm{bs}}},$$

$$a_{\mathrm{bs}}=a_0\left(\frac{m_{\mathrm{ij}}{m}_{\mathrm{k}}}{m_1{m}_2}\right)\left(\frac{a^{\prime }}{a^{\prime }-1}\right),$$

$$a^{\prime }=\frac{a}{a_{\mathrm{c}}},\mathrm{and},a_{\mathrm{c}}=a_0\left(\frac{m_{\mathrm{i}}{m}_{\mathrm{j}}}{m_1{m}_2}\right).$$

$$t_{\mathrm{iso}}=2\pi \frac{a_0^{3/2}}{\sqrt{G{m}_{\mathrm{bs}}}}{\left(\frac{m_{\mathrm{ij}}{m}_{\mathrm{k}}}{m_1{m}_2}\right)}^{3/2}{\left(\frac{a^{\prime }}{a^{\prime }-1}\right)}^{3/2}.$$

$${ϵ}_{I_{\mathrm{ij}}}={\mathcal{E}}^{1/\beta }\mathcal{M}{\left(a_0/\mathcal{R}\right)}^{(1/\beta -1)}\mathcal{G}(a^{\prime },\beta ),$$

$$\mathcal{G}(a^{\prime },\beta )=a^{\prime (1/\beta -1)}{\left(a^{\prime }-1\right)}^{-3/(2\beta )},$$

$${ϵ}_{I_{\mathrm{ij}}}=(1-e_{I_{\mathrm{ij}}}),$$

$$\mathcal{M}=\left(\frac{m_1{m}_2}{m_{\mathrm{i}}{m}_{\mathrm{j}}}\right){\left[{\left(\frac{M}{m_{\mathrm{bs}}}\right)}^2{\left(\frac{m_{\mathrm{bs}}}{{\mu }_{\mathrm{ij}}}\right)}^{3/2}\left(\frac{m_{\mathrm{k}}{m}_{\mathrm{k}}}{m_1{m}_2}\right){\left(\frac{m_{\mathrm{ij}}}{m_{\mathrm{k}}}\right)}^{1/2}\right]}^{1/\beta }.$$

$$F_{\mathrm{tid}}$$

$$F_{\mathrm{bin}}$$

$$a_{\mathrm{l}}^{\prime }$$

$$a_{\mathrm{u}}^{\prime }$$

$$P(I_{\mathrm{ij}}|\mathrm{CI})\approx ⟨N_{\mathrm{IMS}}⟩P({\mathcal{R}}_{\mathrm{I}}),$$

$$P({\mathcal{R}}_{\mathrm{I}})\approx P({\mathcal{R}}_{\mathrm{U}})P({\mathcal{R}}_{\mathrm{I}}|{\mathcal{R}}_{\mathrm{U}}),$$

$$P(I_{\mathrm{ij}}|\mathrm{CI})\approx A_{\mathrm{I}}\left[\frac{⟨N_{\mathrm{IMS}}⟩P({\mathcal{R}}_{\mathrm{U}})}{A_{\mathrm{U}}}\right],$$

$$\mathcal{N}\equiv \left[\frac{⟨N_{\mathrm{IMS}}⟩P({\mathcal{R}}_{\mathrm{U}})}{A_{\mathrm{U}}}\right]$$

$$A_{\mathrm{I}}\approx {\int }_{a_{\mathrm{l}}^{\prime }}^{a_{\mathrm{u}}^{\prime }}{ϵ}_{I_{\mathrm{ij}}}d{a}^{\prime }={\mathcal{E}}^{1/\beta }\mathcal{IM}{\left(a_0/\mathcal{R}\right)}^{(1/\beta -1)},$$

$$\mathcal{I}={\int }_{a_{\mathrm{l}}^{\prime }}^{a_{\mathrm{u}}^{\prime }}\mathcal{G}(a^{\prime },\beta )d{a}^{\prime }.$$

$$\mathcal{I}\approx \frac{1.05}{(1-1.7/\beta )}\frac{\text{ln}(a_{\mathrm{u}}^{\prime })}{(a_{\mathrm{u}}^{\prime }-1{)}^{1/(\beta +1)}}.$$

Full text

11affiliationtext: Department of Astrophysical Sciences, Princeton University, Peyton Hall, 4 Ivy Lane, Princeton, NJ 08544, USA22affiliationtext: School of Natural Sciences, Institute for Advanced Study, 1 Einstein Drive, Princeton, New Jersey 08540, USA33affiliationtext: Department of Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064, USA44affiliationtext: Niels Bohr Institute, University of Copenhagen, Blegdamsvej 17, 2100 Copenhagen, Denmark**affiliationtext: Einstein Fellow

1. Introduction

With the recent detection of gravitational waves (GWs) by LIGO (Abbott et al., 2016a, b, 2017), and upcoming electromagnetic (EM) surveys including LSST (LSST Science Collaboration et al., 2009), JWST (Gardner et al., 2006) and WFIRST (Spergel et al., 2013), a new exciting era in transient astrophysics is eminent. To learn from the population of current detections, that for LIGO sources can be up to 250 Gpc*-3* yr*-1* (The LIGO Scientific Collaboration et al., 2016; Belczynski et al., 2016; de Mink & Mandel, 2016), major effort is now being made to model and understand the range of possible avenues such transient sources can be assembled in the Universe (e.g. Zevin et al., 2017).

One viable avenue to assemble GW sources is through few-body dynamics, where exotic outcomes can form as a result of chaotic or secular exchange of energy and angular momentum. However, not all few-body interactions are equally probable; in high stellar density environments, such as a globular cluster (GC), one finds the most frequent few-body encounters to be binary-single interactions (e.g. Heggie, 1975; Hut & Bahcall, 1983; Hut, 1983, 1993). Such three-body interactions are not only important dynamically (Heggie, 1975), but have also been shown to potentially play a role for EM transients (e.g. Hut & Inagaki, 1985; McMillan, 1986; Fregeau et al., 2004; McMillan & Portegies Zwart, 2007; Rosswog et al., 2009; Lee et al., 2010; Mapelli et al., 2013; Mapelli & Zampieri, 2014; Ramirez-Ruiz et al., 2015; MacLeod et al., 2016; Samsing et al., 2016; Perets et al., 2016), as well as to the assembly of binary black hole (BBH) mergers (Gültekin et al., 2006; Ziosi et al., 2014; Samsing et al., 2014; Spera et al., 2015; Rodriguez et al., 2015; Kimpson et al., 2016; Samsing et al., 2016; Rodriguez et al., 2016a, b; Samsing & Ramirez-Ruiz, 2017).

To understand the evolution of dynamical systems, many $N$-body simulations have been performed, from detailed few-body studies (e.g. Hut & Bahcall, 1983; Samsing et al., 2014; Antognini & Thompson, 2016; Samsing et al., 2016), to full GC simulations (e.g. Wang et al., 2016). Such $N$-body simulations are often evolved using an equation-of-motion (EOM) without terms correcting for possible orbital energy losses and finite size effects. However, recent work have shown that a consistent inclusion of such terms surprisingly seems to play a major role for estimating accurate rates of transients and their corresponding distribution of orbital parameters at merger (e.g. Gültekin et al., 2006; Samsing et al., 2014, 2016; Samsing & Ramirez-Ruiz, 2017).

In this paper we continue our study of binary-single interactions with correction terms that include finite sizes, dynamical tides, and GW emission. The main effect from the inclusion of such terms in the EOM is the formation of high eccentricity compact binaries that assemble during three-body interactions (Gültekin et al., 2006; Samsing et al., 2014, 2016; Samsing & Ramirez-Ruiz, 2017). In order to study the relative importance of this highly eccentric merger population, this paper develops an analytical framework that allows us to robustly estimate their cross section and corresponding rates. This formalism can be used for merging binaries driven by tidal dissipation or gravitational wave emission.

We apply our analytical framework to study the assembly of BBH mergers forming through binary-single interactions as a result of GW emission and its dependence on the initial orbital parameters and the mass hierarchy. This exercise allows us to generalize our previous results, which were derived for equal-mass interactions (Samsing et al., 2014, 2016; Samsing & Ramirez-Ruiz, 2017), and robustly conclude that $>1\%$ of all BBH mergers assembled through binary-single interactions should manifest as two-body GW captures. As the majority of these GW captures form with high eccentricity, our study even suggests that the binary-single interaction channel is likely to dominate the rate of high eccentricity BBH mergers observable by LIGO. Besides being a major motivation for developing non-circular GW templates (e.g. Harry et al., 2016; Huerta et al., 2016), we further conclude that the notable eccentricity associated with GW inspirals might also turn out to play a leading role in helping differentiate between the range of viable BBH merger formation channels. The formation of electromagnetic transients through tidal two-body captures will be explored in a separate paper.

The paper is organized as follows. In Section 2 we present a general analytical framework for calculating the cross section and corresponding rate of two-body captures and collisions that form during binary-single interactions. In Section 3 and 4, we then use this framework to study the absolute and relative formation rates of GW mergers. In Section 5 we use our findings to estimate the total number of BBH mergers assembled throughout the history of a typical GC, as well as current observable rates relevant to Advanced LIGO. Our conclusions are then summarized in Section 6.

2. Close Encounters In Binary-Single Interactions

Binaries interacting with singles in dense stellar systems, such as GCs, often undergo highly chaotic and resonating evolutions before reaching an end-state (e.g. Hut & Bahcall, 1983). During such evolutions, the three objects generally undergo several close pairwise encounters, each of which will either result in a collision or orbital energy losses through mechanisms such as tides and GW emission (e.g. Samsing et al., 2016). While a collision simply leads to a coalescence of the two objects (e.g. Fregeau et al., 2004), the fate of objects driven by orbital energy losses is less understood. However, as illustrated in Gültekin et al. (2006); Samsing et al. (2014, 2016), if the energy loss is significant, the two objects will undergo a dissipative capture resulting in the formation of a compact binary. Depending on the energy loss mechanism, this binary either ends up promptly merging or settling into a tight quasi-stable configuration. In this paper we refer to such two-body dissipative captures as inspirals. An example of an inspiral forming through GW energy losses is shown in Figure 1. If the two objects instead undergo a passage that is smaller than their total radii without inspiraling first, the outcome is here referred to as a collision.

In the sections that follow we calculate the cross section and corresponding rate of inspirals and collisions forming in binary-single interactions. We include the possibility for the three interacting objects to have different masses, which represents a major improvement to our previous work (Samsing et al., 2014, 2016). This allow us to study the role of dissipative effects including tides and GW emission in binary-single systems involving realistic combinations of white dwarfs (WDs), neutron stars (NSs), and BHs.

2.1. Cross Sections and Formation Rates

We start this section by deriving a set of quantities relevant for estimating general cross sections and rates. Below, and in the remaining parts of the paper, we denote an inspiral outcome between object pair $[i,j]$ by $I_{\rm ij}$, and a collision outcome by $C_{\rm ij}$.

2.1.1 Cross Section

Imagine $N_{\rm tot}$ uncorrelated scattering experiments between a randomly orientated binary and an incoming single isotropically sampled across a disk at infinity. In this case the cross section for object pair $[i,j]$ to result in a given outcome $X$, denoted by $X_{\rm ij}$, can be estimated by the following product (Hut & Bahcall, 1983),

[TABLE]

where $N_{X_{\rm ij}}$ is the number of interactions ending as outcome $X_{\rm ij}$, and $b_{\rm max}$ is the radius of the disk at infinity. Here $X_{\rm ij}$ could for example denote an inspiral outcome $I_{\rm ij}$, or a collision outcome $C_{\rm ij}$. The value of $b_{\rm max}$ should be large enough for the triple system to be able to probe all possible pathways that could result in $X_{\rm ij}$. To determine a proper value for $b_{\rm max}$ it is useful to work with the corresponding pericenter distance, $r_{\rm max}$, between the binary center-of-mass (COM) and the incoming single. This distance relates to $b_{\rm max}$ as follows (Samsing et al., 2014),

[TABLE]

where we have assumed the gravitational focusing limit. Here $m_{\rm bs}$ is the total mass of the binary-single system, and $v_{\rm\infty}$ is the initial relative velocity between the binary and the single at infinity. The cross section $\sigma_{X_{\rm ij}}$ can now be expressed in terms of $r_{\rm max}$ as,

[TABLE]

Now, it has been shown that both inspirals and collisions predominantly form in resonant interactions (RIs) due to their long lived chaotic nature that makes it possible for the objects to undergo multiple close passages (Samsing et al., 2014, 2016). A triple system can only enter such a bound resonant state if the initial total orbital energy is negative, which is the case for $v_{\infty}<v_{\rm c}$, where $v_{\rm c}$ is a characteristic velocity given by (Hut & Bahcall, 1983),

[TABLE]

Here $a_{0}$ is the semi-major axis (SMA) of the initial target binary, and the indices ‘1’, ‘2’, and ‘3’ refer to the two objects initially in the target binary and the incoming single object, respectively (see Figure 2). The limit where $v_{\infty}<v_{\rm c}$ is usually referred to as the hard-binary (HB) limit, where $v_{\infty}>v_{\rm c}$ is referred to as the soft-binary (SB) limit. Not all initial conditions (ICs) can lead to a RI in the HB limit, and we will in this work therefore refer to a binary-single interaction which could result in a RI as a close interaction (CI). The value for $r_{\rm max}$ which marks the limit for when interactions with pericenter $r_{\rm p}<r_{\rm max}$ are all CIs is here denoted by $r_{\rm CI}$. Although it can be shown that no exact value of $r_{\rm CI}$ can be defined, one can generally think of $r_{\rm CI}$ as representing the limit dividing democratic and hierarchical RIs (e.g. Hut, 1993; Samsing et al., 2014); a distance that naturally relates to the tidal disruption distance of the binary. In this work we take $r_{\rm CI}=\mathscr{C}a_{0}$, where $\mathscr{C}$ is a dimensionless factor that only depends on the relative mass ratios. With this definition of $r_{\rm CI}$, one can now write the cross section for a CI as,

[TABLE]

From this it follows that the cross section for outcome $X_{\rm ij}$ can be written as the following product,

[TABLE]

where $P(X_{\rm ij}|{\rm CI})$ is the probability for $X_{\rm ij}$ to be an outcome given the interaction is a CI. As seen, in this notation the probability $P(X_{\rm ij}|{\rm CI})$ stands for the factor $N_{X_{\rm ij}}/N_{\rm tot}$ in equation 1. In the following section, we proceed to calculate the associated formation rate.

2.1.2 Formation Rate

The derived cross section in Section 2.1.1 is especially useful for estimating the number of outcomes $X_{\rm ij}$ forming per time interval, referred to as the rate $\Gamma_{X_{\rm ij}}$. Following Samsing et al. (2014), this rate can be written for a given mass hierarchy as,

[TABLE]

Here $\Gamma_{X_{\rm ij}}^{(V)}$ denotes the number of outcomes $X_{\rm ij}$ forming per unit time in a steady-state isotropic $N$-body system with volume $V$ and characteristic velocity dispersion $v_{\rm\infty}$, $N_{\rm bin}$ and $N_{\rm sin}$ are the total number of binaries and singles in the system, respectively, and $P_{\rm bin}(a_{0})$ is the normalized SMA probability distribution. We will later derive formation rates of inspirals and collisions for a few relevant distributions of $P_{\rm bin}(a_{0})$.

2.2. Formation of Inspirals

We now derive analytical expressions for the cross section and rate of inspirals forming in binary-single interactions. For this, we first have to calculate the inspiral probability term $P(I_{\rm ij}|{\rm CI})$, as indicated by Equation (6). The first step in calculating $P(I_{\rm ij}|{\rm CI})$ is to make use of the observation that RIs generally can be described as a series of states characterized by a binary with a bound single (Samsing et al., 2014, 2016). We refer to each of these states as an intermediate state (IMS), and the binary in the state as an IMS binary. Between each IMS the three objects undergo a strong triple interaction, in which they semi-randomly exchange energy and orbital momentum. As a result, a newly formed IMS binary will always have a SMA and an eccentricity that is different from that of the initial target binary. Highly eccentric IMS binaries can therefore form during a RI, even when the target binary is circular. Now, if the IMS binary eccentricity is high enough, or correspondingly if the pericenter distance is small enough, the IMS binary will be able to undergo an inspiral through the loss of orbital energy, while remaining bound to the single in the resonance. In the next section we calculate what the initial orbital parameters of the IMS binary must be for it to undergo such an inspiral.

2.2.1 The Inspiral Boundary

A given IMS can be described as an IMS binary $[i,j]$ with a bound single $k$, where $\{i,j,k\}$ can be any combination of the three objects $\{1,2,3\}$. The IMS binary is formed with a SMA $a$ and eccentricity $e$, according to some distribution. The first question is what combinations of $a,e$ will result in an IMS binary $[i,j]$ undergoing a successful inspiral during the resonance while $k$ is still bound to the system.

For determining the orbital evolution of an IMS binary including orbital energy losses, we assume a model where the binary loses a constant amount of orbital energy, $\Delta{E}_{\rm p}$, during each pericenter passage (Samsing et al., 2016). As a result, in the orbit averaged approximation, the change in orbital energy per unit time can be written as,

[TABLE]

where $m_{\rm ij}=m_{\rm i}+m_{\rm j}$, $\mu_{\rm ij}=m_{\rm i}m_{\rm j}/m_{\rm ij}$, and $E_{\rm orb}$ ($T_{\rm orb}$) denotes the orbital energy (time) of the IMS binary. The solution to this equation leads to a corresponding inspiral time, $t_{\rm insp}$, which we here define as the time it takes for the IMS binary to inspiral from its initial SMA $=a$ to SMA $=0$. The solution for $t_{\rm insp}$ is trivially found by integration of Equation (8),

[TABLE]

To keep our formalism as general as possible, we now assume that $\Delta{E}_{\rm p}$ takes the following generic form,

[TABLE]

where $\mathscr{E}$ is a dimensionless normalization factor, $M$ is the characteristic mass scale, $\mathscr{R}$ is the characteristic length, and $r_{\rm p}$ is the IMS binary pericenter distance at the time of formation. As will be argued, this form can be used to describe orbital energy losses from both tides and GW emission (e.g. Samsing et al., 2016). Using this generic notation one can write the inspiral time as,

[TABLE]

Now, for the IMS binary to uninterrupted undergo an inspiral, its inspiral time, $t_{\rm insp}$, must be shorter than the time it is isolated from the bound single, $t_{\rm iso}$. This isolation time $t_{\rm iso}$ simply equals the Keplerian orbital time of the bound single with respect to the COM of the IMS binary,

[TABLE]

where $a_{\rm bs}$ is the semimajor axis of the single with respect to the COM of the IMS binary. The SMA $a_{\rm bs}$ can be found from energy conservation by assuming the triple system did not lose orbital energy before the formation of the IMS binary. From this assumption it follows that

[TABLE]

from which we find

[TABLE]

where,

[TABLE]

Inserting the relation for $a_{\rm bs}$ into Equation (12), the isolation time can now be expressed as,

[TABLE]

By equating $t_{\rm insp}$ and $t_{\rm iso}$ given by Equation (11) and (16), respectively, we can then find the following conditional relation between the IMS binary orbital parameters $a,e$,

[TABLE]

where

[TABLE]

[TABLE]

$e_{I_{\rm ij}}$ is the eccentricity for which $t_{\rm insp}=t_{\rm iso}$ for a given $a^{\prime}$, and $\mathscr{M}$ is a dimensionless term given by,

[TABLE]

The derived relation between $e_{I_{\rm ij}}$ and $a^{\prime}$ given by Equation (17) defines a boundary in orbital phase space which encloses the region in which inspirals can form, as illustrated in Figure 2. This inspiral region, which for the IMS binary in question we denote $\mathcal{R}_{\rm I}$, is bounded along the $a^{\prime}$-axis by a lower limit $a^{\prime}_{\rm l}$, and an upper limit $a^{\prime}_{\rm u}$. The value111The lower limit should in fact be marginally higher than $1$, since $a^{\prime}_{\rm l}=1$ corresponds to an infinite interaction time. However, uncertainties in how to define the upper limit makes this correction irrelevant for the scales we consider in this work. of $a^{\prime}_{\rm l}$ is $\approx 1$, where the value for $a^{\prime}_{\rm u}$ relates to the limit for when the resonant triple system no longer can be described as an IMS. That is, when $a_{\rm bs}\approx a_{0}$. One way of estimating $a^{\prime}_{\rm u}$ is by considering the force ratio $f_{\rm tid}=F_{\rm tid}/F_{\rm bin}$, where $F_{\rm tid}$ is the tidal force the bound single exerts on the IMS binary, and $F_{\rm bin}$ is the IMS binary’s binding force (Fregeau et al., 2004),

[TABLE]

By the use of Equations (14) and (15), one can solve for $a^{\prime}$ as a function of the mass hierarchy and the value of $f_{\rm tid}$. In this work we take $a^{\prime}_{\rm u}$ to be the value of $a^{\prime}$ for which222Our results depend only weakly on the precise value as $a^{\prime}_{\rm u}-1\propto f_{\rm tid}^{1/3}$. $f_{\rm tid}=0.5$. To summarize, the two boundaries $a^{\prime}_{\rm l}$ and $a^{\prime}_{\rm u}$ approximately take the values,

[TABLE]

In the following section we describe how to relate the inspiral region $\mathcal{R}_{\rm I}$ derived here to the probability term $P(I_{\rm ij}|{\rm CI})$.

2.2.2 Inspiral Probability

The probability $P(I_{\rm ij}|{\rm CI})$ can be written as

[TABLE]

where $\langle N_{\rm IMS}\rangle$ denotes the average number of IMS binaries $[i,j]$ formed in a CI, and $P(\mathcal{R}_{\rm I})$ is the probability for a single IMS binary $[i,j]$ to be formed with [$a^{\prime},e$] within the inspiral region $\mathcal{R}_{\rm I}$. We have here assumed that $P(\mathcal{R}_{\rm I})\ll 1$, and that the sampling of IMS binaries is not correlated.

To proceed we assume that the [$a^{\prime},e$] distribution is uniform at high eccentricity. This essential assumption was first shown to hold for the equal mass case by Samsing et al. (2014), however, as we will show, we find that the distribution also seems to be approximately uniform for the unequal mass case. Figure 5 shows an example of this. From this assumption it follows that the probability $P(\mathcal{R}_{\rm I})$ can be factorized as

[TABLE]

where $\mathcal{R}_{\rm U}$ denotes the region where the [$a^{\prime},e$] sampling is assumed uniform, $P(\mathcal{R}_{\rm U})$ is the probability for the IMS binary to form inside $\mathcal{R}_{\rm U}$, and $P(\mathcal{R}_{\rm I}|\mathcal{R}_{\rm U})$ is the probability for the IMS binary to form inside $\mathcal{R}_{\rm I}$ given that it was formed within $\mathcal{R}_{\rm U}$. This implies that the probability $P(I_{\rm ij}|{\rm CI})$ can be written as

[TABLE]

where $A_{\rm I}$ and $A_{\rm U}$ are the areas of the $\mathcal{R}_{\rm I}$ and $\mathcal{R}_{\rm U}$ regions, respectively. In what follows we define

[TABLE]

and note that the value of $\mathscr{N}$ do not depend to leading order on any finite size effects, orbital energy losses, the absolute mass scale or the initial SMA $a_{\rm 0}$, as further discussed in Section 2.5. Any effects related to orbital energy losses are therefore fully encoded in the area factor $A_{\rm I}$. This factor is given by the integral of $\epsilon_{I_{\rm ij}}$ over $a^{\prime}$, from $a^{\prime}_{\rm l}$ to $a^{\prime}_{\rm u}$,

[TABLE]

where $\mathscr{I}$ denotes the value of the integral,

[TABLE]

This integral has no analytical solution and must be evaluated numerically. We nonetheless find that the solution can be accurately approximated by the following functional form

[TABLE]

A few examples of $\mathscr{I}$ as a function of $a^{\prime}_{\rm u}$ and $\beta$ are shown in Figure 3. By substituting the expression for $A_{\rm I}$, given by Equation (28), into the relation for $P(I_{\rm ij}|{\rm CI})$, given by Equation (26), we find

[TABLE]

We can now make use of this relation to estimate the inspiral cross section and corresponding formation rate.

2.2.3 Inspiral Cross Section

The cross section for object pair [$i,j$] to undergo an inspiral, here referred to as the inspiral cross section $\sigma_{I_{\rm ij}}$, is given by the product of the inspiral probability and the CI cross section, $P(I_{\rm ij}|{\rm CI})\times\sigma_{\rm CI}$, as denoted in Equation (6). Making use of Equation (31) and Equation (5), we can then find

[TABLE]

where we have introduced the cross section term $\sigma_{\rm\mathscr{R}_{\rm ij}}$, and have defined the following quantities

[TABLE]

The term $\sigma_{\rm\mathscr{R}_{ij}}$ equals the cross section for object pair [$i,j$] to pass each other within a fixed distance $\mathscr{R}$ in the limit where no energy loss terms are included in the EOM. As later described in Section 2.3, the value for $\sigma_{\rm\mathscr{R}_{ij}}$ can be derived using Equation (41), and is therefore proportional to the collision cross section. As a result, the terms in the square parenthesis in Equation (32) represent the leading order effect from introducing energy loss terms into the finite-size, unequal-mass binary-single problem.

2.2.4 Inspiral Rate

Having derived an analytical solution to the inspiral cross section, we can then proceed to calculate expressions for the corresponding formation rates. Here we work out rates for two representative cases.

Field Interactions. Binaries in the field are often found to be uniformly distributed in $\log(a_{0})$ (e.g. Chanamé & Gould, 2004; Lépine & Bongiorno, 2007), a result also known as Oepiks law. For target binaries following such a distribution between SMA limits $a_{\rm min}$ and $a_{\rm max}$, the resultant rate of inspirals can be found by direct integration of Equation (7) using the cross section solution derived in Equation (32),

[TABLE]

In systems where the SMA distribution is unknown, one often assumes this distribution for simplicity.

Cluster Interactions. Binaries in dynamical systems, such as GCs, have at late times a SMA distribution which can be described by a Gaussian function in $\log(a_{0})$, as recently pointed out by Rodriguez et al. (2016a). In this case, the inspiral rate can be written as

[TABLE]

where $a_{\rm c}$ is the SMA of the Gaussian peak, and $s$ is the standard deviation. As can be seen, $s$ must be larger than $\beta$ for the assumption of a Gaussian distribution to modify the simple estimate that assumes that all target binaries have the same initial SMA $a_{0}=a_{\rm c}$. As $\log(\sigma_{\rm I})\propto(1/\beta)\log(a_{0})$, the derived rates are only sensitive to the shape of the Gaussian distribution when the inspiral cross section changes significantly over its intrinsic width. We note that since $\beta>3.5$ for both tides and GW emission (Samsing et al., 2016), this condition is in fact rarely met in a typical GC system where $s\approx 0.5$, as inferred from the distributions shown in Rodriguez et al. (2016a).

2.3. Collisions

We now turn our attention to collisions. A newly formed IMS binary [$i,j$] will undergo a collision if its pericenter distance $r_{\rm p}$ is smaller than a characteristic collisional distance $R_{\rm C}$. The maximum pericenter distance from which a collision can happen is therefore per construction given by $r_{\rm p}={R_{\rm C}}$. This limit is met when the orbital parameters of the IMS binary are such that,

[TABLE]

from which it follows that

[TABLE]

where

[TABLE]

and $e_{C_{\rm ij}}$ is the eccentricity that satisfies Equation (38) for a given $a^{\prime}$. The relation between $\epsilon_{C_{\rm ij}}$ and $a^{\prime}$ defines a boundary in orbital phase space, which encloses the region in which IMS binaries will undergo a collision ($r_{\rm p}<{R_{\rm C}}$). The area of this region, denoted $A_{\rm C}$, is given by

[TABLE]

Following the same procedure as we did for inspirals, the probability for a collision outcome ${C_{\rm ij}}$ provided that the interaction is a CI can now be factorized as $P({C}_{\rm ij}|{\rm CI})\approx A_{\rm C}\mathscr{N}$. It then follows that the collision cross section can be written as

[TABLE]

where we have defined $\mathscr{F}\equiv\mathscr{C}\times\mathscr{N}$, a factor that will be discussed in more detail in Section 2.5. The above equation for $\sigma_{{C}_{\rm ij}}$ is similar to Equation (32) but evaluated at $R_{\rm C}=\mathscr{R}$. In fact, our derived form for $\sigma_{{C}_{\rm ij}}$ can be used to estimate the cross section for any fixed distance close encounter, such as, e.g., a tidal disruption event. It is here worth noting that the collision cross section is independent of the initial SMA $a_{0}$, and that the inspiral cross section approaches the collision cross section for $\beta\rightarrow\infty$.

2.3.1 Collision Rate

The rate of collisions can be estimated by integrating Equation (7) using the collision cross section in Equation (41). However, in this case, the solution is particularly simple as the collision cross section does not depend on $a_{0}$. This implies that the corresponding collision rate $\Gamma_{C_{\rm ij}}$ is independent of the SMA distribution. As a result, the rate of collisions takes the following simple form

[TABLE]

2.4. Inspirals Relative to Collisions

With our derived inspiral and collision cross sections, we are now in a position to understand their relative formation rates. Assuming the inspiraling binaries also merge, it is of particular interest to explore whether or not inspirals are expected to significantly contribute to the merger rate. We do so by considering the cross section ratio between inspirals and collisions

[TABLE]

We first note that the factor $\mathscr{F}$ cancels out, which allows us to derive an analytical solution. This is especially important for studying the role of tides, as the tidal inspiral cross section and the collision cross section are often comparable (Samsing et al., 2016).

The above expression shows that the number of inspirals generally increases relative to that of collisions with increasing $a_{\rm 0}$ and decreasing $R_{\rm C}$. The effect of orbital energy losses is thus largest in interactions involving wide binaries and dense objects such as WDs, NSs, and BHs. A similar conclusion was reached for the equal mass case in Samsing et al. (2014, 2016). This counter intuitive scaling also explains why earlier studies on, e.g., tidal interactions not have been able to resolve these notable differences between collisions and inspirals (see e.g. Gaburov et al., 2010). The relevance and applicability of the general analytical formalism derived here is discussed below.

2.5. The Validity of Our Analytical Formalism

Due to the highly chaotic nature of the three-body problem in the HB limit, a simple analytical description, similar to the one that exists in the SB limit (Hut, 1983), is not easily derived. Here we have outlined a simple approach that allows one to calculate the cross section for a variety of outcomes that arise when finite size effects and orbital energy loss term corrections are taken into account. Our solutions make use of the fact that the [$a,e$] distribution is uniform at high eccentricity, which has been validated with the use of numerical experiments (Samsing et al., 2016).

We also note that the scalings in Equation (32) can be derived, as argued in Samsing et al. (2016), by considering the physical attributes of the interaction. For two objects to undergo an inspiral, the orbital energy loss during the first passage, $\Delta{E_{\rm p}}$, must be a notable fraction of the initial orbital energy, $E_{0}$. By equating $\Delta{E_{\rm p}}$ and $E_{0}$ and solving for the corresponding pericenter distance, which we denoted here as $r_{\rm cap}$, one finds

[TABLE]

If the inspiral leads to a merger, one can think of $r_{\rm cap}$ as an effective stellar size. As seen in Equation (44), $r_{\rm cap}$ is not just proportional to the size of the object, but scales with $(a_{0}/\mathscr{R})^{1/\beta}$. This scaling is similar to the one found in Equation (32), giving further credence to our analytical formalism.

We note that we are yet to calculate a precise expression for $\mathscr{F}$, which relates to the difficulties associated with calculating the full IMS binary [$a^{\prime},e$] distribution. However, due to the classical scale-free nature of the three-body problem (Hut & Bahcall, 1983), one can show that $\mathscr{F}$ depends only on the relative mass ratio and not on the absolute mass scale, initial SMA, finite sizes or the precise orbital energy loss mechanism. The results reported in this paper have all been calibrated using numerical experiments, from which we have found, for example, that $\mathscr{F}\approx 6.0$ for the equal mass case. When comparing the relative rate of inspirals to collisions $\mathscr{F}$ cancels out, which allows for a full analytical estimate to be derived.

In general we find that the set of scalings and analytical results derived here agree very well with results from full $N$-body simulations (e.g. Samsing et al., 2014, 2016). In the following sections we will use our analytical model to gain further insight into the dynamical and observational importance of GW inspirals.

3. Gravitational Wave Inspirals

When GR effects are included in the binary-single problem, a close passage between any two of the three objects will lead to GW emission. In this process orbital energy is carried out of the system, and can as a result lead to an inspiral and subsequent merger between the two objects (Gültekin et al., 2006; Samsing et al., 2014, 2016). In the remaining parts of this paper we will discuss several aspects related to these GW mergers that we generally refer to as GW inspirals. For this, we will make use of our analytical framework from Section 2, as well as full numerical $N$-body simulations including PN correction terms.

From an astrophysical perspective, GW inspirals are extremely exciting, as they often enter the LIGO band with non-zero eccentricity (Gültekin et al., 2006; Samsing et al., 2014). In fact, recent studies indicate that they might even dominate the population of observable high eccentricity BBH mergers (Samsing & Ramirez-Ruiz, 2017). Together with the spin distribution (e.g. Rodriguez et al., 2016c; Kushnir et al., 2016; Zaldarriaga et al., 2017), the eccentricity distribution is expected to play a key role in constraining the origin of BBH mergers observable by LIGO and next generation GW observatories (e.g. Chen & Amaro-Seoane, 2017). High eccentricity NS-NS mergers are also promising probes of the NS equation of state (e.g. East & Pretorius, 2012; Lehner & Pretorius, 2014; East et al., 2015; Paschalidis et al., 2015; East et al., 2016a, b).

To provide concrete predictions and scaling solutions relevant for LIGO, we first focus here on the case where all three interacting objects are equal mass BHs. Although this represents an idealized scenario, recent numerical studies do in fact indicate that a significant fraction of the observable population of BBH mergers forming in GCs likely form dynamically through similar mass binary-single interactions (Rodriguez et al., 2016a; Samsing & Ramirez-Ruiz, 2017). The explanation for this is that mass segregation causes heavier BHs to reach the GC center faster than lighter ones, BHs of similar mass are thus more likely to be found in the GC center at similar ages in the evolution of the cluster (e.g. Rodriguez et al., 2016a; Park et al., 2017). Subsequent binary-single interactions also tend to keep BHs of similar mass together (Sigurdsson & Phinney, 1993).

In the sections that follow we start by introducing the GW energy loss term. We then explore how frequent GW inspirals are compared to collisions in encounters involving compact objects. After this we derive a set of solutions to the equal mass case. This is followed by a study on how slight changes to the mass hierarchy affect the GW inspiral cross section. In the last section we use our solutions to estimate the present day GW inspiral rate observable by LIGO.

3.1. GW Energy Loss Term

The energy emitted through GW radiation during a single pericenter passage between object pair $[i,j]$ is at quadrupole order in the high eccentricity limit given by (Peters, 1964; Hansen, 1972),

[TABLE]

In this limit, one finds that $\Delta{E}_{\rm GW}$ can also be expressed by the general form given by Equation (10) by setting,

[TABLE]

With these substitutions the derived relations from Section 2 can now be used to study GW inspirals. Below we proceed by stuyding how frequent GW inspirals are compared to standard collisions.

3.2. GW Inspirals Relative to BH and NS Collisions

The final outcome of both a GW inspiral and a collision is a merger. In this section we explore which of the two channels is expected to dominate the total merger rate. As an example, in main sequence star interactions the collision cross section is typically higher than the tidal inspiral cross section, where for interactions involving high mass WDs the two cross sections are of similar order (Samsing et al., 2016). Here we are interested in calculating the ratio of GW inspirals to collisions for interactions involving BHs and NSs, and how sensitively it depends on the initial SMA, $a_{0}$, and the mass hierarchy. For this study, we use our analytical estimate for ${\sigma_{I_{\rm ij}}}/{\sigma_{C_{\rm ij}}}$ given by Equation (43). Two examples are worked out below.

3.2.1 BH-BH Mergers

We start by considering the following triple BH binary-single interaction,

[TABLE]

where the corresponding masses are $m_{1}=30M_{\odot}$, $m_{2}=20M_{\odot}$, and $m_{3}$ is varied between $10M_{\odot}-40M_{\odot}$. The $\log$ of the cross section ratio of GW inspirals to collisions for object pair $[1,3]$, $\log({\sigma_{I_{13}}}/{\sigma_{C_{13}}})$, as a function of the initial SMA $a_{0}$ and $m_{3}$, is shown in the top panel of Figure 4. As can be seen, the GW inspiral cross section is, for this example, $\approx 10^{2}$ times higher than the collision cross section.

3.2.2 NS-NS Mergers

We now consider a binary-single interaction involving two NSs, and an unspecified binary companion (BC),

[TABLE]

where the corresponding masses are $m_{2}=1.4M_{\odot}$, $m_{3}=1.4M_{\odot}$, and $m_{1}$ is varied between $0.5M_{\odot}-3M_{\odot}$. Our chosen range of $m_{1}$ covers a few astrophysical interesting configurations, from BC1 representing a WD (0.5$M_{\odot}$) to a heavy NS (3$M_{\odot}$). The log of the cross section ratio of [NS, NS] GW inspirals to collisions, $\log({\sigma_{I_{23}}}/{\sigma_{C_{23}}})$, as a function of the initial SMA $a_{0}$ and $m_{1}$, is shown in the bottom panel of Figure 4. In this case, even for a NS, which has a slightly larger physical radius than a BH, GW inspirals clearly dominate over classical collisions.

3.2.3 GW Inspirals versus Collisions for Compact Objects

The number of compact object mergers that dynamically form during resonant three-body interactions is completely dominated by GW inspirals, and not by classical collisions. We find this to be true for both BH and NS mergers. This result is not that surprising when considering Equation (43), which says that the cross section ratio approximately equals the initial SMA $a_{0}$ over the Schwarzschild radius of the two merging BHs to the power of $2/7$. This ratio is a large number for typical astrophysical systems, where $a_{0}$ is about 1 AU ($\approx 10^{8}$ km) and the BH Schwarzschild radius is $\approx 10^{2}$ km. It is thus crucial to have PN terms included in the $N$-body code for estimating a meaningful rate of compact object mergers forming during three-body interactions. As an example, Antonini et al. (2016) did not have PN terms in their few-body code, which likely led to an underestimation of their in-cluster merger rate by a factor of $\approx 10^{1}-10^{2}$ (see Table 1 in Antonini et al. (2016)). Because of this, their rate reported for high eccentricity BBH mergers forming through binary-single interactions is far too low, as we argued in Samsing & Ramirez-Ruiz (2017).

For this reason, in the rest of the paper, we will mainly focus on GW inspirals without making any further inferences for collisions. In fact, the collision population is included in our analytical estimation of GW inspirals, as the inspiral area $A_{\rm I}$ overlaps with the collision area $A_{\rm C}$ in orbital phase space (e.g. Samsing et al., 2016).

3.3. Equal Mass Interactions

In this section we calculate the GW inspiral cross section assuming all three objects have the same mass $m$. As argued earlier, this is an idealized case, but is still likely to provide a reasonable description of the dynamics leading to the majority of BBH mergers forming in GCs observable by LIGO. We note that the equal mass case was also studied in Samsing et al. (2014, 2016), we therefore keep the following section concise. The results presented here will be used in later sections for estimating absolute and relative BBH merger rates.

The inspiral boundary and corresponding inspiral cross section can be calculated using Equations (17) and (32), respectively. By further making use of Equation (46), we find that the GW inspiral boundary can be written as,

[TABLE]

and the corresponding GW inspiral cross section as,

[TABLE]

where $C_{\rm GW}$ is a constant given by,

[TABLE]

The GW inspiral cross section summed over all three possible inspiral pairs, denoted here by ${\sigma}_{I}$, can be written in the more familiar astrophysical units as,

[TABLE]

where we have made use of Equation (50) with $\mathscr{F}=6$ (see Section 2.5).

It is worth noting that the GW inspiral cross section does not scale linearly with mass $m$, as the classical Newtonian outcomes including collisions, but instead as $m^{12/7}$. As a result, the GW cross section for, say, three $30M_{\odot}$ BHs is $\approx 200$ times larger than that for three $1.4M_{\odot}$ NSs, and not just by a factor $\approx 20$ as inferred from gravitational focusing. In the following section we study how slight changes to the mass hierarchy introduces corrections to the equal mass case solution derived here.

3.4. Unequal Mass Interactions

To explore how sensitive our derived GW inspiral cross section is to variations in the mass ratio, we here consider the following unequal mass example,

[TABLE]

where the corresponding masses are denoted by $m_{\rm A}$ and $m_{\rm B}$, respectively. For this interaction we study the formation of [BHB, BHB] GW inspirals. To facilitate comparison, we refer to this case as U (Unequal), and E (Equal) to the case when all three objects have $m=m_{\rm B}$. An example illustrating the formation of a GW inspiral for scenario U, where $m_{\rm A}=10M_{\odot}$ and $m_{\rm B}=20M_{\odot}$, is shown in Figure 1. Below we derive the inspiral boundary and corresponding cross section of case U relative to E, as a function of $m_{\rm A}$.

3.4.1 Dependence on Mass Ratio $m_{\rm A}/m_{\rm B}$

Making use of Equations (17) and (49), the GW inspiral boundary of scenario U relative to E, takes the following form

[TABLE]

where we have introduced the mass ratio $q$ given by

[TABLE]

By the same token, the corresponding GW inspiral cross section for scenario U relative to E can be written as

[TABLE]

where we have used Equations (32) and (50). As can be seen, the above fraction is composed of three different terms. The first term (sampling term) mainly reflects differences in the $[a^{\prime},e]$ distribution and corresponding sampling frequency of IMS binaries, where the second term (kinematic term) relates to changes in the inspiral time and isolation time. The third term arises from the change in gravitational focusing as the total mass is varied.

To explore how sensitive our derived cross section ratio in Equation (56) is to changes in $q$, we now expand it to linear order in $\delta$ where $q=1+\delta$,

[TABLE]

As discussed in Section 2.5, we have not yet been able to write out a precise form for the sampling term and because of this, its derivative has been written explicitly in the above expression. We do have an expression for $\mathscr{I}$, however, when considering variations in $q$ it has to be consistently paired with $\mathscr{F}$. Focusing on the remaining terms, we see that small changes in mass ratio result in small changes of order $(30/21)\delta$. In fact, the resultant change is likely to be even smaller, as our numerical simulations indicate that the sampling term decreases with increasing $\delta$. This is explained by the fact that the heavier objects are generally more prone to form binaries than lighter ones (e.g. Sigurdsson & Phinney, 1993). As a result, the two terms in the square parentheses from Equation (57) almost cancel out, which leads to small differences between scenarios U and E. We show that this is indeed the case in the numerical example explored in the following section.

3.4.2 Formation of GW Inspirals from $m_{\rm A}=10M_{\rm\odot}$ and $m_{\rm B}=20M_{\rm\odot}$

To study a concrete example we here explore a few aspects related to the dynamical formation of [BH(20$M_{\rm\odot}$), BH(20$M_{\rm\odot}$)] GW inspirals for $m_{\rm A}=10M_{\rm\odot}$ and $m_{\rm B}=20M_{\rm\odot}$. We note that $q=0.5$ is actually a rather ‘extreme’ case from a perspective of BBH mergers forming in GCs, as recent simulations indicate that the median mass ratio for the merging population is about $0.9$ (Rodriguez et al., 2016a).

We first consider the [$a^{\prime},e$] distribution of all [BH(20$M_{\rm\odot}$), BH(20$M_{\rm\odot}$)] IMS binaries. Results are shown in Figure 5. As seen, despite the often enormously complex pathway from initial interaction to final GW inspiral (See e.g. Figure 1), we do find excellent agreement between our full $N$-body simulations (large black dots) and our analytical solution for unequal mass given in Equation (54) (black solid line/grey area). Figure 5 also shows that the [$a^{\prime},e$] distribution is indeed approximately uniform at high eccentricity, as we assumed in Section 2.2.2.

The associated GW inspiral cross section is shown in Figure 6 as a function of $a_{0}$. As expected, our generic prediction from Equation (32), which states that the inspiral cross section always scales $\propto a_{\rm 0}^{2/7}$ (for $\beta=7/2$) in the asymptotic limit independently of the mass hierarchy, clearly seems to hold. As a consequence, the GW inspiral cross section increases here by more than an order of magnitude across the considered interval, and will keep increasing until reaching the SB limit, at which it will sharply drop off (Samsing et al., 2014). Considering the normalization, we see that the equal mass solution from Equation (52) with $m=20M_{\odot}$ provides a rather accurate estimate (the equal mass case gives $\approx 1.5\ \text{AU}^{2}$ at $a_{0}=1$ AU), which agrees with that describing slightly unequal mass interactions often can be done assuming the equal mass limit.

3.5. Formation Rate of GW Inspirals

Having derived the cross section of GW inspirals, we are now in a position to estimate the corresponding rate. For this, we consider the population of GW inspirals forming in the cores of GCs. As the BH mass ratios are likely to be near unity at late times (Rodriguez et al., 2016a), we approximate individual scatterings by the equal mass limit – an assumption that also was shown in Samsing & Ramirez-Ruiz (2017) to accurately reproduce several observables related to BBH mergers. We take each GC to have a velocity dispersion $v_{\rm\infty}=10$ km s*-1*, a core volume $V=(0.1\text{pc})^{3}$, a total number of BHs in the core $N_{\rm BH}$ with mass $m_{\rm BH}$, and a corresponding BH binary fraction of $0.5$. The SMA $a_{0}$ distribution is assumed to follow a Gaussian in $\log(a_{0})$ with center value $a_{\rm c}$ and standard deviation $s=0.5$, as has recently been shown to provide a reasonable fit to numerical simulations (Rodriguez et al., 2016a). Adopting a GC number density of $1$ GC Mpc*-3* from (Rodriguez et al., 2015) and making use of Equation (36), we find

[TABLE]

where ${\Gamma_{I}}$ here denotes the rate summed over all three possible inspiral pairs. We thus conclude that if each GC has about 80 interacting BHs, then about $1$ GW inspiral could be observed by LIGO per year. As GW inspirals are likely to enter the LIGO band with non-zero eccentricity (e.g. Samsing & Ramirez-Ruiz, 2017), we note that the possibilities for detecting GW inspirals depend highly on constructing accurate GW templates for varying eccentricity and spin (e.g. Harry et al., 2016; Huerta et al., 2016). Although our formalism is well defined, our estimated for the number of GW inspirals is highly uncertain. This is not a unique problem related to our model, as nearly all proposed BBH merger channels are plagued by the same large uncertainties related to the dynamics and BH demographics of GCs (e.g. Chatterjee et al., 2017).

4. Post-Interaction Mergers and GW Inspirals

The vast majority of binary-single interactions end as a binary with an unbound single, also known as a fly-by or an exchange outcome (e.g. Hut & Bahcall, 1983). Binaries formed in this way, here referred to as post-interaction binaries, will eventually merge through GW emission, however, the merger delay time distribution is generally extremely broad with a tail that often exceeds the Hubble time. This is in stark contrast to GW inspirals that are rare and merge relative promptly.

In this section we first calculate the cross section of such post-interaction binary mergers. We then derive an expression for the cross section ratio between GW inspirals and post-interaction binary mergers. This leads us to the interesting conclusion that the cross section of GW inspirals is always $>1\%$ of the post-interaction binary merger cross section. The corresponding implications for BBH formation and relative rates are studied in Section 5.

4.1. Post-Interaction Binary Merger Cross Section

We start by calculating the cross section for a post-interaction binary $[i,j]$ to merge within a time $\tau$. This cross section, denoted here by $\sigma_{M_{\rm ij}}^{<\tau}$, can be factorized as,

[TABLE]

where $BS_{\rm ij}$ denotes an end-state composed of binary $[i,j]$ and an unbound single, $P(BS_{\rm ij})$ is the probability for outcome $BS_{\rm ij}$, and $P(t_{M_{\rm ij}}<\tau)$ is the probability for that a newly formed $BS_{\rm ij}$ binary has a merger time $t_{M_{\rm ij}}$ that is less than $\tau$. As we will show, the cross section $\sigma_{M_{\rm ij}}^{<\tau}$ first increases with $a_{0}$ until $a_{0}$ approaches a characteristic value $a_{\rm 0}^{<\tau}$, above which $\sigma_{M_{\rm ij}}^{<\tau}$ begins to decrease with $a_{0}$. The transitional SMA $a_{\rm 0}^{<\tau}$ is to leading order the $a_{0}$ for which all corresponding $BS_{\rm ij}$ binaries have $t_{M_{\rm ij}}=\tau$ for eccentricity $e=0$. For $a_{\rm 0}<a_{\rm 0}^{<\tau}$, the cross section $\sigma_{M_{\rm ij}}^{<\tau}$ is therefore simply given by $P(BS_{\rm ij})\times\sigma_{\rm CI}\propto a_{0}$. Our numerical scattering results shown in Figure 6 clearly illustrate this piecewise behavior with $a_{0}$. In the following we estimate $\sigma_{M_{\rm ij}}^{<\tau}$ for $a_{\rm 0}>a_{\rm 0}^{<\tau}$.

For a newly formed $BS_{\rm ij}$ binary to merge within a time $\tau$ for a given $a_{\rm 0}>a_{\rm 0}^{<\tau}$, its eccentricity $e$ must be larger than some characteristic value $e_{\rm\tau}$. Assuming a thermal distribution for the end-state eccentricities (Heggie, 1975), the probability for a $BS_{\rm ij}$ binary to have $e>e_{\rm\tau}$ is simply $1-{e^{2}_{\rm\tau}}$. From which it follows,

[TABLE]

Using the relation between $1-e^{2}$ and the GW merger time of a high eccentricity binary (e.g. Peters, 1964), the cross section $\sigma_{M_{\rm ij}}^{<\tau}$ for $a_{\rm 0}>a_{\rm 0}^{<\tau}$ can now be written as,

[TABLE]

where

[TABLE]

and we have expressed $\sigma_{M_{\rm ij}}^{<\tau}$ in a form that is similar to the inspiral cross section given in Equation (32). Our prediction for $\sigma_{M_{\rm ij}}^{<\tau}$ is shown in Figure 6 as a dash-dotted line. As for the factor $\mathscr{F}$, the term $P(BS_{\rm ij})$ does not depend on either the absolute mass scale, the initial SMA $a_{\rm 0}$, nor on any finite size effects including orbital energy losses. However, a full analytical solution to $P(BS_{\rm ij})$ for a general mass hierarchy is not yet available.

4.1.1 Equal Mass Solution

In the limit where all three objects are compact and have the same mass $m$, we find that the post-interaction binary merger cross section summed over all three potential merger pairs, denoted here as $\sigma_{M}^{<\tau}$, can be written in piecewise form as

[TABLE]

Here the top solution is valid when $a_{0}<a_{\rm 0}^{<\tau}$, while the bottom solution is valid when $a_{0}>a_{\rm 0}^{<\tau}$. For this equal mass case $a_{\rm 0}^{<\tau}$ is approximately given by the SMA for which the initial target binary has a GW merger time that equals $\tau$, from which it follows that

[TABLE]

If we evaluate this for $m=20M_{\odot}$ and $\tau=t_{\rm H}$ we find $a_{\rm 0}^{<\tau}\approx 0.15$ AU. This partly explains why most of the BBH mergers formed in numerical simulation pile up near this SMA (Rodriguez et al., 2016a).

4.2. GW Inspirals Relative to Post-Interaction GW Mergers

The ratio between the GW inspiral cross section and the post-interaction GW merger cross section can be simply written in the equal mass case as

[TABLE]

where we have used Equations (52) and (62), respectively. The top expression here again applies for $a_{0}<a_{\rm 0}^{<\tau}$, while the bottom is for $a_{0}>a_{\rm 0}^{<\tau}$. The value of $\log({\sigma_{I}}/{\sigma_{M}^{<\tau}})$ is illustrated in Figure 7, as a function of initial SMA $a_{0}$ and mass $m$ for $\tau=t_{\rm H}$.

As given by Equation (64), and also seen in both Figures 6 and 7, the ratio ${\sigma_{I}}/{\sigma_{M}^{<\tau}}$ behaves opposite to the scaling of ${\sigma_{M}^{<\tau}}$, as it first decreases with $a_{0}$ until $a_{0}\approx a_{\rm 0}^{<\tau}$, after which it increases with $a_{0}$. This piecewise scaling interestingly implies that the ratio ${\sigma_{I}}/{\sigma_{M}^{<\tau}}$ takes its minimum value at which the initial SMA $a_{0}\approx a_{\rm 0}^{<\tau}$. Evaluating the ratio at $a_{0}=a_{\rm 0}^{<\tau}$, therefore lead us to the following inequality

[TABLE]

which states that the GW inspiral cross section will always be at least $1\%$ of the post-interaction binary merger cross section for interacting compact objects of similar mass $m\gtrsim 10M_{\odot}$.

5. Formation of Binary Black Hole Mergers

Having developed a clear understanding of GW inspirals and post-interaction binary mergers, we are now finally in a position to study their relative contributions to the number of dynamically formed BBHs. Our goal in this section is to examine GW inspirals arising from dynamical interactions in a dense stellar system like a GC. Given the uncertainties surrounding the presence and role of BHs in the evolution of GCs and other dense stellar clusters, we restrict ourselves to two idealized scenarios. In the first Early-Burst Scenario, we study the formation of BBH mergers that form from a population of binaries and singles that interact only over a brief time interval. In the second Steady-State Scenario we study the formation of BBH mergers from an interacting, steady state population. Results from the two scenarios are given below.

5.1. Early-Burst Scenario

Let us first assume a population of binaries and singles undergoing interactions over a short time interval at some initial time $t_{0}$. We speculate that this could represent the dynamical environment of an early collapse associated with the formation of a dense stellar system, such as a GC (e.g. Mapelli et al., 2013; Mapelli & Zampieri, 2014; Ziosi et al., 2014; Spera et al., 2015; Kimpson et al., 2016). For this scenario we study the cumulative number of BBH mergers forming through the GW inspiral channel, $N_{I_{\rm ij}}$, relative to that of the post-interaction binary merger channel, $N_{M_{\rm ij}}$, as a function of time $\tau=t-t_{0}$. We note that only $N_{M_{\rm ij}}$ is time dependent because of its delay time distribution. For this scenario it is straight forward to show that

[TABLE]

The equal mass solution to the Early-Burst Scenario is therefore simply given by Equation (64). This implies that the inequality derived in Equation (65) can be applied to this case, from which we conclude that the cumulative number of GW inspirals in the equal mass limit for this scenario will always be at least $1\%$ of the total number of BBH mergers for $m\gtrsim 10M_{\odot}$.

The solution to Equation (66) also shows that the relative number of GW inspirals scales $\propto\tau^{-2/7}$, which interestingly suggests that the first BBH mergers in the Universe were likely to be GW inspirals. More generally, the solution implies that there exists a finite time interval within which the fraction of BBH mergers is dominated by GW inspirals and not by standard post-interaction binary mergers. In the equal mass case, we can solve for this characteristic time interval, denoted here as $\tau_{\rm IM}$, by setting ${N_{I}}={N_{M}}$, from which we find

[TABLE]

The normalization of the initial SMA $a_{0}$ is here set to $25\ \text{AU}$, which is $\approx 0.1$ times the HB value for $v_{\infty}=10\ \text{km\ s}^{-1}$ and $m=20M_{\odot}$. Despite large uncertainties in what the relevant values of $a_{0}$ and $m$ should on average be, this simple estimate indicates that a characteristic time interval of order $\approx 10^{6}\ \text{years}$ is not unrealistic. We note this is similar to the dynamical time of a typical GC.

These results indicates that GW inspirals might play a role in the formation of a dense stellar system at early times. If we denote the ratio between the number of GW inspirals to the number of binary-single interactions by $f_{\rm IB}$, one finds in the equal mass limit that

[TABLE]

This indicates that for $m\approx 20M_{\odot}$, the initial SMA must be $\lesssim 0.15$ AU for GW inspirals to form in at least $1\%$ of all binary-single interactions. GW inspirals are indeed more likely to contribute to the dynamical evolution of a dense stellar system than post-interaction binary mergers. This is because GW inspirals form and evolve in bound resonant states, which is in contrast to post-interaction binary mergers that always receive a dynamical kick prior to merger. This kick is often high enough to unbound them from the dense stellar system (e.g. Rodriguez et al., 2016a). This could have interesting dynamical consequences, and numerical simulations do in fact show an indication for this to be the case (Gültekin et al., 2004, 2006), yet, more work is needed before firm conclusions can be drawn.

5.2. Steady-State Scenario

As a second representative example, we study the formation of BBH mergers in a GC assuming its distribution of binaries and singles remains constant in time. We imagine this steady-state scenario to approximately describe the late time evolution and formation of BBH mergers in a typical GC. Below we first calculate the cumulative number of BBH mergers and then derive relative rates, from which we find that $\gtrsim 1\%$ of the present day rate of BBH mergers assembled through binary-single interactions is likely to originate from the GW inspiral channel. Finally, in the last section, we estimate the fraction of GW inspirals that is likely to appear in the LIGO band with a particular high eccentricity.

5.2.1 Cumulative Number of BBH Mergers

Here we derive the cumulative number of GW inspirals and post-interaction binary mergers forming as a function of time $\tau=t-t_{0}$, consistently taking into account ongoing interactions as well as post-interaction binary merger delay time distributions. Assuming $a_{0}>a_{\rm 0}^{<\tau}$, one can show that the time dependent cumulative number of post-interaction binary mergers for this scenario can be written as

[TABLE]

and the cumulative number of GW inspirals as,

[TABLE]

where we have assumed that the average inspiral time $\ll\tau$. The cumulative number of GW inspirals relative to that of the post-interaction binary mergers at time $\tau$, is thus given by

[TABLE]

We note that this estimate and the one derived for the early-burst scenario differ by only a factor of $9/7$. This indicates that the relative number of GW inspirals forming through the binary-single channel is not strongly dependent on the encounter history. In the equal mass case, it seems robust to conclude that $>1\%$ of all BBHs forming through the binary-single channel are likely to be in form of GW inspirals.

5.2.2 Relative Rate of GW inspirals

The rate of GW inspirals relative to that of the post-interaction binary mergers evaluated at time $t\gg t_{0}$, can be shown to take the following form

[TABLE]

Assuming the equal mass limit, we can combine this estimate with Equation (65) to obtain

[TABLE]

This again gives further support to our conclusion that the rate of GW inspirals will be $\gtrsim 1\%$ of the total rate for $m\gtrsim 10M_{\odot}$.

From an astrophysical perspective, GW inspirals are of particular interest because of their high eccentricity at formation. However, only a fraction of the GW inspirals will actually appear in the LIGO band with the same high eccentricity they had at formation as a result of circulation during inspiral (e.g. Samsing & Ramirez-Ruiz, 2017). In the section below we derive this fraction and its present day relative rate.

5.2.3 GW Inspirals Forming in the LIGO Band

One can show that for a BBH to appear in the LIGO band with its initial peak eccentricity, it must necessarily form within the band (e.g. Samsing et al., 2014). The fraction of GW inspirals at a high eccentricity is thus given to leading order by the fraction that forms within the observable LIGO band.

By the use of the GW peak frequency fitting formula from Wen (2003), one can solve for the required initial orbital parameters of an inspiraling BBH so that its GW peak frequency is above a certain threshold $f_{\rm GW}$. This provides us with the following relation

[TABLE]

where we have assumed the high eccentricity limit and $r_{f_{\rm GW}}$ is the pericenter distance at which the corresponding GW peak frequency equals $f_{\rm GW}$. If the BBH pericenter distance $r_{\rm p}<r_{f_{\rm GW}}$ then the GW peak frequency will be $\gtrsim f_{\rm GW}$. As $r_{f_{\rm GW}}$ is a fixed distance, we can now conclude that the cross section for a BBH to form with GW peak frequency $\gtrsim f_{\rm GW}$ is to leading order given by Equation (41) with ${R}_{\rm C}=r_{f_{\rm GW}}$, and not by the GW inspiral cross section. From this it follows directly that

[TABLE]

where ${\sigma}_{f_{\rm GW}}$ is the cross section for a BBH to form with an associated GW peak frequency $\gtrsim f_{\rm GW}$ after a binary-single interaction. Comparing this with the GW inspiral cross section from Equation (52), one finds the relative rate in the equal mass limit to be of order

[TABLE]

This illustrates that for the chosen normalizations only $\approx 20\%$ of the GW inspirals will have a GW peak frequency $f_{\rm GW}>10\text{Hz}$ at the time of their formation. In other words, only $\approx 20\%$ of the GW inspirals will appear in the LIGO band with their initial high eccentricity, whereas the remaining $\approx 80\%$ will undergo notable circularization before being detected. Although this seems slightly discouraging in terms of observable rates, we do note that the vast majority of GW inspirals still make it into the LIGO band at $10\text{Hz}$ with significant eccentricity (typically $>0.1$). This was shown by Samsing & Ramirez-Ruiz (2017), who highlighted that if the BBH merger rate is dominated by the GC population, then GW inspirals are likely to dominate the population of eccentric BBH mergers observable by LIGO.

6. Conclusions

The few-body problem with terms correcting for finite sizes, tides, and GR effects, is usually studied using $N$-body codes (e.g. Mardling & Aarseth, 2001; Samsing et al., 2016). However, these generally provide very limited physical insight into the problem. In this paper we have explored analytically the effect from including such correction terms in three-body interactions, which have provided us with a wealth of new insight into this problem. Our study focus on the population of two-body captures that form through tidal and GW energy losses in binary-single interactions – a population that we have been referring to as inspirals. In order to provide a clear physical framework, we have first developed an analytical formalism for calculating the inspiral cross section and corresponding rate. We then applied it to explore the formation of GW inspirals, as we found this population to be highly interesting both dynamically and observationally. Our of findings are summarized below.

6.1. Finite Sizes and Energy Loss Terms

As illustrated in Samsing et al. (2016), the main effect from including finite sizes and orbital energy loss term corrections in binary-single interactions, is an increase in the number of mergers, that either form through standard collisions (by the finite size term) or inspirals (by tides and GW emission). This increase is not related to an overall change in the dynamics, as the addition of the energy loss terms does not significantly modify the number of classical exchange and fly-by outcomes. Instead, the correction terms lead to occasional highly impulsive orbital energy losses that then result in either a collision or an inspiral. As the capture distance is always larger than the corresponding physical size, the inclusion of orbital energy loss terms will always lead to an increase in the number of mergers. Our formalism shows how this effective size, and thereby the number of mergers, varies with the initial orbital energy of the three-body system.

6.2. Inspiral Cross Section and Rates

We found that the cross section of inspirals increases independently of the mass hierarchy as $a_{0}^{1/\beta}$, when the considered orbital energy loss term scales with the pericenter distance as $r_{\rm p}^{-\beta}$. This provides an accurate description for GW emission, but it is only approximate for tides (e.g. Samsing et al., 2016). The collision cross section is found to be independent of $a_{0}$ and, as such, we concluded that the number of mergers relative to that of collisions always increases with $a_{0}$. The relative importance of inspirals is thus partly set by the HB limit; systems with low velocity dispersion are therefore likely to show the largest effect from orbital energy loss term corrections. The corresponding inspiral rate can be written in closed form for a few relevant cases as a result of our analytical solution. As an example, we showed that for a dynamical system where the target binaries follow a Gaussian in $\log(a_{0})$ with standard deviation $s$ and central value $a_{\rm c}$, the Gaussian assumption only plays a role when $s>\beta$. However, this condition is rarely met in typical GC systems333$s$ is often about $\approx 0.5$ as found by Rodriguez et al. (2016a), where $\beta$ is always $>3.5$ as discussed in Samsing et al. (2016), and the rate is, to leading order, simply proportional to the inspiral cross section evaluated at $a_{\rm c}$.

6.3. Formation and Importance of GW Inspirals

The inclusion of GW energy loss terms leads to the formation of GW inspirals. To study this population, we started by investigating the relative contribution of GW inspirals and collisions to the number of GW mergers formed during binary-single interactions. Using our analytical framework, we found that GW inspirals generally are $\approx 100$ times more likely to form than standard collisions for encounters involving BHs and NSs. A consistent inclusion of GW energy loss terms in $N$-body codes is therefore crucial for accurately estimating rates of compact object mergers. We note here that the public code Fewbody (Fregeau et al., 2004; Fregeau, 2012) does not include such terms in its original version, which indeed has led to a long list of studies that significantly underestimate the number of especially highly eccentric GW mergers. This was recently illustrated by Samsing & Ramirez-Ruiz (2017).

6.4. Cross Section and Rate of GW Inspirals

While our formalism directly shows that the GW inspiral cross section varies as $a_{0}^{2/7}$, it is not clear how sensitive it is to changes in the initial mass hierarchy. To explore this, we have performed a controlled analytical experiment, where only one of the masses was varied away from its equal mass value. In this case we found that the change in GW inspiral cross section is only weekly dependent on the fractional mass change. This led us to conclude that the equal mass case actually seems to provide a fairly good description of BBH mergers assembled in a typical GC, as these recently have been found to have a mass ratio close to unity (Rodriguez et al., 2016a). This was also illustrated in Samsing & Ramirez-Ruiz (2017). In this approximation, we have showed that the rate of GW inspirals is of order $\Gamma_{I}\approx 1\ \text{Gpc}^{-3}\ \text{yr}^{-1}$ for a typical GC population, if the number of BHs in each cluster is about $\approx 80$. As the number of BHs in GCs is largely unknown, this rate estimate is of course associated with a large uncertainty. For this reason, we mainly explore relative rates in this paper.

6.5. Formation of Post-Interaction GW Mergers

The majority of GW mergers do not form during the interaction as inspirals, but instead as post-interaction binary mergers. To facilitate comparison, we have calculated the cross section of such post-interaction binary mergers and found it to have a piecewise scaling with $a_{0}$. The cross section is shown to first increases as $\propto a_{0}$ until a characteristic value $a_{\rm 0}^{<\tau}$, after which it decreases as $\propto a_{0}^{-1/7}$. We thus conclude that the post-interaction binary merger cross section reaches its maximum value when $a_{0}=a_{\rm 0}^{<\tau}$. In the equal mass case, the value of $a_{\rm 0}^{<\tau}$ is about the SMA for which the target binary has a GW merger time $\tau$.

6.6. GW Inspirals and Post-Interaction GW Mergers

Relative rates can be determined much more accurately than absolute ones. For this reason, we have explored the cross section of inspirals relative to that of post-interaction binary mergers. We found that the relative number of GW inspirals takes its minimum value when $a_{0}=a_{\rm 0}^{<\tau}$, which led us to the profound conclusion that $>1\%$ of all the BBH mergers assembled through binary-single interactions will be GW inspirals. As a result, if post-interaction binary mergers are as frequent as recently reported by Rodriguez et al. (2016a), GW inspirals are expected to be observed by Advanced LIGO.

6.7. Formation History of GW Inspirals

To gain further insight into the role of GW energy loss term corrections, we have calculated the time dependent number of GW inspirals and post-interaction binary mergers for two scenarios. In the first scenario, we have assumed that binaries and singles interact over a short period. In this case we have found that GW inspirals likely dominate the total GW merger rate within the first $\approx 10^{6}$ years. This interestingly implies that the first mergers in the Universe are likely to be inspirals. In the second scenario, we have explored the rates assuming steady state. For this case, we have found that the present day rate of GW inspirals relative to that of post-interaction binary mergers is approximately given by the ratio of their cross sections evaluated at $\tau=t_{\rm H}$. As a result, GW inspirals are expected to constitute at least $1\%$ of the present day rate of BBHs mergers. We also noted that GW inspirals are more likely to remain in the GC after their formation than standard post-interaction binary mergers, which could have interesting dynamical consequences (e.g. Gültekin et al., 2006).

6.8. Highly Eccentric GW Mergers Observable by LIGO

Although the majority of GW inspirals form with notable eccentricity, many might still experience significant circularization before being observable (e.g. Samsing & Ramirez-Ruiz, 2017). As a result, for a BBH to appear in the LIGO band with an extremely high eccentricity, it must necessarily form within the LIGO band. Using a simple prescription for calculating the GW peak frequency $f_{\rm GW}$, we have showed that a fixed value for $f_{\rm GW}$ corresponds to a fixed pericenter distance $r_{\rm p}(f_{\rm GW})$. The cross section of GW inspirals that form with a GW frequency, say, $>10$Hz, is therefore to leading order given by the collision cross section with $R_{\rm C}=r_{\rm p}(10\text{Hz})$. Making use of this, we derived that $\approx 20\%$ of all GW inspirals are expected to form within the LIGO band and thus appear with a particular high eccentricity. We note that this estimate agrees well with the one derived using $N$-body simulations (Samsing & Ramirez-Ruiz, 2017). Although $\approx 80\%$ all inspirals are not expected to form within the LIGO band, we still find that the majority of GW inspirals will enter the band with an eccentricity $>0.1$ (Samsing & Ramirez-Ruiz, 2017). If the BBH merger rate is dominated by the dynamically assembled population, then about $1\%$ of the observable rate will have an eccentricity $>0.1$.

The results presented here clearly illustrate that GW inspirals are not just a rare curiosity resulting from including GR corrections, but constitute a population with highly interesting dynamical and observational consequences. This study should motivate further work on understanding the role of GR corrections in few-body interactions as well as in full $N$-body GC calculations. As GW inspirals generally merge with notable eccentricity, they are likely to play a key role in differentiating between different BBH merger channels; a test that should become possible with Advanced LIGO.

It is a pleasure to thank C. L. Rodriguez, N. Leigh and T. Ilan for helpful discussions. Support for this work was provided by the David and Lucile Packard Foundation, UCMEXUS (CN-12-578), the Danish National Research Foundation and NASA through an 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.