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

where mi and


where M is the total mass, M* a quantity depending on the masses,


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.
,
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''
and
,
.
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,

where



where


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













With the length
,
for any
given
,
there must be at least one sfor which
.
Therefore,
implies that the system will
eventually disintegrate into at least two permanent subsystems. This
means that
can be used as the aforementioned
characteristic distance describing the escape phenomenon.
By (1), for any given value of ,
the acceleration
of this length can be limited from below as
where




This least upper-bound




Proposition.
(Marchal 1990) If at a given instant
then this inequality is true for all

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,
goes to infinity with
at least as
.
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
![]() |
Figure 1:
Critical x-axes separating different pieces of the piecewise function
|
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
,
with its origin the
barycenter Q and
-plane the plane P.
Obviously, M90 can be applied to the system with kinematic data on
any given axis
.
All such x-axes constitute a
one-parameter family
},
where
denotes the unit vector associated with an axis
or a vector, and they can be parameterized as
with
,
where the operator
rotates an axis or
a vector by the angle
around
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
-dependent quantity
.
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
is chosen
as
,
i.e. the one perpendicularly pointing from Q to
the side
of the triangle. The axis
is defined as the specific x-axis on which the collinear
configuration
satisfies two conditions:
and
,
where
and
are, respectively, the x-coordinates of the mass centers
of
and
.
With subscript permutations, the
axes
and
can be
defined in the same way as
and
,
respectively.
As is clear in Appendix A,
is a piecewise
function with 12 pieces, and it is generically not single-valued at
.
These facts make it
inconvenient to use a generally applied method to find the maximum
value of the function
.
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 .
As mentioned in Sect. 2, for any axis with which the detection is successful,
the characteristic length
will be limited
from below by the increasing function of time
.
Though this does
not necessarily mean that
is itself an
increasing function of time, its value at any
cannot be
less than
.
On the other hand, as a length
characterizing the system size along a given direction,
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.
![]() |
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
|
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
,
which are uniformly and
independently distributed in (0,1). And we take
as the three masses, where
.
Obviously this means that the chosen unit for mass is

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,




where

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
,
which are
uniformly and independently distributed in (0,1). And, we calculate
the following weighted square mean of them
.
A unit for time such that the universal gravitational constant G=1is then

which is simply the mean crossing time defined as



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 (
)
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,
,
the instant at which
,
where
is
the instantaneous semi-major axis of the binary and
the
Jacobian distance measured from the mass center of the binary to the
escaper. In accord with this, the time duration
,
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
.
In this case,
the initial system will be referred to as a post escape system.
![]() |
Figure 3:
Histogram of |
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
,
where
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.
![]() |
Figure 4:
Evolution of the ratio
|
Open with DEXTER |
To verify that M903P detect only escapes that are about to happen,
we draw the histogram in Fig. 3 of
for the
systems in S1. In this figure, the initially post escape systems
(
of the 1000 systems) are included in the leftmost bin, and
the percentages there are the accumulated ones of the initial
systems with
.
From Fig. 3,
we see that an escape detected by M903P generally has small
,
as compared with
.
In Fig 4, the time
evolution of the ratio
is shown for each one of the 6 systems with
.
As can be inferred from this figure,
there will be no close triple approach after t0.
![]() |
Figure 5: Histograms of k0 for the initial systems in S1 and S2. |
Open with DEXTER |
![]() |
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
(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
![]() |
Figure 7:
Percentage of M903P detected escapes as a function of the time elapse from |
Open with DEXTER |
To achieve our second goal, we apply M903P to systems in S2,
using system states at 17 different time instants
.
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
,
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
before
.
But it increases quickly to much higher values, e.g.,
it reaches
at
and
at
.
Recalling that
is the time instant at which
,
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,
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 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.
![]() |
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
|
Open with DEXTER |
![]() |
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
|
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 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 .
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
.
Because of the hierarchical feature of the
system configurations, not all of the typical characteristics of the
function
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
,
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
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 .
![]() |
Figure 10:
Graphs of piecewise functions
|
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 ,
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
,
where E and
are, respectively, the energies
of the whole system and of the ``binary'' formed by the two bodies
with the least mutual distance, and
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
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
,
and the opposite is true if
.
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 ,
to remove the loss in the estimation of
d
/dt2 for different s satisfying
.
In the context of the present work, it should be pointed out that,
for the 3-body case, because of the equivalence between
and
(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 .
Therefore,
this appendix will first study this function. The notations
introduced in Sect. 3 are used without further
explanation.
The piecewise function
depends on

where, for example, ma and








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
Here,





![]() |
Figure A.1:
Open angular regions
|
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 ,
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
.
This implies that
is a
critical axis with respect to the ordering of the bodies when they
are projected on
.
In fact, the order changes
from
xa<xc<xb to
xa<xb<xc as
rotates
anticlockwise across
.
This order will remain unchanged
until this rotating
reaches
,
on which
we have xa=xb.
Under permutation
,
the above
description of the angular region from
to
also applies to the one from
to
(from
to
,
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
.
To understand the role of the
direction
in Fig. A.1, we calculate the
difference between
and
for an x-axis
inside the angular region going from
anticlockwise to
,
where









Based on the above arguments, we have for
(closure of I1), i.e.
,
where








![\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}](/articles/aa/full_html/2009/34/aa09392-08/img189.png)
It is then obvious that, on



![$[\alpha_{\vec{r}_1}-\alpha_{\vec{j}_2},\alpha_{\vec{r}_1}]$](/articles/aa/full_html/2009/34/aa09392-08/img193.png)








The equation in (A.4) is solved for the angle


![$[-\gamma_b,\gamma_b]$](/articles/aa/full_html/2009/34/aa09392-08/img204.png)

![$[\pi-\gamma_b,\pi+\gamma_b]$](/articles/aa/full_html/2009/34/aa09392-08/img206.png)











![$\beta_c
~\overline{\in}~
[\alpha_{\vec{r}_1}-\alpha_{\vec{j}_2},\alpha_{\vec{r}_1}]$](/articles/aa/full_html/2009/34/aa09392-08/img213.png)
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
(see Fig. A.1),
i.e.
.
Here, we
have to again include the axis
,
which was considered as a boundary axis of
.
This is
because, though we have
for this axis, q=q1assuming
is generally different from
assuming
.
To understand this, one simply notes that
only depends on the position of the 3-body system
in the configuration space, while q also depends on the
velocities.
We now consider
.
For this, it is
helpful to recall the equivalence between
and
.
In terms of the notations used in
Fig. A.1, this equivalence can be interpreted as that
between
and
,
which means that we can
forget
as long as
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
applies exactly to
.
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
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
, assuming
; calculate
(if it is in the required interval
, then the corresponding critical axis
is a local minimizer of p, or equivalently, a local maximizer of q1). The required formulae for this step are
where k is an integer. It can be checked thatis a necessary condition for
. This means that we only need to try these two values of k to decide whether
is in the required interval.
- 3.
- Calculate q1 according to (A.3) with
, and
(if
).
- 4.
- Using the maximum q1 obtained at step 3, apply test
(4).

The above discussion does not apply to the measure zero class of Dcharacterized by the collinearity of the three bodies. This is
because
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
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 (
or
)
or there are two superposed pieces
(
)
of
associated
with each given ordering.
![]() |
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
that contains the
known motion states by choosing any one of the two opposite axes
perpendicular to
as the
-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
.
Then, denoting the two
exterior bodies by ma and mc, we set
if the
resulting
,
where
.
If it happens that the
so-defined
,
we use two different coordinate
systems, referred to as non-primed and primed ones. For the
non-primed one, we define
,
and
for the primed one
.
Finally, we complete the coordinate system(s) by setting
(and
).
By the previously mentioned equivalence between any two opposite
x-axes, we only need to consider the x-axes with
(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
,
one easily verifies directly or by adapting (A.2) to the
present setting,
From (A.7), it is clear that we either have only one piece of




In the case of ,
we have q=q1. And, (A.3) with the
present specific
is reduced to
which can also be checked directly with the present setting. According to (A.8),






Now we come to the case of d=0, for which we have to consider both
and
in searching for the
maximum value of q. Again, the trick is to consider instead
and
,
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
.
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.
as illustrated in the right panel of
Fig. A.2, where the direction of
-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
remains true for all x-axes directing
upward. And, changing from
to
makes no difference to this discussion for the
case where there is no superposed body and
.
References
- Aarseth, S. J. 2003a, Gravitational N-body Simulations: Tools and Algorithms (Cambridge University Press) (In the text)
- Aarseth, S. J. 2003b, Regularization tools for binary interactions, in Astrophysical Supercomputing Using Particle Simulations, ed. J. Makino, & P. Hut (ASP), 295 (In the text)
- Aarseth, S. J. 2004, in The Environment and Evolution of Double and Multiple Stars, IAU Colloq. 191, RevMexAA, Ser. Conf., 21, 156
- Aarseth, S. J., & Zare, K. 1974, Celest. Mech., 10, 185 [NASA ADS] [CrossRef] (In the text)
- Aarseth, S. J., Anosova, J. P., Orlov, V. V., & Szebehely, V. G. 1994, Celest. Mech. & Dyn. Astron., 60, 131 [NASA ADS] [CrossRef] (In the text)
- Agekian, T. A., & Anosova, J. P. 1967, Astron. Zh., 44, 1261 [NASA ADS]; Soviet Astron., 2, 1006 (In the text)
- Agekian, T. A., & Anosova, J. P. 1968, Astrofizika, 4, 31
- Agekian, T. A., & Anosova, J. P. 1990, Celest. Mech. & Dyn. Astron., 49, 145 [NASA ADS] [CrossRef] (In the text)
- Anosova, J. P. 1986, Ap&SS, 124, 217 [NASA ADS] [CrossRef] (In the text)
- Anosova, J. P.. 1990, Celest. Mech. Dyn. Astron., 48, 357 [NASA ADS] [CrossRef]
- Anosova, J. P., & Orlov, V. V. 1992, A&A, 260, 473 [NASA ADS] (In the text)
- Anosova, J. P., & Orlov, V. V. 1994, Celest. Mech. & Dyn. Astron., 59, 327 [NASA ADS] [CrossRef] (In the text)
- Anosova, J. P., Orlov, V.V., & Aarseth, S. J. 1994, Celest. Mech. Dyn. Astron., 60, 365 [NASA ADS] [CrossRef]
- Birkhof, G. D. 1927, Dynamical Systems, American Mathematical Society, 247 (In the text)
- Chenciner, A. 2007, Three body problem, Scholarpedia, http://www.scholarpedia.org/article/Three_Body_Problem (In the text)
- ESA 1997, The Hipparcos and Tycho Catalogues, 1239, 0 (In the text)
- Georgakarakos, N. 2008, Celest. Mech. & Dyn. Astron., 100, 151 [NASA ADS] [CrossRef] (In the text)
- Griffith, J. S., & North, R. D. 1974, Celest. Mech., 8, 473 [NASA ADS] [CrossRef] (In the text)
- Henry, T. J., & McCarthy, D. W., Jr. 1993, AJ, 106, 773 [NASA ADS] [CrossRef] (In the text)
- Kiseleva, L. G., Eggleton, P. P., & Orlov, V. V. 1994, MNRAS, 270, 936 [NASA ADS] (In the text)
- Kuznetsova, E. A., & Orlov, V. V. 1983, Trudy XII Stud. Conf., Leninfgrad Univ. (In the text)
- Laskar, J., & Marchal, C. 1984, Celest. Mech., 32, 15 [NASA ADS] [CrossRef] (In the text)
- Lattanzi, M. G., Jancart, S., Morbidelli, R., et al. 2004, in The Three-Dimensional Universe with Gaia, ed. C. Turon, K. S. OFlaherty, & M. A. C. Perryman, 251 (In the text)
- Marchal, C. 1974, Celest. Mech., 9, 381 [NASA ADS] [CrossRef] (In the text)
- Marchal, C., Yoshida, J., & Sun, Y. 1984a, Celest. Mech., 33, 193 [NASA ADS] [CrossRef] (In the text)
- Marchal, C., Yoshida, J., & Sun, Y. 1984b, Celest. Mech., 34, 65 [NASA ADS] [CrossRef]
- Marchal, C. 1990, The Three-body Problem (Elsevier) (In the text)
- Mikkola, S. 1983, MNRAS, 203, 1107 [NASA ADS] (In the text)
- Orlov, V. V., & Zhuchkov, R. Ya. 2005, Astron. Rep., 49, 3 [CrossRef] (In the text)
- Orlov, V. V., Petrova, A. V., & Martynova, A. I. 2002, MNRAS, 333, 495 [NASA ADS] [CrossRef]
- Saslaw, W. C., Valtonen, M. J., & Aarseth, S. J. 1974, ApJ, 190, 253 [NASA ADS] [CrossRef] (In the text)
- Standish, E. M. 1971, Celest. Mech., 4, 44 [NASA ADS] [CrossRef] (In the text)
- Standish, E. M. 1972, A&A, 21, 185 [NASA ADS]
- Szebehely, V. 1972, Celest. Mech., 6, 84 [NASA ADS] [CrossRef]
- Tevzadze, G. A. 1962, Izv. Akad. Nauk Armyan. S. S. R., 15, No. 5, 67 (In the text)
- Tokovinin, A. A. 1997, A&AS, 124, 75 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
- Valtonen, M., Mylläri, Orlov, V., & Rubinov, A. 2005, MNRAS, 364, 91 [NASA ADS] [CrossRef] (In the text)
- Yoshida, J. 1972, PASJ, 24, 391 [NASA ADS] (In the text)
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
![]() |
Figure 1:
Critical x-axes separating different pieces of the piecewise function
|
Open with DEXTER | |
In the text |
![]() |
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
|
Open with DEXTER | |
In the text |
![]() |
Figure 3:
Histogram of |
Open with DEXTER | |
In the text |
![]() |
Figure 4:
Evolution of the ratio
|
Open with DEXTER | |
In the text |
![]() |
Figure 5: Histograms of k0 for the initial systems in S1 and S2. |
Open with DEXTER | |
In the text |
![]() |
Figure 6: Respective distributions of systems in S1 and S2 on the mass ratio plane. |
Open with DEXTER | |
In the text |
![]() |
Figure 7:
Percentage of M903P detected escapes as a function of the time elapse from |
Open with DEXTER | |
In the text |
![]() |
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
|
Open with DEXTER | |
In the text |
![]() |
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
|
Open with DEXTER | |
In the text |
![]() |
Figure 10:
Graphs of piecewise functions
|
Open with DEXTER | |
In the text |
![]() |
Figure A.1:
Open angular regions
|
Open with DEXTER | |
In the text |
![]() |
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.