Free Access
Issue
A&A
Volume 504, Number 1, September II 2009
Page(s) 277 - 289
Section Celestial mechanics and astrometry
DOI https://doi.org/10.1051/0004-6361/200809392
Published online 02 July 2009

Adapting Marchal's test of escape to real triple stars

P. J. Li1 - Y. N. Fu2 - Y. S. Sun1

1 - Department of Astronomy, Nanjing University, 22 Hankou Road, Nanjing 210093, PR China
2 - Purple Mountain Observatory, Chinese Academy of Sciences, 2 West Beijing Road, Nanjing 210008, PR China

Received 14 January 2008 / Accepted 11 May 2009

Abstract
Context. For a general N-body system, Marchal constructed an analytical test of escape, which uses only a one-dimensional projected motion state of the system at any given instant. This test is well adapted to identifying real, disintegrating small stellar systems, of which the full motion states are generally unavailable. However, to our knowledge, there has been no practical application of this test until the present-day.
Aims. In this paper, we aim at adapting the above test to visual triple stars with estimable component masses and known kinematic data on the plane perpendicular to the line-of-sight. As illustrating examples, our goal is to identify disintegrating Hipparcos linear triple systems.
Methods. The fundamental techniques of analytical geometry were used to adapt the test of escape to practical applications, and the Monte Carlo method used to cope with the unavoidable observational errors, so that the confidence probability of a real triple star disintegrating could be obtained.
Results. A practical algorithm was designed to make full use of the two-dimensional kinematic data in testing usual visual triple stars. This algorithm is then applied to 24 Hipparcos linear triple systems with estimable component masses and the disintegration probability given.

Key words: celestial mechanics - instabilities - stellar dynamics - binaries: general - methods: analytical

1 Introduction

By the well-known Lagrange-Jacobi equation, we know that escape is unavoidable for an N-body system with non-negative energy. When the energy is negative, bounded motions exist, e.g. homographic and choreographic solutions (e.g. Chenciner 2007). However, extensive numerical experiments show that escape, usually associated with a close triple approach (e.g. Agekian & Anosova 1990; Aarseth et al. 1994), remains a dominant phenomenon (e.g., Anosova & Orlov 1994; Kiseleva et al. 1994; Valtonen et al. 2005, and references therein).

For the case of N=3, several analytical tests of escape exist that are valid for systems with negative energy (Birkhoff 1927; Tevzadze 1962; Standish 1971; Griffith & North 1974; Yoshida 1972; Marchal 1974; Kuznetsova & Orlov 1983; Marchal et al. 1984a,b; for a recent review, see Georgakarakos 2008). Underlying these tests, the common idea is to confine two bodies to a certain region so that their joint attraction on the third body can be bounded from above.

Apparently, it is not easy to extend the above idea to much more complicated cases where N>3. As far as we know, the only kind of analytical test of escape for an N-body system, N>3, with negative energy stems from a different idea. By analyzing the projected motion of the general N-body system on a given axis, Marchal (1990) obtained several different versions of a sufficient condition for escape, each of which corresponds to a balance between simplicity and efficiency.

The condition mentioned above involves only a one-dimensional projected motion state of the system at any given instant. This makes it well-adapted to studying real stellar systems, multiples in particular, of which the full motion states are generally unavailable. Indeed, directly determining the projected mutual distances along the line-of-sight between component stars of multiples remains difficult for present-day astrometry. Though there is a common practice for determining the motion state of a hierarchical multiple by ``binary'' orbit fitting, the result is questionable if the considered multiple is actually in the disintegration phase. Moreover, for systems displaying almost no orbital motion during the observation periods, such as the existing Hipparcos or the coming Gaia linear systems, the observational data can only be used to determine the linear motion states of the component stars, which are effectively the instant ones at the middle of the observational time intervals. This would make Marchal's test a preferable, if not indispensable, tool for identifying escape phenomena in the above systems. In fact, the overall precisions of the kinematic data will otherwise be significantly reduced by incorporating much less precise long-term observations.

In his 1990 book entitled The Three-Body Problem, Marchal mentions the possible application of his test of escape to the Galaxy ``for all directions of the galactic plane''. On the other hand, however, we have found no practical application of the test until now. It would be too ambitious for the time being to apply the test to the whole Galaxy, but it is practical if we consider multiple stars. At this point, it should also be noted that an escape event is much more significant for a multiple star than for the Galaxy. In particular, such an event occurring in a triple star can be equally considered as the disintegration of the system. Restricting ourselves to this simple yet interesting case, we aim at adapting the above test to real applications. And as examples, our goal is to identify disintegrating Hipparcos linear triple systems with estimable masses and known kinematic data on the tangent plane, that is, the plane perpendicular to the line-of-sight.

In Sect. 2, we briefly recall Marchal's test. And in Sect. 3, we initiate the problem of how to adapt this test to the study of visual triple stars. This problem is solved in Appendix A, where a special-purpose algorithm is presented. Section 4 discusses the relevance of the algorithm with the help of 3-body experiments. Section 5 is devoted to applications to Hipparcos linear triple systems, and discussions are given in Sect. 6.

2 Marchal's test of escape

Consider the Newtonian N-body problem in the barycentric coordinate system. The equations of motion are

\begin{displaymath}m_i \frac{{\rm d}^2 \vec{r}_i}{{\rm d}t^2}=\sum\limits_{j=1,j...
...r}_i)}
{\vert\vec{r}_j-\vec{r}_i\vert^3},\quad
(i=1,2,...,N),
\end{displaymath}

where mi and $\vec{r}_i$ are, respectively, the mass and the position vector of the ith body, and t the time variable, Gthe gravitational constant. Together with the the above notations, the following usual ones are used without further explanation

\begin{displaymath}\begin{array}{ll}
M=\sum\limits_{i=1}^Nm_i,\\
M^*=\sum\lim...
...tyle{\frac{m_im_j}{\vert\vec{r}_i-\vec{r}_j\vert}},
\end{array}\end{displaymath}

where M is the total mass, M* a quantity depending on the masses, $\rho$ the mean quadratic distance, I the moment of inertia, $\nu$ the mean harmonic distance, and U the force function.

For the purpose of constructing a test of escape for a general N-body system, Marchal analyzed the motion of the system along an arbitrarily chosen x-axis. We follow Marchal by relabeling the bodies such that their x-coordinates are in a non-decreasing order. This implies that (m1,...,mN) is no longer a fixed N-uple. In fact, the subscript associated with a particular body may be changed when the considered body has the same x-coordinate with one or more other bodies. At this point, it should be noted that, in the following, we do not attempt to test solutions that admit x1=xN, i.e. $\vec{x} = (x_1,x_2,...,x_N)=0$, at a given future instant. Bearing in mind that this kind of solution, if any, must be removed from the ones passing the Marchal's test of escape introduced below, we exclude them from the following consideration; that is, we assume xN-x1>0 for future time.

To describe an escape phenomenon of a subsystem with respect to another subsystem, one needs a characteristic distance between these two subsystems. Considering such a phenomenon along the x-axis, one divides the whole system into two ``subsystems'' $\{m_{1},...,m_{s}\}$ and $\{m_{s+1},...,m_{N}\}$, $s \in
\{1,2,...,N-1\}$. Even for a given s, this division is generally temporary in the sense that such a subsystem may be composed of different bodies at different instants. As a matter of fact, using this temporary division is the key idea in avoiding the difficulty of confining a permanent subsystem of a multi-body system to a certain region. Since one actually cannot predefine s such that the corresponding two subsystems become permanent and finally escape from one another, one first introduces a characteristic distance for any possible s and any possible temporary ordering. This is simply the distance between the mass centers of the two temporary subsystems,

\begin{displaymath}\lambda_{s}(\vec{x})=\frac{1}{M-M_{s}}\sum\limits_{j=s+1}^Nm_{j}x_j-\frac{1}{M_{s}}\sum\limits_{i=1}^sm_{i}x_i,
\end{displaymath}

where $M_{s}=\sum_{i=1}^sm_{i}$ is the total mass of the temporary subsystem $\{m_{1},...,m_{s}\}$. The assumption xN-x1>0 ensures that $\lambda_{s}>0$ for any s, and one has

\begin{displaymath}
\frac{{\rm d}^2\lambda_s}{{\rm d}t^2}\ge-\frac{GM}{\lambda_s^2}A_s(\vec{\kappa}(\vec{x})),
\end{displaymath} (1)

where $\vec{\kappa}=(\kappa_1,\kappa_2,...,\kappa_N)=\vec{x}/\lambda_s$, and

\begin{displaymath}A_{s}(\vec{\kappa})=\frac{1}{M_{s}(M-M_{s})}\sum_{l=1}^{s}\sum_{k=s+1}^{N}m_{l}m_{k}(\kappa_{k}-\kappa_{l})^{-2}.
\end{displaymath}

We note that (1) becomes an equality if and only if all bodies are on the x-axis. This implies that, in order to bound d $^2\lambda_s$/dt2 from below, one can always assume that all bodies are on the x-axis. Since $\lambda_s$ is linearly dependent on $\vec{x}$, $\vec{\kappa}$ is a scale-independent N-uple. Therefore, As can be considered as a function defined on the quotient collinear configuration space, $\mathcal{C}$, obtained by identifying collinear configurations with the same ordering and shape. For any given ordering, a configuration in $\mathcal{C}$ is specified by the homogeneous coordinate $\vec{x}$. In the following, we denote the set of all possible orderings by O, and o an arbitrary element of O. This means that a general element of $\mathcal{C}$ can be specified by the pair $(o,\vec{x})$. Now it is clear that the quantity A depends on the triple $(o,\vec{x},s)$. An ordering and its reverse, though geometrically the same, are not equivalent to each other in the present text, because the direction of x-axis has been fixed. Nevertheless, the triple $(o,\vec{x},s)$is equivalent to $(-o,\vec{\tilde{x}},N-s)$, where $
\vec{\tilde{x}}=(\tilde{x}_1,...,\tilde{x}_N)=(-x_N,...,-x_1)$ and -o represents the reverse of o.

With the length $\lambda=\sup\{\lambda_1,\lambda_2,\ldots,\lambda_{N-1}\}$, for any given $(o,\vec{x}) \in \mathcal{C}$, there must be at least one sfor which $\lambda=\lambda_s$. Therefore, $\lim_{t \rightarrow
\infty}\lambda \rightarrow \infty$ implies that the system will eventually disintegrate into at least two permanent subsystems. This means that $\lambda$ can be used as the aforementioned characteristic distance describing the escape phenomenon.

By (1), for any given value of $\lambda$, the acceleration of this length can be limited from below as

\begin{displaymath}
\frac{{\rm d}^2\lambda}{{\rm d}t^2} \geq -\frac{GM}{\lambda^2}\bar{A},
\end{displaymath} (2)

where $\bar{A} \equiv \sup_{(o,\vec{x},s)~\in
\Omega}A(o,\vec{x},s)$, with $\Omega=\mathcal{C} \times
S(o,\vec{x})$. Here, $S(o,\vec{x}) \subset \{1,...,N-1\}$ collects all s satisfying $\lambda_s=\lambda$. It can be verified that

\begin{displaymath}
\bar{A}=\frac{M^2}{M-m_1}\sum_{k=2}^N\frac{m_k}{[m_1+2(m_2+...+m_{k-1})+m_k]^2}\cdot
\end{displaymath} (3)

This least upper-bound $\bar{A}$ is attained at $(o,\vec{x},s)=(\hat{o},\hat{\vec{x}},1)$, where the configuration $(\hat{o},\hat{\vec{x}})$ is specified by the following two conditions: the bodies are arranged in a non-decreasing mass order along the x-axis, and all $\lambda_s$ are equal. The above discussions lead to the following conclusion,

Proposition.

(Marchal 1990) If at a given instant $\tau$

\begin{displaymath}
{\rm d}\lambda/{\rm d}t\ge\sqrt{2GM\bar{A}/\lambda},
\end{displaymath} (4)

then this inequality is true for all $t>\tau$.

Based on this proposition, one can divide, for any given instant, the phase space into two complementary regions according to whether (4) is satisfied or not. In the former region, hereafter escape region, $\lambda$ goes to infinity with $t\rightarrow\infty$ at least as $(9GM\bar{A}/2)^{1/3}(t-\tau)^{2/3}$. This makes (4) a test of escape (hereafter M90).

To complete, we point out that the previous assumption, i.e. xN-x1>0 for all future time, is true if the considered solution starts from the escape region. This is because xN-x1>0 is true in the escape region, and it is impossible for a solution to go from this region to its complementary region.

3 Applying M90 to a spatial 3-body system with motion state available on a plane

 \begin{figure}
\par\includegraphics[width=9cm,clip]{09392fg1.eps}
\end{figure} Figure 1:

Critical x-axes separating different pieces of the piecewise function $q(\hat{\vec{x}})$. The symbol $\odot $, situated at the barycenter Q, indicates that $\vec{\zeta}$ is perpendicular to the paper and is directed toward the reader.

Open with DEXTER

We consider a spatial 3-body system for which the position and velocity vectors are only known by their projections onto a given fixed plane P, which is the case of a visual triple star with known motion state on the plane of the sky. To be more specific, we consider the system in a right handed rectangular coordinate system $\{Q,\vec{\xi}\vec{\eta}\vec{\zeta}\}$, with its origin the barycenter Q and $\vec{\xi}\vec{\eta}$-plane the plane P. Obviously, M90 can be applied to the system with kinematic data on any given axis $\vec{x} \in P$. All such x-axes constitute a one-parameter family $X=\{\hat{\vec{x}}: \hat{\vec{x}} \in P$}, where $\hat{\cdot}$ denotes the unit vector associated with an axis or a vector, and they can be parameterized as $\hat{\vec{x}}=R_{\alpha}(\hat{\vec{\xi}})$ with $\alpha \in
[0,2\pi)$, where the operator $R_{\alpha}(\cdot)$ rotates an axis or a vector by the angle $\alpha $ around $\vec{\zeta}$ according to the righthand rule. As is obvious, to make the full use of the known information in applying M90, it is sufficient for us to find the optimal upper bound of the $\alpha $-dependent quantity $q(\alpha)=\sqrt{\lambda}\frac{{\rm d}\lambda}{{\rm d}t}$.

In Fig. 1, the projected triangle on P of the 3-body configuration is shown, together with a representative x-axis and some specific ones on P. The reference axis $\vec{\xi}$ is chosen as $\vec{i}_{\rm a}$, i.e. the one perpendicularly pointing from Q to the side $\overline{m_bm_c}$ of the triangle. The axis $\vec{j}_{\rm b}$is defined as the specific x-axis on which the collinear configuration $\{x_{\rm a},x_{\rm b},x_{\rm c}\}$ satisfies two conditions: $x_{\rm a}<x_{\rm b}<x_{\rm c}$ and $x_{\rm bc}-x_{\rm a}=x_{\rm c}-x_{\rm ab}$, where $x_{\rm bc}$ and $x_{\rm ab}$ are, respectively, the x-coordinates of the mass centers of $\{m_{\rm b},m_{\rm c}\}$ and $\{m_{\rm a},m_{\rm b}\}$. With subscript permutations, the axes $\vec{i}_{\rm b}\ (\vec{i}_{\rm c})$ and $\vec{j}_{\rm c}\ (\vec{j}_{\rm a})$ can be defined in the same way as $\vec{i}_{\rm a}$ and $\vec{j}_{\rm b}$, respectively.

As is clear in Appendix A, $q(\alpha)$ is a piecewise function with 12 pieces, and it is generically not single-valued at $\pm\vec{j}_a,\pm\vec{j}_b,\pm\vec{j}_c$. These facts make it inconvenient to use a generally applied method to find the maximum value of the function $q(\alpha)$. Therefore, we decide to discuss in detail this function, so that a highly efficient special-purpose algorithm[*] can be designed for implementing M90 with kinematic data known on a plane. This rather technical discussion can be found in Appendix A.

It should be noticed that this algorithm will not implement M90 using kinematic data on any axis outside the prescribed plane P. Therefore, when the full 3-dimensional motion state of a general spatial 3-body system is known, a more sophisticated algorithm is necessary in order to make full use of the known data. On the other hand, however, in cases of both a planar 3-body system and a general spatial 3-body system with motion state available only on a plane, implementations of M90 and of our algorithm are equivalent. Since our main concern in this work is the real application of M90 to very distant three celestial body systems, we find ourselves in the latter case. To emphasize that M90 in this case will be applied to a spatial 3-body system by only using kinematic data on a plane, the above-mentioned algorithm will be abbreviated as M903P.

4 Relevance of M903P

We start this section with a few general remarks on M90 in the context of a general spatial 3-body system with a full motion state being available. Consider an escape that is detected by M90 using the system state at a certain instant $\tau$. As mentioned in Sect. 2, for any axis with which the detection is successful, the characteristic length $\lambda(t)\ (t\geq\tau)$ will be limited from below by the increasing function of time $\lambda(\tau)+(9GM\bar{A}/2)^{1/3}(t-\tau)^{2/3}$. Though this does not necessarily mean that $\lambda(t)\ (t\geq\tau)$ is itself an increasing function of time, its value at any $t>\tau$ cannot be less than $\lambda(\tau)$. On the other hand, as a length characterizing the system size along a given direction, $\lambda(t)$must decrease dramatically with the system size in the process of the system tending to a close triple approach, including of course the final ones shortly before all escapes seen in numerical experiments on systems with negative total energy and responsible for them (e.g., Anosova 1986,1990; Aarseth et al. 1994). This implies that the time needed for the escape to be initiated has to be quite short. In the last section, we argue that this should also be the case for other analytical tests of escape.

We now go back to the case of our immediate interest, i.e., applying M90 to general spatial 3-body systems with motion state available only on a plane. To our knowledge, M90 is the only analytical test that can be applied to this case, where it is equivalent to M903P. But we do not know if M903P is efficient enough to be really useful. Therefore, in this section we discuss the relevance of M903P with the help of 3-body experiments. The first goal is to verify that M903P detect only escapes that are about to happen and to show some characteristics of the systems with such an escape. The second goal is to show the efficiency of M903P.

 \begin{figure}
\par\includegraphics[width=9cm,clip]{09392fg2.eps}
\end{figure} Figure 2:

Normalized initial configurations of randomly constructed 3-body systems in both S1 and S2. The configuration is normalized by putting the bodies m1 and m2 at (-0.5,0) and (0.5,0), respectively, and select the position of the third body m3 such that $r_{12}\geq r_{13}\geq r_{23}$ and $y\geq 0$, where rij is the distance between mi and mj and ythe y-coordinate of m3. This means that m3 must lie on or inside the boundary curve (x+0.5)2+y2=1 due to $1=r_{12}\geq r_{13}$, on or to the right of the boundary line x=0 due to $r_{13}\geq r_{23}$, and, on or above the boundary line y=0 due to $y\geq 0$. Note that the obvious overdensity of the chosen systems in S1 near the point (0.5,0) reveals the fact that M903P is significantly more efficient in detecting escape systems with hierarchical configuration than those with non-hierarchical configuration.

Open with DEXTER

4.1 Sampling systems and Aarseth's 3-body code

To achieve both our goals, we first need to construct sampling initial 3-body systems in a unit system that allows us to quantitatively compare them. All of the systems mentioned in this section are randomly constructed using the process and the unit system described below (Agekian & Anosova 1967,1968; Anosova 1986,1990; Anosova & Orlov 1992).

To generate the three masses, we first generate three random numbers $(\tilde{m}_1,\tilde{m}_2,\tilde{m}_3)$, which are uniformly and independently distributed in (0,1). And we take $m_i=\tilde{m}_i/\bar{m}$ as the three masses, where $\bar{m}=(\tilde{m}_1\tilde{m}_2+\tilde{m}_2\tilde{m}_3+\tilde{m}_3\tilde{m}_1)/(\tilde{m}_1+\tilde{m}_2+\tilde{m}_3)$. Obviously this means that the chosen unit for mass is

\begin{displaymath}u_M=\frac{M^*}{M}\equiv\frac{m_1m_2+m_2m_3+m_3m_1}{m_1+m_2+m_3}\cdot
\end{displaymath}

In generating initial 3-body configurations, we always start with the ones uniformly normalized by taking their largest mutual distance as the unit for length. To generate such a normalized configuration, we consider a plane with two axes, $\vec{x}$ and $\vec{y}$, forming a cartesian rectangular coordinate system. As shown in Fig. 2, we put the body with mass m1(resp. m2) at (-0.5,0) (resp. (0.5,0)). Since all masses have been generated independently, we can assume that the mutual distances satisfy $r_{12}\geq r_{13}\geq r_{23}$ without any loss of generality. Under this assumption, all possible 3-body (normalized) configurations can be generated with (x3,y3), the position of m3, uniformly distributed in the area enclosed by the two positive half axes and the unit circle centered at (-0.5,0), shown in Fig. 2 as the area enclosed by the heavy solid boundaries. Then we change the unit for length. To obtain the target unit, we introduce k0=T0/(-U0)>0, where T0 and U0 are the initial kinematic and potential energies, respectively. In order for the total energy E<0, k0 must be less than 1. We now generate randomly k0, assumed to be uniformly distributed in (0,1). The unit for length is then fixed as

\begin{displaymath}u_L=\frac{\bar{d}}{2(1-k_0)},
\end{displaymath}

where $\bar{d}=M^*/(m_1m_2/r_{12}+m_2m_3/r_{23}+m_3m_1/r_{31})$ is a weighted harmonic mean of mutual distances. At this stage, we have enough information to work with a barycentric coordinate system.

Finally, we arrive at the selection of the three velocities in the barycentric system. Their directions can be selected by generating, randomly and independently, three unit vectors in the 3-dimensional space, e.g., with ending point uniformly distributed on the unit sphere centered at the starting point. To randomly select the magnitudes of the velocities (v1,v2,v3), we generate three random numbers $(\tilde{v}_1,\tilde{v}_2,\tilde{v}_3)$, which are uniformly and independently distributed in (0,1). And, we calculate the following weighted square mean of them $\bar{v}=\sqrt{(m_1\tilde{v}_1^2+m_2\tilde{v}_2^2+m_3\tilde{v}_3^2)/M}$. A unit for time such that the universal gravitational constant G=1is then

\begin{displaymath}u_T=\frac{\sqrt{k_0}\ \bar{d}}{2(1-k_0)^{3/2}\bar{v}},
\end{displaymath}

which is simply the mean crossing time defined as $\tau_{\rm c}=G\sqrt{M}M^*/(-2E)^{3/2}$ (e.g. Anosova 1990). Together with the unit for length, this gives the unit for velocity $u_V=u_L/u_T=\sqrt{(1-k_0)/k_0}~\bar{v}$. And in this unit, $(v_1,v_2,v_3)=(\tilde{v}_1/u_V,\tilde{v}_2/u_V,\tilde{v}_3/u_V)$.

To numerically integrate a sampling 3-body system, we make use of the 3-body code Triple by Aarseth (ftp:// ftp.ast.cam.ac.uk/pub/sverre/triple/triple.tar.Z, Aarseth & Zare 1974; Aarseth 2003a), incorporating it with the escape test or criterion given by Marchal (formula (38) of Marchal 1974, hereafter M74). M74 is one of the existing analytical criteria predicting an escape by using full system state at a given instant (Birkhoff 1927; Tevzadze 1962; Standish 1971; Griffith & North 1974; Yoshida 1972; Marchal 1974; Kuznetsova & Orlov 1983; Marchal et al. 1984a,b; for some related analytical estimations, see also Laskar & Marchal 1984; and for a review on the practical usage of these criteria, see Anosova 1986). According to Anosova (1986), there is no essential difference between these criteria when used with 3-body simulation, provided that they are applied, for example, at each integration step. Of course, a system with an escape is identified at different evolution stages of the system by different criteria. Therefore, the instant ( $t_{\rm M74}$) at which an escape is detected by M74 is not a suitable reference time instant for the escape. In this context, it should be mentioned that the starting and ending time instants of an escape process are actually not easy to be defined such that they are mathematically precise while, at the same time, useful practically. Here, we simply take as a single reference time instant for escape, $t_{\rm e}$, the instant at which $D_{\rm J}/a_{\rm B}=20$, where $a_{\rm B}$ is the instantaneous semi-major axis of the binary and $D_{\rm J}$ the Jacobian distance measured from the mass center of the binary to the escaper. In accord with this, the time duration $T_{\rm e}=t_{\rm e}-t_0$, where t0 is the initial instant, will be referred to as the escape time, or the time needed for an escape to happen. A randomly constructed initial system can have negative $T_{\rm e}$. In this case, the initial system will be referred to as a post escape system.

 \begin{figure}
\par\includegraphics[width=9cm,clip]{09392fg3.eps}
\end{figure} Figure 3:

Histogram of $T_{\rm e}$ for the initial systems in S1.

Open with DEXTER

We use two sets (S1 and S2) of initial triple systems with negative total energy. S1 contains 1000 systems with an escape detected by M903P. To obtain a system in S1, we proceed as follows. We generate an initial system and apply M903P using the kinematic data on a randomly selected coordinate plane. This process is repeated until an escape is detected; that is, a system in S1is obtained. The systems in S2 are obtained as follows. We first construct a set of 1000 initial systems. Then, by numerical integration and M74, we select those systems with an escape detected before $T_{\rm L}=100\tau_{\rm c}$, where $T_{\rm L}$ can be thought of as the mean lifetime of unbounded triple systems with negative total energy (Anosova 1990). This gives us S2, which contains 604 systems with an escape.

4.2 Characteristics of initial systems with an escape detected by M903P

By performing numerical experiments, we first check and verify that all systems in S1 are indeed systems with an escape. This implies that our computer implementation of M903P gives no false detection of escape.

 \begin{figure}
\par\includegraphics[width=9cm,clip]{09392fg4.eps}
\end{figure} Figure 4:

Evolution of the ratio $D_{\rm J}/a_{\rm B}$ for the initial systems in S1 with $T_{\rm e}>5\tau _{\rm c}$.

Open with DEXTER

To verify that M903P detect only escapes that are about to happen, we draw the histogram in Fig. 3 of $T_{\rm e}$ for the systems in S1. In this figure, the initially post escape systems ($28.5\%$ of the 1000 systems) are included in the leftmost bin, and the percentages there are the accumulated ones of the initial systems with $T_{\rm e}<i\tau_{\rm c}\ (i=1,...,6)$. From Fig. 3, we see that an escape detected by M903P generally has small $T_{\rm e}$, as compared with $T_{\rm L}$. In Fig 4, the time evolution of the ratio $D_{\rm J}/a_{\rm B}$ is shown for each one of the 6 systems with $T_{\rm e}>5\tau _{\rm c}$. As can be inferred from this figure, there will be no close triple approach after t0.

 \begin{figure}
\par\includegraphics[width=9cm,clip]{09392fg5.eps}
\end{figure} Figure 5:

Histograms of k0 for the initial systems in S1 and S2.

Open with DEXTER

 \begin{figure}
\par\includegraphics[width=9cm,clip]{09392fg6.eps}
\end{figure} Figure 6:

Respective distributions of systems in S1 and S2 on the mass ratio plane.

Open with DEXTER

The above discussion implies the importance of the knowledge about the initial systems in which an escape can be detected by M903P. Without this kind of knowledge, the results obtained by applying M903P to a set of systems could be statistically interpreted with significant bias. Therefore, we show some characteristics of the initial systems in S1 relative to those in S2. For systems in both S1 and S2, the normalized configurations, the histograms of k0, and the relative masses are shown in Figs. 2, 5, and 6, respectively.

From Fig. 2, we observe that many systems in S2 are in a non or weakly hierarchical configuration. This is qualitatively the same as the results presented in Fig. 1 of Anosova & Orlov (1992). However, the systems in S1 generally have a significantly hierarchical configuration, as is manifested by the obvious overdensity of the systems in S1 near the point (0.5,0). This reveals the fact that M903P is significantly more efficient in detecting escape systems with hierarchical configuration than those with non-hierarchical configuration. We know that, as an easy consequence of the energy integral, such a configuration is inevitable for any unbounded 3-body system with negative total energy, but, as is well known, it is also typical among systems with very long lifetime $T_{\rm L} \gtrsim 1000\tau_{\rm c}$ (e.g., Anosova & Orlov 1992). This is of practical interest because hierarchical configuration is typical among real triple stars (see the next section for further discussion).

From Fig. 5, we observe that k0 of systems in S2 is distributed rather uniformly over (0,1), though a shallow dip around k0=0.5 is observable (the appearance of this dip is understandable since a system with k0=0.5 is initially in an instant virial equilibrium). The situation is very different for S1, and it is evident that M903P is more efficient in detecting escapes in systems with larger k0, that is, systems in a more significant expansion phase.

From Fig. 6, we observe that the distribution of systems over the shown mass ratio plane has a similar tendency for S1 and S2. Concretely speaking, the number of systems increases with the middle mass, and fewer systems distributed in the lower left corner (one mass dominant over the other two) than the upper right corner (two compatible masses dominant over the remaining one). This is expected, since systems with one dominant mass can usually be divided into two weakly perturbed 2-body systems, and the perturbation becomes more and more significant if one or both lower masses become higher and higher. On the other hand, however, the tendency is obviously stronger for S1 than for S2, which means that M903P is most efficient in detecting escapes in systems with two compatible dominant masses.

4.3 Efficiency of M903P

 \begin{figure}
\par\includegraphics[width=9cm,clip]{09392fg7.eps}
\end{figure} Figure 7:

Percentage of M903P detected escapes as a function of the time elapse from $t_{\rm e}$.

Open with DEXTER

To achieve our second goal, we apply M903P to systems in S2, using system states at 17 different time instants $t_{\rm i}=t_{\rm e}+i\tau_{\rm c}\
(i=-8,...,0,...8)$. We denote the set of system states at ti by S2i. To obtain these 17 sets, we first integrate the initial systems in S2 to $t_{\rm e}$, respectively. This gives S20. By further backward or forward integration, we obtain the 16 other sets.

Figure 7 shows the time dependence of the percentage of the escapes detected by M903P. From this figure, we see again that the percentage is very low if t is several $\tau_{\rm c}$before $t_{\rm e}$. But it increases quickly to much higher values, e.g., it reaches $\sim $$30\%$ at $t=t_{\rm e}-\tau_{\rm c}$ and $\sim $$65\%$ at $t=t_{\rm e}+8\tau_{\rm c}$. Recalling that $t_{\rm e}$ is the time instant at which $D_{\rm J}/a_{\rm B}=20$, we conclude from the above results that M903P is relevant in detecting escapes after the escaper is already rather far away from the remaining binary. This makes M903P a possibly useful tool for identifying real triple stars with an escape, for the reason that many triple stars have a significantly hierarchical configuration.

5 Applications to Hipparcos linear triple systems

It would be reasonable to believe that most frequently observed objects or structures are practically stable or have very long lifetimes. But one should be cautious with the particular case of real triple stars. The rate of occurrence of a hierarchical configuration among these stars is dramatically higher than randomly constructed triple systems. On the one hand, this can be explained as the result of the fact that, dynamically speaking, systems with non-hierarchical configurations are generically unstable. On the other hand, however, it does not necessarily mean that the observed hierarchical triple stars are generically stable. Indeed, three additional facts have to be taken into consideration. One is disclosed by extensive simulations on triple systems with negative total energy; that is, $\sim $$95\%$ of such systems are systems with an escape (e.g. Anosova 1990). The second comes from numerical experiments on binary-binary encounters (Mikkola 1983) and dynamical evolution of star clusters (e.g. Aarseth 2003b,2004); that is, there are triples formed via binary-binary encounters, which can be destroyed by the slingshot mechanism (Saslaw et al. 1971). The remaining one is simply that, if a triple system with negative total energy is a system with an escape, then a hierarchical configuration of the system is inevitable, and the hierarchical nature of the configuration becomes persistent shortly after the final close triple approach. These facts imply the possibility that non-negligible number of hierarchical triple stars are or will soon be in their respective disintegration phase. As can be inferred from the results presented in Figs. 2, 3, and 7, M903P should be a useful tool in identifying these triple stars. And when the initial data needed for numerically tracking the dynamical evolution of a tripe star is unavailable, it is probably an irreplaceable tool.

As explained in the introduction, these initial data can only be calculated from a full orbital solution of the considered triple star, usually obtained by fitting observations to the double two-body model in terms of Jacobian vectors. But there are many triple stars for which such a solution is not available. In particular, this happens when the period, if it exists, of the ``binary - third body'' system is much longer than the time duration spanned by available observations. Indeed, even with observational data as precise as those from the coming five-year long Gaia mission, orbital determination can only be used to obtain a reliable orbital solution for a system with orbital period less than or compatible to the time duration spanned by observations (Lattanzi et al. 2004). Actually, orbital coverage of observations is one of the major factors affecting the grade of an orbit solution, and if the observed arc is much shorter than the full elliptical orbit, then the resulting solution must be of grade ``Indeterminate'' meaning that ``the orbital elements may not even be approximately correct'' (Hartkopf & Mason, Sixth Catalog of Orbits of Visual Binary Stars, http://ad.usno.navy.mil/wds/orb6/orb6text.html). This excludes the possibility of using 3-body integration to detect an escape in a real triple star if the period of its ``binary - third body'' system is much longer than the time duration spanned by all available observations. This kind of duration should be no more than $\sim $100 years for almost every real triple star. A natural question is then ``is M903P really helpful or efficient enough to detect an escape in such kind of a real triple star?'' In this section, we show that the answer is affirmative.

 \begin{figure}
\par\includegraphics[width=9cm,clip]{09392fg8.eps}
\end{figure} Figure 8:

a) Motion state of the Hipparcos linear triple system HIP 2173(2175) on the tangent plane at the epoch J1991.25(TT). The positions of the component stars are marked by stars, to each of which the attached directed line segment represents the velocity of the star relative to the system mass center Q. The two intersect line segments parallel to the axes +RA and +Dec respectively indicate the used scales for the shown velocity components in the respective directions. The printed value of velocity applies to both directions. b) Graph of the piecewise function $q(\hat{\vec{x}})$. The circles indicate the boundaries between different pieces of q, and the dash-dotted horizontal line indicates the value of the constant $\sqrt{2GM\bar{A}}$.

Open with DEXTER

 \begin{figure}
\includegraphics[width=9cm,clip]{09392fg9.eps}
\end{figure} Figure 9:

a) Motion state of the Hipparcos linear triple system HIP 14193 on the tangent plane at the epoch J1991.25(TT). The radial line segments starting from the mass center Q are half of the boundary axes separating the different pieces of the function $q(\hat{\vec{x}})$ (the other half can be obtained by reversing each of the shown axes). In this plot, the scales of both axes +RA and +Dec are exactly the same, so the shape of the system configuration, as well as the angular distances between the above-mentioned boundary axes, is real. b) Graph of the piecewise function $q(\hat{\vec{x}})$. c) Locally enlargement of the panel b). All of the unmentioned parts of the legends for the present figure can be found in the caption of Fig. 8.

Open with DEXTER

As illustrating examples, we apply M903P to the Hipparcos linear triple stars found in the Double and Multiple Systems Annex (ESA 1997). In this Annex, the component solutions are categorized into 3 types: F, I, and L. An F-type system has only data describing the integral motion of the system as a whole, and an I-type system is thought to be optical. This left only type L systems for us to consider here.

An L-type system was identified as physical mainly by pre-Hipparcos long-term observations, but there were no reliable parameters fully describing its kinematics. According to the Multiple Star Catalog (MSC, Tokovinin 1997), which is now an extensive and dynamic data source available on the world wide web (http://www.ctio.noao.edu/~atokovin/stars/intro.html), this situation remain unchanged until the present day. This is understandable, since the expected period of ``binary - third body'' system ranges from $\sim $500 years to more than 10 000 years (MSC).

There are 35 L-type triple systems. For each system, the data describing the motion states of all component stars on the tangent plane at J1991.25 are provided with their respective standard deviations. As expected, the individual parallaxes of the component stars are not given. But the parallax of each system as a whole is given, which allows us to transform angular distances to linear ones on the tangent plane. Unfortunately, however, only 22 out of these L-type systems have all their respective component masses found in MSC. And, in the remaining 13 systems, two have all their component visual magnitudes individually determined by Hipparcos consortia, which allows us to use the empirical stellar mass-luminosity relation (Henry & McCarthy 1994) to estimate the component masses. MSC does not provide standard errors for the estimated masses, which reflects various difficulties and complexities in estimating stellar masses. Following Orlov & Zhuchkov (2005), we here assume that the relative error of mass estimation is $10\%$.

We then have 24 real triple systems having data exactly as required by M903P. These systems are listed in Table 1. For the purpose of illustrating, we first concentrate on the mathematically expected systems (that is, taking the mathematically expected values of all of the above-mentioned data, ignoring their errors). We show in Figs. 8 and 9, respectively, two concrete examples of the system geometry and the piecewise function $q(\hat{\vec{x}})$. Because of the hierarchical feature of the system configurations, not all of the typical characteristics of the function $q(\hat{\vec{x}})$ are discernible from the (b) panels of these two figures. By locally enlarging the graph in the (b) panel of Fig. 9, the (c) panel of the same figure clearly shows those characteristics and their correspondences to the geometry depicted in the panel (a).

According to the (b) panels of Figs. 8 and 9, in which the horizontal dash-dotted lines indicate the mass-dependent quantity $\sqrt{2GM\bar{A}}$, M903P tells us that the triple system HIP 2173(2175) will eventually be disintegrated, while no conclusion can be made for HIP 14193. Figure 10 schematically shows the respective graphs of $q(\hat{\vec{x}})$ for all 24 mathematically expected systems. By this figure, it is easy to identify those systems in which an escape will definitely happen.

We check that our code gives the same results as those given by Fig. 10. These results are presented as +1 or +0 in the last column of Table 1. In the first case, the mathematically expected system will be disintegrated, and in the second case, no conclusion can be drawn from M903P. In Table 1, the disintegration probabilities of the considered systems are also shown. For obtaining these statistical results, we generate for each system a set of 1000 sampling systems, assuming that the error of each used datum follows the standard normal distribution. It is important to be aware that the shown values of the probability are actually their respective lower bounds. This is made explicit by using the sign $\geq$.

 \begin{figure}
\par\includegraphics[width=16.5cm,clip]{09392f10.eps}
\end{figure} Figure 10:

Graphs of piecewise functions $q(\hat{\vec{x}})$ for the Hipparcos L-type triple stars. The legends are the same as panel b) of Fig. 8. The numbers in the panels are the Hipparcos ID numbers, and the tag ``E'' indicates that an escape event is unavoidable and the system will eventually be disintegrated.

Open with DEXTER

Table 1:   A sample of triple stars and their disintegration probability.

6 Discussion

In the 3-body problem, two kinds of initial phase points, leading to bounded and unbounded orbits, are entangled in certain parts with negative total energy of the phase space. Therefore, the ``escape region'' composed of only initial conditions leading to an escape, as defined by an analytical criterion (a certain inequality), is expected to be well away from these parts of the phase space. For example, the analytical test by Marchal et al. (MYS, Marchal et al. 1984a), which is believed to be the most efficient one known so far (Anosova 1986; Georgakarakos 2008), requires the precondition that there is an isolated ``third body''. In this test, the distance used to characterize the system size is $D_{\rm J}$, i.e. the length of the position vector of this third body relative to the mass center of the other two. This distance will have at most one minimum after an escape is detected (Marchal et al. 1984a), which implies that the system should have already gone through the complex process of the final close triple approach, and the third body is not far from being ejected without return. In this respect, as we argued in the first paragraph of Sect. 4, M90 should be no exception.

The main advantage of M90 over other existing analytical tests is that it can be applied to a general spatial N-body system with only partial information on the system state. In particular, M90 allows us to design the algorithm M903P and apply it successfully to some real triple stars. M903P is quick because only 18 (at most) function evaluations are required (for example, it takes only several seconds of CPU time on a usual personal computer to test all our 24 000 sampling systems of Hipparcos linear triple stars). Therefore, it should be a useful tool to detect triple stars in disintegration phase among linear ones from, for example, the coming Gaia mission (http://www.esa.int/science/gaia). These linear systems should be characterized by the same features as the Hipparcos ones we considered in this work. First, the observational data can be used to determine the linear motion states of the component stars on the plane of sky. And, to apply M903P, it is also required that the component star masses can be reasonably estimated, e.g., by using a stellar mass-luminosity relation. Here, it should be emphasized again that M903P is efficient only for triple stars in hierachical configuration; in other words, the non-hierachical category of triple stars is basically irrelevant for the use of M903P. As for applying M90 to multiple stars with more than 3 components, one has to first extend M903P to systems with the same number of bodies.

Theoretically speaking, M903P and its extension could be of interest, but mainly in the case of planar systems. In this case, applying them is equivalent to applying M90. For spatial N-body systems, however, it would be necessary to develop an algorithm implementing M90 with all axes in 3-dimensional physical space. Such an algorithm is certainly useful in the study of final fate of general spatial N-body systems (N>3) with negative total energy. If N=3, there are other existing criteria for escape, which are often used with numerical integration. In fact, analytical criteria and numerical integrations are both important in this kind of study. On the one hand, an analytical criterion is indispensable for one to know for sure whether there will be an escape. On the other hand, numerical integration can provide system states at different evolution stages to which an analytical test can be applied. This is crucial because a single implementation of an analytical test is not efficient enough in the study. According to Anosova (1986), all analytical tests of escape known to her at that time (but already include all that mentioned in Georgakarakos 2008) are practically of no significant difference when used with numerical integration. A natural question is then whether M90, when well programmed, can make any difference?

A hint to this question can be obtained by comparing the efficiencies of M90 and M74 in the case of the planar 3-body problem, where M90 is equivalent to M903P. A simple bidirectional comparison is made as explained below. Among a set of 104randomly selected initial planar 3-body systems with an escape detected by M74 (resp. M90), there are 3741 (resp. 4467) systems that can be identified as system with an escape by M90 (resp. M74). It can be inferred from this comparison that M74 is more efficient than M90, but this is not in the sense that an escape detected by M90 can also be detected by M74. In fact, it is reasonable for M74 and M90 to be complementary to each other, since, as we already mentioned in the introduction, M90 is constructed using an idea different from the one used to derive other existing tests of escape. This mutual complementarity would be interesting because it implies that alternatively applying M74 and M90 in the process of numerical integration could be helpful for saving computing time. However, this may not be very interesting in practice, because the final triple approach usually produces a very energetic escaping body, and the time instants at which such an escape can be detected by different criteria, respectively, are not expected to be very different.

Similar comparison between M90 and an energetic criterion is also made. The energetic criterion is expressed, following Aarseth, as $E_{\rm e}=E-E_{\rm B}>0$, where E and $E_{\rm B}$ are, respectively, the energies of the whole system and of the ``binary'' formed by the two bodies with the least mutual distance, and $E_{\rm e}$ can be understood as the orbital energy of the distant body orbiting around the binary. This criterion is not analytical, and so, it can lead to a false detection of escape, especially when the system configuration is not sufficiently hierarchical. Therefore, a precondition such as $D_{\rm J}/a_{\rm B}>r$ where r is a free parameter, has to be checked before applying the energetic criterion. The comparison results show that M90 is more efficient than the energetic criterion with $r\geqslant
10.0$, and the opposite is true if $r\lesssim 9.5$. Therefore, the main advantage of M90 over the energetic criterion is mainly that M90 is analytical and so its results is more determinate.

As already stated, M90 applies to systems with more than three bodies. Therefore, it is interesting to realize M90 by an algorithm making use of a full motion state of a general spatial 4 or more body system. Improving the efficiency of M90 is also a worthy goal (Marchal, private communication). In his 1990 book, Marchal presented a simple way to improve the test for an equal mass system to, in a certain sense, its most efficient version. The simple idea behind this is that one can use a different characteristic distance, replacing $\lambda$, to remove the loss in the estimation of d$^2\lambda$/dt2 for different s satisfying $\lambda=\lambda_s$. In the context of the present work, it should be pointed out that, for the 3-body case, because of the equivalence between $(o,\vec{x},1)$ and $(-o,\vec{\tilde{x}},2)$ (see the second section), the original test is already the most efficient one in this sense.

Acknowledgements
We are grateful to Christian Marchal suggesting the present research topic. And we are deeply indebted to an anonymous referee and Dr. Beust for their instructive and constructive comments, which significantly improved the overall presentation. Many thanks to Fang Xia and Shulin Ren for various helpful discussions on the mass and kinematic data of real triple stars. Alain Chenciner is thanked for instructive discussions of the initial conditions leading respectively to bounded and unbounded motions of a 3-body system. Alain Albouy, Christian Marchal again, and other colleagues at IMCCE are thanked for discussions after the second author had reported this work at the IMCCE. Thanks go to Sverre Aarseth for his very informative comments on numerical study of N-body systems. This work is supported by National Basic Research Program of China (973 Program) 2007CB814800 and NSFC 10473025 & 10833001.

Appendix A: An algorithm implementing M903P

According to Sect. 3, an algorithm implementing M90 with all axes on a plane, i.e. M903P, is basically the one finding the maximum value of the piecewise function $q(\alpha)$. Therefore, this appendix will first study this function. The notations introduced in Sect. 3 are used without further explanation.

The piecewise function $q(\alpha)$ depends on

\begin{displaymath}D=\{(m_a,\vec{r}_a,\vec{v}_a);(m_b,\vec{r}_b,\vec{v}_b);(m_c,\vec{r}_c,\vec{v}_c)\}
\end{displaymath}

where, for example, ma and $(\vec{r}_a,\vec{v}_a)$ are respectively the mass and the projected motion state on the coordinate plane P of the body a moving in the 3-dimensional space. And D is supposed to be all our known data. Of course, the concrete form of $q(\alpha)$ also depends on how we parameterize $X=\{\hat{\vec{x}}: \hat{\vec{x}} \in P$}, that is, on the chosen reference axes $\vec{\xi} \in P$ and $\vec{\zeta}$ perpendicular to P. But this is not an intrinsic aspect of the problem, so we not limited to the choice shown in Fig. 1. Indeed, we are free to use different reference axes for different pieces of q, and, as we see later, all pieces can be treated in the same way in conveniently chosen piece-dependent coordinate systems. Therefore, we first proceed geometrically by considering the mapping $q : X \mapsto
\mathcal{R}$ without referring to any particular coordinate system. The instant goal is to categorize, for a given D, the axes in Xsuch that each category corresponds to a piece of $q=q(\hat{\vec{x}})$, $\hat{\vec{x}} \in X$.

The critical axes separating two adjacent categories of X are shown in Fig. 1, where the projection point of a body on the coordinate plane is simply referred to, for brievity, as the body itself. The shown directions are

\begin{displaymath}
\begin{array}{l}
\vec{\vec{\zeta}}=\vec{\vec{\zeta}}_{abc}=...
...c=R_{\frac{\pi}{2}}(\vec{r}_c-\vec{r}_a-\vec{r}_b).
\end{array}\end{displaymath} (A.1)

Here, $\vec{\zeta}_{abc}$ specifies the opposite direction the reader looks at the geometrical object, and technically speaking, it defines an orientation of the plane and makes the rotation $R_{\alpha}$ a well-defined operator. The geometry depicted in Fig. 1 remains qualitatively the same for any given D, except for the degenerate ones that involve three collinear bodies and unoriented $\vec{\zeta}_{abc}=0$. We consider this practically unimportant case later for theoretically completeness, but, until then, we exclude it from our consideration and refer to the remaining D as noncollinear D. This implies, in particular, that each formula in the last three rows of (A.1) defines a non-zero vector, so it can be used to define two axes in X, e.g. $\hat{\vec{x}}=\hat{\vec{i}}_a$ and $\hat{\vec{x}}=-\hat{\vec{i}}_a$.

 \begin{figure}
\par\includegraphics[width=9cm,clip]{09392fA1.eps}
\end{figure} Figure A.1:

Open angular regions $I_1,I_2,I_1\hspace {0.2mm}'$, and $I_2\hspace {0.2mm}'$ collecting x-axes on which a given body is projected strictly as the middle one.

Open with DEXTER

To discuss the geometrical meaning of each shown axis and the correspondence of the angular regions bounded by any two adjacent axes to the pieces of the function q defined on X, we start with the axis $\vec{i}_a$, which points from Q perpendicularly to the line passing through mb and mc. Obviously, these two bodies have the same x-coordinate if we set $\hat{\vec{x}}=\hat{\vec{i}}_a$. This implies that $\vec{i}_a$ is a critical axis with respect to the ordering of the bodies when they are projected on $\hat{\vec{x}} \in X$. In fact, the order changes from xa<xc<xb to xa<xb<xc as $\hat{\vec{x}}$ rotates anticlockwise across $\vec{i}_a$. This order will remain unchanged until this rotating $\hat{\vec{x}}$ reaches $-\vec{i}_c$, on which we have xa=xb.

Under permutation $\{abc\} \mapsto \{bca\}~(\{cab\})$, the above description of the angular region from $\vec{i}_a$ to $-\vec{i}_c$also applies to the one from $\vec{i}_b$ to $-\vec{i}_a$ (from $\vec{i}_c$ to $-\vec{i}_b$, respectively). In the way Marchal labels the bodies (see the second section), the aforementioned three angular regions can be uniformly illustrated as in Fig. A.1, where we set an order-dependent $\hat{\vec{\xi}}=\hat{\vec{i}}_1$. To understand the role of the direction $\vec{j}_2$ in Fig. A.1, we calculate the difference between $\lambda_1$ and $\lambda_2$ for an x-axis inside the angular region going from $\vec{i}_1$ anticlockwise to $-\vec{i}_3$,

\begin{displaymath}
\left\{
\begin{array}{l}
\lambda_1 = -\frac{Mx_1}{m_2+m_3} ...
...ha_{\vec{j}_2}-\alpha)}{(m_1+m_2)(m_2+m_3)}
\end{array}\right.
\end{displaymath} (A.2)

where $\langle \cdot, \cdot \rangle$ denotes the usual scalar product, and $\alpha_{\vec{v}}$ the angle sweeping positively from $\hat{\vec{\xi}}=\hat{\vec{i}}_1$ to a vector or an axis $\vec{v}$, e.g. $\alpha_{\hat{\vec{x}}}=\alpha$, as shown in Fig. A.1. From this calculation, we know that $\lambda_1=\lambda_2$ if and only if $\hat{\vec{x}}=\hat{\vec{j}}_2$, and that $\lambda=\lambda_1
> \lambda_2$ when $\hat{\vec{x}} \in I_1$ (see Fig. A.1).

Based on the above arguments, we have for $\hat{\vec{x}} \in
\overline{I}_1$ (closure of I1), i.e. $\alpha \in
[0,\alpha_{\vec{j}_2}]$,

\begin{displaymath}
\begin{array}{ll}
q & = q_1(\alpha) = \sqrt{\lambda_1} \fra...
...c{r}_1\vert\vert\vec{v}_1\vert^2
}{(m_2+m_3)^3}}},
\end{array}\end{displaymath} (A.3)

where $\phi_1$ is the position angle of $\vec{v}_1$ measured from $\vec{r}_1$, $\beta=\alpha_{\vec{r}_1}-\alpha$, the position angle of $\hat{\vec{r}}_1$ measured from $\hat{\vec{x}}$, and $p(\beta)=\sqrt{-\cos\beta}\cos(\beta+\phi_1)$. It should be noted that p is a well-defined real quantity. Actually, we have $\cos\beta<0$ because

\begin{displaymath}\beta=\alpha_{\vec{r}_1}-\alpha \in
[\alpha_{\vec{r}_1}-\alph...
...\vec{r}_1}]
\subset \left(\frac{\pi}{2},\frac{3}{2}\pi\right).
\end{displaymath}

It is then obvious that, on $\bar{I}_1 \ni \hat{x}$, one piece of the function q is defined, associated with a given ordering of bodies and the condition $\lambda=\lambda_1\geq\lambda_2$. When restricted to this piece, our goal of finding maximum q for a given D is reduced to minimizing the continuous function $p(\beta)$ over $[\alpha_{\vec{r}_1}-\alpha_{\vec{j}_2},\alpha_{\vec{r}_1}]$. The minimum is attained at one of the following values of $\beta$: the two extreme values $\alpha_{\vec{r}_1}-\alpha_{\vec{j}_2}$ ( corresponding to $\hat{\vec{x}}=\hat{\vec{j}}_2$) and $\alpha_{\vec{r}_1}$ (corresponding to $\hat{\vec{x}}=\hat{\vec{i_1}}$), and potential local minimizer(s) $\beta_c\in
(\alpha_{\vec{r}_1}-\alpha_{\vec{j}_2},\alpha_{\vec{r}_1})$satisfying $\frac{{\rm d}p}{{\rm d}\beta}(\beta_c)=0$ and $\frac{{\rm d}p^2}{{\rm d}\beta^2}(\beta_c) \geq 0$; that is,

\begin{displaymath}
\left\{
\begin{array}{l}
\sin \phi_1+3 \sin (2 \beta_c +\phi_1)=0\\
\cos (2\beta_c +\phi_1) \geq 0.
\end{array} \right.
\end{displaymath} (A.4)

The equation in (A.4) is solved for the angle $\gamma=2
\beta_c +\phi_1 \in (\phi_1+\pi,\phi_1+3\pi)$. This gives only two roots. One is congruent to an angle (modulo $2\pi$) in $[-\gamma_b,\gamma_b]$, where $\gamma_b={\rm Arcsin}(1/3)$, and the other in $[\pi-\gamma_b,\pi+\gamma_b]$. Therefore, the inequality in (A.4) is actually strict and tells us that there is one and only one local minimum of $p(\beta)$ attained at $\beta_c$ corresponding to the first root. But it depends on $\phi_1
\in [0,2\pi)$, as well as on the considered piece, whether this $\beta_c\in
(\alpha_{\vec{r}_1}-\alpha_{\vec{j}_2},\alpha_{\vec{r}_1})$. And similarly, because $\phi_1$ is independent of the above-mentioned extreme values of $\beta$, we do not have enough information to predict exactly at which $\beta \in
\{\alpha_{\vec{r}_1}-\alpha_{\vec{j}_2},\alpha_{\vec{r}_1},\beta_c\}$, $p(\beta)$ reaches its minimum value. Nonetheless, it is already clear from the above arguments that $q_1(\alpha)$ has an attainable maximum value $\sup(q_1(0),q_1(\alpha_{\vec{j}_2}),q_1(\alpha_{\vec{i}_1}-\beta_c))$, or $\sup(q_1(0),q_1(\alpha_{\vec{j}_2}))$ if the obtained $\beta_c
~\overline{\in}~
[\alpha_{\vec{r}_1}-\alpha_{\vec{j}_2},\alpha_{\vec{r}_1}]$. This makes it easy to design an algorithm that is valid for any non-collinear D, of finding the maximum of q1.

To complete the discussion for the axes on which the projections of the bodies are in a given order, it is necessary for us to also consider $\hat{\vec{x}} \in \overline{I}_2$ (see Fig. A.1), i.e. $\alpha \in [\alpha_{\vec{j}_2},\alpha_{-\vec{i}_3}]$. Here, we have to again include the axis $\hat{\vec{x}}=\hat{\vec{j}}_2$, which was considered as a boundary axis of $\overline{I}_1$. This is because, though we have $\lambda_1=\lambda_2$ for this axis, q=q1assuming $\lambda=\lambda_1$ is generally different from $q=q_2=\sqrt{\lambda_2} \frac{d\lambda_2}{dt}$ assuming $\lambda=\lambda_2$. To understand this, one simply notes that $\hat{\vec{j}}_2$ only depends on the position of the 3-body system in the configuration space, while q also depends on the velocities.

We now consider $\hat{\vec{x}} \in \overline{I}_2$. For this, it is helpful to recall the equivalence between $(o,\vec{x},s)$ and $(-o,\vec{\tilde{x}},N-s)$. In terms of the notations used in Fig. A.1, this equivalence can be interpreted as that between $(\hat{\vec{x}},\{x_1,x_2,x_3\},2)$ and $(-\hat{\vec{x}},\{x_1\hspace{0.2mm}',x_2
\hspace{0.2mm}',x_3\hspace{0.2mm}'\},1)$, which means that we can forget $\hat{\vec{x}} \in \overline{I}_2$ as long as $\hat{\vec{x}}
\in \overline{I}_1\hspace{0.2mm}'$ is also considered. Noting that the plane orientation is reversed if one goes from the non-primed to the primed situations, he or she easily checks that the above discussion about $\overline{I}_1$ applies exactly to $\overline{I}_1\hspace{0.2mm}'$. This gives us all the knowledge we need to design a concise algorithm for finding maximum q and implementing M90 based on a noncollinear D.

The algorithm can be described as follows. We perform a loop over all six permutations of the three bodies; that is, we associate each time the triple label 123 with a different body order abc,bca,cab,acb,cba, or bac. For each permutation, we go through the following steps.

1.
Calculate the necessary directions

\begin{displaymath}
\left\{\begin{array}{l}
\vec{\vec{\zeta}}_{123}=(\vec{r}_2-...
...ac{\pi}{2}}(\vec{r}_2-\vec{r}_3-\vec{r}_1).
\end{array}\right.
\end{displaymath} (A.5)

Since the involved vectorial calculations can be easily realized in any fixed rectangular coordinate system with origin Q, it is not necessary in practice for us to be bothered with the ones we introduced in the above discussion. This notice applies also to similar cases in the subsequent steps.

2.
Solve the equation in (A.4) for the angle $\gamma=2\beta_c +\phi_1$, assuming $\gamma \epsilon
[-\gamma_b,\gamma_b]$; calculate $\beta_c$ (if it is in the required interval $(\alpha_{\vec{r}_1}-\alpha_{\vec{j}_2},\alpha_{\vec{r}_1})$, then the corresponding critical axis $\vec{x}_{c}=R_{-\beta_c}(\vec{r}_1)$ is a local minimizer of p, or equivalently, a local maximizer of q1). The required formulae for this step are

\begin{displaymath}
\left\{\begin{array}{l}
\gamma = {\rm Arcsin}(-\frac{\sin \p...
...i_1}{2} \in
(\frac{\pi}{2},\frac{3\pi}{2}),
\end{array}\right.
\end{displaymath} (A.6)

where k is an integer. It can be checked that $k \in \{1,2\}$ is a necessary condition for $\beta_c \in (\frac{\pi}{2},\frac{3\pi}{2})
\supset (\alpha_{\vec{r}_1}-\alpha_{\vec{j}_2},\alpha_{\vec{r}_1})$. This means that we only need to try these two values of k to decide whether $\beta_c$ is in the required interval.

3.
Calculate q1 according to (A.3) with $\beta=\alpha_{\vec{r}_1}-\alpha_{\vec{j}_2}, \alpha_{\vec{r}_1}$, and $\beta_c$ (if $\beta_c\in
(\alpha_{\vec{r}_1}-\alpha_{\vec{j}_2},\alpha_{\vec{r}_1})$).

4.
Using the maximum q1 obtained at step 3, apply test (4).

If step 4 ends up (4) being satisfied, we stop the program and conclude that the system will eventually be disrupted. Otherwise, we continue with the next permutation until all permutations are tested. $\blacksquare$

The above discussion does not apply to the measure zero class of Dcharacterized by the collinearity of the three bodies. This is because $\hat{\vec{\zeta}}$ in (A.1) or (A.5) is no longer a well-defined direction. To see some intrinsic specialities pertinent to this class, we first concentrate on its subclass, for which there are no superposed bodies on the plane. For D in this subclass, it is the same body that lies strictly in the middle of the projected configuration on almost all x-axes. The exceptional x-axes are the two opposite ones orthogonal to the line $\ell$ containing all three bodies. On these two axes, all three projected bodies are superposed at Q, which allows no conclusion to be drawn from M903P, so both of them are excluded from the following discussion. We can then say that the number of possible orderings of the projected bodies is reduced to 2, or, intrinsically speaking, 1, if we take into account the previously mentioned equivalence. Moreover, the collinearity ensures that all projected configurations, regardless of the chosen x-axis, have the same intrinsic shape. Therefore, there is only one piece ( $\lambda=\lambda_1
> \lambda_2$ or $\lambda=\lambda_2>\lambda_1$) or there are two superposed pieces ( $\lambda=\lambda_1=\lambda_2$) of $q(\hat{\vec{x}})$ associated with each given ordering.

 \begin{figure}
\par\includegraphics[width=9cm,clip]{09392fA2.eps}
\end{figure} Figure A.2:

Coordinate systems and notations in the case of collinear configuration. a) Subcase where there is no overlapped body. b) Subcase where there are a pair of overlapped bodies.

Open with DEXTER

All of the above-mentioned geometrical facts will become clearer if we interpret them using a conveniently fixed coordinate system. To construct such a system, shown in the left panel of Fig. A.2, we first fix the orientation of the plane $\wp$ that contains the known motion states by choosing any one of the two opposite axes perpendicular to $\wp$ as the $\vec{\zeta}$-axis. Here, it should be noted that, this plane, e.g. the one orthogonal to the line-of-sight connecting the observer and Q, can be well-specified even if all three velocities are parallel to $\ell$. Then, denoting the two exterior bodies by ma and mc, we set $\hat{\vec{\eta}}=(\vec{r}_b-\vec{r}_c-\vec{r}_a)/d$ if the resulting $\hat{\vec{\eta}}\neq 0$, where $d=\vert\vec{r}_b-\vec{r}_c-\vec{r}_a\vert$. If it happens that the so-defined $\hat{\vec{\eta}}=0$, we use two different coordinate systems, referred to as non-primed and primed ones. For the non-primed one, we define $\hat{\vec{\eta}}=(\vec{r}_c-\vec{r}_a)/\vert\vec{r}_c-\vec{r}_a\vert$, and for the primed one $\hat{\vec{\eta}}\hspace{0.2mm}'=(\vec{r}_a-\vec{r}_c)/\vert\vec{r}_a-\vec{r}_c\vert$. Finally, we complete the coordinate system(s) by setting $\hat{\vec{\xi}}=\hat{\vec{\eta}} \times \hat{\vec{\zeta}}$ (and $\hat{\vec{\xi}}\hspace{0.3mm}'=\hat{\vec{\eta}}\hspace{0.2mm}'
\times \hat{\vec{\zeta}}$).

By the previously mentioned equivalence between any two opposite x-axes, we only need to consider the x-axes with $\alpha \in
(0,\pi)$ (see Fig. A.2). We first focus on the non-primed notations in Fig. A.2. We have m1=ma,m2=mb, and m3=mcin the strict sense, that is xa<xb<xc. Noting that $\alpha_{\vec{r}_1}=\frac{3\pi}{2}~{\rm and}~\alpha_{\vec{r}_3}=\frac{\pi}{2}$, one easily verifies directly or by adapting (A.2) to the present setting,

\begin{displaymath}
\left\{
\begin{array}{l}
\lambda_1 = \frac{M\vert\vec{r}_1\...
...eq 0\\ =0~~~{\rm if}~d=0.\end{array}\right.
\end{array}\right.
\end{displaymath} (A.7)

From (A.7), it is clear that we either have only one piece of $q(\hat{\vec{x}})$ ( $\lambda=\lambda_1$) if $d\neq 0$, or two superposed pieces ( $\lambda=\lambda_1=\lambda_2$) if d=0.

In the case of $d\neq 0$, we have q=q1. And, (A.3) with the present specific $\beta=\frac{3\pi}{2}-\alpha$ is reduced to

\begin{displaymath}
\begin{array}{lll}
q_1=\sqrt{\sin\alpha }\sin(\alpha-\phi_1...
...vec{r}_1\vert\vert\vec{v}_1\vert^2}{(m_2+m_3)^3}}},
\end{array}\end{displaymath} (A.8)

which can also be checked directly with the present setting. According to (A.8), $\lim_{\alpha\rightarrow 0}
q_1=\lim_{\alpha\rightarrow \pi} q_1=0$. On the one hand, this confirms that no conclusion can be drawn by choosing $\hat{\vec{x}}=\pm\hat{\vec{\eta}}$, because q1>0 is needed for M903P to be effective. On the other hand, it implies that if there is an $\alpha \in
(0,\pi)$ such that q1>0, then a positive maximum of the continuous function $q_1(\alpha)\ (\alpha \in
(0,\pi))$ must be attained at some $\alpha \in
(0,\pi)$. For finding this maximum, one can either follow the previous procedure through the intermediate variable $\beta$ or formulate an equivalent one directly from (A.8).

Now we come to the case of d=0, for which we have to consider both $\lambda=\lambda_1$ and $\lambda=\lambda_2$ in searching for the maximum value of q. Again, the trick is to consider instead $\lambda=\lambda_1$ and $\lambda=\lambda_1\hspace{0.2mm}'$, using respectively the primed and non-primed settings shown in Fig. A.2. The discussion of seeking the maximum of q1 (or q1') is exactly the same as the case where $d\neq 0$.

We end this appendix by pointing out that no further speciality exists with the case of one body being superposed over another one, e.g. $\vec{r}_b=\vec{r}_c$ as illustrated in the right panel of Fig. A.2, where the direction of $\vec{\eta}$-axis is fixed as from the non-superposed body (ma) to the superposed ones (mband mc). Indeed, disregarding the case of all three bodies being superposed at Q to which M903P does not apply, we now have the simplest case. This is because d=0 is now impossible, and $\lambda_1>\lambda_2$ remains true for all x-axes directing upward. And, changing from $\{m_2=m_b,m_3=m_c\}$ to $\{m_2=m_c,m_3=m_b\}$ makes no difference to this discussion for the case where there is no superposed body and $d\neq 0$.

References

Footnotes

... algorithm[*]
A Fortran 90 code implementing this algorithm can be found at http://159.226.71.170/M903P.f90.

All Tables

Table 1:   A sample of triple stars and their disintegration probability.

All Figures

  \begin{figure}
\par\includegraphics[width=9cm,clip]{09392fg1.eps}
\end{figure} Figure 1:

Critical x-axes separating different pieces of the piecewise function $q(\hat{\vec{x}})$. The symbol $\odot $, situated at the barycenter Q, indicates that $\vec{\zeta}$ is perpendicular to the paper and is directed toward the reader.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=9cm,clip]{09392fg2.eps}
\end{figure} Figure 2:

Normalized initial configurations of randomly constructed 3-body systems in both S1 and S2. The configuration is normalized by putting the bodies m1 and m2 at (-0.5,0) and (0.5,0), respectively, and select the position of the third body m3 such that $r_{12}\geq r_{13}\geq r_{23}$ and $y\geq 0$, where rij is the distance between mi and mj and ythe y-coordinate of m3. This means that m3 must lie on or inside the boundary curve (x+0.5)2+y2=1 due to $1=r_{12}\geq r_{13}$, on or to the right of the boundary line x=0 due to $r_{13}\geq r_{23}$, and, on or above the boundary line y=0 due to $y\geq 0$. Note that the obvious overdensity of the chosen systems in S1 near the point (0.5,0) reveals the fact that M903P is significantly more efficient in detecting escape systems with hierarchical configuration than those with non-hierarchical configuration.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=9cm,clip]{09392fg3.eps}
\end{figure} Figure 3:

Histogram of $T_{\rm e}$ for the initial systems in S1.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=9cm,clip]{09392fg4.eps}
\end{figure} Figure 4:

Evolution of the ratio $D_{\rm J}/a_{\rm B}$ for the initial systems in S1 with $T_{\rm e}>5\tau _{\rm c}$.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=9cm,clip]{09392fg5.eps}
\end{figure} Figure 5:

Histograms of k0 for the initial systems in S1 and S2.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=9cm,clip]{09392fg6.eps}
\end{figure} Figure 6:

Respective distributions of systems in S1 and S2 on the mass ratio plane.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=9cm,clip]{09392fg7.eps}
\end{figure} Figure 7:

Percentage of M903P detected escapes as a function of the time elapse from $t_{\rm e}$.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=9cm,clip]{09392fg8.eps}
\end{figure} Figure 8:

a) Motion state of the Hipparcos linear triple system HIP 2173(2175) on the tangent plane at the epoch J1991.25(TT). The positions of the component stars are marked by stars, to each of which the attached directed line segment represents the velocity of the star relative to the system mass center Q. The two intersect line segments parallel to the axes +RA and +Dec respectively indicate the used scales for the shown velocity components in the respective directions. The printed value of velocity applies to both directions. b) Graph of the piecewise function $q(\hat{\vec{x}})$. The circles indicate the boundaries between different pieces of q, and the dash-dotted horizontal line indicates the value of the constant $\sqrt{2GM\bar{A}}$.

Open with DEXTER
In the text

  \begin{figure}
\includegraphics[width=9cm,clip]{09392fg9.eps}
\end{figure} Figure 9:

a) Motion state of the Hipparcos linear triple system HIP 14193 on the tangent plane at the epoch J1991.25(TT). The radial line segments starting from the mass center Q are half of the boundary axes separating the different pieces of the function $q(\hat{\vec{x}})$ (the other half can be obtained by reversing each of the shown axes). In this plot, the scales of both axes +RA and +Dec are exactly the same, so the shape of the system configuration, as well as the angular distances between the above-mentioned boundary axes, is real. b) Graph of the piecewise function $q(\hat{\vec{x}})$. c) Locally enlargement of the panel b). All of the unmentioned parts of the legends for the present figure can be found in the caption of Fig. 8.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=16.5cm,clip]{09392f10.eps}
\end{figure} Figure 10:

Graphs of piecewise functions $q(\hat{\vec{x}})$ for the Hipparcos L-type triple stars. The legends are the same as panel b) of Fig. 8. The numbers in the panels are the Hipparcos ID numbers, and the tag ``E'' indicates that an escape event is unavoidable and the system will eventually be disintegrated.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=9cm,clip]{09392fA1.eps}
\end{figure} Figure A.1:

Open angular regions $I_1,I_2,I_1\hspace {0.2mm}'$, and $I_2\hspace {0.2mm}'$ collecting x-axes on which a given body is projected strictly as the middle one.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=9cm,clip]{09392fA2.eps}
\end{figure} Figure A.2:

Coordinate systems and notations in the case of collinear configuration. a) Subcase where there is no overlapped body. b) Subcase where there are a pair of overlapped bodies.

Open with DEXTER
In the text


Copyright ESO 2009

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.