An analytic method for identifying dynamically-formed runaway stars

TL;DR

This paper develops an analytic method to identify dynamically-formed runaway stars by relating their velocity vectors to binary properties, validated through simulations, with applications to GAIA data and specific runaway stars.

Contribution

It introduces a simple analytic formula linking ejected stars' velocities to binary separation, aiding the identification of dynamical runaways in observational data.

Findings

Analytic formula matches simulation results with high accuracy.

Method can be applied to GAIA data to identify runaway stars.

Application to specific stars constrains background potential history.

Abstract

In this paper, we study the 3-body products (two single stars and a binary) of binary-binary (2+2) scattering interactions. This is done using a combination of analytic methods and numerical simulations of 2+2 scattering interactions, both in isolation and in a homogeneous background potential. We derive analytically a simple formula relating the angle between the velocity vectors of the two ejected single stars and the orbital separation of the remaining binary. We compare our analytic formulation to numerical scattering simulations, and illustrate that the agreement is excellent, both in isolation and in a homogeneous background potential. Our results are ideally suited for application to the GAIA database, which is expected to identify many hundred runaway stars. The analytic relation presented here has the potential to identify runaway stars formed dynamically with high confidence. Finally, by applying our method to the runaways AE Aur and $\mu$ Col, we illustrate that it can be used to constrain the history of the background potential, which was denser than the presently observed density in the case of the Trapezium cluster.

Figures

Tables

Table 1. Table 1: p s , 1 subscript 𝑝 s 1 p_{\rm s,1} and p s , 2 subscript 𝑝 s 2 p_{\rm s,2} are the momenta of the first and second ejected star, respectively, i.e., p s , i = m i ​ v s , i subscript 𝑝 s 𝑖 subscript 𝑚 i subscript 𝑣 s 𝑖 p_{{\rm s},i}=m_{\rm i}v_{{\rm s},i} ( i = 1 , 2 𝑖 1 2 i=1,~{}2 ). B = ( m t ​ v s , 2 ) / ( m s , 1 ​ v s , 1 ) 𝐵 subscript 𝑚 t subscript 𝑣 s 2 subscript 𝑚 s 1 subscript 𝑣 s 1 {B}=(m_{\rm t}v_{\rm s,2})/(m_{\rm s,1}v_{\rm s,1}) where m t subscript 𝑚 t m_{\rm t} is the mass of the substellar system which the first ejected star left behind, i.e., m t = m s , 2 + m b subscript 𝑚 t subscript 𝑚 s 2 subscript 𝑚 b m_{\rm t}=m_{\rm s,2}+m_{\rm b} .

case	$v_b^{\prime }$	$\xi$
case 1	$\frac{1}{m_b}\sqrt{p_{s,1}^2+p_{s,2}^2+2p_{s,1}{p}_{s,2}\cos \xi }$	${\cos }^{-1}\left(\frac{-1}{B}\right)\le \xi \le \text{uppi}$ for $B>1$.
case 2 (++)	No solutions	-
case 2 (+-)	$\frac{1}{m_b}\sqrt{p_{s,1}^2+p_{s,2}^2+2p_{s,1}{p}_{s,2}\cos \xi }$	for .
case 2 (-+)	$\frac{1}{m_b}\sqrt{p_{s,1}^2+p_{s,2}^2+2p_{s,1}{p}_{s,2}\cos \xi }$	for $B\ge 1$ and for .
case 2 (–)	No solutions	-
case 3	No solutions	-
case 4 (++)	No solutions	-
case 4 (+-)	No solutions	-
case 4 (-+)	$\frac{1}{m_b}\sqrt{p_{s,1}^2+p_{s,2}^2+2p_{s,1}{p}_{s,2}\cos \xi }$	$0\le \xi \le \frac{\text{uppi}}{2}$ for any $B$
case 4 (–)	No solutions	-

Equations

$${\text{p}}_{\mathrm{s},1}+{\text{p}}_{\mathrm{s},2}+{\text{p}}_{\mathrm{b}}=\text{0}.$$

$$v_{\mathrm{b}}$$

$$\cos \xi$$

$$p_{\mathrm{s},2}=\sqrt{p_{\mathrm{s},1}^2+p_{\mathrm{b}}^2+2p_{\mathrm{s},1}{p}_{\mathrm{b}}\cos \Psi }.$$

$$\cos \Psi =-\frac{p_{\mathrm{s},1}+p_{\mathrm{s},2}\cos \xi }{\sqrt{p_{\mathrm{s},1}^2+p_{\mathrm{s},2}^2+2p_{\mathrm{s},1}{p}_{\mathrm{s},2}\cos \xi }}.$$

$$v$$

$$M_{\mathrm{en},\mathrm{bg}}(r)$$

$${\text{f}}_{\mathrm{bg}}(r)$$

$$\displaystyle=\begin{cases}-\frac{4}{3}\uppi Gm\rho\;\textbf{{r}}\hfill\hskip 21.68121ptr\leq r_{\rm bg}\,;\\ -\frac{4}{3}\uppi Gm\rho\Big{(}\frac{r_{\rm bg}}{r}\Big{)}^{3}\textbf{{r}}\hfill\hskip 21.68121ptr>r_{\rm bg}\,,\\ \end{cases}$$

$$V_{\mathrm{bg}}(r)$$

$$v_{\mathrm{esc}}=\sqrt{\frac{2G{M}_{\mathrm{bg}}}{r_{\mathrm{bg}}}}=\sqrt{2}{r}_{\mathrm{bg}}w.$$

$$r(t)$$

$$v(t)$$

$$[(r_{\mathrm{bg}}-r_{\mathrm{ej}})w{]}^2+v_{\mathrm{bg}}^2=v_{\mathrm{ej}}^2,$$

$$\displaystyle v_{\rm ej}=\sqrt{[(r_{\rm bg}-r_{\rm ej})w]^{2}+v_{\rm bg}^{2}}=\sqrt{\frac{1}{2}v_{\rm esc}^{2}\Big{[}1-\frac{\sqrt{2}r_{\rm ej}w}{v_{\rm esc}}\Big{]}^{2}+v_{\rm bg}^{2}}.$$

$$\frac{\Delta v}{v_{\mathrm{ej}}}$$

$$\displaystyle\frac{\Delta v}{v_{\rm ej}}\sim\Big{|}1-\sqrt{1-\frac{1}{2}\Big{(}\frac{v_{\rm esc}}{v_{\rm ej}}\Big{)}^{2}}\Big{|}\sim\Big{(}\frac{v_{\rm esc}}{v_{\rm ej}}\Big{)}^{2}+O(\frac{v_{\rm esc}^{4}}{v_{\rm ej}^{4}}).$$

$$e_{\mathrm{bg}}=\frac{1}{2}{v}_{\mathrm{bg}}^2-\frac{G{M}_{\mathrm{bg}}}{r_{\mathrm{bg}}},$$

$$e(t)=\frac{1}{2}v(t{)}^2-\frac{G{M}_{\mathrm{bg}}}{r(t)}.$$

$$v(t)=\frac{dr(t)}{dt}=\sqrt{2e_{\mathrm{bg}}}\sqrt{1+\frac{G{M}_{\mathrm{bg}}}{e_{\mathrm{bg}}}\frac{1}{r}}.$$

$${\int }_{t_{\mathrm{bg}}}^{t}dt\sqrt{2e_{\mathrm{bg}}}$$

$$\sqrt{2e_{\mathrm{bg}}}(t-t_{\mathrm{bg}})$$

$$\frac{r_{\mathrm{bg}}}{R}$$

$$\displaystyle=\Big{(}\frac{v_{\rm bg}}{v_{\rm esc}}\Big{)}^{2}-1\simeq\Big{(}\frac{v_{\rm bg}}{v_{\rm esc}}\Big{)}^{2}\simeq 10^{2}\Big{(}\frac{v_{\rm bg}}{70\;\mathrm{km}/\;\mathrm{s}}\Big{)}^{2}\Big{(}\frac{v_{\rm esc}}{10\;\mathrm{km}/\;\mathrm{s}}\Big{)}^{-2}\,,$$

$$\displaystyle r(t)=\frac{1}{2}{R}\Bigg{[}-1\pm\sqrt{1+4\Big{[}\sqrt{\frac{r_{\rm bg}}{{R}}+\Big{(}\frac{r_{\rm bg}}{{R}}\Big{)}^{2}}+\sqrt{\frac{2GM_{\rm bg}}{{R}^{3}}}(t-t_{\rm bg})\Big{]}^{2}}\Bigg{]}$$

$$\displaystyle v(t)=\frac{dr(t)}{dt}=\frac{2\sqrt{\frac{2GM_{\rm bg}}{{R}}}\Big{[}\sqrt{\frac{r_{\rm bg}}{{R}}+\Big{(}\frac{r_{\rm bg}}{{R}}\Big{)}^{2}}+\sqrt{\frac{2GM_{\rm bg}}{{R}^{3}}}(t-t_{\rm bg})\Big{]}}{\sqrt{1+4\Big{[}\sqrt{\frac{r_{\rm bg}}{{R}}+\Big{(}\frac{r_{\rm bg}}{{R}}\Big{)}^{2}}+\sqrt{\frac{2GM_{\rm bg}}{{R}^{3}}}(t-t_{\rm bg})\Big{]}^{2}}}$$

$$v(t=t_{\mathrm{bg}})$$

$$\displaystyle=2v_{\rm bg}\frac{\Big{(}\frac{v_{\rm bg}}{v_{\rm esc}}\Big{)}^{2}-1}{2\Big{(}\frac{v_{\rm bg}}{v_{\rm esc}}\Big{)}^{2}-1}\simeq v_{\rm bg}$$

$$v(t\to \infty )$$

$$v(t)$$

$$r(t)$$

$$v(t)=$$

Full text

††pubyear: 2016††pagerange: An analytic method for identifying dynamically-formed runaway stars–9

1 introduction

In the general gravitational $N$-body problem, all masses are non-zero and their initial positions and velocities are not arranged in any specific way. The simplest subset of this more general problem, called the three-body problem, involves only three particles. And yet, the problem is analytically intractable without simplifying assumptions. This is in contrast to the two-body problem, in which the positions and velocities of both particles are known exactly at any time given any initial configuration. This complication stems from the fact that, for $N>$ 2, there are no coordinate transformations that simplify the problem sufficiently (e.g., the centre of mass coordinate system). Consequently, it is not possible to express the positions and velocities of all particles at any future time given a set of initial conditions using analytic theory alone. What’s more, the long-term evolution of these systems are very sensitive to the initial conditions. If one waits long enough, however, a final stable state will always emerge. This usually occurs after the system has disrupted and all particles can be regarded as having escaped to spatial infinity. But the final state cannot be calculated directly from the initial conditions without the need for computer simulations. The problem is inherently chaotic.

Few-body systems have been extensively investigated over the last few decades (e.g Heggie, 1975; Mikkola, 1983, 1984). These studies typically adopted numerical approaches. However, analytic theory remains a powerful tool in the study of gravitational dynamics. For example, if the outcome of the interaction is known a priori, aspects of the problem can become deterministic. That is, energy and momentum conservation can be invoked to directly relate the initial conditions to the final outcome properties.

Here, we apply this logic to the four-body problem and binary-binary scattering. In particular, for (nearly) identical initial binaries and large viral ratios, most four-body interactions decay to produce two single stars and a binary (Leigh et al., 2016). As we will show, knowledge of the initial interaction energy and the momentum are sufficient to directly relate the properties of the ejected single stars to the remaining binary. This offers a clear prediction for observations of dynamically-formed runaway stars; young O/B stars that have been ejected from their host star cluster at high velocity ($\gtrsim$ 30 km s*-1*) (e.g Fujii & Portegies Zwart, 2011; Gualandris et al., 2004). Consequently, the method presented here can be used to unambiguously distinguish runaways formed during binary-binary interactions from those formed from supernova explosions in a binary system (e.g Stone, 1979, 1982), by making specific predictions for the observable properties of the left-over binary. Our results are particularly useful in light of the expected hundreds of detections of runaways stars with the GAIA satellite (i.e. Kenyon et al. 2014).

Our paper is organized as follows: Section 2 presents the derivation of the analytical formulation of the 2+1+1 problem, first in the purely stellar dynamical case (Sec.2.1), and then in the presence of a homogeneous background potential (Sec.2.2). A comparison of the predictions of the analytical formulae and the results of numerical scattering experiments is presented in Section 3. Astrophysical implications of our findings are discussed in Section 4, with a particular emphasis on the application to runaway stars. We summarize in Section 5.

2 Analytic formulation of 2+1+1 case

2.1 Analytic formula : Purely stellar dynamics

When two binaries collide, a chaotic interaction ensues until two single stars are ejected sequentially, leaving behind a binary pair. After the first ejection event, a single star of mass $m_{\rm s,1}$ and the remaining system of mass $m_{\rm t}$ recede in opposite directions (in the CM frame of all four stars). We assume a velocity of $\textbf{{v}}_{\rm s,1}$ (so the momentum $\textbf{{p}}_{\rm s,1}=m_{\rm s,1}\textbf{{v}}_{\rm s,1}$) for the first ejected single star. The remaining system (with three stars) recedes with velocity $\textbf{{v}}_{\rm t}=(m_{\rm s,1}/m_{\rm t})\textbf{{v}}_{\rm s,1}$. At the time of the second ejection event, another single star is ejected with velocity $\textbf{{v}}_{\rm s,2}$ (the momentum $\textbf{{p}}_{\rm s,2}$) and leaves behind the binary with a recoil velocity of $\textbf{{v}}_{\rm b}$ (the momentum $\textbf{{p}}_{\rm b}$). We present a schematic diagram in Figure 1 showing the two ejection events.

Now, in the CM of the four stars, the three momentum vectors satisfy the following relation,

[TABLE]

Using Equation 1, we find an expression for $v_{\rm b}$ as follows,

[TABLE]

where $p_{\rm s,1}$, $p_{\rm s,2}$ and $p_{\rm b}$ are the momenta of the first, second ejected stars and the binary, respectively, i.e., $p_{i}=m_{\rm i}v_{i}$. And $\xi$ is the relative angle between the two ejected single stars, defined as,

[TABLE]

We also provide the derivations for Equation 2 using their vector components in Appendix A. In this derivation, we consider the two ejection events in their respective center of mass frames separately. This allows us to find explicit relations for the relative angles between the two ejected stars in the different reference frames. This approach could be useful for identifying the first ejected single star, when applying our method to observations of runaways where only one single star and an associated binary have been identified.

Using Equation 2, we can find an expression for the relative angle between the first ejected star and the binary as a function of $p_{\rm s,1}$, $p_{\rm s,2}$ and $\xi$. We denote this angle as $\Psi$ (see Figure 1). Then the angle between the second ejected single star and the binary (denoted by $\Phi$) is automatically determined since $\xi+\Psi+\Phi=2\pi$. Equation 1 also gives the following equation,

[TABLE]

Plugging Equation 2 into Equation 4 yields

[TABLE]

Therefore, given the velocities and the masses of the two ejected stars, with Equations 2 and 5 we can determine the momenta of the binary (and hence the binary speed, given the binary mass) and the relative angles between the binary and each of the two ejected stars ($\Psi$ and $\Phi$).

There are several points worth emphasizing here.

We assume only conservation of momentum to derive Equations 2, 3 and 5 without additional conditions. 2. 2.

These formulas can be applied independently of the details of the initial binary-binary encounter (i.e., head-on, or with non-zero impact parameter) since we consider only the velocities at the time of the last ejection event. The only condition is that two binaries interact to produce a 2+1+1 outcome. As an example, in Section 3.2, we apply our analytic formula to the 2+1+1 outcome from scatterings of two binaries with equal/unequal energies. 3. 3.

In Equation 2, the two ejected stars are distinguished in terms of which star is ejected first. However, Equation 2 is symmetric upon exchanging the two single star momenta. Clearly, $v_{\rm b}$ at $\xi=0$ depends only on the total sum of the momenta, whereas $v_{\rm b}$ at $\xi=\pi$ depends on the difference between the two momenta. Hence the order of ejection events is irrelevant in terms of estimating the final velocity of the binary (see the green line in the left panel of Figure 3). In other words, the subscripts 1 and 2 do not necessarily correspond to the order of ejection events. This is of critical value for applying our method to observational data, since in practice it is very hard to distinguish the order of ejection events. 4. 4.

We can generalize Equation 2 to a stellar system going through any number of ejection events. Consider a stellar system consisting of $N_{\rm tot}$ stars. In the CM of the whole system, the speed of a substellar group left behind after the $N$th ejection event can be expressed as

[TABLE]

where $p_{i}$ represents the momentum of the $i_{\rm th}$ ejected object. 5. 5.

If scatterings between two binaries occur in a background potential (eg., a star cluster), an escape condition can be derived by comparing the escape velocity of the potential with $v_{\rm b}$, as estimated from our formula. We will return to this in the next section.

2.2 Analytic formula: Homogeneous background potential

In this section we generalize our analytic formulation to include a homogeneous background potential.

When we observe at time $t=t_{\rm obs}$ the velocities of two field stars (outside the potential) suspected of being a causally-related pair of runaway stars, we can trace their velocities backward in time. In other words, we can estimate the velocities at any time $t~{}(<t_{\rm obs})$ post-interaction (i.e., after the second ejection event). Starting with the velocities of the ejected single stars, $\textbf{{v}}_{\rm s,1}(t=t_{\rm obs})$ and $\textbf{{v}}_{\rm s,2}(t=t_{\rm obs})$ as observed outside of the potential ($r>r_{\rm bg}$) at $t=t_{\rm obs}$, we integrate the equations of motion backward in time in order to estimate the stars’ initial ejection velocities inside the potential. Then, we can use equation 2 to calculate the velocity of the final binary at the time of the last ejection event. Integrating the equation of motion for the binary forward in time gives us the speed of the binary $v_{\rm b}(t=t_{\rm obs})$ at $t=t_{\rm obs}$.

We consider a spherically symmetric potential with constant density $\rho$ and outer boundary $r_{\rm bg}$. For a given total background mass $M_{\rm bg}$, the outer boundary is automatically set. The mass enclosed within a spherical volume of radius $r$ can be written,

[TABLE]

The gravitational force imparted by the background mass on a given star particle at $r$ is given by:

[TABLE]

where $m$ is the mass of the star and r is the vector pointing from the system CM to the star. The term $\sqrt{(4\uppi G\rho)/3}$ corresponds to the frequency of the resulting harmonic motion, denoted by $w$. Accordingly, the background potential has the following form:

[TABLE]

and the escape velocity of the potential $v_{\rm esc}$ at $r=r_{\rm bg}$ is written as,

[TABLE]

This is the same background potential as the one adopted by Ryu et al. (2017b). See Ryu et al. (2017b) for more details.

Given the adopted background potential, we solve for the equations of motion of the stars inside and outside the potential individually. To do this, we make the following assumptions: At each ejection event, the ejected star reaches a sufficiently high velocity to escape the underlying potential well (without becoming trapped and undergoing subsequent oscillations within the background potential). Therefore the angles between stars do not significantly change while they escape, which we will show in Figure 6.

2.2.1 Inside the potential, $r\leq r_{\rm bg}$

Given the radial components of the velocities of stars, we can then solve the equations of motion in 1-D. For $r\leq r_{\rm bg}$, the solution for Equation 2.2 is,

[TABLE]

The speed is

[TABLE]

where $v_{\rm ej}$ and $r_{\rm ej}$ denote the speed and radial distance from the system CM (or the center of the background potential) at the time of an ejection event $t=t_{\rm ej}$.

When a star crosses the outer boundary $r_{\rm bg}$ with speed $v_{\rm bg}$ at time $t=t_{\rm bg}$, Equations 11 and 12 yield the following relations,

[TABLE]

or,

[TABLE]

It is useful to consider how much the speeds of the stars change inside the potential from the moment of ejection until the time of escape. For this, we introduce a fractional difference between $v_{\rm ej}$ and $v_{\rm bg}$,

[TABLE]

Especially for stars ejected near centre ($r_{\rm ej}w\ll v_{\rm esc}$) at sufficiently high velocities to escape,

[TABLE]

This implies that: 1) for a given $v_{\rm ej}$ and $r_{\rm ej}$, $v_{\rm esc}$ is a good diagnostic (rather than either $M_{\rm bg}$ or $\rho$ individually) to gauge the overall decrease in the speeds of the escaping stars; 2) the fractional difference ($\Delta v/v_{\rm ej}$) or absolute difference ($\Delta v$) in $v_{\rm ej}$ depends on $v_{\rm ej}$ itself.

2.2.2 Outside the potential, $r>r_{\rm bg}$

For $r>r_{\rm bg}$, we apply conservation of energy instead of directly solving the second-order differential equation in Equation 2.2. At $r=r_{\rm bg}$, the total specific energy $e_{\rm bg}$ is defined as,

[TABLE]

Similarly, at $r>r_{\rm bg}$ and $t>t_{\rm bg}$,

[TABLE]

Equating Equations 17 and 18 and rearranging terms,

[TABLE]

Defining $GM_{\rm bg}/e_{\rm bg}\equiv{R}$ and integrating both sides of the above gives,

[TABLE]

We are unable to solve further for a simple expression for $r(t)$. However, we note that the second term (with $\log$) on the right hand side can be neglected when $r/{R}$ is sufficiently large.

It useful to further consider the range in ${R}/r_{\rm bg}$ and its physical meaning. $r_{\rm bg}/{R}$ measures the relative difference between the absolute magnitude of the kinetic energy and that of the potential energy at the outer boundary of the potential. Assuming $v_{\rm bg}>v_{\rm esc}$, we can express ${R}$ in terms of $r_{\rm bg}$ and $v_{\rm esc}$,

[TABLE]

where $v_{\rm bg}\sim 60-70\;\mathrm{km}\;\mathrm{s}^{-1}$ are typical speeds of runaway stars formed during binary-binary numerical scattering experiments (Ryu et al. 2017a, “simulation paper” in the following). Typically, $r_{\rm bg},~{}r(t)>{R}$. Such large values for $r_{\rm bg}/{R}$ imply that the kinetic energy is dominant over the background potential energy when stars are escaping from the potential. Moreover, this may mean that the decrease in speed outside of the potential is insignificant as long as $v_{\rm esc}$ is not too high (see also the $v_{\rm ej}$ dependence in Equation 15).

Given a high value for $r/{R}$ (i.e., for rapidly moving stars, hence a suitable value for runaway stars), we can ignore the log-term in Equation (2.2.2). Then we can find a relatively simple expression for $r(t)$ at $t\geq t_{\rm bg}$, which is,

[TABLE]

where we must choose the plus sign since $r(t=t_{\rm bg})=r_{\rm bg}$. Differentiating with respect to $t$ gives,

[TABLE]

Since we made an approximation in Equation 2.2.2 in order to derive the above, we must also ensure that $v(t)$ is continuous at $t=t_{\rm bg}$ and well-defined at $t\rightarrow\infty$.

$v(t)$ at $t=t_{\rm bg}$.

[TABLE]

for $v_{\rm bg}/v_{\rm esc}\gg 1$, which is the same assumption that we made to solve Equation 2.2.2. 2. 2.

$v(t)$ at $t\rightarrow\infty$.

[TABLE]

which is the velocity at $t=\infty$ (see equation 18).

We can simplify the solution further for stars moving at high speeds (i.e., large $v_{\rm bg}$). We can approximate their speeds as being constant since $v(t)$ is a monotonically decreasing function with the same two boundary values ($v_{\rm bg}$ at $t=t_{\rm bg}$ and $t\rightarrow\infty$). Note that for such high speeds, $\sqrt{2e_{\rm bg}}\simeq v_{\rm bg}$. This is confirmed in the simulations of our simulation paper. Figure 2 shows the time evolution of the speeds of the binary and the two single stars formed during binary-binary encounters (Model 1 and Model 2 of the simulation paper) at $t=$1, 2, 3 and 4$\;\mathrm{Myr}$. We differentiate the two single stars by their speeds: $S_{\rm fast}$ ($S_{\rm slow}$) refers to the single star with the higher (slower) speed, for a given run. At each $t$, the three vertically distributed dots correspond to (from left to right) the binary (black dots), $S_{\rm fast}$ (red dots) and $S_{\rm slow}$ (blue dots), respectively. The horizontal dotted line indicates the escape velocity. To avoid overcrowding, we offset the distributions for $S_{\rm fast}$ and $S_{\rm slow}$ by $+0.1\;\mathrm{Myr}$ and $+0.2\;\mathrm{Myr}$ for a given $t$ (i.e., shift to the right). In the simulations, more than 95% of all outcomes escape from the background potential before $t=1\;\mathrm{Myr}$. Figure 2 shows that the speeds of rapidly-moving stars ($v/v_{\rm esc}\geq 5-7$ for both single stars) barely change once they have escaped from the potential. We also refer to the cumulative velocity distributions for the two single stars and the binary shown in Appendix A of the simulation paper.

However, this approximation may lose its validity for stars whose speeds at $r=r_{\rm bg}$ are not sufficiently large compared to $v_{\rm esc}$ (i.e., typically more massive escapers). In our numerical experiments (simulation paper), we find that for most final binaries, $v_{\rm bg}/v_{\rm esc}\simeq 1-3$. Therefore, in the following (Section 3), we will use Equation 23 to estimate the final speeds of the binaries along with a constant $v(t)$ for ejected single stars moving at high speeds. Using the following equations of motions for stars at $t$: for $r(t)<r_{\rm bg}$,

[TABLE]

and for $r(t)\geq r_{\rm bg}$ (Equation 22)

[TABLE]

where $v_{\rm ej}\cos w~{}t_{\rm bg}=v_{\rm bg}$, we can calculate $v_{\rm ej}$ of the two ejected stars given their observed values. After getting $v_{\rm b}(t=t_{\rm ej})$ from Equation 2, we can calculate what we would observe for the binary at time $t=t_{\rm obs}$, i.e., $v_{\rm b}(t=t_{\rm obs})=v_{\rm b,obs}$. In order to do this, we first use Equation 14 to find $v_{\rm bg}$, then Equation 23 to find $v_{\rm b,obs}$.

Finally, we summarize the above procedure as follows. For a background potential with an escape velocity not too high compared to the typical ejection velocities of the stars (a valid assumption for runaway stars), when we observe two runaway stars with speeds $v_{\rm s,1,obs}$ and $v_{\rm s,2,obs}$ with a relative angle $\xi$ at $t=t_{\rm obs}$ 111If their radial distances $r_{\rm s,1,obs}$ and $r_{\rm s,2,obs}$ (in the CM of four stars) are known, $t_{\rm obs}(\gg t_{\rm bg})$ can be approximated for the two fast runaway stars as,

$\displaystyle t_{\rm obs}\simeq r_{\rm s,1,obs}/v_{\rm s,1,obs}\simeq r_{\rm s,2,obs}/v_{\rm s,2,obs}.$

(26)

, the velocity of the binary at $r=r_{\rm bg}$ and $t=t_{\rm bg}~{}(<t_{\rm obs})$ left behind by the two runaway stars can be estimated as,

[TABLE]

where

[TABLE]

where $r_{\rm ej,1}$ ($r_{\rm ej,2}$) represents the radial distance between the CM of the background potential and the CM of the system made of the ejected star and the remaining $N=3$ ($N=2$) stars at the first (second) ejection event.

Finally, the speed and radial distance of the final binary outside the potential ($r\geq r_{\rm bg}$) at any $t~{}(\geq t_{\rm bg})$ are described as,

[TABLE]

[TABLE]

where $\mathfrak{R}=r_{\rm bg}/[(v_{\rm b,bg}/v_{\rm esc})^{2}-1]$ 222As a sanity check, we explore two extreme limits, 1 $t\rightarrow\infty$ and 2 $\rho,~{}M_{\rm bg}\rightarrow 0$ .

$t\rightarrow\infty$

$\displaystyle v_{\rm b}(t>t_{\rm bg})\sim\sqrt{\frac{2GM_{\rm bg}}{\mathfrak{R}}}\sim\sqrt{v_{\rm b,bg}^{2}-v_{\rm esc}^{2}}\sim\begin{cases}v_{\rm bg}\hskip 7.22743pt\text{for}\,\,\,\ v_{\rm b,bg}\gg v_{\rm esc};\\ 0\hskip 14.45377pt\text{for}\,\,\,\ v_{\rm b,bg}\simeq~{}v_{\rm esc}.\end{cases}$

(32)

Both of the two limiting values are consistent with the speeds at spatial infinity for an $\sim(1/r)$ gravitational potential.

$\rho,~{}M_{\rm bg}\rightarrow 0$

In this limit,

$\displaystyle\mathfrak{R}$ $\displaystyle=\frac{r_{\rm bg}v_{\rm esc}}{v_{\rm b,bg}^{2}-v_{\rm esc}^{2}}\sim 0,$

(33)

$\displaystyle\sqrt{\frac{2GM_{\rm bg}}{\mathfrak{R}}}$ \displaystyle=\sqrt{\frac{2GM_{\rm bg}}{r_{\rm bg}}\Big{[}\frac{v_{\rm b,bg}}{v_{\rm esc}}-1\Big{]}^{2}}=\sqrt{v_{\rm b}^{2}-v_{\rm esc}^{2}}\sim v_{\rm b}.

(34)

Now, as the terms with $\mathfrak{R}^{k}$ ($k>0$) converge toward zero,

$\displaystyle r_{\rm b}(t>t_{\rm ej})\sim\sqrt{\frac{2GM_{\rm bg}}{\mathfrak{R}}}t\sim v_{\rm b}t$

(35)

which is expected for ejected stars after being completely isolated from their parent stellar group.

.

Given the equations above, we find that, when the following relation is satisfied,

[TABLE]

it is more likely that the final binary will have remained bound to the potential. Using Equation 36, for two runaway stars and a binary which are believed to form from a 2+2 encounter, we can estimate a rough lower limit for $v_{\rm esc}$ of the parent cluster at the moment of the last ejection event. We will apply the equation above for the Trapezium cluster in Section 4.

3 Results

In this section, we apply our analytic formulae to simulated 2+1+1 outcomes formed during binary-binary encounters in the presence/absence of a homogeneous background potential. The numerical data are taken from the simulation paper. There, we performed suites of numerical scattering experiments between two binaries (which serve as proxies for O-type stars), both in isolation and in a homogeneous background potential. However, since we considered binary-binary encounters occurring near the system CM and $r_{\rm ej}\leq r_{\rm bg}$, we assume that $r_{\rm ej}$ is negligible in the following analysis, i.e., $r_{\rm ej}\simeq 0$. We compare our analytic formulae to the results of these simulations, and confirm their validity. We further study here the effects of the background potential on the formation of runaway stars, in terms of affecting the properties of the runaway stars and final binary. For more details regarding the setup of the simulations, see the simulation paper.

3.1 Results for scatterings in isolation

In Figure 3, we show the speeds of the final binaries $v_{\rm b}$ (upper panel), the corresponding semimajor axes $a$ (middle panel) and the relative angle $\Psi$ between the binary and the faster single star for each simulation (bottom panel), given by our analytic formulae (Equations 2 and 5, solid/dotted lines) for different combinations of the ejected single star speeds ($v_{\rm s,1}$, $v_{\rm s,2}$) in unit of$\;\mathrm{km}\;\mathrm{s}^{-1}$. We plot equal values for both $v_{\rm s}$ parameters in the left column and unequal values in the right column. To compare, we mark via the data points the results of the numerical experiments of binary-binary scatterings (2+1+1 outcome) performed in the simulation paper. For the simulated data points, the final speeds of the ejected stars are constrained to be within 5% of the speeds shown in the plots. The lines for $v_{\rm b}$, $a$ and $\Psi$ (also in all plots below) are estimated using the median values of the speeds of the ejected stars. Hence the actual combinations of the speeds for the estimates are slightly different from those indicated in the plots. For the analytic calculations, we assume the same total initial energy and masses for the stars as adopted in the simulated data. For each combination of the single star ejection speeds, the lines and dots share the same color. For comparison, we over-draw the lines and dots for ($v_{\rm s,1}$, $v_{\rm s,2}$)=(100, 100) in the right panel in faint blue.

The results of our numerical scattering experiments are in excellent agreement with our analytic formulae. As shown in the upper panel of Figure 3, the speed of the final binary is the smallest at $\xi=180^{\circ}$, and rises steadily towards $\xi=0^{\circ}$. This is easily understood via conservation of momentum. Finally, in order to confirm the symmetric nature of Equation 2, we mark the data points with $(v_{\rm s,1},v_{\rm s,2})=(65,20)$ (right column) by green hollow circles, and those with $(20,65)$ by green solid circles. As expected, both sets of dots follow the theoretical line, which shows that the order of ejection events is irrelevant in terms of estimating the final binary speed using Equation 2.

In the middle panel of Figure 3, we estimate the semimajor axes $a$ of the final binaries assuming conservation of energy. The simulated data from are recorded when all stars (single stars and binaries) are sufficiently far from the system CM that no further (significant) gravitational interaction occurs between them. At this point, the mutual gravitational potentials between the stars are negligible. Hence, the total energy $E_{\rm tot}$ is the sum of the total kinetic energy for all objects ($KE$) and the binding energy of the binary ($E_{\rm bind}$), namely,

[TABLE]

where $m_{\rm b,1}$ and $m_{\rm b,2}$ are the masses of the component stars of the final binary (i.e., $m_{\rm b}=m_{\rm b,1}+m_{\rm b,2}$). Combining the results for $v_{\rm b}$ and $a$ in both panels, for a given combination of speeds for the two ejected stars (i.e., the same $KE_{\rm s,1}+KE_{\rm s,2}$ or same $KE_{\rm b}+E_{\rm binding}$), when two single stars are ejected with a smaller angle ($\xi$) and at high velocities, conservation of momentum requires that the final binary moves at a higher velocity (i.e., larger $KE_{\rm b}$ for a given $KE_{\rm b}+E_{\rm binding}$). Consequently, a more compact binary is formed.

In the bottom panel, the simulated data and the lines from our analytic formula (Equation 5) show good agreement. However, we see a discrepancy between the results for equal and unequal values of the two single star velocities (left and right panels, respectively). For equal speeds, independently of the values for the speeds, the angles $\Psi$ and $\xi$ have an unique relation. This is because, as the speeds of the two single stars of same mass become equal, $\Psi$ loses the dependence on the speeds themselves. For encounters with the exactly same values of the single star masses $m_{\;\mathrm{s},1}$ and $m_{\rm s,2}$, Equation 5 reduces to the following equation,

[TABLE]

We depict this case in the bottom left panel of Figure 3 using the cyan-dotted line, running diagonally from $\Psi=180^{\circ}$ (left upper corner) to $\xi=180^{\circ}$ (right bottom corner). However, the solid lines in the figure for equal single star velocities (left panel) increase again up to $\Psi=180^{\circ}$ as $\xi$ increases to $\xi=180^{\circ}$, similar to the lines for unequal velocities (right panel). This is because the single star velocities actually used for these estimates are not exactly identical (we noted above that the lines are drawn using the median values of the ejected stars’ speeds). We further note that the simulated data for both $(v_{\rm s,1},v_{\rm s,2})=(65,~{}20)$ (green hollow dots) and (20, 60) (green solid dots) nicely satisfy Equation 5. This is because we define $\Psi$ as the angle between the velocity vector of the binary and that of the more rapidly-moving single star.

3.2 Results for scatterings in a background potential

In Figure 4 and 5, we show $v_{\rm b}$ (left panels) and $a$ (right panels) for the final binaries which have escaped from the background potential using Equation 30 for equal and unequal single star velocities, respectively. We calculate the semimajor axes as before, assuming conservation of energy (Equation 37)333Note that for such high speeds of the single stars and large distances from the potential, $V_{\rm bg}$ has insignificant contribution to the total energy.. However, we have used $\mathcal{P}_{\rm s,2}/m_{\rm s,2}$ (Equation 28) and $\mathcal{P}_{\rm s,2}/m_{\rm s,2}$ (Equation 29) for the speeds of the two single stars and $1/m_{\rm b}\sqrt{\mathcal{P}_{\rm s,1}^{2}+\mathcal{P}_{\rm s,2}^{2}+2\mathcal{P}_{\rm s,1}\mathcal{P}_{\rm s,2}\cos\xi}$ (Equation 27) for that of the binary. For comparison, we also show the isolated (i.e., no background potential) case (solid lines with the same colors as used in Figure 3). Note that the final speeds of the ejected stars are constrained to be within 7% of the speeds shown in the plots. Even when the background potential is taken into account, our analytic formulae show good agreement with the results of the simulations444Note the increased scatter about the analytic predictions in the right panels of Figures 4 and 5. We emphasize that this is due to the simplifying assumptions made for the velocity in deriving the equations of motion. Therefore, a more accurate calculation can easily be done for the binary semi-major axis, provided this assumption is not made and the higher-order terms are included.. As the background potential gets deeper (i.e., higher escape velocities), the binary speed is reduced more. Moreover, given the adopted potentials and corresponding escape velocities, $\xi$ has an upper limit. For example, for $(v_{s,1},v_{s,2})=(70,~{}70)$, the line corresponding to $v_{\rm esc}=20\;\mathrm{km}\;\mathrm{s}^{-1}$ (purple line) drops quickly to zero at $\xi\simeq 110^{\circ}$. This maximum angle refers to the angle required to give a sufficiently high recoil kick to the binary that it can escape to infinity, as imposed by Equation 36. The maximum angle decreases as $v_{\rm esc}$ increases (see also Section 3.4 in the simulation paper).

In Figure 6, we show the relative angles $\Psi$ between the escaped binaries and the faster single stars formed in a background potential. The angles still nicely fulfill the relation with $\xi$ expected from Equation 5 for both equal and unequal single star velocities. In order to avoid overcrowding, we do not mark the simulated data for $v_{\rm esc}=0$ (see the bottom panels in Figure 3 for the case without a background potential). We use different dot types to distinguish the potential models with different depths : $v_{\rm esc}=6.3\;\mathrm{km}\;\mathrm{s}^{-1}$ (square), $v_{\rm esc}=9.2\;\mathrm{km}\;\mathrm{s}^{-1}$ (triangle) and $v_{\rm esc}=20\;\mathrm{km}\;\mathrm{s}^{-1}$ (diamond). The dots and the lines share the same color. For example, the blue square dots refer to the simulated data for $(v_{s,1},v_{s,2})=(100,~{}100)$ for the background potential model with $v_{\rm esc}=6.3\;\mathrm{km}\;\mathrm{s}^{-1}$.

As an example that our formulae apply to any 2+2 encounters, we compare in Figure 7 the final speeds and semimajor axes of the binaries formed from two different 2+2 scenarios, namely between binaries with equal and unequal orbital energies. We show the simulation results for Model 1 and Model 1-1 from the simulation paper. The final speeds of the simulated ejected stars are within 7-10% of the analytic predictions. In both models, we assume the same homogenous potential, with $M_{\rm bg}=3\times 10^{3}$ and $v_{\rm esc}=6.3\;\mathrm{km}\;\mathrm{s}^{-1}$. In Model 1, the two initial binaries have the same binding energy whereas in Model 1-1 the more massive binary is more tightly bound than the less massive binary by a factor of two. For the given combinations of the final speeds for the two ejected single stars, we estimate $v_{\rm b}$ (solid lines) using Equation 31. The hollow circles correspond to the case with equal binding energies ($E_{\rm b,1}/E_{\rm b,2}=1$) and the solid circles correspond to the case with unequal binding energies ($E_{\rm b,1}/E_{\rm b,2}=2$). Regardless of the orbital properties of the initial binaries, our formulae reproduce very well the results of the simulations.

In Figure 9 in Appendix B, we also show the radial distance $r_{\rm b}$ from the core at $t=4\;\mathrm{Myr}$ for both equal $v$ (left column) and unequal $v$ (right column) with/without the background potential using Equation 31.

4 Discussion: Astrophysical applications

Given the ubiquity of a background potential in the astrophysical sites of dynamical interactions, proper theoretical tools that account for its presence are needed. Our analytical method derived above offers one such tool, and is ideal for application to observations of runaway stars.

Given an observed pair of runaway stars with a common origin, our method can be used to identify the left-over binary. First, using Equation 36, we can estimate the probability of observing the left-over binary in its host cluster potential (i.e., to evaluate whether or not it should have escaped). Next, we can infer the history of the background potential of the host cluster. For example, consider the two runaway stars AE Aur and $\mu$ Col, in conjunction with the $\iota$ Orionis binary (all thought to have formed during the same binary-binary interaction; Blaauw & Morgan 1954; Gies & Bolton 1986; Hoogerwerf et al. 2001). Given their observed speeds ($\sim 100\;\mathrm{km}\;\mathrm{s}^{-1}$), their locations on the sky ($\simeq 250\;\mathrm{pc}$ from the CM of all four stars) and the relative angle between the velocity vectors of the two ejected single stars ($\xi\simeq 140-150^{\circ}$) (Gualandris et al., 2004), Equation 36 (assuming $r_{\rm ej}\simeq 0$) gives $v_{\rm esc}>15.5\;\mathrm{km}\;\mathrm{s}^{-1}$ at the time of the last ejection event. Assuming the background potential has a total mass $M_{\rm bg}=10^{4}\;\mathrm{M}_{\odot}$ (Tielens & Hollenbach, 1985; Hillenbrand & Hartmann, 1998) and $\rho=1800-3000\;\mathrm{M}_{\odot}/\;\mathrm{pc}^{3}$ (Herbig & Terndrup, 1986) for the Trapezium cluster, the escape velocity of the cluster at present (assuming a homogeneous medium) is,

[TABLE]

Since we know that the $\iota$ Orionis binary resides in the Trapezium cluster, this suggests that $2.4\;\mathrm{Myr}$ ago (i.e., when the two runaway stars formed) the Trapezium cluster was denser and/or more massive than it is now. However, we note that this prediction is based on a number of simplified assumptions, such as a homogeneous potential in the core region. For a more precise estimation, we have to take into account a more realistic (non-homogeneous) potential in the core in addition to the subsequent dynamical evolution of the cluster, including the gas-depletion process (Boily & Kroupa, 2003), evolution in an external tidal field (Baumgardt & Makino, 2003) and mass segregation (Bonnell & Davies, 1998; Portegies Zwart et al., 2010).

Second, using the temporal information extracted from the observed positions and velocities of the runaway stars, our method can be used to constrain the lifetimes of massive stars. What’s more, in the event that one of the components of the left-over binary has ended its life to form a black hole (BH) via direct collapse, then our method can be used to constrain the initial-final mass relation for massive stellar remnants. In this case, our method provides a relation between the initial mass of the BH progenitor and the initial binary semi-major axis. The most probable combination of these two parameters can then be constrained using further numerical scattering simulations, by evaluating the probability of forming the observed single star velocities for every allowed combination of the initial progenitor mass and binary semi-major axis.

Going one step further, if both binary components form BHs and merge (or will BHs in the first place), their inspiral could be detectable via the associated GW emission. Our method could then be applied to identify an associated pair of runaway stars, and further constrain the properties of the BH-BH binary before merger.

In the derivations above, we consider a particular outcome consisting of two single stars and a binary. However, our analytical treatment is applicable to any kind of three-body outcome (two ejected systems and one left-over system), such as triple-binary encounters that produce two binaries and one single star (e.g. Leigh & Geller, 2013).

We also note that, in reality, it may be more likely that two objects (single, binary or higher-order hierarchical system) encounter each other with non-zero centre of mass velocities. However, when it comes to runaway stars, it is often still valid to treat the ejected objects in their CM frame and to use the analytic formula derived above. For encounters in clusters, the center of mass velocities of the objects (with respect to the CM of the cluster) should be roughly comparable to the stellar velocity dispersion. Hence, for clusters with low velocity dispersions relative to the observed speeds of the runaway stars, the centre of mass velocities can be ignored.

5 Summary

In this paper, we have studied binary-binary interactions, focusing on the formation of one binary and two single stars (i.e., the 2+1+1 outcome). Assuming only conservation of momentum in the inertial frame of each ejection event, we have derived an analytic relation between the orbital properties of the final binary (speed and semimajor axis) and the relative angles between the two ejected stars. This is first done assuming that the interactions occur in isolation (see Section 2.1). Going one step further, we have also found an analytic formula to account for the effects of a homogeneous background potential on the subsequent motions of the stars (see Section 2.2). Then we apply our formula to compare with the results of numerical simulations taken from a companion simulation paper. We have found an overall good agreement with the numerical simulations.

We summarize the results of this analytic study as follows:

For purely stellar interactions, when a binary and two single stars form, the speeds of the binary $v_{\rm b}$ and the relative angle $\xi$ of two single stars are described as follows,

[TABLE]

or, even more simply,

[TABLE]

where $p_{\rm s,1}$, $p_{\rm s,2}$ and $p_{\rm b}$ are the momenta of the first, second ejected stars and the binary, respectively, i.e., $p_{i}=m_{i}v_{i}$. The relative angle $\Psi$ between the binary and the first ejected star can be written as,

[TABLE] 2. 2.

Since we consider the velocities at the time of the last ejection event, these formulae can be applied independently of the details of the initial binary-binary encounter. In Section 3.2, we compared $v_{\rm b}$ as estimated from the analytic formula to the results of simulations of binary-binary scatterings assuming equal/unequal binding energies (Figure 7). We showed that our formula matches the simulated data in all cases. 3. 3.

The formula is symmetric upon exchanging the momenta of the two ejected stars. In our derivations, we differentiate the two single stars in terms of the order of ejection events (i.e., first/second ejected star). But the order of ejection events is irrelevant for estimating the final binary velocity. We have shown this in Figure 3 (green line). This symmetric property of the 2+1+1 outcome makes comparisons to observations straightforward, since the order of ejection events does not need to be known. 4. 4.

For a background potential of total mass $M_{\rm bg}$ with outer boundary $r_{\rm bg}$ and escape velocity $v_{\rm esc}$, given two observed runaway stars at $t=t_{\rm obs}$ with velocities measured in the CM frame of all four stars (i.e., their speeds are $v_{\rm obs,1}$, $v_{\rm obs,2}$ and the relative angle is $\xi$), the speed of the final binary at the outer boundary of the potential (i.e., $r=r_{\rm bg}$ and $t=t_{\rm bg}$) can be expressed as follows,

[TABLE]

where

[TABLE]

Here $m_{\rm s,1},m_{\rm s,2}$ and $m_{\rm b}$ are, respectively, the masses of the two single stars and the binary. Finally, the speed and the location of a binary which has escaped from the potential at $t~{}(>t_{\rm bg}$) can be estimated as,

[TABLE]

[TABLE]

where $\mathfrak{R}=r_{\rm bg}/[(v_{\rm b,bg}/v_{\rm esc})^{2}-1]$. We have also derived the equations of motion for the stars while they are still moving outward in the potential (see Equation 2.2.2). 5. 5.

Given the derived formulae, we have found a condition for the complete escape of a binary from the potential,

[TABLE] 6. 6.

In summary, we have derived analytic formulae which can be applied to three-body outcomes (two ejected systems and one left-over system) with/without a background potential. We expect that these formulae can be applied directly to observations of runway stars, and hence to be particularly useful with the expected detection of hundreds of runaway stars with the GAIA satellite.

Acknowledgements

The authors are grateful to the referee, Joseph Antognini, for his very meaningful comments and suggestions which helped improve the presentation of our paper. Results in this paper were obtained using the high-performance LIred computing system at the Institute for Advanced Computational Science at Stony Brook University, which was obtained through the Empire State Development grant NYS #28451.

Appendix A Derivation for $v_{\rm b}$ with vector components-purely stellar dynamics

When two binaries collide, a chaotic interaction ensues until two single stars are ejected sequentially, leaving behind a binary pair. For each ejection event, we adopt the CM frame of the ejected star and the remaining system (i.e., excluding any previously ejected stars). This derivation could prove useful for identifying the first ejected single star formed during a binary-binary interaction, when applying our method to observations of runaways where only one single star and an associated binary have been identified.

After the first ejection event, a single star of mass $m_{\rm s,1}$ and the remaining system of mass $m_{\rm t}$ recede in opposite directions (in the CM frame of all four stars). We assume a velocity of $\textbf{{v}}_{\rm s,1}$ for the first ejected single star, directed along the $(-z)$ axis from the CM of all four stars (denoted frame $S$). The remaining $N=3$ system recedes with velocity $\textbf{{v}}_{\rm t}=(m_{\rm s,1}/m_{\rm t})\textbf{{v}}_{\rm s,1}$ along the $z$-axis. At the time of the second ejection event, another single star is ejected with velocity $\textbf{{v}}^{\prime}_{\rm s,2}$ in the $x-z$ plane in the CM frame of all three remaining stars (denoted frame $S^{\prime}$, with its $z^{\prime}-$axis oriented such that it coincides with the $z-$axis in frame $S$). The final binary recedes in the opposite direction with velocity $\textbf{{v}}^{\prime}_{\rm b}=-(m_{\rm s,2}/m_{\rm b})\textbf{{v}}^{\prime}_{\rm s,2}$ relative to the second ejected single star in frame S. Hereafter, the velocities defined in frame S’ will be denoted with a prime symbol (’). We present a schematic diagram in Figure 8 showing the two ejection events in the two different reference frames. The velocities of both single stars and the final binary can be expressed in frame $S$ using polar coordinates as follows,

[TABLE]

where $\theta$ is the polar angle relative to the $z^{\prime}-$axis. Accordingly, their speeds in frame $S$ are,

[TABLE]

where $\alpha=\cos\theta$. Given the preceding velocities, we can estimate the relative angle $\xi$ between the two ejected single stars,

[TABLE]

where ${A}=v_{\rm t}/v^{\prime}_{\rm s,2}$.

Next, we motivate our definitions for the angles $\xi$ and $\alpha$. We have introduced two inertial frames of reference ($S$ and $S^{\prime}$). In frame $S$ (i.e., the CM frame of all four stars), $\xi$ corresponds to the relative angle between the velocity vectors of the two ejected single stars. The variable $\alpha$, on the other hand, is defined in the $S^{\prime}$ frame (i.e., the CM frame of the remaining system of three stars left-over after the first ejection event). It measures the relative angle between the velocity vector of the second ejected star and the CM velocity vector of the remaining (temporary) $N=$ 3 system (i.e., frame $S^{\prime}$ itself, or the direction opposite to the direction of motion of the first ejected single star). In what follows, we will identify every allowed combination of these two angles.

In order to find an expression for $v_{\rm b}$ as a function of $v_{\rm s,1},~{}v_{\rm s,2}$, $\xi$ and the masses $m$ of the stars, i.e., $v_{\rm b}=v_{\rm b}(v_{\rm s,1},~{}v_{\rm s,2},~{}\xi,~{}m)$, we combine Equations 52b and A and solve for $\alpha$ as a function of $\xi$. Hereafter, we will use the vector $m$ to refer to the masses of two ejected single stars and one binary collectively, or $(m_{\rm s,1},~{}m_{\rm s,2},~{}m_{\rm b})$ unless it is necessary to denote each of them separately. Plugging $\alpha(\xi)$ back into equation A, we obtain ${A}=$ ${A}(\xi,v_{\rm s,2})$. Finally, using the relation $v^{\prime}_{\rm b}=(m_{\rm s,2}/m_{\rm b})v^{\prime}_{\rm s,2}$ in conjunction with Equation 52c, we obtain $v_{\rm b}=v_{\rm b}(v_{\rm s,1},~{}v_{\rm s,2},~{}\xi,~{}m)$. For tractability of the newly defined parameters, we use Greek letters (e.g., $\alpha$ and $\xi$) for angle-related parameters and capital Roman letters (e.g., ${A}$ and ${B}$) for velocity- (momentum-) or distance-related parameters.

Now, we start with calculating an expression for ${A}$ from Equation A as follows,

[TABLE]

We can get another expression for ${A}$ from Equation 52b. That is, solving for $v^{\prime}_{\rm s,2}$ in Equation 52b,

[TABLE]

where we define ${B}=v_{\rm s,2}/v_{\rm t}$. Using the relation ${A}=v_{\rm t}/v^{\prime}_{\rm s,2}$, we find,

[TABLE]

Equating Equations 54 and 56,

[TABLE]

Given the equation above, we can find an expression for $\alpha$ in terms of $\xi$ and ${B}$. There are four possible solutions for $\alpha(\xi,{B})$ (or two, adopting a different range for $\xi$ and ${B}$):

[TABLE]

The final relative orientation between the two ejected single stars, characterized by the parameters $\xi$ and $\alpha$, allow or forbid different combinations of the plus and minus signs (in each closed bracket) in Equation 57. This, in turn, specifies the solution for $\alpha(\xi,{B})$, of the four possible solutions in Equation 58. Hereafter, we consider two separate ranges for $\xi$, i.e., $0<\xi<\uppi/2$ and $\uppi/2<\xi<\uppi$. For each range of $\xi$, we must further distinguish between all possible cases using $\alpha$ (which can be positive/negative). We will henceforth use the "same (opposite) direction" to refer to cases where the relative angle between the final velocity vectors of the two single stars is smaller (larger) than $\uppi/2$.

A.1 Case 1. $\frac{\uppi}{2}<\xi\leq\uppi$ and positive $\alpha$ ($0\leq\theta\leq\frac{\pi}{2}$)

In Mode 1, a single star is ejected in the same direction as the CM velocity of frame $S$ ($0\leq\theta\leq\frac{\pi}{2}$) and in the opposite direction with respect to the first ejected star ($\textbf{{v}}_{\rm s,2}\cdot\textbf{{e}}_{\rm z}>0$ for $\frac{\uppi}{2}<\xi\leq\uppi$ where $\textbf{{e}}_{z}$ is the unit vector in the $z$ direction). For this case, there are two possible directions of motion for the final binary, depending on the relative magnitudes of the speed of the binary and the CM speed of frame $S^{\prime}$. In other words, if the binary is kicked in the $-z^{\prime}$ direction at a speed smaller (higher) than the CM speed of frame $S^{\prime}$, then in frame $S$ the binary is seen as moving in the $z(-z)$ direction. These geometry-based arguments allow for only one sign for both ${A}$ and $v^{\prime}_{\rm s,2}$. If ${A}>0$, $v^{\prime}_{\rm s,2}>0$, $\cot\xi<0$ and $\alpha\geq 0$, then the only allowed expressions for ${A}$ and $v^{\prime}_{\rm s,2}$ are,

[TABLE]

[TABLE]

Given a real number for $v^{\prime}_{\rm s,2}$, we impose the condition that $\alpha^{2}-1+{B}^{2}>0$. We will return to this condition later (Condition 1 below). Equation 57 now reduces to,

[TABLE]

The corresponding solutions for $\alpha(\xi,{B})$ are:

For ${B}\geq 1$,

[TABLE]

For ${B}<1$,

[TABLE]

Since $\alpha(\xi,{B})\geq 0$, all solutions with a negative sign in front are ruled out. However, before we can solve for $v_{\rm b}$, we must consider a couple of additional conditions. These are designed to ensure the validity of our expression for $\alpha$, and to find a proper range for $\xi$ that is allowed in this case. For this and every subsequent case, we will refer to the following conditions :

$\alpha^{2}-1+{B}^{2}\geq 0$.

This is to ensure that $v^{\prime}_{\rm s,2}$ is a real number. In case 1, this condition is always satisfied,

[TABLE]

This is satisfied for all solutions (Equation 58)555In general, for $\alpha\sim\frac{|1\pm{B}~{}\cos\xi|}{\sqrt{1\pm{B}^{2}+2{B}\cos\xi}}$,

\displaystyle\alpha^{2}-1+{B}^{2}={B}^{2}\Big{[}\frac{\cos^{2}\xi+{B}^{2}\pm 2{B}\cos\xi}{1+{B}^{2}\pm 2{B}\cos\xi}\Big{]}>0

(66)

. Therefore, we assume $v^{\prime}_{\rm s,2}$ is a real number for all allowed expressions for $\alpha$. 2. 2.

${A}>0$.

For the solution for $\alpha(\xi,{B})$ (Equation 63), it is enough to use either of two different expressions for ${A}$ (Equations 59 and 60) to show that this condition is satisfied. Using Equation 60, it is easy to see that $v^{\prime}_{\rm s,2}$ is positive for ${B}>1$. For ${B}<1$,

[TABLE] 3. 3.

$\textbf{{v}}_{s,2}\cdot\textbf{{e}}_{z}>0$.

This condition is clearly satisfied for positive $\alpha$.

Finally, with Equations 60, 52c and $v_{\rm b}=(m_{\rm s,2}/m_{\rm b})v_{\rm s,2}$, we are able to find a solution for $v^{\prime}_{\rm b}$ for Case 1,

[TABLE]

for $\cos\xi\leq\frac{-1}{{B}}$ and ${B}>1$.

A.2 Case 2. $\frac{\uppi}{2}\leq\xi\leq\uppi$ and negative $\alpha$ ($\frac{\pi}{2}<\theta\leq\pi$)

In this case, both the final binary ($-1\leq\alpha<0$) and the second ejected star ($\frac{\uppi}{2}\leq\xi\leq\uppi$) are moving in the opposite direction (or $\textbf{{v}}_{s,2}\cdot\textbf{{e}}_{z}>0$) to the first ejected star, in frame $S$. This geometric property allows for both signs for both ${A}$ and $v^{\prime}_{\rm s,2}$. Hence, it is necessary to consider each combination of the signs individually. We refer to each combination as (+,+), (+,-), (-,+) and (-,-), respectively.

A.2.1 (+,+) case

For negative $\alpha$, we know that both ${A}$ and $v^{\prime}_{\rm s,2}$ are positive,

[TABLE]

Accordingly, Equation 57 becomes

[TABLE]

There is only one solution for the equation above, namely

[TABLE]

where we have used $\cos\xi\leq 0$ in the second equality. As mentioned in Case 1, there are a couple of conditions for $\xi$ which must be satisfied in order to ensure physically meaningful solutions. Condition 1 below needs to be confirmed since in this case the solutions for ${A}$ and $v_{\rm s,2}$ are not manually chosen to be positive. The conditions are:

${A}>0$.

As mentioned before, we can use either of the two expressions for ${A}$ to show that this condition is satisfied. Clearly, ${A}>0$ from Equation 70 for $\alpha<0$. 2. 2.

$\textbf{{v}}_{s,2}\cdot\textbf{{e}}_{z}\geq 0$.

We write this relation as follows,

[TABLE]

We find from the relation above that

[TABLE]

However, for any ${B}$ (and given an expression for $\alpha$; Equation A.2.1), $\alpha$ is always $<-\frac{1}{\sqrt{1+{B}^{2}}}$ 666For $\cos\xi<0$,

$\displaystyle\alpha^{2}-\frac{1}{1+{B}^{2}}$ $\displaystyle\sim(1-{B})^{2}(1+{B}^{2})-(1+{B}^{2}-2{B}\cos\xi)$

$\displaystyle\sim{B}^{2}\cos\xi({B}^{2}\cos\xi+\cos\xi-2{B})\geq 0$

(76)

. Therefore, we can not find solutions for $\alpha$ which satisfy all the required conditions in this case.

A.2.2 (+,-) case

Given these signs for ${A}$ and $v^{\prime}_{\rm s,2}$, we have

[TABLE]

and Equation 57 becomes

[TABLE]

We have different solutions for $\alpha(\xi,{B})$ for ${B}>1$ and ${B}<1$:

For ${B}>1$,

[TABLE]

For ${B}<1$,

[TABLE]

Since the solutions 81 and 83 are positive definite, they can be ruled out due to the fact that $\alpha<0$. In addition, for ${B}>1$ (i.e., where the solution 80 is valid), we exclude solution 80 because $v^{\prime}_{\rm s,2}$ becomes negative. Only solution 83 remains. However, we need to check the following conditions for the remaining solution to ensure it is physically allowed:

${A}>0$.

Starting with Equation 78,

[TABLE]

where in the second equality, we have used the conditions for ${B}$ that ${B}<1$ and ${B}+\cos\xi>0$. For these conditions, $1+{B}\cos\xi>0$. 2. 2.

$\textbf{{v}}_{s,2}\cdot\textbf{{e}}_{z}>0$.

We write this relation as follows,

[TABLE]

Re-arranging terms, we find that

[TABLE]

Inserting Equation 82 into Equation 86, it turns out that this condition is always satisfied since

[TABLE]

Hence, we have only one solution for this case,

[TABLE]

[TABLE]

for $0<{B}<1$ and $-1<\cos\xi<-{B}$.

A.2.3 (-,+) case

Given these signs for ${A}$ and $v_{\rm s,2}$,

[TABLE]

Equation 57 becomes

[TABLE]

We have obtained the same solutions as for the (+,-) case, but with different ranges for $\xi$:

For ${B}\geq 1$,

[TABLE]

For ${B}<1$,

[TABLE]

Solutions 94 and 96 (positive definite) can be ruled out if $\alpha<0$. Note that, for ${B}<1$, $1+{B}\cos\xi>0$ (the numerator of solution 95). As before, the conditions that must be satisfied to ensure the validity of the solutions are:

${A}>0$.

This condition is already fulfilled given the choice of signs in Equations 90 and 91. 2. 2.

$\textbf{{v}}_{s,2}\cdot\textbf{{e}}_{x}\geq 0$.

Given the same expressions as for Equations 2, 74 and 75,

the necessary conditions are described in terms of $\xi$,

[TABLE]

Combining the allowed ranges for $\xi$ (see Equations 93 and 95), we have different conditions for $\xi$ depending on whether or not ${B}$ is larger than 1, i.e.,

[TABLE]

Note that provided $-(2{B})/(1+{B}^{2})<-{B}$ for $0\leq{B}<1$, the allowed range for $\cos\xi$ is stringently determined by ${-B}$, i.e., $-B\leq\cos\xi$777. A hierarchical comparison of the magnitudes of $\frac{2{B}}{1+{B}^{2}}$, $\frac{1}{{B}}$ and ${B}$ is as follows,

$\displaystyle\begin{cases}\frac{1}{{B}}>\frac{2{B}}{1+{B}^{2}}>{B}\hskip 18.06749pt\text{for}\,\,\,\ 0<{B}<1,\\ {B}\geq\frac{2{B}}{1+{B}^{2}}\geq\frac{1}{{B}}\hskip 18.06749pt\text{for}\,\,\,\,{B}\geq 1,\\ \end{cases}$

(100)

. Finally, the solutions can be written,

[TABLE]

[TABLE]

for

[TABLE]

A.2.4 (-,-) case

In this case, there is no solution for negative values of $\alpha$ that satisfy Equation 57.

A.3 Case 3. $0\leq\xi<\frac{\uppi}{2}$ and positive $\alpha$ ($0\leq\theta<\frac{\pi}{2}$)

In this case, the ranges for $\xi$ and $\alpha$ are incompatible. For positive $\alpha$, the second ejected star is moving in the same direction as frame $S$ itself (or, the positive $z$ component of $\textbf{{v}}_{\rm s,2}$). Therefore, $\textbf{{v}}_{s,2}\cdot\textbf{{e}}_{z}$ should be positive, which forbids the range $0\leq\xi<\frac{\uppi}{2}$. Therefore, this case is not physically allowed.

A.4 Case 4. $0\leq\xi\leq\frac{\uppi}{2}$ and negative $\alpha$ ($\frac{\pi}{2}\leq\theta\leq\pi$)

In this case, the second ejected star is moving (in frame $S^{\prime}$) in the direction opposite to the CM motion of frame $S^{\prime}$ ($\frac{\pi}{2}\leq\theta\leq\pi$). In frame S, the two ejected stars recede in the same direction (-z direction, $0\leq\xi<\frac{\uppi}{2}$). The binary, however, is recoiled in the opposite direction (z direction, $\textbf{{v}}_{s,2}\cdot\textbf{{e}}_{z}>0$). In this case, it is more likely for the binary to move at high speeds relative to the other cases. As for case 2, since it is possible for ${A}$ and $r_{\rm s,2}$ to have either sign, we have to consider each case individually.

A.4.1 (+,+) case

We take the positive signs in both Equations 54 and A,

[TABLE]

[TABLE]

We have two sets of solutions for different ranges of ${B}$,

For ${B}\geq 1$,

[TABLE]

and for ${B}<1$,

[TABLE]

The condition $-1\leq\alpha<0$ excludes the solutions 106 and 108. As before, we must check the following conditions:

${A}>0$.

This condition is already satisfied given the choice of signs in Equations 103 and 104. 2. 2.

$\textbf{{v}}_{s,2}\cdot\textbf{{e}}_{z}\leq 0$.

We already derived an expression for $\textbf{{v}}_{s,2}\cdot\textbf{{e}}_{z}$ (see Equations 2, 74 and 75). The corresponding condition can be written as,

[TABLE]

We find that

[TABLE]

Here, for ${B}<1$, Equation 109 imposes one other constraint such that $\alpha<-\sqrt{1-{B}^{2}}$. This condition is less strict compared to the above since $-\frac{1}{\sqrt{1+{B}^{2}}}<-\sqrt{1-{B}^{2}}$. The conditions above are described in terms of $\xi$,

[TABLE]

Here, we find that Equation 111 is in direct conflict with the conditions for $\xi$ identified by 105 and 107 (see footnote 7). Therefore, no solutions exist in this case.

A.4.2 (+,-) case

In this case, solutions with negative $\alpha$ do not exist.

A.4.3 (-,+) case

Given the choice of the signs in the expressions for ${A}$ and $v^{\prime}_{\rm s,2}$,

[TABLE]

[TABLE]

there exists a solution with negative $\alpha$,

[TABLE]

Accordingly, the conditions for ${A}$ and $\xi$ are,

${A}>0$

Clearly, $v^{\prime}_{\rm s,2}$ (Equation 113) is positive for $\alpha<0$. 2. 2.

$\textbf{{v}}_{s,2}\cdot\textbf{{e}}_{z}\leq 0$.

This is satisfied for any range of $\xi$. This can be understood as follows. When we compare $\alpha(\xi,{B})$ in Equations 114 and $-1/\sqrt{1+{B}^{2}}$ in Equation 110 (still valid in this case), they are the same at $\cos\xi=0$. Since $\alpha(\xi,{B})$ is a monotonically decreasing function of $\cos\xi$ for fixed ${B}$ (as $\cos\xi$ increases from 0 to 1, $\alpha$ decreases), $\alpha$ is always smaller than $-1/\sqrt{1+{B}^{2}}$, which means Equation 110 is fulfilled for any $\xi$.

Overall, in this case, we have one solution that is valid for any value of ${B}$, which is,

[TABLE]

[TABLE]

A.4.4 (-,-) case

With negative signs in both of the expressions for ${A}$ and $v^{\prime}_{\rm s,2}$,

[TABLE]

[TABLE]

we obtain the following solutions for $\alpha(\xi,{B})$:

For ${B}\geq 1$,

[TABLE]

and for ${B}<1$,

[TABLE]

Equation 119 only satisfies the required condition if $\alpha<0$. For the same reason as in Case 2 (+,-), for ${B}>1$, $v^{\prime}_{\rm s,1}$ becomes negative, such that no solutions exist in this case.

After exploring all cases, we have found one simple expression for $v_{\rm b}$ which is valid for all ranges of ${B}$ and $\xi$,

[TABLE]

Or even more simply,

[TABLE]

where $p_{\rm s,1}$, $p_{\rm s,2}$ and $p_{\rm b}$ are the momenta of the first, second ejected stars and the binary, respectively, i.e., $p_{i}=m_{i}v_{i}$. We summarize the physically-allowed solutions for each case in Table 1.

Appendix B The radial distance $r_{\rm b}$ as a function of $\xi$