EJECTION AND CAPTURE DYNAMICS IN RESTRICTED THREE-BODY ENCOUNTERS

The disruption of a star by a massive black hole (BH) is one of the most spectacular examples of the tidal phenomena (Komossa & Bade 1999; Donley et al. 2002; Grupe et al. 1995). A star that wanders too close to a massive BH is torn apart by gravitational forces. Almost half the debris would escape on hyperbolic orbits, while the other half would traverse elliptic orbits and return to periapsis before producing a conspicuous flare (e.g., Rees 1988). The disruption process has been numerically investigated in detail (Evans & Kochanek 1989; Laguna et al. 1993; Ayal et al. 2000; Kobayashi et al. 2004; Guillochon et al. 2009), and the new generation of all-sky surveys are expected to detect many tidal flares (Strubbe & Quataert 2009; Lodato & Rossi 2011). Recently, a possible discovery of the onset of rapid BH accretion has been reported (Burrows et al. 2011; Zauderer et al. 2011; Levan et al. 2011; Bloom et al. 2011).

Once a star gets deeply inside the tidal radius of a BH, the tidal force dominates over the self-gravity and thermal pressure of the star. A very simplified description of the disruption process could be the encounter between a star cluster (or a cluster of point masses) and a massive BH. The simplest case consists of a binary and a massive BH in which after the tidal disruption, one star would escape to infinity, while the other could be captured by the BH. This is actually one of the leading models for the formation of hypervelocity stars (Hills 1988; Yu & Tremaine 2003). The captured stars may explain the S-stars in the Galactic center (Gould & Quillen 2003; Ginsburg & Loeb 2006; Ghez et al. 2005; Genzel et al. 2010).

Hypervelocity stars are stars with a high velocity exceeding the escape velocity of the Galaxy. After the discovery of such stars in a survey of blue stars within the Galactic halo (Brown et al. 2005; Hirsch et al. 2005; Edelmann et al. 2005), many authors have predicted the properties of hypervelocity stars (Gualandris et al. 2005; Bromley et al. 2006; Sesana et al. 2007; Perets et al. 2007; Kenyon et al. 2008; Tutukov & Fedorova 2009; Antonini et al. 2011; Zhang et al. 2010). These investigations so far have used three-body simulations or analytic methods that relied on results from three-body simulations.

The six orders of magnitude mass ratio between the Galactic center BH and the binary stars allows us to formulate the problem in the restricted three-body approximation. In a previous paper (Sari et al. 2010, hereafter SKR), we have shown that the approximation is efficient and useful to understand how binary stars behave at the tidal breakup when the binary's center of mass approaches the BH in a parabolic orbit. In this paper, we generalize the approximation for orbits with arbitrary eccentricity. This enables us to give a complete picture of the ejection and capture process. We also provide the ejection and capture probability distributions that can be simply rescaled in terms of binary masses, their initial separation, and the binary-to-black hole mass ratio when applied to a specific system. Our method is computationally more efficient than full three-body simulations, and it is easier to grasp the nature of the tidal interaction.

In Section 2, we outline the restricted three-body approximation. In Section 3, we evaluate how much energy each member gains or loses at the tidal encounter and we discuss how the energy change evaluated for a parabolic orbit can be used to study non-parabolic orbit cases, and in Section 4, we give qualitative discussion on the ejection and capture preferences. In Section 5, we study high-velocity encounters. In Section 6, the numerical results are discussed. In Section 7, we use our results to describe the capture process of Triton around Neptune. Finally, in Section 8, we summarize the results.

The equation of motion for each of the binary members is given by

$$\begin{equation} \ddot{\bf r}_1=-{{\it G M} \over r_1^3}{\bf r}_1 + {Gm_2 \over |{\bf r}_1- {\bf r}_2|^3} ({\bf r}_2- {\bf r}_1), \end{equation} \tag{ 1 }$$

$$\begin{equation} \ddot{\bf r}_2=-{{\it G M} \over r_2^3}{\bf r}_2 - {Gm_1 \over |{\bf r}_1-{\bf r}_2|^3} ({\bf r}_2- {\bf r}_1), \end{equation} \tag{ 2 }$$

where r₁ and r₂ are the respective distance from the massive point mass with M. We will call the point mass the BH, though the binary is assumed to travel well outside the event horizon and our results can be applied to any systems which include a Newtonian massive point mass. The equation for the distance between the two r ≡ r₂ − r₁ is

$$\begin{equation} \ddot{\bf r}=-{{\it G M} \over r_2^3} {\bf r}_2+ {{\it G M} \over r_1^3}{\bf r}_1 - {Gm \over r^3} {\bf r}, \end{equation} \tag{ 3 }$$

where m = m₁ + m₂ ≪ M. We assume that the two masses are much closer to each other, and to the trajectory of the center of mass of the binary r_m, than each of them to the BH. Both energy and orbit obtained under the approximation are fairly accurate except for a part of the orbit just around the periapsis passage (see SKR for the details).

Linearizing the first two terms of Equation (3) around the center of mass orbit r_m, we find that the zero orders cancel out. Then, rescaling the distance between the bodies by (m/M)^1/3r_p and the time by $\sqrt{\vphantom{A^A}\smash{\hbox{${r_{\rm p}^3/{\it GM}}$}}}$, where r_p is the distance of the closest approach between the center of mass of the binary and the BH, we can re-write Equation (3) in terms of the dimensionless variables: ${\boldsymbol{\eta }}\equiv (M/m)^{1/3} ({\bf r}/r_{\rm p})$ and t:

$$\begin{equation} \ddot{\boldsymbol{\eta }}=\left(r_{\rm p} \over r_{\rm m}\right)^3 [ - {\boldsymbol{\eta }}+3 ({\boldsymbol{\eta }} {\hat{\bf r}}_{\rm m}){\hat{\bf r}}_{\rm m}] - { {\boldsymbol{\eta }} \over |{\boldsymbol{\eta }}|^3}, \end{equation} \tag{ 4 }$$

where ${\hat{\bf r}}_{\rm m}$ is a unit vector pointing the center of mass of the binary. We define the orbit of the center of mass to be a conic orbit r_m/r_p = (1 + e)/(1 + ecos f), where e is the eccentricity and the true anomaly f  is the angle from the point of closest approach. Since ${\hat{\bf r}}_{\rm m}=(\cos f, \sin f, 0)$, and we set ${\boldsymbol{\eta }}=(x,y,z)$, explicit equations in terms of dimensionless Cartesian coordinates read

$$\begin{eqnarray} \ddot{x} &=& \frac{(1+e \cos f)^3}{ (1+e)^3} \left[ - x+3 (x \cos f+y \sin f)\cos f \right]\nonumber\\ && - {x \over (x^2+y^2+z^2)^{3/2} }\,, \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} \ddot{y} &=& \frac{(1+e \cos f)^3}{ (1+e)^3} \left[ -y+3 (x \cos f+y \sin f)\sin f \right]\nonumber\\ && - {y \over (x^2+y^2+z^2)^{3/2} }\,, \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} \ddot{z} &=& - \frac{(1+e \cos f)^3}{ (1+e)^3} z - {z \over (x^2+y^2+z^2)^{3/2} }\,,\qquad\ \ \ \end{eqnarray} \tag{ 7 }$$

where the eccentricity

$$\begin{equation} e=1+\frac{2r_{\rm p}E}{\it GMm} \end{equation} \tag{ 8 }$$

is related to the energy of the center of mass which is given by

$$\begin{eqnarray} E &=& \frac{m }{2}|\dot{\bf r}_{\rm m}|^2-\frac{\it GMm}{r_{\rm m}}. \end{eqnarray} \tag{ 9 }$$

Using the dimensionless time, the conservation of the angular momentum can be expressed as

$$\begin{equation} \dot{f}=\left(1+e \right)^{-3/2}\left(1+e\cos f \right)^2. \end{equation} \tag{ 10 }$$

Analytically, one has relations between r_m and t through a parameter which are given by (e.g., Landau & Lifshitz 1976)

$$\begin{eqnarray} E<0 \quad r_{\rm m}/r_{\rm p}&=&(1-e)^{-1}\left(1-e \cos \xi \right), \nonumber\\ t&=&(1-e)^{-3/2} \left(\xi - e \sin \xi \right),\quad\ \ \ \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} \ E=0 \quad r_{\rm m}/r_{\rm p}&=&(1+\xi ^2 ), \quad t=\sqrt{2} (\xi +\xi ^3/3), \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} E>0 \quad r_{\rm m}/r_{\rm p}&=& (e-1)^{-1} \left(e \cosh \xi -1 \right), \nonumber\\ t&=&(e-1)^{-3/2} \left(e \sinh \xi -\xi \right),\quad\ \end{eqnarray} \tag{ 13 }$$

where the closest approach r_m = r_p happens at t = 0.

We are interested in the fate of stars in a binary, following its encounter with a massive BH. In order to study the ejection and capture process, we evaluate the energies of the stars as functions of time. When the binary is at a large distance from the BH, the binary members rotate around their center of mass which gradually accelerates toward the BH. The specific self-gravity energy of the binary is about −v²₀ ≡ −Gm/a. Analytic arguments (SKR) suggest that at the tidal breakup one member gets additional energy of the order of v_mv₀, where v_m is the velocity of the center of mass at the tidal radius r_t = (M/m)^1/3a. If the binary approaches the BH with negligible orbital energy, the velocity is v_m = (GM/r_t)^1/2 = v₀(M/m)^1/3. The additional energy is larger than the self-gravity energy by a factor of (M/m)^1/3 ≫ 1. Therefore, we will neglect the self-gravity term in the following energy estimates. This treatment is valid as long as the binary is injected into the orbit r_m at a radius much larger than the tidal radius.

The energy of one binary member m_i is given by

$$\begin{equation} E_i=\frac{m_i}{2}|\dot{\bf r}_i|^2-\frac{\it GMm_i}{r_i}. \end{equation} \tag{ 14 }$$

Linearizing the kinetic and potential energy terms around the orbit of the center of mass r_m and using the initial energy I_i ≡ (m_i/m)E, we obtain

$$\begin{eqnarray} E_i&=& I_i+ \Delta E_i, \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} \Delta E_i &\equiv & m_i \dot{\bf r}_{\rm m}\left(\dot{\bf r}_i-\dot{\bf r}_{\rm m} \right) +\frac{\it GMm_i}{r_{\rm m}^3}{\bf r}_{\rm m}\left({\bf r}_i-{\bf r}_{\rm m} \right), \end{eqnarray} \tag{ 16 }$$

Since in our limit the total energy of the system is E, considering ΔE₂ = −ΔE₁, we get

$$\begin{equation} \Delta E_2= \frac{m_1m_2}{m} \left(\dot{\bf r}_{\rm m} \dot{\bf r} +\frac{\it GM}{r_{\rm m}^3}{\bf r}_{\rm m} {\bf r} \right). \end{equation} \tag{ 17 }$$

Using our rescaled variables, the additional energy is given by

$$\begin{eqnarray} \Delta E_2 &=& -\Delta E_1= \frac{Gm_1m_2}{a}\left(\frac{M}{m} \right)^{1/3} \Delta \bar{E},\quad\ \ \ \ \end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray} \Delta \bar{E} &\equiv & D^{-1} \left[ \frac{-\dot{x}\sin f+\dot{y}(e+\cos f)}{\sqrt{1+e}}\right. \nonumber\\ &&\left.+\,\frac{(1+e\cos f)^2}{(1+e)^2}\left(x\cos f+y\sin f \right) \right], \end{eqnarray} \tag{ 19 }$$

where D = r_p/r_t is the penetration factor which is useful to characterize the tidal encounter. Once the binary dissolves, $\Delta \bar{E}$ becomes a constant because the body is eventually moving only under the conservative force of the BH. Hereafter, the energy change $\Delta \bar{E}$ means the constant value after the disruption, otherwise we specify it. The equation of motion (4) indicates that the negative of a solution r = r(t) is also a solution. The energies ΔE_i are also linear in the coordinates. Therefore, another binary starting with a phase difference π will have the same additional energy in absolute value but opposite in sign. A uniform distribution in the binary phase implies that, when the binary is disrupted, each body has a 50% chance of gaining energy (and a 50% chance of losing energy).

As we have discussed, the typical energy change is larger than the self-gravity energy by a factor of ∼(M/m)^1/3, it is of the order of (Gm₁m₂/a)(M/m)^1/3. Then, the dimensionless quantity $\Delta \bar{E}$ is an order-of-unity constant after the disruption. Its exact value depends on orbital parameters, but for qualitative discussion we just need to know that $\Delta \bar{E}$ is about unity. Later, we will numerically show that $\Delta \bar{E}$ is an order of unity in the relevant parameter regime,⁶ and numerical values will be used to estimate the ejection and capture probabilities.

Rescaling energies by the typical value of the energy change, the energies of the binary members after the disruption are given by

$$\begin{equation} \bar{E}_1 = \bar{I}_1-\Delta \bar{E}, \ \ \bar{E}_2 = \bar{I}_2+\Delta \bar{E}, \end{equation} \tag{ 20 }$$

where bar denotes energy scaled by (Gm₁m₂/a)(M/m)^1/3. An interesting outcome of the encounter between a binary system and a massive BH is the "three-body exchange reaction" (Heggie 1975; Hills 1975) where one member of the binary is expelled and its place is taken by the BH, i.e., one binary member is captured by the BH and the other is ejected to infinity. In order for a member m_i to escape from the BH, the initial binding energy should be smaller than the energy gain: $|\bar{I}_i| < |\Delta \bar{E}| \sim 1$. The same condition is required when a member of the binary in a hyperbolic orbit loses energy and is bound around the BH. Therefore, when we discuss the ejection or capture process associated with a massive BH, the absolute value of the initial energy should be comparable or less than unit: $|\bar{I}_i|\lesssim 1$.

Since the energy, penetration factor (periapsis radius), and eccentricity are related by Equation (8) or equivalently

$$\begin{equation} e=1+2D\bar{E}\left(\frac{m}{M} \right)^{1/3} \left(\frac{m_1}{m} \right) \left(\frac{m_2}{m} \right), \end{equation} \tag{ 21 }$$

only two of them are independent parameters used to describe the binary orbit. Considering $|\bar{I}_i|\lesssim 1$ together with the mass ratio (m/M)^1/3 ≪ 1 and the tidal disruption condition D ≲ 1, the eccentricity should be almost unity |1 − e| ≲ D(m/M)^1/3(m_par/m) for a member m_i to be ejected or captured where m_par is the mass of the partner (m_par = m₂ for i = 1 and m₁ for i = 2). If we use the semimajor axis r_a ≡ r_tD/(1 − e), the condition can be rewritten as |r_a/r_t| ≳ (M/m)^1/3(m/m_par) where r_a is negative for hyperbolic orbits.

Such orbits differ very little from parabolic orbits with the same periapsis distance, especially around the tidal radius and inside it. Therefore, the energy change $\Delta \bar{E}$ is expected to be almost identical to that for the parabolic case. As long as we study the exchange reaction, we can approximate $\Delta \bar{E}$ by the parabolic results $\Delta \bar{E}_{e=1}$. However, e ∼ 1 does not necessarily mean $|\bar{E}|\ll 1$. In general, we need to take into account the offset of the final energy due to the non-zero initial energy, which would affect the ejection and capture probabilities. The final energies are approximately given by

$$\begin{equation} \bar{E}_1 = \bar{I}_1-\Delta \bar{E}_{e=1}, \ \ \bar{E}_2 = \bar{I}_2+\Delta \bar{E}_{e=1}. \end{equation} \tag{ 22 }$$

We here consider a simple question: which member is ejected or captured if an unequal-mass binary is tidally disrupted by a massive BH? If $\bar{E}>0$ (hyperbolic orbits), then one binary member could lose energy and get captured by the BH, while the other flies away with a larger energy. Assuming a uniform distribution in the binary phase, each member has a 50% chance of losing energy (and gaining energy). However, since the lighter one (the secondary) has a smaller initial energy, it is preferentially captured and the heavier one (the primary) has more chance to be ejected.

For elliptical orbits, by considering a plausible semimajor axis r_a, we can obtain tighter constraints on the eccentricity and energy, compared to the requirements from the exchange reaction. This is particularly relevant for studies of hypervelocity stars. If r_a is around the radius of influence of the BH r_h ∼ GM/σ² where σ is the local stellar velocity dispersion, for the Galactic center, it is about r_a ∼ a few parsecs ∼10⁵r_t for a ∼ several solar radii. Then, we get 1 − e = D(r_t/r_a) ≲ 10⁻⁵ and $|\bar{E}| \sim (r_t/r_a) (M/m)^{1/3}(m/m_1)(m/m_2)\sim 10^{-3}$. Our previous estimates based on parabolic orbits are appropriate to study the production of hypervelocity stars for which an equal ejection chance is expected (SKR). When the semimajor axis is as small as r_a ∼ (M/m)^1/3(m/m_par)r_t, the initial energy $|\bar{I}_i|$ would be of the order of unity, as we have discussed, and affect the ejection preference. Since the secondary has less negative initial energy, it is preferentially ejected.

Recently, Antonini et al. (2011) performed N-body simulations of unequal-mass binaries with m₁ = 6 M_☉, m₂ = 1 or 3 M_☉, and a = 0.1 AU in elliptical orbits around a supermassive BH M = 4 × 10⁶ M_☉. They find that the initial distance of the binary from the central BH plays a fundamental role in determining which member is ejected: for a large initial distance d = 0.1 pc, or equivalently r_a ∼ 3 × 10³r_t, the ejection probability is almost independent of the stellar mass, while for d = 0.01 pc or r_a ∼ 3 × 10²r_t, the lighter star is preferentially ejected. Considering that the ejection probability significantly decreases if r_a becomes smaller than ∼(M/m)^1/3(m/m_par)r_t ∼ 80(m/m_par)r_t, these results are consistent with our analysis.

These ejection preferences for hyperbolic and elliptic orbits are naturally understood if we consider a large mass ratio for the binary members. The energy of the primary practically does not change at the tidal encounter. Whether it is ejected or captured after the tidal breakup simply depends on the initial energy ∼E, while the secondary might have a chance to make a transition between bound and unbound orbits around the BH (Bromley et al. 2006). In the large mass ratio limit, the exchange reaction condition (i.e., the transition condition for the secondary) is |E| ≲ (Gm²/a)(M/m)^1/3 or equivalently |r_a| ≳ (M/m)^2/3a.

If a binary has a large orbital energy $\bar{E}\gg 1$, then both members are ejected after the BH encounter as a binary system or two independent objects.⁷ Although the high-energy regime is not important in the context of the three-body exchange reaction, we discuss the regime to clarify the parameter dependence of the numerical results in the next section. A high orbital energy $\bar{E}\gg (M/m)^{1/3}(m/m_1)(m/m_2)$ affects the velocity of their center of mass at the encounter $v_m \sim (E/m+{\it GM}/r_{\rm m})^{1/2}=(M/m)^{1/3}v_0\sqrt{(e-1)/2D+r_t/r_{\rm m}}$. Then, the tidal disruption radius (i.e., where a binary is disrupted) can be defined in three different ways. We here order them from a large to small radius. (1) Relative acceleration: the radius at which the BH tidal force becomes comparable to the mutual gravity of the binary. This is r_t. (2) Relative velocity: the radius at which the tidal force induces the relative velocity between the binary members comparable to the binary escape velocity v₀. (3) Relative position: the radius at which the difference in position increases by more than the initial binary separation.

The duration that the center of mass is around r_m is of the order of Δt ∼ r_m/v_m. During this period, the tidal acceleration of the relative motion of the binary members by the BH is of the order of A ∼ GMa/r³_m. The two radii (2) and (3) can be estimated from two conditions: Δv = AΔt ∼ v₀ and Δx = AΔt² ∼ a, provided that the duration of the encounter is comparable to or shorter than the binary rotation timescale: Δt ≲ a/v₀. If the energy is high $(e-1)/D=2\bar{E}(m/M)^{1/3}(m_1/m)(m_2/m)\gg 1$, then these conditions give r_m = r_tD^1/4/(e − 1)^1/4 and r_tD/(e − 1), respectively. Since they should be larger than the periapsis distance, only the cases that satisfy D ≲ (e − 1)^−1/3 for the radius (2) or e ≲ 2 for the radius (3) lead to the disruption. The radius (3) is basically the place at which the orbit of the center of mass makes its turn (i.e., the periapsis). If the energy is low (e − 1)/D ≪ 1, then all the estimates give the original tidal radius r_t.

When we discuss the energy change ΔE at the tidal encounter, there are two important points which we should emphasize. First, the energy of each of the binary members in the BH frame changes only due to the mutual force between the binary members. Second, most of the work done by one member on the other, which is ΔE, is done outside the tidal radius r_t. The mutual force is of the order of Gm₁m₂/a². During the binary rotation timescale a/v₀, the force acts over a length ∼(v_m/v₀)a in the BH frame. Therefore, the work is W ∼ (m₁m₂/m)v_mv₀. Since the direction of the mutual force changes with the binary rotation, ΔE(t) oscillates with the amplitude of W. When the binary is disrupted, ΔE becomes a constant value which is basically determined by the binary phase at the disruption. Then, we might expect that the final value of ΔE is a sinusoidal function of the binary phase for circular binaries. As we will see later, this is actually the case for the high-energy encounters. Even with the largest estimate of the tidal radius (i.e., r_t), the duration of the encounter Δt ∼ r_t/v_m is shorter than the binary rotation timescale by a factor of ${\sim} \sqrt{D/(e-1)}\ll 1$. The work during the encounter is negligible compared to the work W which has been done outside r_t. On the other hand, in the low-energy regime, the duration of the encounter is comparable to the binary rotation timescale. Considering that at the encounter the orbits of the members in the comoving frame of their center of mass should be significantly deformed from the original orbits (e.g., circular orbits) before they finally break up, the work during the encounter could induce deviation of ΔE(ϕ) from a simple sinusoidal function. However, the typical value is still expected to be about ΔE ∼ (m₁m₂/m)v₀v_m.

In both the low- and high-energy regime, the typical energy change is given in a dimensionless form by

$$\begin{eqnarray} \Delta \bar{E} &\sim & \sqrt{1+\bar{E}\left(\frac{m}{M} \right)^{1/3}\left(\frac{m_1}{m} \right)\left(\frac{m_2}{m} \right)}\nonumber\\ &=& \sqrt{1+\frac{e-1}{2D}}=\sqrt{1-\frac{r_t}{2r_a}}, \end{eqnarray} \tag{ 23 }$$

where we have assumed r_m = r_t to estimate v_m. In the high-energy regime, the disruption might happen at a smaller radius, but v_m is determined by the orbital energy and it is insensitive to the choice of r_m. When (e − 1)/D ≫ 1, the energy change becomes much larger than unity. However, the energy gain is not significant compared to the original energy $\bar{E}$, and one finds that the tidal encounter is not an efficient acceleration process anymore.

In the high-energy regime, the energy change Equation (19) can be evaluated by assuming that a binary is abruptly disrupted at the tidal radius r_t, since the work during the tidal encounter is negligible. For a circular coplanar binary: (x, y) = D⁻¹(cos ϕ_t, sin ϕ_t) and $(\dot{x},\dot{y})=\pm D^{1/2}(-\sin \phi _t,\cos \phi _t)$, we obtain

$$\begin{eqnarray} \Delta \bar{E}&=&\left(1\pm \frac{1}{\sqrt{D(1+e)}} \right)\sin f_t \sin \phi _t \nonumber\\ &&+\left(1\pm \frac{1+e/\cos f_t}{\sqrt{D(1+e)}} \right)\cos f_t \cos \phi _t, \end{eqnarray} \tag{ 24 }$$

where ϕ_t is the binary phase at the tidal radius, f_t is the negative value solution of 1 + ecos f_t = (1 + e)D, and the signature indicates a prograde (+) or retrograde (−) orbit. This is a sinusoidal function of the binary phase as we expected, and the square of its amplitude is $3-r_t/r_a\pm 2\sqrt{(1+e)D}$, which is larger for prograde orbits and the difference between prograde and retrograde orbits becomes smaller for deep penetrators D ≪ 1, because in this limit the binary center of mass approaches the BH in an almost radial fashion.

If the disruption is abrupt, then the ejection and capture preference could be roughly illustrated in terms of velocity (e.g., Morbidelli 2006; Agnor & Hamilton 2006). The binary members rotate around their center of mass, such that their own motion is half of the time with and half of the time against, the motion of the center of mass $\dot{\bf r}_{\rm m}$. The net velocity of the members relative to the BH is accordingly increased or reduced. Since the secondary has a higher rotation velocity, it has more of a chance that the net velocity exceeds or drops below the escape velocity from the BH. Then, it is preferentially ejected to infinity or captured in a bound orbit. However, for the full discussion of the process, we also need to take into account the variation in the escape velocity or the variation in the potential. The displacement of order a in the position of each member of the binary, at a distance of about r_t from the BH, results in a change in gravitational energy of GMa/r²_t ∼ v₀²(M/m)^1/3; this is comparable to the variation in the kinetic energy. As we have done, it should be easier to discuss the overall effect in the energy domain. In our formula, the energy change (17) includes both the variation of kinetic energy and potential energy. For prograde orbits, the kinetic and potential terms cooperate and the net energy change is larger, the member on the "outside track" is expected to be ejected and its partner is captured (the "outside track" could be well defined, especially when the orbital energy is large because the duration of the encounter is much shorter than the binary rotation timescale). On the other hand, for retrograde orbits the variation in the gravitational energy would counteract that in the kinetic energy.

In this paper, we focus on results for circular coplanar binaries, though our formulae can be used to study the evolution of a binary with arbitrary orbital parameters. The orbit of a binary is assumed to be initially circular in the comoving frame of the binary center of mass. The center of mass of the binary is in a prograde or retrograde orbit around the BH (see SKR for the details of the numerical setup).

For M/m ≫ 1, the problem can be reduced to the motion of a single particle in a time-dependent potential ("the restricted three-body approximation") described by Equations (5)–(7) and (10). The energy change, Equation (19), depends only on the penetration factor D, the eccentricity e, the initial binary phase ϕ, and the binary rotation direction. As we have shown, when a binary member is captured by the BH and the other is ejected, we can further reduce the number of the parameters by assuming e = 1 to approximate the additional energy. The effect of the eccentricity e ≠ 1 is taken into account through the non-zero initial orbital energy of the center of mass. This method (22) will be called "the parabolic approximation." For a large orbital energy $\bar{E}\gg (M/m)^{1/3}(m/m_1)(m/m_2)$, "the sudden disruption approximation" (24) would become valid, but the three-body exchange reaction does not take place in this regime. In all the numerical codes, the time evolution of objects are evaluated by using a fourth-order Runge–Kutta integration scheme.

6.1. Energy Change and Probability Distributions

We have tested the restricted three-body approximation against the full three-body simulations of a binary evolving around a massive object (SKR). The full three-body orbit is accurately reproduced by the approximation equations, and the energy change, for example, differs at a 0.1% level when the exchange reaction happens. The comparison of the energy change between the full three-body and restricted three-body results is shown in Figure 1. Since they are in excellent agreement, in the following discussion, we will use the restricted three-body approximation to test the parabolic and the sudden disruption approximation.

Zoom In Zoom Out Reset image size

In Figure 1, one could notice the symmetry $\Delta \bar{E}(\phi +\pi)=- \Delta \bar{E}(\phi)$. The top panel shows that the energy change for |e − 1|/D = 0.1 is very similar to the parabolic case with the same D, especially if we take into account the phase shifts.⁸ When the three-body exchange reaction happens, the value |e − 1|/D ≲ (m/M)^1/3(m_par/m) ≪ 1 should be very small, in such a case non-parabolic results tend toward a perfect overlapping with the parabolic results. In the bottom panel, the energy change in the high-energy regime, (e − 1)/D = 100, is shown. The sudden disruption approximation well reproduces the three-body results, except for some notable spikes.

Most binaries are disrupted at the tidal encounter with a BH. However, there are always finite-phase regions where binaries survive (SKR). The narrow gaps in the top panel of Figure 1 correspond to such regions. Although it is not evident in the bottom panel, very narrow gaps exist just in the middle of the spikes. The high-energy regime is usually realized with a deep penetrator D ≪ 1. In such a case, we can ignore the self-gravity term r/r³ in Equations (5)–(7) and free solutions are obtained. Actually, one of the free solutions, (x, y)∝(− sin f, e + cos f), which corresponds to the case that binary members have the same trajectory but are slightly separated in time, dominates around the periapsis passage (SKR). Then, the binary is once disrupted at the tidal encounter, but after the periapsis passage at t = 0, they come close to each other. If we fine tune the initial binary phase, they form a binary again. This produces the narrow gaps at the spikes. If the binary phase is slightly different from the fine-tuned values, they almost form a binary, but they eventually break up. The additional work at t > 0 due to the mutual forces between the binary members produces the notable spikes.

The top panel of Figure 2 shows the phase-averaged absolute value of the energy change $\langle |\Delta \bar{E}| \rangle$ as a function of the penetration factor D. The average is taken over the phase space where binaries are disrupted. We alternatively fix the orbital eccentricity or the orbital energy $\bar{E} \propto (e-1)/D$. For a given eccentricity, a smaller D corresponds to a larger orbital energy, while for a given orbital energy, it corresponds to e → 1. For (e − 1)/D ≲ 1, the energy change remains of the order of unity as we expect from Equation (23) (see the black and green lines and the green crosses for the whole disruption range, and the red lines and red crosses for D ≳ 10⁻²). For (e − 1)/D ≳ 1, instead, $\left\langle |\Delta \bar{E}| \right\rangle$ increases toward smaller D (the red lines and red crosses for D ≲ 10⁻²). Correspondingly, the disruption probability (bottom panel) becomes almost 100% for (e − 1)/D ≳ 1. On the other hand, for the fixed energies (e − 1)/D = 0 and 0.1, the disruption chance is about 80% even in the deep penetration limit D ≪ 1 (the black and green lines). Note that in the sudden disruption approximation the disruption probability is 100% by definition, and we do not show it in the bottom panel. For our choice of (e − 1)/D = 0.1, the behavior of $\left\langle |\Delta \bar{E}| \right\rangle$ always remains very similar to the parabolic orbit case, even when we use the sudden disruption approximation especially for small D (the green crosses in the top panel). The same applies to the disruption probability function, but the sudden disruption approximation overestimates the disruption probability. As it is well known, the retrograde binaries tend to be more stable against the tidal encounter.

Zoom In Zoom Out Reset image size

For a given initial energy of a binary member $\bar{I}=(m_i/m)\bar{E}$ (in the following discussion we drop the subscript i of the initial energy for simplicity), considering the symmetry of $\Delta \bar{E}(\phi)$, the probability $P_{{\rm eje}}(\bar{I})$ that the member is ejected after the tidal encounter is determined by the fraction of the binary phase region [0, 2π] that satisfies $\Delta \bar{E}(\phi) < \bar{I}$, while the capture probability $P_{{\rm cap}}(\bar{I})$ is determined by the fraction satisfying $\Delta \bar{E}(\phi) > \bar{I}$. The sum is equal to the disruption probability: P_dis = P_eje + P_cap. Since $P_{{\rm eje}}(-\bar{I})=P_{{\rm cap}}(\bar{I})$, once we evaluate the ejection probability $P_{{\rm eje}}(\bar{I})$ for negative and zero energy, the distribution for positive energy $P_{{\rm eje}}(\bar{I})=2P_{{\rm eje}}(0)-P_{{\rm eje}}(-\bar{I})$ and the capture probability for any energy $P_{{\rm cap}}(\bar{I})=P_{{\rm eje}}(-\bar{I})$ are also obtained. The disruption probability is P_dis = 2P_eje(0).

We show in Figure 3 how the ejection probability of a binary member behaves for different penetration factors and binary orientations. Since the orbital energy and the eccentricity are not independent parameters when D is fixed, $e-1 \propto \bar{E}$, we need to assume a smaller eccentricity, in principle, for a larger negative orbital energy $\bar{I}$ (or larger negative $\bar{E}$) to estimate the energy change $\Delta \bar{E}(\phi)$, the disruption probability, and the ejection probability. However, the energy change and the disruption probability are not so sensitive to the eccentricity, we thus evaluate them by using parabolic orbit results (the parabolic approximation) when the three-body exchange reaction happens. In the figure, the parabolic approximation calculations (the lines) are in a good agreement with the restricted three-body calculations (the circles). For prograde orbits, the energy change $\Delta \bar{E} (\phi)$ is more sensitive to the phase ϕ around zero points, only very narrow phase regions satisfy $|\Delta \bar{E}|\lesssim 1$ (see Figure 1, top panel), then the ejection probabilities have a plateau around $\bar{I}=0$. The long low energy tail of the prograde D = 0.1 case (the green solid line) reflects the fact that the energy gain would be quite large for a narrow binary phase region for this case. The ejection (capture) probability for $-\infty < \bar{I} < \infty$ is a monotonically increasing (decreasing) function of the energy $\bar{I}$. The probabilities rapidly change around $|\bar{I}| \sim 1$ and their function form depends mainly on D and the binary rotation direction.

Zoom In Zoom Out Reset image size

6.2. Probability Distributions for Various Semimajor Axes

When we discuss an actual astrophysical system, it is more physically intuitive to use the semimajor axis rather than the energy $\bar{I}$ or $\bar{E}$. By specifying the semimajor axis and the mass ratios, the ejection (or capture) probability and the energy after the disruption can be evaluated. The orbital energy $\bar{E}$ is given by

$$\begin{eqnarray} &&\bar{E}=-2.7 \left(\frac{r_a/r_t}{100} \right)^{-1} \left(\frac{M/m}{10^6} \right)^{1/3} \left(\frac{m_1/m}{3/4} \right)^{-1} \left(\frac{m_2/m}{1/4} \right)^{-1},\nonumber\\ \end{eqnarray} \tag{ 25 }$$

where M/m = 10⁶ and m₁/m₂ = 3 will be assumed in this section.

Figure 4 shows the characteristics of the ejection process when a binary approaches the BH in elliptic orbits. To evaluate the ejection probability (the top panel), we have not distinguished which member is ejected. Since both members are never ejected together when $\bar{E}< 0$, it is just the sum of the ejection probabilities of the two members. For parabolic orbits, one of the members is always ejected if the binary is disrupted, then the black solid line indicates the disruption probability. As r_a/r_t becomes smaller, the deviation of the orbit from the parabolic one becomes larger. The effect of the non-zero orbital energy is expected to become significant around r_a/r_t ∼ (M/m)^1/3(m/m₂) ∼ 400, and the ejection probability rapidly decreases. The middle panel indicates which member in a binary is ejected more frequently. There is no preference in the parabolic case. However, the secondary is preferentially ejected for small r_a/r_t. The phase-averaged ejection energy is shown in the bottom panel. For parabolic orbits this is equivalent to the phase-averaged energy change $\left\langle |\Delta \bar{E}| \right\rangle$.

Zoom In Zoom Out Reset image size

If we consider $P_{{\rm cap}}(\bar{I})=P_{{\rm eje}}(-\bar{I})$, it is possible to interpret Figure 4 as the capture probability and capture preference for the hyperbolic orbit cases with the same semimajor axes in absolute value. The phase-averaged energy (the bottom panel) gives the absolute values of the averaged capture energies.

Over 150 satellites are orbiting the giant planets in the solar system. About one-third of these are classified as regular, with nearly circular and planar orbits. The majority of the satellites, however, are irregular ones which are more distant from their planet and typically have larger eccentricities and/or inclination. Interestingly, a large fraction of the irregular satellites orbit their planet in the retrograde direction (Jewitt & Haghighipour 2007). Triton, Neptune's largest moon, is among them. Because of the retrograde orbit and composition similar to Pluto's, Triton is thought to have been captured from the Kuiper belt. Recently, Agnor & Hamilton (2006) demonstrated that the gravitational encounter between Neptune and a binary system, which had included Triton as a member, is an effective mechanism to capture Triton. Since the mass ratio between Neptune m_N ∼ 17.2 M_⊕ and Triton m_T ∼ 3.58 × 10⁻³ M_⊕ is reasonably large m_N/m_T ∼ 4800, we here revisit the capture process from the point of view of the restricted three-body problem.

We consider an encounter between Neptune and a binary system with m₁ = m_T, m₂ = 0.1m_T, and a = 20R_T. The binary is assumed to be in a prograde hyperbolic orbit with the periapsis r_p = 8R_N where R_T ∼ 1.35 × 10³ km and R_N ∼ 2.46 × 10⁴ km are the radii of Triton and Neptune, respectively. This set of the parameters is identical to what Agnor & Hamilton (2006) have assumed to obtain their Figure 2. Since the tidal radius is about r_t ∼ 18R_N, the penetration factor for the orbit is D ∼ 0.45. For a given encounter velocity at infinity v_∞, we evaluate the binary disruption chance and the capture probability of the binary members. The velocity at infinity v_∞ = 1 km s⁻¹ is related to the semimajor axis of the hyperbolic orbit as

$$\begin{eqnarray} |r_a/r_t| &=& 16 \left(\frac{v_\infty }{1 \,{\rm km\,s^{-1}}} \right)^{-2} \left(\frac{a}{20R_T} \right)^{-1} \left(\frac{m}{1.1 m_T} \right) \nonumber\\ &&\times \left(\frac{M/m}{4.37 \times 10^3} \right)^{2/3}, \end{eqnarray} \tag{ 26 }$$

which is about the critical value |r_a/r_t| ≳ (M/m)^1/3 ∼ 16 for the three-body exchange reaction. The numerical results are shown in Figure 5, and the relevant case (the solid lines in the top panel) is in good agreement with Figure 2 in Agnor & Hamilton (2006). Since the disruption chance depends mainly on D, for the fixed D, it is constant ∼97%. For high velocities, the probability for the primary drops sharply, while it is still ∼50% for the secondary because of its smaller initial energy $\bar{I}$. If we assume a larger secondary mass m₂ = m_T/2, then the energy change ΔE at the encounter would be larger. The 50% capture probability of the primary now extends to higher velocities, while the capture rate of the secondary is similar because the primary mass is fixed (the dashed lines). The black dash-dotted lines indicate the capture probabilities obtained by using the parabolic approximation for m₂ = 0.1m_T. Even at the high velocity v_∞ ∼ 1.6 km s⁻¹ or equivalently high eccentricity e − 1 = D(m/M)^2/3(av²_∞/Gm) ∼ 7 × 10⁻², it reasonably agrees with the restricted three-body results. In the bottom panel, the capture probability is shown as a function of a scaled velocity $(m/m_T)^{1/6}(m_{{\rm par}}/m_T)^{-1/2}v_\infty =(m_N/m_T)^{1/6} (2\bar{I} Gm_T/a)^{1/2}$. Although the difference among them is small, if the parabolic approximation is used for all the cases, then the solid and dashed lines should perfectly overlap with the black dash-dotted line.

Zoom In Zoom Out Reset image size

At the disruption, each member is captured with 50% probability if the initial energy $\bar{I}$ is close to zero. Using the relation $P_{{\rm cap}}(\bar{I})=P_{{\rm eje}}(-\bar{I})$, Figure 3 also shows that at high energies the capture is very rare. A transition occurs at an intermediate energy as shown in Figure 3. For D = 0.45 prograde orbits, the capture probability sharply drops from ∼50% around $\bar{I} \sim 1.4$. The critical value of v_∞ at which the capture probability of a binary member m₁ drops from ∼50% is given by

$$\begin{eqnarray} v_{\infty, {\rm crit}} &=&\sqrt{\frac{2\bar{I} Gm_2}{a}}\left(\frac{M}{m} \right)^{1/6} \sim 0.49 \left(\frac{\bar{I}}{1.4} \right)^{1/2} \left(\frac{m_2}{0.1m_T} \right)^{1/2}\nonumber\\ &&\times \left(\frac{a}{20R_T} \right)^{-1/2} \left(\frac{M/m}{4.37\times 10^3} \right)^{1/6} \,{\rm km\,s^{-1}}. \end{eqnarray} \tag{ 27 }$$

This estimate well explains the critical velocities in Figure 5. If Triton has a heavier companion and/or a member of a harder binary, then the critical velocity could be higher, provided the binary is disrupted by Neptune: D ≲ 1. The semimajor axis of the captured member m₁ is

$$\begin{eqnarray} r_a&=& \frac{a}{2|\bar{E}_1|}\left(\frac{M}{m_2} \right)\left(\frac{m}{M} \right)^{1/3} \sim 1.6\times 10^3R_N |\bar{E}_1|^{-1} \nonumber\\ &&\times \left(\frac{a}{20R_T} \right) \left(\frac{M/m_2}{4.80\times 10^4} \right) \left(\frac{M/m}{4.37\times 10^3} \right)^{-1/3},\quad \end{eqnarray} \tag{ 28 }$$

where $|\bar{E}_1|$ is about $|\Delta \bar{E}|\sim 1$. After capture, the orbit of Triton needs to shrink to the current observed r_a ∼ 14R_N and e ∼ 10⁻⁵, either through tides or other means (Correia 2009; Nogueira et al. 2011). The critical velocity and the semimajor axis for the secondary m₂ are obtained by exchanging the subscript 1 and 2 in Equations (27) and (28).

Although a detailed modeling of the solar system is beyond the scope of this paper, the Sun cannot be generally ignored in the capture process (e.g., Philpott et al. 2010; Gaspar et al. 2011). Since the tidal radius r_t ∼ 18R_N(m/m_T)⁻¹(a/20R_T) is much smaller than the radius of Neptune's Hill sphere r_H ∼ 4700R_N, Neptune completely dominates the attraction of a binary during the tidal encounter. Tidal effects from the Sun are negligible for the disruption process itself, which is the main focus of the paper. However, the permanent capture of Triton requires an additional condition that after the tidal breakup a capture orbit should be well inside the Hill sphere. Detailed numerical simulations show that retrograde satellites (those orbiting in the opposite sense as Neptune orbits the Sun⁹) of Neptune are more stable than prograde ones, and they are stable to distance of r_stable ∼ 0.4r_H (Nesvorný et al. 2003; Holman et al. 2004). Capture orbits are more eccentric than assumed in the numerical simulations, the stability region might be slightly smaller. Considering that the apocenter distance is about twice the semimajor axis (28) for highly eccentric orbits, the binary which had included Triton as a member should satisfy a relation

$$\begin{equation} \left(\frac{a}{20 r_T} \right) \left(\frac{m_2}{0.1 m_T} \right)^{-1} \left(\frac{m}{1.1 m_T} \right)^{1/3} \lesssim 0.4 \left(\frac{r_{{\rm stable}}}{0.3r_H} \right), \end{equation} \tag{ 29 }$$

where we have used $|\Delta \bar{E}| \sim 1$. A smaller value of the ratio a/m₂ gives a smaller apocenter distance and a higher critical velocity (27). As shown in the top panel of Figure 2 (the black solid and dashed lines), the phase-averaged value $\left\langle |\Delta \bar{E}| \right\rangle$ is slightly larger than unity for prograde binaries, and slightly smaller for retrograde binaries when the disruption probability is high. For a given encounter velocity, a retrograde orbit of the binary center of mass around Neptune with a prograde binary rotation is optimal to produce a stable capture orbit.

The ellipticity _N ∼ 1.7 × 10⁻² of Neptune induces a small deviation in its own gravitational potential from the point-mass estimate, especially in a non-planar configuration. However, the deviation in the attraction force ΔF/F and the tidal acceleration Δa_t/a_t is of the order of (R_N/r_t)²_N ∼ 5 × 10⁻⁵ at the tidal radius. The effect is clearly negligible.

We have discussed how the members of a binary are ejected or captured after a tidal encounter with a massive object. We have shown that the ejection and capture dynamics can be well understood in the framework of the restricted three-body approximation. When the three-body exchange reaction happens, the orbit of the center of mass of the binary should be almost parabolic |1 − e| ≲ D(m/M)^1/3(m_par/m) ≪ 1 or equivalently the semimajor axis is large |r_a/r_t| ≳ (M/m)^1/3(m/m_par). The essential quantity to characterize the disruption process is the energy change $\Delta \bar{E}$ at the encounter, which practically depends only on three parameters: the penetration factor D = r_p/r_t, the binary rotation direction, and the binary phase ϕ. Except for the phase shift of π, the energy change is exactly the same for the two members with arbitrary mass ratio.

In principle, for any positive (negative) orbital energy E of the center of mass of a binary, the heavier member has more chance to be ejected (captured), because it carries a larger fraction of the orbital energy. However, if the orbital energy is close to zero, the difference between their ejection (capture) probabilities becomes small, and there is practically no ejection and capture preference. For a parabolic orbit, each member is ejected in exactly 50% of the cases. The preference becomes significant when the absolute value of the energy |E| is comparable to the typical energy change (Gm₁m₂/a)(M/m)^1/3. On the other hand, if |E| is much larger than the typical energy change, they are both ejected for E > 0 or captured for E < 0, and there is no ejection or capture preference.

Corresponding to the typical energy change, we can define a critical encounter velocity $v_{\infty, {\rm crit}}=(2 \bar{I} Gm_{{\rm par}}/a)^{1/2}(M/m)^{1/6}$ for the capture process where the critical initial energy $\bar{I}\sim 0.5\hbox{--}1.5$ mainly depends on the penetration factor and the orientation of the binary. Since the distribution of the energy change is the same for the two members, the secondary star has a higher critical velocity. The tidal encounter might be responsible for the capture of not only Triton but also other irregular satellites in the solar system (Morbidelli 2006, but see also Vokrouhlický et al. 2008; Philpott et al. 2010. The stability of capture orbits for the solar perturbations also should be taken into account). However, the other irregular satellites are much lighter than Triton. The capture mechanism v_{∞, crit}∝m^1/2_par requires that these irregular satellites have been in a binary with a very massive partner, and that the less massive member has been predominantly captured.

If a binary has a large orbital energy E ≫ (M/m)^2/3(Gm²/a), then the disruption would happen slightly inside the usual tidal radius at which the BH tidal forces become comparable to the binary mutual gravity forces. The abrupt disruption approximation provides a good estimate of the energy change, and the typical energy change in the high-energy regime is much larger than in the low-energy regime. However, the change is not significant compared to the initial energy E, and the tidal encounter is not an efficient acceleration process anymore. In the high-energy regime, both members are ejected after the BH encounter.

In recent years, observations have identified a remarkable number of hypervelocity stars (Brown et al. 2009; Tillich et al. 2011). The distribution of line-of-sight velocities of hypervelocity star candidates shows a long tail in the high-velocity region (v > 275 km s⁻¹) which includes comparable numbers of unbound and bound stars assuming that the escape velocity of the Galaxy at 50 kpc is ∼360 km s⁻¹ (Brown et al. 2009; Kenyon et al. 2008). Since the initial orbital energy is negligible when the three-body exchange reaction happens, the velocity distribution should reflect the distribution of the energy change ΔE at the tidal encounter (and the Galactic potential) with its dependence on the periapsis distance, the binary orientation, and phase.

Binaries are supplied to the BH at the Galactic center predominantly from its radius of influence r_h ∼ a few pc or even beyond it (e.g., Perets et al. 2007; another interesting possibility is that they might arise from the nuclear stellar disk, e.g., Madigan et al. 2011). Those that come from about the radius of influence are on elliptical orbits. Since the radius of influence is much larger than the tidal radius, the preferential ejection for high- and low-mass members (m₁ > m₂) is therefore irrelevant for stellar binaries (m₁/m₂ ≲ 10) if the binary semimajor axis is smaller than ∼1 AU. However, stars with planets (m₁/m₂ ≳ 10³) of a semimajor axis of 1 AU will predominantly eject the planets.

We thank the referee, Hagai Perets, for his constructive comments, and Mostafa Ahmadi, Phil Armitage, Witold Maciejewski, and Ruth Murray Clay for useful discussion. This research was supported by STFC, ERC, and IRG grants, and Packard, Guggenheim, and Radcliffe fellowships.