An Analytical Portrait of Binary Mergers in Hierarchical Triple Systems

We assume three bodies with masses m₀, m₁, and m₂. In an inertial frame, the positions of the three bodies are described by three vectors, r₀, r₁, and r₂, respectively. With Newtonian gravitation, the Hamiltonian of the system is given by

$$\begin{eqnarray}H & = & \displaystyle \frac{1}{2{m}_{0}}| {{\boldsymbol{p}}}_{0}{| }^{2}+\displaystyle \frac{1}{2{m}_{1}}| {{\boldsymbol{p}}}_{1}{| }^{2}+\displaystyle \frac{1}{2{m}_{2}}| {{\boldsymbol{p}}}_{2}{| }^{2}-\displaystyle \frac{{{Gm}}_{0}{m}_{1}}{| {{\boldsymbol{r}}}_{0}-{{\boldsymbol{r}}}_{1}| }\\ & & -\displaystyle \frac{{{Gm}}_{0}{m}_{2}}{| {{\boldsymbol{r}}}_{0}-{{\boldsymbol{r}}}_{2}| }-\displaystyle \frac{{{Gm}}_{1}{m}_{2}}{| {{\boldsymbol{r}}}_{1}-{{\boldsymbol{r}}}_{2}| },\end{eqnarray} \tag{ 1 }$$

where ${{\boldsymbol{p}}}_{i}={m}_{i}{\dot{{\boldsymbol{r}}}}_{i}\,(i=0,1,2)$ is the momentum conjugate to ${{\boldsymbol{r}}}_{i}$, and G is the gravitational constant, which equals $4{\pi }^{2}$ in astronomical units (au, yr, M_⊙).

It is well known that the Keplerian two-body problem can be recast into two independent motions, namely, the (trivial) motion of the mass center and the motion of the reduced mass relative to the mass center. The three-body problem considered here can be treated similarly. We group m₀ and m₁ as a binary system, called the "inner" binary. The inner binary and m₂ then form an effective two-body system, called the "outer" binary. Thus, we introduce the following set of variables:

$$\begin{eqnarray}m & \equiv & {m}_{0}+{m}_{1},\ M\equiv {m}_{0}+{m}_{1}+{m}_{2},\\ {\mu }_{1} & \equiv & \displaystyle \frac{{m}_{0}{m}_{1}}{{m}_{0}+{m}_{1}},\ {\mu }_{2}\equiv \displaystyle \frac{({m}_{0}+{m}_{1}){m}_{2}}{{m}_{0}+{m}_{1}+{m}_{2}},\end{eqnarray} \tag{ 2 }$$

which correspond to the total masses of the inner binary and the outer binary and the reduced masses of the inner binary and the outer binary, respectively. We also introduce the following position vectors:

$$\begin{eqnarray}{\boldsymbol{R}} & = & \displaystyle \frac{1}{M}({m}_{0}{{\boldsymbol{r}}}_{0}+{m}_{1}{{\boldsymbol{r}}}_{1}+{m}_{2}{{\boldsymbol{r}}}_{2}),\ {{\boldsymbol{R}}}_{1}={{\boldsymbol{r}}}_{1}-{{\boldsymbol{r}}}_{0},\\ {{\boldsymbol{R}}}_{2} & = & {{\boldsymbol{r}}}_{2}-\displaystyle \frac{1}{m}({m}_{0}{{\boldsymbol{r}}}_{0}+{m}_{1}{{\boldsymbol{r}}}_{1}),\end{eqnarray} \tag{ 3 }$$

and the corresponding momenta,

$$\begin{eqnarray}&&{\boldsymbol{P}}=M\dot{{\boldsymbol{R}}},\ {{\boldsymbol{\Pi }}}_{1}={\mu }_{1}{\dot{{\boldsymbol{R}}}}_{1},\ {{\boldsymbol{\Pi }}}_{2}={\mu }_{2}{\dot{{\boldsymbol{R}}}}_{2}.\end{eqnarray} \tag{ 4 }$$

Then, it is easy to see that the Hamiltonian (1) can be rewritten as

$$\begin{eqnarray}H & = & \displaystyle \frac{1}{2M}| {\boldsymbol{P}}{| }^{2}+\left(\displaystyle \frac{1}{2{\mu }_{1}}| {{\boldsymbol{\Pi }}}_{1}{| }^{2}-\displaystyle \frac{{Gm}{\mu }_{1}}{| {{\boldsymbol{R}}}_{1}| }\right)\\ & & +\ \left(\displaystyle \frac{1}{2{\mu }_{2}}| {{\boldsymbol{\Pi }}}_{2}{| }^{2}-\displaystyle \frac{{GM}{\mu }_{2}}{| {{\boldsymbol{R}}}_{2}| }\right)+H^{\prime} ,\end{eqnarray} \tag{ 5 }$$

where the system breaks down to three "independent" motions, namely, the motion of the mass center of the triple (which is trivial and will be neglected), the motion of the inner and outer binary, plus the "interaction" between the inner and outer binary described by $H^{\prime}$,

$$\begin{eqnarray}&&H^{\prime} =\displaystyle \frac{{{Gmm}}_{2}}{| {{\boldsymbol{R}}}_{2}| }-\displaystyle \frac{{{Gm}}_{0}{m}_{2}}{\left|{{\boldsymbol{R}}}_{2},+,\tfrac{{m}_{1}}{m},{{\boldsymbol{R}}}_{1}\right|}-\displaystyle \frac{{{Gm}}_{1}{m}_{2}}{\left|{{\boldsymbol{R}}}_{2},-,\tfrac{{m}_{0}}{m},{{\boldsymbol{R}}}_{1}\right|}.\end{eqnarray} \tag{ 6 }$$

Assuming that the inner and outer binaries are "weakly coupled," we can perform a perturbative multipole expansion of $H^{\prime}$ at position R₂ as follows:

$$\begin{eqnarray}H^{\prime} & = & \displaystyle \sum _{{\ell }=2}^{\infty }{H}^{({\ell })},\\ {H}^{({\ell })} & = & -\displaystyle \frac{{{Gm}}_{0}{m}_{1}{m}_{2}}{{m}^{{\ell }}}[{m}_{0}^{{\ell }-1}+{(-1)}^{{\ell }}{m}_{1}^{{\ell }-1}]\displaystyle \frac{{R}_{1}^{{\ell }}}{{R}_{2}^{{\ell }+1}}{P}_{{\ell }}(\cos \varphi ),\end{eqnarray} \tag{ 7 }$$

where φ is the angle between R₁ and R₂, and ${P}_{{\ell }}(z)$ is the Legendre polynomial. The leading nonvanishing term is the quadrupole interaction with ℓ = 2,

$$\begin{eqnarray}&&{H}^{(2)}=-\displaystyle \frac{{{Gm}}_{0}{m}_{1}{m}_{2}}{2m}\displaystyle \frac{{R}_{1}^{2}}{{R}_{2}^{3}}(3{\cos }^{2}\varphi -1),\end{eqnarray} \tag{ 8 }$$

which will be our main focus. To quantify the meaning of the weak coupling, we compare the magnitude of H⁽²⁾ with the Hamiltonian of the two binaries, i.e., the terms in parentheses in Equation (5). It is easy to see that the conditions for weak coupling between the inner and outer binaries are

$$\begin{eqnarray}&&\displaystyle \frac{{m}_{2}}{m}{\left(\displaystyle \frac{{R}_{1}}{{R}_{2}}\right)}^{3}\ll 1,\ \displaystyle \frac{{\mu }_{1}}{m}{\left(\displaystyle \frac{{R}_{1}}{{R}_{2}}\right)}^{2}\ll 1.\end{eqnarray} \tag{ 9 }$$

Therefore, in the weak coupling regime, we can retain the quadrupole term (8) only and study its perturbation on the motion of the two binaries. Before doing so, in the next subsection we will introduce orbital parameters that are more suitable for perturbation theory.

In the Hamiltonian (5) the motion of the triple breaks down to two Keplerian two-body orbits, plus the perturbative interaction between them. It is thus useful to review the standard orbital parameters describing the configuration and orientation of the two-body orbit. In this subsection we focus on one orbit only and thus drop the subscript (1, 2) distinguishing the inner and outer orbits.

Since we are concerned only with bound orbits, which are always elliptical, the size and the shape of the orbit can always be described by two parameters, namely, the semimajor axis a and the eccentricity e. For example, the semiminor axis b is given by $b=a\sqrt{1-{e}^{2}}$, and the distance from the periapsis to the focus is $a(1-e)$. Then, the location of the rotating body in the orbital plane can be determined by one additional parameter ψ, called the true anomaly, which is defined to be the angle from the periapsis to the rotating body in the orbital plane, as shown in Figure 1. In addition to $(a,e,\psi )$, we need three Euler angles to characterize the orientation of the orbit, relative to some reference plane, also shown in Figure 1. Conventionally, the three angles are chosen to be $(\vartheta ,I,\gamma )$, where ϑ is the angle from the reference direction to the ascending node on the reference plane (called the longitude of ascending node), I is the angle between the reference plane and the orbital plane (called the inclination), and γ is the angle from the ascending node to the periapsis in the orbital plane (called the argument of periapsis). The six-parameter set $(a,e,\psi ,\theta ,I,\gamma )$ then completely characterizes the position of the rotating body.

Zoom In Zoom Out Reset image size

In addition to these, we also use alternative widely used parameters in the following. First, the mean anomaly β and the eccentric anomaly u are related to the true anomaly ψ through the following relations:

$$\begin{eqnarray}&&\beta =u-e\sin u,\quad \cos \psi =\displaystyle \frac{\cos u-e}{1-e\cos u}.\end{eqnarray} \tag{ 10 }$$

The mean anomaly β undergoes periodic motion and has the same period as the true anomaly ψ, but it is defined such that the motion is uniform in time, namely, β = ωt, with the orbital frequency ${\omega }^{2}={Gm}/{a}^{3}$. The eccentric anomaly u makes the relation between the true anomaly ψ and the mean anomaly β explicit. Furthermore, since $1-{e}^{2}$ is a frequently appearing combination, it is useful to define it as a new parameter ,

$$\begin{eqnarray}&&\epsilon \equiv 1-{e}^{2}.\end{eqnarray} \tag{ 11 }$$

The Hamiltonian (5) for the triple system is expressed in terms of the coordinates of the three bodies and their conjugate momenta. While this set of coordinates is intuitively simple, they are less suitable for our study here since eventually we want to trace the slow variation of orbital parameters rather than the exact locations of the three bodies. For systems involving librations or periodic motions like ours, there are well-known angle-action variables that are in place for our purpose. In celestial mechanics, a widely used set of angle-action variables are Delauney variables. For a Keplerian two-body system, the three angle variables (β, γ, ϑ) are given by the mean anomaly β, the argument of the periapsis γ, and the longitude of ascending node ϑ, respectively. The three conjugate action variables $({J}_{\beta },{J}_{\gamma },{J}_{\vartheta })$ are given by

$$\begin{eqnarray}{J}_{\beta } & = & \mu \sqrt{{Gma}},\quad {J}_{\gamma }=\mu \sqrt{{Gma}(1-{e}^{2})},\\ {J}_{\vartheta } & = & \mu \sqrt{{Gma}(1-{e}^{2})}\cos I,\end{eqnarray} \tag{ 12 }$$

where m is the total mass of the binary, a is the semimajor axis of the orbit, e is the eccentricity, and I is the inclination. We note that the grouping of angle and action variables mixes up the position and orientation variables. The three angle variables include one position variable (the mean anomaly) and two Euler angles (γ, θ), while the third Euler angle I is in the action variable J_ϑ. We also note that ${J}_{\gamma }$ coincides with the magnitude of the angular momentum, while J_ϑ is the component of the angular momentum orthogonal to the reference plane. In terms of the Delauney variables, the Hamiltonian of an isolated binary becomes

$$\begin{eqnarray}&&{H}_{\mathrm{binary}}=-{\left(\displaystyle \frac{{Gm}}{{J}_{\beta }}\right)}^{2}.\end{eqnarray} \tag{ 13 }$$

The advantage of the Delauney variables can now be understood by noting that the Hamiltonian has explicit dependence on J_β only, which means that all Delauney variables but β are conserved for an isolated binary system, while β is under periodic motion for an elliptic orbit. After introducing a perturbation, the rest of the variables should undergo only slow motion, and we shall derive the corresponding equations in the next subsection.

The equations governing the long-term evolution of various orbital elements have three types of perturbations, namely, the tidal perturbation at quadrupole order, the PN precession of the orbit, and the back-reaction of the GW radiation. We consider these three pieces in turn.

2.3.1. Tidal Perturbation at Quadrupole Order

As mentioned earlier, an isolated Keplerian binary system has an elliptical orbit with fixed orbital parameters (a, e) and orientation. After turning on a perturbation, those parameters could undergo slow time variation. Specifically, for the hierarchical triple described by the Hamiltonian (5), we can rewrite the two binary terms in big parentheses in Delauney variables as in Equation (13). It remains to recast the interaction term $H^{\prime}$ similarly, which can be replaced by H⁽²⁾ in Equation (8) in quadrupole approximation. To this end, for the moment we choose the reference plane to be the plane of the outer orbit, and we choose the reference direction (the x-axis in Figure 1) to be the periapsis of the outer orbit. Then, we use $({a}_{i},{e}_{i},{\psi }_{i}),\,(i=1,2)$ to denote the orbital parameters of the inner and outer orbits, respectively. Since the three Euler angles $(\gamma ,I,\vartheta )$ define the orientation of the inner orbit relative to the outer orbit (but not the orientations of both orbits with respect to a third reference plane), we do not assign a subscript to them. With this notation, the leading-order binary motions of the two orbits can be represented by their position vectors ${{\boldsymbol{R}}}_{\mathrm{1,2}}$ as

$$\begin{eqnarray}{{\boldsymbol{R}}}_{1} & = & {R}_{1}\left(\begin{array}{ccc}\cos \vartheta & -\sin \vartheta & 0\\ \sin \vartheta & \cos \vartheta & 0\\ 0 & 0 & 1\end{array}\right)\left(\begin{array}{ccc}1 & 0 & 0\\ 0 & \cos I & -\sin I\\ 0 & \sin I & \cos I\end{array}\right)\\ & & \times \ \left(\begin{array}{ccc}\cos \gamma & -\sin \gamma & 0\\ \sin \gamma & \cos \gamma & 0\\ 0 & 0 & 1\end{array}\right)\left(\begin{array}{c}\cos {\psi }_{1}\\ \sin {\psi }_{1}\\ 0\end{array}\right),\end{eqnarray} \tag{ 14 }$$

and

$$\begin{eqnarray}&&{{\boldsymbol{R}}}_{2}={R}_{2}\left(\begin{array}{c}\cos {\psi }_{2}\\ \sin {\psi }_{2}\\ 0\end{array}\right),\end{eqnarray} \tag{ 15 }$$

where ${R}_{i}=| {{\boldsymbol{R}}}_{i}| ={a}_{i}(1-{e}_{i}^{2})/(1+{e}_{i}\cos {\psi }_{i}),\,(i=1,2)$ describe the familiar elliptical motions. Then the quadrupole Hamiltonian (8) can be written as

$$\begin{eqnarray}{H}^{(2)} & = & -\displaystyle \frac{{{Gm}}_{0}{m}_{1}{m}_{2}}{2m}\displaystyle \frac{{a}_{1}^{2}{(1-{e}_{1}^{2})}^{2}{(1+{e}_{2}\cos {\psi }_{2})}^{3}}{{a}_{2}^{3}{(1-{e}_{2}^{2})}^{3}{(1+{e}_{1}\cos {\psi }_{1})}^{2}}\\ & & \times \,\{3[\cos ({\psi }_{1}+\gamma )\cos ({\psi }_{2}-\vartheta )\\ & & +\,{\sin ({\psi }_{1}+\gamma )\sin ({\psi }_{2}-\vartheta )\cos I]}^{2}-1\}.\end{eqnarray} \tag{ 16 }$$

To trace the slow motion of the orbital parameters under the perturbation of ${H}^{(2)}$, as is standard we "integrate out" the fast periodic motions of both the inner and outer orbits, which yields an "effective Hamiltonian." In the literature this is called the secular approximation, and the resulting effective Hamiltonian is called the doubly averaged Hamiltonian (Lidov & Ziglin 1976). Conventionally, this is done by averaging Equation (16) over both mean anomalies β_1,2, since mean anomalies are defined to have uniform motion in time. From Equation (10) we know that the mean anomaly is related to the true anomaly by

$$\begin{eqnarray}&&d{\beta }_{i}=\displaystyle \frac{{(1-{e}_{i}^{2})}^{3/2}}{{(1+{e}_{i}\cos {\psi }_{i})}^{2}}d{\psi }_{i}.\,\,\,\,(i=1,2).\end{eqnarray} \tag{ 17 }$$

Therefore, the averaged Hamiltonian ${\overline{H}}^{(2)}$ can be worked out as

$$\begin{eqnarray}{\overline{H}}^{(2)} & \equiv & \,\displaystyle \frac{1}{4{\pi }^{2}}{\displaystyle \int }_{0}^{2\pi }d{\beta }_{1}d{\beta }_{2}\,{H}^{(2)}\\ & = & -\displaystyle \frac{{{Gm}}_{0}{m}_{1}{m}_{2}}{32m}\displaystyle \frac{{a}_{1}^{2}}{{a}_{2}^{3}{(1-{e}_{2}^{2})}^{3/2}}\\ & & \times \,[(2+3{e}_{1}^{2})(1+3\cos 2I)+30{e}_{1}^{2}\cos 2\gamma {\sin }^{2}I].\end{eqnarray} \tag{ 18 }$$

It is convenient to rewrite ${\overline{H}}^{(2)}$ in terms of a conserved and dimensionless function W, following Lidov & Ziglin (1976), as

$$\begin{eqnarray}&&{\overline{H}}^{(2)}=-K\left(W+\displaystyle \frac{5}{3}\right),\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&K\equiv \displaystyle \frac{3{{Gm}}_{0}{m}_{1}{m}_{2}}{8m}\displaystyle \frac{{a}_{1}^{2}}{{a}_{2}^{3}{(1-{e}_{2}^{2})}^{3/2}},\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&W\equiv (-2+{\cos }^{2}I)(1-{e}_{1}^{2})+5{e}_{1}^{2}({\cos }^{2}I-1){\sin }^{2}\gamma .\end{eqnarray} \tag{ 21 }$$

Some remarks about the perturbed Hamiltonian are as follows.

• 1.  
  Alternatively, it is possible to apply a canonical transformation, known as a Von Zeipel transformation, to the original Hamiltonian (5) to eliminate the fast modes corresponding to the mean anomalies ${\beta }_{\mathrm{1,2}}$ at the quadrupole level (Naoz et al. 2013a). In other words, the two mean anomalies can be eliminated without spoiling the canonical structure of the Hamiltonian. Thus, we conclude that the action variables conjugate to ${\beta }_{\mathrm{1,2}}$, namely, ${J}_{\beta 1}$ and J_β2, are conserved quantities, and it follows immediately that the semimajor axes and energies of both orbits are separately conserved by the doubly averaged Hamiltonian.
• 2.  
  The Hamiltonian ${\overline{H}}^{(2)}$ is obtained by substituting the leading-order solutions (i.e., Keplerian orbits), and this result resembles the expectation value of the perturbed Hamiltonian calculated in the time-independent perturbation theory of quantum mechanics. In this sense we can say that Equation (19) is an "on-shell" expression for ${\overline{H}}^{(2)}$. Consequently, while Equation (19) can be very useful when applying "on-shell" arguments such as energy conservation, the (off-shell) functional dependence on canonical variables in Equation (19) has been obscured. Therefore, one cannot immediately derive the equations of motion from Equation (19). In particular, it would be wrong to conclude that ${J}_{\vartheta 1}$ and ${J}_{\vartheta 2}$ are separately conserved from the fact that Equation (19) is apparently independent of ${\vartheta }_{\mathrm{1,2}}$.

To derive the equations for the secular evolution of the orbital parameters, we rewrite the Hamiltonian in terms of Delauney variables of both the inner and the outer orbits. To define these Delauney variables, it is helpful to switch the choice of reference plane to the one perpendicular to the total angular momentum of the triple system, known as the invariant plane. We then have six Delauney variables for each of the two orbits, but not all of them are independent dynamical variables. First, we note that ${J}_{\gamma i}=| {{\boldsymbol{J}}}_{i}|$ and ${J}_{\vartheta i}={{\boldsymbol{J}}}_{i}\cdot {\boldsymbol{J}}/| {\boldsymbol{J}}| \,(i=1,2)$, where ${{\boldsymbol{J}}}_{\mathrm{1,2}}$ are angular momenta of the inner and outer orbits, and ${\boldsymbol{J}}={{\boldsymbol{J}}}_{1}+{{\boldsymbol{J}}}_{2}$ is the total angular momentum. Therefore, these variables are related to each other via

$$\begin{eqnarray}&&{J}_{\gamma 1}^{2}-{J}_{\vartheta 1}^{2}={J}_{\gamma 2}^{2}-{J}_{\vartheta 2}^{2},\quad {J}_{\vartheta 1}+{J}_{\vartheta 2}=J.\end{eqnarray} \tag{ 22 }$$

Thus, we can eliminate ${J}_{\vartheta \mathrm{1,2}}$ in favor of ${J}_{\gamma \mathrm{1,2}}$,

$$\begin{eqnarray}&&{J}_{\vartheta 1}=\displaystyle \frac{{J}^{2}+{J}_{\gamma 1}^{2}-{J}_{\gamma 2}^{2}}{2J},\quad {J}_{\vartheta 2}=\displaystyle \frac{{J}^{2}-{J}_{\gamma 1}^{2}+{J}_{\gamma 2}^{2}}{2J}.\end{eqnarray} \tag{ 23 }$$

Note that the quadrupole Hamiltonian (19) is independent of ${\vartheta }_{\mathrm{1,2}}$, and thus we have removed the conjugate pair $({J}_{\vartheta ,i},{\vartheta }_{i})$. Second, we have noted before that ${J}_{\beta \mathrm{1,2}}$ are conserved quantities. Finally, one can see that ${\gamma }_{1}=\gamma$, where γ is defined above with respect to the outer orbital plane, since the angular momenta of the two orbits are coplanar with the total angular momentum. Therefore, we see that the quadrupole Hamiltonian (19) is also independent of γ₂. These arguments reduce the triple system into an integrable system with only one pair of conjugate variables γ₁ and ${J}_{\gamma 1}$. We further note that the total inclination $I={I}_{1}+{I}_{2}$, i.e., the inclination of the inner orbit relative to the outer orbit, is related to various angular momenta via

$$\begin{eqnarray}&&\cos I=\displaystyle \frac{{J}^{2}-{J}_{\gamma 1}^{2}-{J}_{\gamma 2}^{2}}{2{J}_{\gamma 1}{J}_{\gamma 2}}.\end{eqnarray} \tag{ 24 }$$

Now we can rewrite the quadratic Hamiltonian (19) in terms of γ₁ and ${J}_{\gamma \mathrm{1,2}}$ as

$$\begin{eqnarray}{\overline{H}}^{(2)} & = & -\displaystyle \frac{3G{\mu }_{1}{m}_{2}}{8}\displaystyle \frac{{a}_{1}^{2}{J}_{\beta 2}^{3}}{{a}_{2}^{3}{J}_{\gamma 2}^{3}}\left\{\displaystyle \frac{{J}_{\gamma 1}^{2}}{{J}_{\beta 1}^{2}}\left[-2+{\left(\displaystyle \frac{{J}_{\gamma 1}^{2}+{J}_{\gamma 2}^{2}-{J}^{2}}{2{J}_{\gamma 1}{J}_{\gamma 2}}\right)}^{2}\right]\right.\\ & & +\,\left.5\left(1-\displaystyle \frac{{J}_{\gamma 1}^{2}}{{J}_{\beta 1}^{2}}\right)\left[{\left(\displaystyle \frac{{J}_{\gamma 1}^{2}+{J}_{\gamma 2}^{2}-{J}^{2}}{2{J}_{\gamma 1}{J}_{\gamma 2}}\right)}^{2}-1\right]{\sin }^{2}{\gamma }_{1}+\displaystyle \frac{5}{3}\right\}.\end{eqnarray} \tag{ 25 }$$

Then, the equations for the orbital parameters can be derived from the canonical equations ${\dot{\gamma }}_{1}=\partial {\overline{H}}^{(2)}/\partial {J}_{\gamma 1}$ and ${\dot{J}}_{\gamma 1}=-\partial {\overline{H}}^{(2)}/\partial {\gamma }_{1}$:

$$\begin{eqnarray}&&{\left.\displaystyle \frac{{{de}}_{1}}{{dt}}\right|}_{\mathrm{KL}}=\,\displaystyle \frac{5K}{{J}_{\gamma 1}}{e}_{1}(1-{e}_{1}^{2})(1-{\cos }^{2}I)\sin 2{\gamma }_{1},\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}{\left.\displaystyle \frac{d{\gamma }_{1}}{{dt}}\right|}_{\mathrm{KL}} & = & \,2K\left[\displaystyle \frac{1}{{J}_{\gamma 1}}(2(1-{e}_{1}^{2})-5(1-{e}_{1}^{2}-{\cos }^{2}I)\right.\\ & & \times \,\left.{\sin }^{2}{\gamma }_{1})+\displaystyle \frac{1}{{J}_{\gamma 2}}(1-{e}_{1}^{2}+5{e}_{1}^{2}{\sin }^{2}{\gamma }_{1})\cos I\right].\end{eqnarray} \tag{ 27 }$$

We now have the equations of motion governing the secular evolution of the inner binary. When applied to inspiraling BBHs, they should be supplemented by the equations from the PN correction and GW back-reactions, which we elaborate in the following.

2.3.2. Post-Newtonian Correction

As is well known, the first nontrivial order of the PN correction to the Keplerian potential is a trivial constant shift plus a new term proportional to the inverse square of the distance, for which the net effect is to generate a precession for the periapsis of an elliptical orbit. The precession rate is given by (see, e.g., Weinberg 1972)

$$\begin{eqnarray}&&{\left.\displaystyle \frac{d{\gamma }_{1}}{{dt}}\right|}_{\mathrm{PN}}=\displaystyle \frac{3}{{c}^{2}{a}_{1}(1-{e}_{1}^{2})}{\left(\displaystyle \frac{{Gm}}{{a}_{1}}\right)}^{3/2}.\end{eqnarray} \tag{ 28 }$$

While we can just include this term in the evolution equation for γ₁ and study its effect by solving the equation, it will be useful to reconstruct a corresponding term in the Hamiltonian, which will allow us to apply an energy conservation argument to estimate the merger time in Equation (41). Since a $1/{r}^{2}$ potential results from a conservative force, the corresponding Hamiltonian term must also be conserved. In principle, one can derive the desired Hamiltonian again from the original PN Hamiltonian by double-averaging over ${\beta }_{\mathrm{1,2}}$. Here we adopt a simpler poor man's derivation by rewriting ${\dot{\gamma }}_{1}{| }_{\mathrm{PN}}$ in terms of Delauney's variables and integrating:

$$\begin{eqnarray}\displaystyle \frac{\partial {H}_{\mathrm{PN}}}{\partial {J}_{\gamma 1}} & = & {\left.\displaystyle \frac{d{\gamma }_{1}}{{dt}}\right|}_{\mathrm{PN}}=\displaystyle \frac{3}{{c}^{2}}\displaystyle \frac{{Gm}{\mu }^{2}}{{J}_{\gamma 1}^{2}}{\left(\displaystyle \frac{{Gm}\mu }{{J}_{\beta 1}}\right)}^{3}\,\,\Rightarrow \,\,{H}_{\mathrm{PN}}\\ & = & \displaystyle \int {{dJ}}_{\gamma 1}{\left.\displaystyle \frac{d{\gamma }_{1}}{{dt}}\right|}_{\mathrm{PN}}=-\displaystyle \frac{3{G}^{2}{m}^{2}\mu }{{c}^{2}{a}_{1}^{2}\sqrt{1-{e}_{1}^{2}}}.\end{eqnarray} \tag{ 29 }$$

It is convenient to rewrite ${H}_{\mathrm{PN}}$ as

$$\begin{eqnarray}&&{H}_{\mathrm{PN}}=-{{KW}}_{\mathrm{PN}},\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{W}_{\mathrm{PN}}=\displaystyle \frac{8{{Gm}}^{2}{a}_{2}^{3}{(1-{e}_{2}^{2})}^{3/2}}{{c}^{2}{m}_{2}{a}_{1}^{4}{(1-{e}_{1}^{2})}^{1/2}}.\end{eqnarray} \tag{ 31 }$$

where K is defined in Equation (20). In the following, we shall also use the parameter ${{\rm{\Theta }}}_{\mathrm{PN}}\equiv {W}_{\mathrm{PN}}\sqrt{1-{e}_{1}^{2}}$, which is conventional in the literature. We note that Equation (30), as with Equation (19), is an "on-shell" expression that can be used when applying energy conservation but is not suitable for deriving the equation of motion.

There are additional PN corrections that we have not included. The correction to the quadrupole coupling between the orbits is subdominant in that it is smaller than the Newtonian quadrupole, which is already a perturbation. The PN correction to the outer orbit gives nonzero ${\dot{\gamma }}_{2}$, which is irrelevant to the KL oscillations at quadrupole order since it is independent of γ₂. This is no longer valid in systems such as GCs when octupole interactions play an essential role (Naoz et al. 2013b). An N-body code of the triple system with PN corrections was presented in Bonetti et al. (2016), and further subtleties in the calculation of GW radiation from a hierarchical triple were discussed in Bonetti et al. (2017).

2.3.3. GW Radiation

The radiated GWs from the inner binary carry energy and angular momentum away from the system, leading to a reduction of both the semimajor axis a and the eccentricity e. Throughout the paper we assume that the semimajor axis of the outer orbit is much greater than that of the inner orbit, so that the GW back-reaction on the outer orbit can be neglected and we can consider the influence of GWs only on a₁ and e₁. The equations governing ${\dot{a}}_{1}$ and ${\dot{e}}_{1}$, known as Peters's equations (Peters 1964), have been reviewed in Randall & Xianyu (2018), and we only quote the results here:

$$\begin{eqnarray}&&{\left.\displaystyle \frac{{{da}}_{1}}{{dt}}\right|}_{\mathrm{GW}}=-\displaystyle \frac{64}{5}\displaystyle \frac{{G}^{3}{\mu }_{1}{m}^{2}}{{c}^{5}{a}_{1}^{3}}\displaystyle \frac{1}{{(1-{e}_{1}^{2})}^{7/2}}\\ &&\,\times \left(1+\displaystyle \frac{73}{24}{e}_{1}^{2}+\displaystyle \frac{37}{96}{e}_{1}^{4}\right),\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&{\left.\displaystyle \frac{{{de}}_{1}}{{dt}}\right|}_{\mathrm{GW}}=-\displaystyle \frac{304}{15}\displaystyle \frac{{G}^{3}{\mu }_{1}{m}^{2}}{{c}^{5}{a}_{1}^{4}}\displaystyle \frac{{e}_{1}}{{(1-{e}_{1}^{2})}^{5/2}}\left(1+\displaystyle \frac{121}{304}{e}_{1}^{2}\right).\end{eqnarray} \tag{ 33 }$$

We solve a₁ as a function of e₁ from the above two equations as ${a}_{1}=g({e}_{1})$, with g(e) defined as

$$\begin{eqnarray}g(e)=\displaystyle \frac{{e}^{12/19}}{1-{e}^{2}}{\left(1+\displaystyle \frac{121}{304}{e}^{2}\right)}^{870/2299}\simeq \left\{\begin{array}{cc}{e}^{12/19} & (e\ll 1),\\ \displaystyle \frac{1.1352}{1-{e}^{2}} & (e\lesssim 1).\end{array}\right.\end{eqnarray} \tag{ 34 }$$

The GW radiation affects the evolution of a₁ and e₁ not only through Peters's equation above but also implicitly through Equations (26) and (27) since the right-hand side of these two equations depends on the inclination, which in turn depends on the total angular momentum J through Equation (24). Since the total angular momentum is no longer conserved after including GWs, we should add another equation describing the loss of it. Remember that we consider the GW back-reaction to the inner orbit only, and thus the loss of angular momentum happens only to $| {{\boldsymbol{J}}}_{1}| ={J}_{\gamma 1}$ but not to $| {{\boldsymbol{J}}}_{2}|$. This loss of angular momentum has been reviewed in Randall & Xianyu (2018) and can be written as

$$\begin{eqnarray}&&\displaystyle \frac{{{dJ}}_{\gamma 1}}{{dt}}=-\displaystyle \frac{32}{5}\displaystyle \frac{{G}^{7/2}{\mu }_{1}^{2}{m}^{5/2}}{{c}^{5}{a}_{1}^{7/2}}\displaystyle \frac{1}{{(1-{e}_{1}^{2})}^{2}}\left(1+\displaystyle \frac{7}{8}{e}_{1}^{2}\right).\end{eqnarray} \tag{ 35 }$$

Then, differentiating Equation (24), we get $\dot{J}\,=({J}_{\gamma 1}+{J}_{\gamma 2}\cos I){\dot{J}}_{\gamma 1}/J$, where we have used the fact that the GW does not affect the inclination I since it only reduces the magnitude of the angular momentum but does not influence its direction. Consequently,

$$\begin{eqnarray}{\left.\displaystyle \frac{{dJ}}{{dt}}\right|}_{\mathrm{GW}} & = & -\displaystyle \frac{32}{5}\displaystyle \frac{{G}^{7/2}{\mu }_{1}^{2}{m}^{5/2}}{{c}^{5}{a}_{1}^{7/2}}\\ & & \times \,\displaystyle \frac{{J}_{\gamma 1}+{J}_{\gamma 2}\cos I}{J}\displaystyle \frac{1}{{(1-{e}_{1}^{2})}^{2}}\left(1+\displaystyle \frac{7}{8}{e}_{1}^{2}\right).\end{eqnarray} \tag{ 36 }$$

2.3.4. Summary

All three pieces—the quadrupole interaction between the inner and outer orbits, the PN precession of the inner orbit, and the GW back-reaction—can now be assembled to give the following set of equations:

$$\begin{eqnarray}&&\displaystyle \frac{{{da}}_{1}}{{dt}}=-\displaystyle \frac{64}{5}\displaystyle \frac{{G}^{3}{\mu }_{1}{m}^{2}}{{c}^{5}{a}_{1}^{3}}\displaystyle \frac{1}{{(1-{e}_{1}^{2})}^{7/2}}\left(1+\displaystyle \frac{73}{24}{e}_{1}^{2}+\displaystyle \frac{37}{96}{e}_{1}^{4}\right),\end{eqnarray} \tag{ 37a }$$

$$\begin{eqnarray}\displaystyle \frac{{{de}}_{1}}{{dt}} & = & \,\displaystyle \frac{5K}{{J}_{\gamma 1}}{e}_{1}(1-{e}_{1}^{2})(1-{\cos }^{2}I)\sin 2{\gamma }_{1}\\ & & -\,\displaystyle \frac{304}{15}\displaystyle \frac{{G}^{3}{\mu }_{1}{m}^{2}}{{c}^{5}{a}_{1}^{4}}\displaystyle \frac{{e}_{1}}{{(1-{e}_{1}^{2})}^{5/2}}\left(1+\displaystyle \frac{121}{304}{e}_{1}^{2}\right),\end{eqnarray} \tag{ 37b }$$

$$\begin{eqnarray}\displaystyle \frac{d{\gamma }_{1}}{{dt}} & = & 2K\left[\displaystyle \frac{1}{{J}_{\gamma 1}}(2(1-{e}_{1}^{2})-5(1-{e}_{1}^{2}-{\cos }^{2}I){\sin }^{2}{\gamma }_{1})\right.\\ & & +\left.\,\displaystyle \frac{1}{{J}_{\gamma 2}}(1-{e}_{1}^{2}+5{e}_{1}^{2}{\sin }^{2}{\gamma }_{1})\cos I\right]\\ & & +\,\displaystyle \frac{3}{{c}^{2}{a}_{1}(1-{e}_{1}^{2})}{\left(\displaystyle \frac{{Gm}}{{a}_{1}}\right)}^{3/2},\end{eqnarray} \tag{ 37c }$$

$$\begin{eqnarray}\displaystyle \frac{{dJ}}{{dt}} & = & -\displaystyle \frac{32}{5}\displaystyle \frac{{G}^{7/2}{\mu }_{1}^{2}{m}^{5/2}}{{c}^{5}{a}_{1}^{7/2}}\\ & & \times \,\displaystyle \frac{{J}_{\gamma 1}+{J}_{\gamma 2}\cos I}{J}\displaystyle \frac{1}{{(1-{e}_{1}^{2})}^{2}}\left(1+\displaystyle \frac{7}{8}{e}_{1}^{2}\right).\end{eqnarray} \tag{ 37d }$$

This set of equations is to be supplemented by the definition of K in Equation (20), the formulae for ${J}_{\gamma \mathrm{1,2}}$ in Equation (12), and the relation given by Equation (24) between the inclination I and the total angular momenta J. Then, everything can be expressed in terms of the four independent variables $({a}_{1},{e}_{1},{\gamma }_{1},J)$.

At this point, in principle we can sample the initial parameters following a given distribution and solve this set of equations to get the final distribution of eccentricity at the LIGO threshold. But it is far preferable to find a direct analytical relation between the initial parameters and the final eccentricity. Such a solution would presumably elucidate the dependence of the final answer on the input parameters, and the calculation would also be much faster than directly solving the equations case by case—hopefully yielding some insight into the environments in which observable BBH mergers occur. The challenge, of course, is the large number of parameters, but an analytical solution will help pinpoint what might be possible to extract. In the next section we set out to find this analytical mapping between initial parameters and final eccentricity. We can then get the final eccentricity distribution by integrating over the initial parameter space, given any initial distribution of binaries.

This problem involves a large number of parameters, but some should be measurable by LIGO, and only some of the others significantly impact the resultant eccentricity. In principle, this makes it feasible to learn about the environment of the merger. Specifically, the solutions depend on the following parameters: the three masses $({m}_{0},{m}_{1},{m}_{2})$, two constant parameters of the outer orbit $({a}_{2},{e}_{2})$, and four initial conditions, which we can choose to be $({a}_{1},{e}_{1},{\gamma }_{1},I)$. Here we replace the magnitude of the total angular momentum J in Equation (37) by the inclination I by means of Equation (24). Given these parameters, we can follow the evolution of the triple system until the inner binary separation a₁ reaches the observational band. Among these parameters, two masses $({m}_{0},{m}_{1})$, together with the inner orbital parameter $({a}_{1},{e}_{1})$, are measurable after entering the LIGO band. The other parameters are not directly accessible from LIGO observation, but they can make an impact on the observables e₁. In the next section we will see that for slow mergers $({a}_{2},{e}_{2})$ and I have the strongest impact on the final e₁, while the impact of $({a}_{1},{e}_{1},{\gamma }_{1})$ is relatively mild.

In this subsection we present several qualitatively distinct solutions for the inner binary merger, all with initially highly inclined orbits. As we explore below, the difference in relative strengths of the KL oscillation and the PN precession is critical to the merger, as the former increases eccentricity while the latter suppresses it. The type of numerical solutions we present are not new, but we review them to gain some intuition that will be helpful when constructing our analytical method. Our numerical samples are generated using the secular Equations (37), which are free of numerical subtleties that could occur in a direct N-body simulation. A nice agreement between integrating secular equations and the direct N-body simulation has been shown in Antonini & Perets (2012) using the ARCHAIN integrator. Similar numerical solutions have been presented in Wen (2003) but for different parameters.

Throughout this subsection, we set the binary masses ${m}_{0}={m}_{1}=10\,{M}_{\odot }$ and the tertiary mass ${m}_{2}=4\times {10}^{6}\,{M}_{\odot }$. We assume that the initial inner semimajor axis is ${a}_{10}=0.1$ au and the initial inner and outer eccentricities are ${e}_{10}\,={e}_{20}=0.1$. Furthermore, we choose the initial inclination to be ${I}_{0}=89\buildrel{\circ}\over{.} 9$ and the initial argument of inner periapsis to be ${\gamma }_{10}=0^\circ$. By changing the initial value of the outer semimajor axis a₂₀, we can adjust the strength of the KL oscillation relative to PN corrections.

For our first example, we take ${a}_{20}=80$ au, which means that the inner BBHs are close enough to the central SMBH that the KL resonance is extremely effective. The maximal value of the eccentricity in this case can be so large that the binary merges within one KL cycle. As we can see from Figure 2, the merger happens when ${\epsilon }_{1}=1-{e}_{1}^{2}$ reaches its minimum ${\epsilon }_{1\min }\sim { \mathcal O }(0.001)$ on the timescale of the KL oscillation, which is about ≃5 yr in this example.

Zoom In Zoom Out Reset image size

In our second example, we take ${a}_{20}=150$ au, which is about twice the distance of the first example. The KL timescale is correspondingly longer. In several initial KL cycles, the KL timescale $\sim {(K/{J}_{\gamma 1})}^{-1}$ is about 50 yr, whereas the PN timescale ${(d{\gamma }_{1}/{dt}{| }_{\mathrm{PN}})}^{-1}$ is about 1 yr at ${\epsilon }_{1\min }$, as can be found from Equation (28). Hence, the PN precession is more effective in suppressing the maximal value of e₁ than it is in the first example. As a result, the reduction of a₁ in each KL cycle is less significant, and the merger of the inner binary takes a large number of cycles.

From Figure 3 we see that the merger time in this example is ∼7000 yr. In its early stage, ₁ can reach a minimum of ∼0.006. This minimum does not change significantly in the first half of the binary's life, but it becomes larger in the later stage owing to the stronger effect of GW emission. Were there no GWs, the frequency of the oscillation and the value of ${\epsilon }_{1\min }$ would stay constant as one can see from the gray curve in the top panel of Figure 3, since PN precession conserves energy and angular momentum. Accordingly, the semimajor axis ${a}_{1}(t)$ changes more slowly than in the first example. However, if we zoom in to examine the function ${a}_{1}(t)$ on smaller timescales, we see that, at least during the early stage of the merger, the reduction of a₁ occurs chiefly at the minimum of ₁, so that ${a}_{1}(t)$ has a stair-wise behavior. In addition, it can be seen from Figure 3 that the period of the KL oscillation decreases over time and that the minimal/maximal e₁ (i.e., maximal/minimal ₁) in each KL cycle increases/decreases monotonically with time until the KL oscillation is almost fully suppressed. All these phenomena can be understood by noting that the period of the KL oscillation in e₁ is effectively determined by γ₁ as one can see from Equation (26). The PN correction always advances the phase of γ₁ relative to the phases in its absence and thus effectively reduces the period and amplitude of the KL oscillation and hence eccentricity. The strength of PN precession, when expressed in terms of ${W}_{\mathrm{PN}}$ in Equation (31), is proportional to ${a}_{1}^{-4}$. Therefore, the PN correction gets stronger at later stages when a₁ is smaller. Consequently, the amplitude of the KL oscillation, which is in any case smaller with smaller a₁, is even smaller at later times.

Zoom In Zoom Out Reset image size

For our third example, we choose ${a}_{2}=420$ au. In this case, the PN precession has such a strong effect on the KL cycle and the maximum of e₁ (and thus the minimum of ₁) that the KL oscillation is strongly suppressed. From Figure 4 we see that ₁ in this case never becomes small enough to boost GW radiation. Consequently, the stair-wise function ${a}_{1}(t)$ in the second example gets "melted" into a smoother function of time. Therefore, the merger time in such cases is well approximated by the merger time of an isolated binary with initial eccentricity given by ${e}_{1\max }$, where ${e}_{1\max }$ is the maximal value of eccentricity during its tiny oscillations and can be extremely long.

Zoom In Zoom Out Reset image size

Now we estimate the merger time τ of the inner binary, which is defined to be the time the binary takes to coalesce with given initial parameters. Recall that this parameter is critical to determining whether a binary will merge or evaporate. As the previous numerical solutions show, there are three distinct regions of parameter space depending on the maximal value of eccentricity that the inner binary can reach: (1) For highly inclined binary orbits that are very close to a central BH, the eccentricity can be close to 1, in which case the merger happens within one KL cycle. We call this scenario a fast merger. (2) For the intermediate case, which we call KL-boosted, the merger happens after many KL cycles, and the reduction of a₁ happens chiefly at those times within a cycle when the eccentricity reaches the maximum. (3) For binaries very far from the central BH, no effective KL cycles exist and the reduction of a₁ is essentially a smooth function of t. We call this case the isolation limit. We call both KL-boosted binaries and binaries in the isolation limit slow mergers.

Due to the qualitatively different behavior in the three cases, we need to estimate the merger time for each of them independently. We assume arbitrary initial values of the triple parameters, including ${a}_{\mathrm{1,2}}$, ${e}_{\mathrm{1,2}}$, and I. However, we assume ${\gamma }_{1}=0$ because we could always let the inner binary evolve to ${\gamma }_{1}=0$ were it not so initially. Clearly, this assumption breaks down for the fast merger because the binary can merge before γ₁ evolves to 0. However, even in this case it is still sensible to assume ${\gamma }_{10}=0$ since the merger time estimated in this way is correct within an order of magnitude. Furthermore, the initial value of e₁₀ will change if we shift the initial time to $\gamma =0$. Therefore, assuming ${\gamma }_{1}=0$ initially will affect the initial distribution of e₁₀. But we will show later on that the final eccentricity at the LIGO threshold is insensitive to the initial distribution of e₁₀, so this correction is unimportant when integrating over initial distributions.

Now we estimate the merger time for the three cases: first the isolation limit, then the KL-boosted case, and finally the fast mergers.

Isolation limit. When tidal forces are small, the merger time ${\tau }_{\mathrm{iso}}$ can be well approximated by the merger time of an isolated binary of the same initial a₁₀ and e₁₀. Therefore, we can get ${\tau }_{\mathrm{iso}}$ by integrating Equation (32) from ${a}_{1}={a}_{10}$ to the coalescence ${a}_{1}=0$. Practically, it is more convenient to integrate Equation (33) instead, from ${e}_{1}={e}_{10}$ to ${e}_{1}=0$, since we know from Equation (76) that e₁ is small at small a₁ at coalescence. See Maggiore (2008) for more details. The resulting ${\tau }_{\mathrm{iso}}$ is

$$\begin{eqnarray}&&{\tau }_{\mathrm{iso}}=\displaystyle \frac{5}{256}\displaystyle \frac{{c}^{5}{a}_{10}^{4}}{{G}^{3}{m}_{0}{m}_{1}m}G({e}_{10}){(1-{e}_{10}^{2})}^{7/2},\end{eqnarray} \tag{ 38 }$$

with $G({e}_{10})$ defined by the following integral:

$$\begin{eqnarray}&&G({e}_{10})=\displaystyle \frac{48}{19{g}^{4}({e}_{10}){(1-{e}_{10}^{2})}^{7/2}}{\int }_{0}^{{e}_{10}}{de}\,\displaystyle \frac{{g}^{4}(e){(1-{e}^{2})}^{5/2}}{e\left(1+\tfrac{121}{304}{e}^{2}\right)}.\end{eqnarray} \tag{ 39 }$$

Using $g(e)\sim {e}^{19/12}$ for small e in Equation (76), we can see that the function $G({e}_{10})$ equals 1 when ${e}_{10}=0$ and increases monotonically with e₁₀, but it remains very close to 1 for most of ${e}_{10}\in (0,1)$, until it rapidly rises to 768/425 ≃ 1.80 when ${e}_{10}\to 1$. Therefore, it is a good approximation to neglect G₁₀ and just use the following expression for our estimate:

$$\begin{eqnarray}&&{\tau }_{\mathrm{iso}}\simeq \displaystyle \frac{5}{256}\displaystyle \frac{{c}^{5}{a}_{10}^{4}}{{G}^{3}{m}_{0}{m}_{1}m}{(1-{e}_{10}^{2})}^{7/2}.\end{eqnarray} \tag{ 40 }$$

KL-boosted. For KL-boosted mergers, the reduction of a₁ happens mostly when ₁ is around its minimum (i.e., the eccentricity e₁ reaches its maximum) during each KL cycle, since we learn from Equation (32) that ${\dot{a}}_{1}\propto {\epsilon }_{1}^{-7/2}$. More explicitly, when ₁ increases from its minimum ${\epsilon }_{1\min }$ by a factor of 2, ${\dot{a}}_{1}\propto$ drops to 9% of its maximal value at ${\epsilon }_{1\min }$. We now show how to estimate the merger time by considering an imagined and isolated binary, with initial separation given by a₁₀ and initial eccentricity set to ${e}_{1\max }$. Here ${e}_{1\max }$ is taken to be the maximal value of e₁ in the first several KL cycles, where ${e}_{1\max }$ does not change significantly. To proceed along these lines, we need to determine the maximal eccentricity of the first KL cycle ${e}_{1\max }$, or equivalently ${\epsilon }_{1\min }\equiv 1-{e}_{1\max }^{2}$.

For those mergers where the GW back-reaction is negligible in the first several KL cycles, the value ${\epsilon }_{1\min }$ can be estimated by the conservation of energy, that is, the sum ${\overline{H}}^{(2)}+{H}_{\mathrm{PN}}$, or more conveniently $W+{W}_{\mathrm{PN}}$, is a constant. Here we have three time-dependent parameters e₁, I, and γ₁. We input the initial values ${e}_{10}$ and I₀, while, as observed earlier, the initial value of γ₁ can always be set to zero without loss of generality because otherwise we can just wait until γ₁ returns to 0. On the other hand, we see from either Equation (21) or Equation (26) that the eccentricity e₁ reaches its maximum when ${\gamma }_{1}=\pi /2$. Therefore, we can solve ${e}_{1\max }$ from the following equation:

$$\begin{eqnarray}&&{[W+{W}_{\mathrm{PN}}]}_{{e}_{1}={e}_{10},I={I}_{0},{\gamma }_{1}=0}={[W+{W}_{\mathrm{PN}}]}_{{e}_{1}={e}_{1\max },{\gamma }_{1}=\pi /2}.\end{eqnarray} \tag{ 41 }$$

It is possible to derive an analytical approximation for ${\epsilon }_{1\min }$ when ${\epsilon }_{1\min }\ll 1$ (Wen 2003). To see this, we spell out Equation (41) explicitly as follows:

$$\begin{eqnarray}&&(-2+{\cos }^{2}{I}_{0}){\epsilon }_{10}+\displaystyle \frac{{{\rm{\Theta }}}_{\mathrm{PN}}}{\sqrt{{\epsilon }_{10}}}=(-2+{\cos }^{2}{I}_{\star }){\epsilon }_{1\min }\\ &&\quad +\,5(1-{\epsilon }_{1\min })({\cos }^{2}{I}_{\star }-1)+\displaystyle \frac{{{\rm{\Theta }}}_{\mathrm{PN}}}{\sqrt{{\epsilon }_{1\min }}},\end{eqnarray} \tag{ 42 }$$

where ${{\rm{\Theta }}}_{\mathrm{PN}}={W}_{\mathrm{PN}}\sqrt{1-{e}_{1}^{2}}$, as we defined earlier below Equation (31), and ${I}_{\star }$ is the inclination I evaluated at ${e}_{1}={e}_{1\max }$ and ${\gamma }_{1}=\pi /2$. The value of ${I}_{\star }$ is not independent and can be determined by the conservation of angular momentum. To see this, we rewrite Equation (24) as

$$\begin{eqnarray}&&{J}_{\gamma 1}\left(\cos I+\displaystyle \frac{{J}_{\gamma 1}}{2{J}_{\gamma 2}}\right)=\displaystyle \frac{{J}^{2}-{J}_{\gamma 2}^{2}}{2{J}_{\gamma 2}}.\end{eqnarray} \tag{ 43 }$$

Both J and ${J}_{\gamma 2}$ are conserved quantities; therefore, the combination on the right-hand side (and hence the left-hand side) of the above formula must be conserved as well. We evaluate this combination with both $({\epsilon }_{10},{I}_{0},{\gamma }_{10}=0)$ and $({\epsilon }_{1\min },{I}_{\star },{\gamma }_{1}=\pi /2)$ and equate them, so that ${I}_{\star }$ can be solved to be

$$\begin{eqnarray}&&{\cos }^{2}{I}_{\star }\simeq \displaystyle \frac{{\epsilon }_{10}}{{\epsilon }_{1\min }}{\left(\cos {I}_{0}+\displaystyle \frac{\sqrt{{\epsilon }_{10}}}{{J}_{\gamma 2}/{J}_{\beta 1}}\right)}^{2},\end{eqnarray} \tag{ 44 }$$

where we have used the condition $\sqrt{{\epsilon }_{1\min }}/({J}_{\gamma 2}/{J}_{\beta 1})\ll \cos {I}_{\star }$, which is true because both factors on the left-hand side are much smaller than 1 while the right-hand side is of ${ \mathcal O }(1)$. Substituting this solution back into Equation (42) and keeping the leading terms in the ${\epsilon }_{1\min }\ll 1$ limit, we get a quadratic equation for $\sqrt{{\epsilon }_{1\min }}$, which can be readily solved as

$$\begin{eqnarray}&&\sqrt{{\epsilon }_{1\min }}=\displaystyle \frac{{{\rm{\Theta }}}_{\mathrm{PN}}+\sqrt{{{\rm{\Theta }}}_{\mathrm{PN}}^{2}+{AC}}}{2A},\end{eqnarray} \tag{ 45 }$$

$$\begin{eqnarray}A & = & (-2+{\cos }^{2}{I}_{0}){\epsilon }_{10}+\displaystyle \frac{{{\rm{\Theta }}}_{\mathrm{PN}}}{\sqrt{{\epsilon }_{10}}}\\ & & +\,4{\epsilon }_{10}{\left(\cos {I}_{0}+\displaystyle \frac{\sqrt{{\epsilon }_{10}}}{{J}_{\gamma 2}/{J}_{\beta 1}}\right)}^{2}+5,\end{eqnarray} \tag{ 46 }$$

$$\begin{eqnarray}&&C=20{\epsilon }_{10}{\left(\cos {I}_{0}+\displaystyle \frac{\sqrt{{\epsilon }_{10}}}{{J}_{\gamma 2}/{J}_{\beta 1}}\right)}^{2}.\end{eqnarray} \tag{ 47 }$$

We may want to identify the lifetime of KL-boosted merger ${\tau }_{\mathrm{KL}}$ with the lifetime ${\tau }_{\mathrm{iso}}$ in Equation (40) of an imaginary isolated binary with initial ellipticity ${\epsilon }_{1\min }$. But this is not quite right because the inner binary spends only a small portion of its time around ${\epsilon }_{1\min }$ in each KL cycle, and thus a₁ is essentially inert the rest of the time. The proportion of time around ${\epsilon }_{1\min }$ in each KL cycle can be estimated by asking how long it takes to increase ₁ from ${\epsilon }_{1\min }$ by ${\rm{\Delta }}{\epsilon }_{1}\sim {\epsilon }_{1\min }$. From Equation (27), we see that the KL timescale is roughly ${\dot{\gamma }}_{1}^{-1}\sim {J}_{\gamma 1}/K$. To find the time duration that ₁ stays within $({\epsilon }_{1\min },2{\epsilon }_{1\min })$, we Taylor-expand the function ${\epsilon }_{1}(t)$ around its minimum ${\epsilon }_{1\min }={\epsilon }_{1}(0)$ to quadratic order, ${\epsilon }_{1}(t)\simeq {\epsilon }_{1\min }+{\ddot{\epsilon }}_{1}(0){t}^{2}$, where ${\ddot{\epsilon }}_{1}\sim {\epsilon }_{1\min }{\dot{\gamma }}_{1}\,\sim {\epsilon }_{1\min }K/{J}_{\gamma 1}$ by Equations (26) and (27). In this way we see that ₁ stays within $({\epsilon }_{1\min },2{\epsilon }_{1\min })$ with the time duration $\sim \sqrt{{\epsilon }_{1\min }}{J}_{\gamma 1}/K$ in each KL cycle, or, the inner binary spends only a proportion of $\sqrt{{\epsilon }_{1\min }}$ of the whole KL cycle around ${\epsilon }_{1\min }$. As a result, the merger time of KL-boosted mergers can be estimated as

$$\begin{eqnarray}&&{\tau }_{\mathrm{slow}}\simeq \displaystyle \frac{{\tau }_{\mathrm{iso}}}{\sqrt{{\epsilon }_{1\min }}}\simeq \displaystyle \frac{5}{256}\displaystyle \frac{{c}^{5}{a}_{10}^{4}}{{G}^{3}{m}_{0}{m}_{1}m}{\epsilon }_{1\min }^{3}.\end{eqnarray} \tag{ 48 }$$

We note that this expression also works well for binaries in the isolation limit so long as e₁₀ is not very close to 1, for the reasons that ${\epsilon }_{1\min }\simeq 1-{e}_{10}^{2}$ and that the ratio between Equations (40) and (48) is $\sqrt{1-{e}_{10}^{2}}\sim { \mathcal O }(1)$, so the formula for the merger time in Equation (48) covers both Kozai-boosted mergers and the merger isolation limit.

Fast mergers. This is the simplest case, as the merger time is simply the timescale of a KL oscillation, which is the time γ₁ takes to evolve from 0° to 90° as can be seen from Equation (26). Therefore, we can read the merger time directly from Equation (27) for $d{\gamma }_{1}/{dt}$ as

$$\begin{eqnarray}&&{\tau }_{\mathrm{fast}}\simeq {\dot{\gamma }}^{-1},\end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray}&&\dot{\gamma }\simeq \displaystyle \frac{K}{{J}_{\gamma 1}}+\displaystyle \frac{3}{{c}^{2}{a}_{10}{\epsilon }_{10}}{\left(\displaystyle \frac{{Gm}}{{a}_{10}}\right)}^{3/2}.\end{eqnarray} \tag{ 50 }$$

Both ${\tau }_{\mathrm{fast}}$ and ${\tau }_{\mathrm{slow}}$ underestimate the merger time if extrapolated beyond their respective validity ranges, and therefore the best estimate of the merger time is simply the maximum of the two,

$$\begin{eqnarray}&&\tau =\max \{{\tau }_{\mathrm{fast}},{\tau }_{\mathrm{slow}}\}.\end{eqnarray} \tag{ 51 }$$

In Figure 5 we plot the merger time τ computed from Equation (51) as a function of a₂ with other parameters fixed. In the same figure we also quote the merger time computed from directly integrating Equations (37) as was done in Antonini & Perets (2012). We see the nice agreement between the analytical estimate (51) and the direct integration.

Zoom In Zoom Out Reset image size

KL-boosted, revisited. Here we provide a complementary estimate of the merger time for KL-boosted binaries. The method used here is technically more involved, but we present it because it has a clear physical picture. This estimation uses two assumptions: (1) the reduction of a₁ happens only when e₁ reaches its maximum in each KL cycle, and (2) the maximal e₁ does not change very much over much of the lifetime of the binary. The second assumption is not very good at late stages, so the result derived in this way tends to underestimate the merger time, but it is still reasonably good.

To proceed, we determine the amount of a₁ reduced in each KL cycle, namely, the height of each step in the stair-like function ${a}_{1}(t)$ in Figure 3. Since this happens only in a very narrow range around the moment ${t}_{\star }$ where e₁ reaches the maximum, namely, ${e}_{1}({t}_{\star })={e}_{1\max }$, we can approximate e(t) around each ${t}_{\star }$ by

$$\begin{eqnarray}&&{e}_{1}(t)\simeq {e}_{1\max }+\displaystyle \frac{1}{2}{\ddot{e}}_{1}({t}_{\star }){(t-{t}_{\star })}^{2}.\end{eqnarray} \tag{ 52 }$$

Here ${\ddot{e}}_{1}({t}_{\star })$ can be easily obtained from the evolution equations by noting that ${\dot{e}}_{1}({t}_{\star })=0$,

$$\begin{eqnarray}&&{\ddot{e}}_{1}({t}_{\star })=-\displaystyle \frac{10K}{{J}_{\gamma 1}}{e}_{1\max }(1-{e}_{1\max }^{2})\dot{\gamma }({t}_{\star }),\end{eqnarray} \tag{ 53 }$$

$$\begin{eqnarray}&&\dot{\gamma }({t}_{\star })\simeq \displaystyle \frac{K}{{J}_{\gamma 1}}+\displaystyle \frac{3}{{c}^{2}{a}_{1}(1-{e}_{1\max }^{2})}{\left(\displaystyle \frac{{Gm}}{{a}_{1}}\right)}^{3/2}.\end{eqnarray} \tag{ 54 }$$

Then the reduced a₁ in each KL cycle can be found by integrating $\dot{a}(t)$ with e(t) chosen as Equation (52). The result is

$$\begin{eqnarray}| {\rm{\Delta }}{a}_{1}| & \simeq & \displaystyle \frac{64{G}^{3}{m}_{0}{m}_{1}m}{5{c}^{5}{a}_{1}^{3}({t}_{\star })}\left(1+\displaystyle \frac{73}{24}{e}_{1\max }^{2}+\displaystyle \frac{37}{96}{e}_{1\max }^{4}\right)\\ & & \times \,{\displaystyle \int }_{-\infty }^{\infty }\displaystyle \frac{{dt}}{{(1-{e}_{1}^{2}(t))}^{7/2}}\\ & \simeq & \displaystyle \frac{1024{G}^{3}{m}_{0}{m}_{1}m}{75{c}^{5}{a}_{1}^{3}({t}_{\star }){(1-{e}_{1\max }^{2})}^{3}}\left(1+\displaystyle \frac{73}{24}{e}_{1\max }^{2}+\displaystyle \frac{37}{96}{e}_{1\max }^{4}\right)\\ & & \times \,\displaystyle \frac{1}{\sqrt{{e}_{1\max }| {\ddot{e}}_{1}({t}_{\star })| }}.\end{eqnarray} \tag{ 55 }$$

Then it remains to estimate the width of the step, which is simply the timescale of KL oscillation, i.e.,

$$\begin{eqnarray}&&{({\rm{\Delta }}t)}^{-1}\simeq \dot{\gamma }\simeq \displaystyle \frac{K}{{J}_{\gamma 1}}+\displaystyle \frac{3}{{c}^{2}{a}_{1}(1-{e}_{10}^{2})}{\left(\displaystyle \frac{{Gm}}{{a}_{1}}\right)}^{3/2},\end{eqnarray} \tag{ 56 }$$

which has the same form as Equation (54) except that it is evaluated with initial e₁₀ rather than the maximum ${e}_{1\max }$.

With both the height and the width of the steps known, we can now write the merger time as

$$\begin{eqnarray}&&\tau ={\int }_{0}^{{a}_{10}}{{da}}_{1}\,\displaystyle \frac{{\rm{\Delta }}t}{| {\rm{\Delta }}{a}_{1}| },\end{eqnarray} \tag{ 57 }$$

which essentially agrees with Equation (48) when ${\epsilon }_{1\min }$ is not too large, which we have checked numerically. Indeed, we see that both Equations (57) and (48) have the same dependence on a₁₀ and ${e}_{1\max }$ when ${e}_{1\max }\sim 1$.

We are now in a position to figure out the eccentricity distribution of the BBHs at the time of entering the LIGO window. This amounts to mapping the initial parameter space to the values of eccentricity ${e}_{\mathrm{LIGO}}$ at the LIGO threshold, which we take to be 10 Hz. The ultimate goal is, given an observed distribution of ${e}_{\mathrm{LIGO}}$, to use this map as a way of probing the initial distribution of orbital parameters, and hence the BH mergers' environments.

The relevant initial orbital parameters include initial values of the inner binary separation a₁₀, of the eccentricity e₁₀, of the inclination I₀, and also of the outer orbital parameters a₂ and e₂, which are approximately constant, together with the tertiary mass m₂. When we identify the tertiary body to be the SMBH in the galactic center, m₂ will be fixed for a given galaxy. For slow mergers, we can always choose the initial value ${\gamma }_{10}=0$ without loss of generality since we can always wait until γ₁ becomes zero.

To establish the map from the initial parameter space to the eccentricity at the LIGO threshold, it is again helpful to consider fast mergers and slow mergers separately. Once again, by fast mergers we mean those binaries merging within ${ \mathcal O }(1)$ KL cycles, while slow mergers undergo many KL oscillations. For fast mergers, there is no obvious analytical method to estimate ${e}_{\mathrm{LIGO}}$, as the three effects, KL, PN, and GW, are all important in each KL cycle, and none can be neglected when estimating ${e}_{\mathrm{LIGO}}$. On the other hand, as will be elaborated in the following, we can estimate ${e}_{\mathrm{LIGO}}$ for slow mergers quite well using simple analytical methods, since there are then two distinct stages of binary evolution, dominated first by KL oscillations and then by GW circularization. The second phase is controlled by the well-known Peters Equations (32) and (33). We show below that the first state is also analytically tractable to good approximation. Though slow mergers tend to give rise to small ${e}_{\mathrm{LIGO}}$, we show in the Appendix that LIGO could in principle be sensitive to such small eccentricities of ${ \mathcal O }(0.01)$. By measuring smaller eccentricities, we can probe more of the parameter space. Therefore, we shall take ${e}_{\star }=0.01$ as the sensitivity threshold for ${e}_{\mathrm{LIGO}}$. More importantly, future GW detectors like LISA could observe such BBHs at a much earlier stage when their binary separations are several orders of magnitude larger than at the LIGO threshold. The eccentricity will then be significantly larger than it is in LIGO since it will have undergone less GW circularization (Breivik et al. 2016; Michaely & Perets 2018; Nishizawa et al. 2017).

We now estimate the eccentricity ${e}_{\mathrm{LIGO}}$ for slow mergers. We take an imaginary and isolated binary with the same masses ${m}_{\mathrm{0,1}}$ and initial separation a₁₀ as the inner binary in question. Then we set the eccentricity of this imagined isolated binary, denoted by ${\widehat{e}}_{1}$, such that it has the same lifetime as the inner binary. Then, the eccentricity of the inner binary when entering the LIGO band can be approximated by the eccentricity of the isolated binary at the time of entering the LIGO band. The reason behind this identification is that the lifetime of a slow merger is much longer than the KL timescale, and therefore we can think of the imagined isolated binary as the inner binary with fast KL oscillations averaged away. In the same way, we can think of ${\widehat{e}}_{1}$ as an averaged eccentricity of the inner binary over the initial several KL cycles. By the time that the inner binary enters the LIGO band, its KL oscillations have long ceased owing to both PN and GW effects, and therefore we can identify at this moment the eccentricity of the inner binary by the corresponding value of the imagined and isolated binary.

Putting the above words into equations, we relate ${\widehat{e}}_{1}$ to the lifetime τ of the inner binary by Equation (40),

$$\begin{eqnarray}&&\tau =\displaystyle \frac{5}{256}\displaystyle \frac{{c}^{5}{a}_{10}^{4}}{{G}^{3}{m}_{0}{m}_{1}m}{\widehat{\epsilon }}_{1}^{\ 7/2},\end{eqnarray} \tag{ 58 }$$

where ${\widehat{\epsilon }}_{1}\equiv 1-{\widehat{e}}_{1}^{\,2}$. On the other hand, we know from Section 3.2 that the merger time of the inner binary is well approximated by Equation (48),

$$\begin{eqnarray}&&\tau =\displaystyle \frac{5}{256}\displaystyle \frac{{c}^{5}{a}_{10}^{4}}{{G}^{3}{m}_{0}{m}_{1}m}{\epsilon }_{1\min }^{3},\end{eqnarray} \tag{ 59 }$$

and therefore ${\widehat{\epsilon }}_{1}={\epsilon }_{1\min }^{6/7}$, where ${\epsilon }_{1\min }$ is the minimal value of ₁ in the first several KL cycles. Then, at a later time when the KL oscillation is totally suppressed, the eccentricity can just be read from the binary separation a₁ via Equation (76), namely,

$$\begin{eqnarray}&&{e}_{1}={g}^{-1}\left[\displaystyle \frac{{a}_{1}}{{a}_{10}}g(\sqrt{1-{\epsilon }_{1\min }^{6/7}})\right].\end{eqnarray} \tag{ 60 }$$

In particular, the eccentricity of the binary at the LIGO threshold, ${e}_{\mathrm{eLIGO}}$, is

$$\begin{eqnarray}&&{e}_{\mathrm{LIGO}}={g}^{-1}\left[\displaystyle \frac{{a}_{\mathrm{LIGO}}}{{a}_{10}}g(\sqrt{1-{\epsilon }_{1\min }^{6/7}})\right],\end{eqnarray} \tag{ 61 }$$

where ${a}_{\mathrm{LIGO}}\simeq 513\,\mathrm{km}\times {(m/{M}_{\odot })}^{1/3}$ if the lower end of the LIGO frequency band is taken to be 10 Hz. Though this expression looks quite simple, it should be noted that most of the information about the initial condition of the binary is encoded in the value of ${\epsilon }_{1\min }$ in a rather complicated way.

We now can also find the eccentricity distribution f(e) when the binary leaves the tidal sphere of influence, which we define as the region over which KL dominates GWs. In Randall & Xianyu (2018), this distribution was one of two unknowns required to find the final eccentricity distribution. Here we explicitly included the PN correction so we could calculate the final ellipticity distribution directly. Nonetheless, now that we have the expression for the eccentricity at all values of a₁ for which the KL oscillation is suppressed, we immediately know that the desired eccentricity can be solved from Equation (60) with a₁ replaced by its value at the time of exiting the tidal sphere of influence (Randall & Xianyu 2018),

$$\begin{eqnarray}&&{a}_{1}={\left[{a}_{2}^{6}{R}_{m}^{5}{\left(\displaystyle \frac{m}{M}\right)}^{2}\displaystyle \frac{1}{{(1-{e}_{1}^{2})}^{6}}\right]}^{1/11},\end{eqnarray} \tag{ 62 }$$

where ${R}_{m}=2{Gm}/{c}^{2}$ is the Schwarzschild radius associated with binary mass m.

In Figure 6, we show the approximate eccentricity (60) in the $({e}_{1},{a}_{1})$-plane for our second example in Section 3.1 (namely, the slow merger in Figure 3) compared to the numerical solution (gray curve). It can be seen clearly that the numerical solution undergoes a number of KL oscillations at an early stage while the binary separation stays essentially constant (upper gray belt). Once the KL oscillation is suppressed, the standard GW circularization takes place and drives the binary to smaller eccentricity. It can be seen that our analytical estimate is a good approximation to e₁ for the circularization stage where KL oscillations are suppressed, which is essential to analytically estimate the value of ${e}_{\mathrm{LIGO}}$ when the binary enters the LIGO window.

Zoom In Zoom Out Reset image size

In the same plot we also show the observational bands of LIGO ($10\,\mathrm{Hz}\leqslant {f}_{\mathrm{GW}}\leqslant 2000$ Hz) and LISA ($0.001\,\mathrm{Hz}\,\leqslant {f}_{\mathrm{GW}}\leqslant 0.1\,\mathrm{Hz}$). Both bands are tilted toward large a₁ when e₁ gets large because the peak frequency ${f}_{\mathrm{peak}}$ of emitted GWs moves to larger harmonics as e₁ get larger (Wen 2003),

$$\begin{eqnarray}&&{f}_{\mathrm{peak}}({a}_{1},{e}_{1})\simeq \displaystyle \frac{\sqrt{{Gm}}}{\pi {[{a}_{1}(1-{e}_{1}^{2})]}^{3/2}}{(1+{e}_{1})}^{1.1954}.\end{eqnarray} \tag{ 63 }$$

Since LISA is able to probe stellar BBHs with much greater separation, it can in principle capture more information stored in eccentricity before the circularization washes it away. In the future, it is possible that a compact binary enters both LISA and LIGO windows at different stages of inspiral, and a joint analysis with both ground and space GW detectors can be very powerful in measuring the properties of inspiral binaries, and a LIGO event without a LISA counterpart can also provide us with useful information about the inspiral history (Breivik et al. 2016; Michaely & Perets 2018; Nishizawa et al. 2017).

In Figure 7 we compare the numerical solution and analytical estimate (61) for the eccentricity ${e}_{\mathrm{LIGO}}$ of the BBH at the LIGO threshold. We vary a₂ to show different strength KL oscillations, while all other parameters are taken to be the same as in the three examples in Section 3.1. As expected, estimate (61) works quite well for small eccentricity owing to a large number of KL oscillations and a clear distinction between KL domination and GW domination. At large eccentricity estimate (61) does not agree with the numerical solution as well, but it still serves as a good indication in the order-of-magnitude sense. In general, large eccentricities correspond to fast mergers, and such events are expected to be rarer than ones with smaller eccentricity in the galactic center since they happen only in the small-volume inner region that is only slowly replenished.

Zoom In Zoom Out Reset image size

As we showed in Section 3, the final eccentricity ${e}_{\mathrm{LIGO}}$ depends on a number of parameters, including the two masses of the binary $({m}_{0},{m}_{1})$, the mass of the central SMBH m₂, the semimajor axis and eccentricity of the outer orbit $({a}_{2},{e}_{2})$, the initial semimajor axis and the eccentricity of the inner binary $({a}_{10},{e}_{10})$, and the inclination I₀ of the inner orbit relative to the outer orbit.

Among these parameters, only the binary masses $({m}_{0},{m}_{1})$ are directly measurable by LIGO, and thus we do not know most of them for each event. However, these parameters can leave their impact on the final eccentricity ${e}_{\mathrm{LIGO}}$, and the hope is that we will obtain useful information about these parameters with more statistics of observations. To this end, it is important to understand how the parameters affect the distribution of ${e}_{\mathrm{LIGO}}$ and how well we know about their initial distribution. We comment on these dependencies now.

In general, most parameters enter Equation (61) for ${e}_{\mathrm{LIGO}}$ through the value of ${e}_{1\max }$, the maximal eccentricity in the first several KL cycles. Therefore, their importance to the final ${e}_{\mathrm{LIGO}}$ basically depends on how much they affect the value of ${e}_{1\max }$.

The final distribution of ${e}_{\mathrm{LIGO}}$ does not depend strongly on e₁₀ and has very mild dependence on a₁₀. For e₁₀, this is because the information about initial eccentricity is largely washed out by many KL cycles, unless e₁₀ is extremely large. However, binaries with extremely large e₁₀ are rare, because we expect the distribution of e₁₀ in a very dense environment to be flat in ${e}_{10}^{2}$, which follows from an equipartition argument first given by Jeans (1919).

On the other hand, we know very little about the a₁₀ distribution in an NC. But we do not need detailed knowledge about it either because binaries with very large a₁₀ will evaporate anyway, and very small a₁₀ is hard to form. In typical NCs, this means that we can assume that a₁₀ follows some distribution ranging from ${ \mathcal O }(0.01\,\mathrm{au})$ to ${ \mathcal O }(10\,\mathrm{au})$, and we find that the distribution of a₁₀ in this range has only a tiny impact on the final answer for ${e}_{\mathrm{LIGO}}$. Therefore, we will assume a flat distribution for a₁₀ in the following.

The most influential parameter affecting the final eccentricity distribution is the inclination I₀ because it largely determines ${e}_{1\max }$ in the absence of the PN correction. We will assume that the orientation of the inner binary follows a random uniform distribution in all possible directions—that is, we assume a flat distribution in cos I₀. This is also the assumption adopted in Antonini & Perets (2012). Eventually, binaries with low inclination will evaporate since they receive little KL oscillations and thus do not merge fast enough. We will take account of the evaporation constraint explicitly shortly.

After including the PN correction, the outer orbital parameters $({a}_{2},{e}_{2})$ become important because they affect the rate of KL oscillations, and this rate competes with the PN correction in determining ${e}_{1\max }$. The distribution of e₂ is important, because a relatively large e₂ can reduce the distance of periapsis for fixed a₂, and the binary may receive larger perturbations there. In our examples we follow Antonini & Perets (2012) and take the distribution of e₂ to be thermal, i.e., $\propto {{de}}_{2}^{2}$, although ultimately this is an assumption that should be observationally checked.

We know little about the a₂ distribution in an NC, and finding a method for observationally determining this is potentially one of the most important goals of the approach we describe in this paper since at this point it is extremely difficult to probe this distribution close to the galactic center in other ways. Examples of possible distributions as described in the literate are as follows: For less relaxed systems with a core profile for the stellar distribution, one possibility is that the distribution of a₂ follows the background. As in Antonini & Perets (2012), we take the profile to be ${a}_{2}^{2-\beta }{{da}}_{2}$ with β = 0.5 for the core model. For the fully relaxed case with a Bahcall–Wolf cusp, it is expected that the BH binaries would be more cuspy owing to mass segregation if the BH masses are much greater than the background stellar mass when the BBH is not the dominant component of the NC (Alexander & Hopman 2009). In this case we shall take the mass-segregated distribution with β = 2. But for lighter binaries with $m\sim {M}_{\odot }$, there is no segregation at work and we shall take β = 7/4 to be in accordance with the Bahcall–Wolf cusp. Throughout we restrict a₂ within the inner 0.1 pc from the central SMBH.

This leaves the most important remaining parameter as I. As discussed, we take a flat distribution initially. However, mergers compete with evaporation in determining the fate of BH binaries, which favors large inclinations for successful mergers. Hence, the final consideration we apply is comparing merger and evaporation rates, the latter of which is already discussed in the literature, and which we review in the next section.

In dense environments such as NCs, BH binaries can be destroyed by background stellar mass objects through binary–single interactions and thus "evaporate" rather than merge (Leigh et al. 2018). Therefore, to be observed by LIGO or other detectors, it is important that the merger time is shorter than the evaporation timescale. This introduces a constraint on the parameter space of the initial distribution of binaries. This constraint was taken into account numerically in previous studies (Wen 2003; Antonini & Perets 2012). To impose the constraint analytically, we use the cross section of the process $\sigma (\mathrm{binary}+\mathrm{single}\to \mathrm{three}\,\mathrm{singles})$ in conjunction with the analytical merger calculation we did above.

Intuitively, the evaporation cross section can at most be the geometrical cross section, $\sigma \sim \pi {a}_{1}^{2}$, according to which binaries with large separation evaporate more readily. The cross section also depends crucially on the velocity of the scattered object ${v}_{\star }$ relative to the binary mass center. In fact, the geometrical cross section $\pi {a}_{1}$ can be reached only when ${v}_{\star }$ is mildly greater than the rotational velocity of the binary ${v}_{b}\sim \sqrt{{Gm}/{a}_{1}}$, because the time duration of the close encounter is too short to destroy the binary if ${v}_{\star }\gg {v}_{b}$, and there is not enough energy to evaporate the binary if ${v}_{\star }\ll {v}_{b}$.

A binary is conventionally said to be "hard" if ${v}_{b}\gt {v}_{\star }$ and "soft" if ${v}_{b}\lt {v}_{\star }$. In galactic nuclei with an SMBH, the typical velocity of field stars is ${v}_{\star }\sim \sqrt{{{Gm}}_{2}/{a}_{2}}$, with m₂ and a₂ the mass of and the distance to the central SMBH, respectively. Therefore, a binary is hard if ${a}_{1}\lesssim (m/{m}_{2}){a}_{2}$. For typical masses $m\sim 10\,{M}_{\odot }$ and ${m}_{2}\sim {10}^{6}\,{M}_{\odot }$, a binary in the central region ${a}_{2}\lt {10}^{4}$ au of an NC will in general be soft if ${a}_{1}\gt 0.1$ au. Hard binaries do not evaporate. For very soft binaries, an analytical expression for the cross section of evaporation is known (Hut 1983),

$$\begin{eqnarray}&&{\sigma }_{\mathrm{evap}}=\displaystyle \frac{40\pi G}{3}\displaystyle \frac{{m}_{\star }^{2}}{m}\displaystyle \frac{a}{{v}_{\star }^{2}},\end{eqnarray} \tag{ 64 }$$

where ${m}_{\star }$ is the mass of the field star and ${v}_{\star }$ is its velocity relative to the binary. From this expression we can find the evaporation timescale ${\tau }_{\mathrm{evap}}$ of the binary,

$$\begin{eqnarray}&&{\tau }_{\mathrm{evap}}^{-1}={n}_{\star }\langle {\sigma }_{\mathrm{evap}}{v}_{\star }\rangle =\displaystyle \frac{40\pi {Ga}}{3}\displaystyle \frac{{\rho }_{\star }{m}_{\star }}{m}\langle {v}_{\star }^{-1}\rangle ,\end{eqnarray} \tag{ 65 }$$

where ${n}_{\star }$ and ${\rho }_{\star }$ are number density and mass density of the background stars, respectively, and the average $\langle \cdots \rangle$ is over the velocity distributions with a cutoff at v_b. We assume that the velocity of field stars follows the Maxwell distribution, $f(v)={e}^{-{v}^{2}/(2{\bar{v}}^{2})}{v}^{2}{dv}$. Then, the variance of the velocity distribution $\bar{v}$ can be determined in terms of densities by solving the Jeans equation. At distances a₂ close to the SMBH its mass dominates, so the velocity dispersion is approximately

$$\begin{eqnarray}&&{\bar{v}}^{2}\sim \displaystyle \frac{{{GM}}_{\mathrm{SMBH}}}{{a}_{2}}.\end{eqnarray} \tag{ 66 }$$

Additional velocity contributions (Leigh et al. 2016) due to stellar mass are subdominant for the parameters of interest here. Then, the evaporation timescale is given by

$$\begin{eqnarray}&&{\tau }_{\mathrm{evap}}=\displaystyle \frac{3m\bar{v}}{40\sqrt{2\pi }G{\rho }_{\star }{m}_{\star }{a}_{10}}.\end{eqnarray} \tag{ 67 }$$

We now impose the constraint ${\tau }_{\mathrm{evap}}\gt \tau$, where τ is given in Equation (48). In the situation considered in Randall & Xianyu (2018), where the initial binary separation is fixed and the inclination is assumed to be large enough, this evaporation constraint just reduces to a cutoff on a₂, which was one of two unknowns in Randall & Xianyu (2018), but which we can now obtain analytically.

The background density profile ${\rho }_{\star }$ makes an impact on the final result through the evaporation constraint, which means that the distribution of ${e}_{\mathrm{LIGO}}$ can be sensitive to the background profile. In the central region of an NC, we assume that the density has a power-law profile,

$$\begin{eqnarray}&&{\rho }_{\star }({a}_{2})={\rho }_{0}{\left(\displaystyle \frac{{a}_{2}}{{a}_{20}}\right)}^{-\alpha },\end{eqnarray} \tag{ 68 }$$

where ρ₀ is the density at a benchmark distance a₂₀. We use models for the density profile in the literature, though our analysis would apply to any proposed profile. As mentioned earlier, we take the Bahcall–Wolf cusp profile with $\alpha =7/4$ for fully relaxed galaxies and a core model with $\alpha =1/2$ for less relaxed galaxies.

To conclude this subsection, let us consider some limits of expression (61) to gain a bit of insight into the eccentricity ${e}_{\mathrm{LIGO}}$. Here we consider instances where the final eccentricity is small but not unobservable, i.e., ${e}_{\mathrm{LIGO}}\sim 0.01\mbox{--}0.1$. We assume that ${a}_{10}\gg {a}_{\mathrm{LIGO}}$ so that ${\epsilon }_{1\min }$ is small. With these assumptions, we can use the $e\simeq 1$ limit of g(e) in Equation (76) for the inner layer of Equation (61) and use the $e\ll 1$ limit for ${g}^{-1}$. As a result, we have

$$\begin{eqnarray}&&{e}_{\mathrm{LIGO}}\simeq 1.22{\left(\displaystyle \frac{{a}_{\mathrm{LIGO}}}{{a}_{10}}\right)}^{19/12}{\epsilon }_{1\min }^{-19/14}.\end{eqnarray} \tag{ 69 }$$

Let us now consider the limit of a highly inclined inner binary with $I\simeq 90^\circ$ so that ${\epsilon }_{1\min }$ in Equation (45) is dominated by ${{\rm{\Theta }}}_{\mathrm{PN}}$. In this case we have ${{\rm{\Theta }}}_{\mathrm{PN}}^{2}\gg {AC}$ and $A\simeq 3$ in Equation (45), and thus ${\epsilon }_{1\min }\simeq {{\rm{\Theta }}}_{\mathrm{PN}}^{2}/9$. We also use ${a}_{\mathrm{LIGO}}\simeq 513\mathrm{km}\times {(m/{M}_{\odot })}^{1/3}$, which is true for small ${e}_{\mathrm{LIGO}}$. Putting these two expressions into Equation (69), we get

$$\begin{eqnarray}{e}_{\mathrm{LIGO}} & \simeq & 0.01{\left(\displaystyle \frac{20{M}_{\odot }}{m}\right)}^{1235/252}{\left(\displaystyle \frac{{m}_{2}}{4\times {10}^{6}{M}_{\odot }}\right)}^{19/7}\\ & & \times \,{\left(\displaystyle \frac{{a}_{10}}{0.1\mathrm{au}}\right)}^{779/84}{\left(\displaystyle \frac{100\mathrm{au}}{{b}_{2}}\right)}^{57/7},\end{eqnarray} \tag{ 70 }$$

where ${b}_{2}={a}_{2}\sqrt{1-{e}_{2}^{2}}$ appears from ${W}_{\mathrm{PN}}$ in Equation (31). We learn from this expression that the binary gains more eccentricity with smaller binary mass, larger central SMBH, and shorter distance to the central SMBH. The eccentricity increases with a₁₀ because tidal force is stronger for larger binary separation.

As another interesting limit, we take the inner binary to be very close to the tertiary body but with mildly large inclination, which means that the dominant contribution to ${\epsilon }_{1\min }$ in Equation (45) is the not-so-large inclination angle. From Equation (31) we see that ${{\rm{\Theta }}}_{\mathrm{PN}}\propto {a}_{2}^{3}$, and thus the PN effect is suppressed for very small a₂, so we have ${AC}\gg {{\rm{\Theta }}}_{\mathrm{PN}}$, and we can neglect ${{\rm{\Theta }}}_{\mathrm{PN}}$ in Equation (45) and use the approximation $A\simeq 3+5{\cos }^{2}{I}_{0}$ and $C\simeq 20{\cos }^{2}{I}_{0}$, and thus ${\epsilon }_{1\min }\,\simeq 5{\cos }^{2}{I}_{0}/(3+5{\cos }^{2}{I}_{0})$. Using Equation (69) again, we get

$$\begin{eqnarray}&&{e}_{\mathrm{LIGO}}\simeq 0.01{\left(\displaystyle \frac{m}{20{M}_{\odot }}\right)}^{19/36}{\left(\displaystyle \frac{0.1\mathrm{au}}{{a}_{10}}\right)}^{19/12}{\left(\displaystyle \frac{\cos 88.8^\circ }{\cos {I}_{0}}\right)}^{19/7}.\end{eqnarray} \tag{ 71 }$$

The m and a₁₀ dependence in Equation (71) is from the ${a}_{\mathrm{LIGO}}$ dependence in Equation (69) and is in the opposite trend to Equation (70). As we will see later, the final distribution of ${e}_{\mathrm{LIGO}}$ will be anticorrelated with m because Equation (70) has stronger m-dependence than Equation (71) and also because more elliptical events are from the high-inclination limit (70) than from the small-a₂ limit (71). Clearly, the m-dependence depends on where the events arise and merits further study, as it can ultimately be very interesting in studying the BH merger environments.

We show in Figure 8 the eccentricity ${e}_{\mathrm{LIGO}}$ at the LIGO threshold in various sections of the initial parameter space, using the analytical estimate (61). From the three panels in the first row, we can see that larger-${e}_{\mathrm{LIGO}}$ regions (with darker shade) move from the left side to the top boundary as a₁₀ increases from a smaller value to a larger value. This means that harder binaries get large eccentricities when they are closer to the central SMBH compared with softer ones, and that softer binaries could reach large eccentricities when they are far away but have very large inclinations. The lower left sides of these contours correspond to the limit described by Equation (71), while the upper right sides correspond to the limit of Equation (70).

Zoom In Zoom Out Reset image size

In the second row, the left panel shows that ${e}_{\mathrm{LIGO}}$ is not sensitive to the initial eccentricity e₁₀ unless e₁₀ is very large, and the right panel shows that ${e}_{\mathrm{LIGO}}$ is not sensitive to initial separation a₁₀ either, unless a₁₀ is extremely small. Finally, the right panel shows that highly inclined binaries can be very eccentric.

We are now ready to map binaries prepared with a given distribution of initial parameters to their distribution of eccentricity ${f}_{\mathrm{LIGO}}(e)$ at the LIGO threshold. For fixed binary masses $({m}_{0},{m}_{1})$ and mass of the central SMBH m₂, the initial parameters include the binary separation a₁₀ and eccentricity e₁₀, the semimajor axis a₂ and eccentricity e₂ of the outer orbit, and the mutual inclination angle I₀. Given the analytical expression (61), we are in principle free to explore any distribution of initial parameters.

In our examples here, we fix ${m}_{2}=4\times {10}^{6}\,{M}_{\odot }$ and consider two cases for the inner binary masses: one with ${m}_{0}={m}_{1}=10\,{M}_{\odot }$ and the other with ${m}_{0}={m}_{1}={M}_{\odot }$. For the inner orbital parameters, we take ${e}_{10}=0$ since the final answer is not sensitive to this value, and we take a uniform distribution $\propto {{da}}_{10}$ ranging from 0.1 to 5 au. For the outer orbital parameters, we will take a thermal distribution for the eccentricity $\propto {{de}}_{2}^{2}$ and take the semimajor axis distribution to be ${a}_{2}^{2-\beta }{{da}}_{2}$ with $\beta =2$ for the cusp model and $\beta =0.5$ for the core model.

We can now integrate over the parameter space spanned by $({a}_{10},{a}_{2},{e}_{2},{I}_{0})$ for fixed masses $({m}_{0},{m}_{1},{m}_{2})$, where we assume that the mass parameters m₀ and m₁ can be independently determined and m₂ is a fixed parameter in a given galaxy. To find the final eccentricity distribution at the LIGO threshold, we note that our analytical Equation (61) cuts this parameter space into equal-${e}_{\mathrm{LIGO}}$ slices, and the probability of finding a binary with eccentricity between e and $e+{\rm{\Delta }}e$ is proportional to the integral over the slice between the e-surface and the $(e+{\rm{\Delta }}e)$-surface, weighted by the distribution of initial parameters as elaborated above.

On top of this, we also impose the constraint of nonevaporation, $\tau \lt {\tau }_{\mathrm{evap}}$ as discussed before, and the constraint of binary stability (Randall & Xianyu 2018),

$$\begin{eqnarray}&&{a}_{2}(1-{e}_{2})\gt {\left(\displaystyle \frac{4{m}_{2}}{m}\right)}^{1/3}{a}_{10},\end{eqnarray} \tag{ 72 }$$

which says that the binary should not be too close to the central SMBH when it reaches the periapsis, or the tidal force would overwhelm the self-gravity of the binary and disintegrate it. Comparing with Equation (9), we see that this condition ensures the weak coupling between the inner and outer orbits, and thus the stability of the hierarchical triple system.

In summary, the probability $P(e,{\rm{\Delta }}e)$ that a BH binary near an SMBH enters the LIGO window with eccentricity between e and $e+{\rm{\Delta }}e$ is given by

$$\begin{eqnarray}&&P(e,{\rm{\Delta }}e)=\displaystyle \frac{{\int }_{{e}_{\mathrm{LIGO}}\in (e,e+{\rm{\Delta }}e)}{dx}\,{f}_{\mathrm{ini}}(x)C(x)}{\int {dx}\,{f}_{\mathrm{ini}}(x)C(x)},\end{eqnarray} \tag{ 73 }$$

where $x=({a}_{10},{e}_{10},{a}_{2},{e}_{2},{I}_{0})$ denotes compactly the initial parameters, ${f}_{\mathrm{ini}}(x)$ represents the initial distributions of x, while C(x) contains the two constraints that can be expressed in terms of Heaviside step function $\theta (x)$,

$$\begin{eqnarray}&&C(x)\equiv \theta ({\tau }_{\mathrm{evap}}-\tau )\theta [{a}_{2}(1-{e}_{2})-{(4m/{m}_{2})}^{1/3}{a}_{10}].\end{eqnarray} \tag{ 74 }$$

In Figure 9 we take ${m}_{0}={m}_{1}=10\,{M}_{\odot }$, ${m}_{2}=4\times {10}^{6}\,{M}_{\odot }$ and show two layers of equal-${e}_{\mathrm{LIGO}}$ slices (in teal blue) for ${e}_{\mathrm{LIGO}}=0.01$ (upper and light) and ${e}_{\mathrm{LIGO}}=0.1$ (lower and dark) in the parameter space of $({a}_{10},{a}_{2},\cos {I}_{0})$ with ${e}_{20}=0$. Also shown are regions excluded by the evaporation constraint (yellow) and tidal disruption (gray). The evaporation constraint is computed from the cusp model (68) with ${m}_{\star }={M}_{\odot }$, $\alpha =7/4$, ${\rho }_{0}={10}^{6}\,{M}_{\odot }\,{\mathrm{pc}}^{-3}$, and ${a}_{20}=0.1$ pc.

Zoom In Zoom Out Reset image size

We perform the integration (73) for several sets of initial distributions and show the resultant probability P(e) in Figure 10. For the cusp model, we take the same background profile as the one we take for Figure 9, while the core profile corresponds to replacing $\alpha =7/4$ by $\alpha =1/2$. It is clear from the figure that the cusp profile tends to produce a more elliptic distribution than a core model when other parameters are fixed. More interestingly, lighter binaries tend to gain more eccentricity in NC than heavier binaries, which means that there is an anticorrelation between the binary mass and eccentricity in this formation channel.

Zoom In Zoom Out Reset image size

This is different from the binaries in GCs, where the mass has little impact on eccentricity distribution (Wen 2003), and is in contrast to what is claimed in Gondán et al. (2018), who considered an alternative eccentric BBH formation channel with direct two-body encounter and found that the binary mass and the eccentricity are positively correlated. Though we have yet to analyze such situations, such parameter dependence might ultimately be used to distinguish different formation scenarios.

In general, it is clear that most binaries in NCs will have small eccentricities. A careful measurement of eccentricity distribution in this range could be very important in revealing the formation of binaries. It is likely that LIGO can search only for the tail of this distribution in the not-too-small e region ($e\geqslant 0.01$; see the Appendix), which contributes 5% (24%) of all mergers in NCs in the cusp model with $m=20\,{M}_{\odot }$ ($m=2\,{M}_{\odot }$). Further, a joint observation with future GW detectors such as LISA would be a lot more powerful in measuring the distribution in the region of tiny (LIGO) eccentricities, and it is even possible to reveal the peak in the distributions if this formation channel contributes significantly to total merger events.