A&A 378, 569586 (2001)
DOI: 10.1051/00046361:20011189
K. Gozdziewski^{1,2}  E. Bois^{1}  A. J. Maciejewski^{3}  L. KiselevaEggleton^{4,5}
1  Observatoire de Bordeaux, UMR/CNRS/INSU 5804,
BP 89, 33270 Floirac, France
2 
Torun Centre for Astronomy,
N. Copernicus University,
Gagarina 11, 87100 Torun, Poland
3 
Institute of Astronomy, Zielona Góra University,
Lubuska 2, 65265 Zielona Góra, Poland
4 
IGPP, LLNL,
L413, 7000 East Ave, Livermore, CA 94550, USA
5 
UC Davis, Dept. of Physics, Davis, CA 95616, USA
Received 11 June 2001 / Accepted 22 August 2001
Abstract
In this paper we apply a new technique alternative to the numerically
computed Lyapunov Characteristic Number (LCN) for studying the dynamical
behaviour of planetary systems in the framework of the gravitational Nbody
problem. The method invented by P. Cincotta and C. Simó is called the Mean
Exponential Growth of Nearby Orbits (MEGNO). It provides an efficient way for
investigation of the fine structure of the phase space and its regular and
chaotic components in any conservative Hamiltonian system. In this work we use
it to study the dynamical behaviour of the multidimensional planetary
systems. We investigate the recently discovered And planetary
system, which consists of a star of
and three Jupitersize
planets. The two outermost planets have eccentric orbits. This system appears
to be one of the best candidates for dynamical studies. The mutual
gravitational interaction between the two outermost planets is strong.
Moreover the system can survive on a stellar evolutionary time scale as it is
claimed by some authors (e.g., Rivera & Lissauer 2000b). As the main
methodological result of this work, we confirm important properties of the
MEGNO criterion such as its fast convergence, and short motion times (of the
order of 10^{4} times the longest orbital period) required to distinguish
between regular and chaotic behaviors. Using the MEGNO technique we found
that the presence of the innermost planet may cause the whole system to become
chaotic with the Lyapunov time scale of the order of 10^{3}10^{4} yr only.
Chaos does not induce in this case visible irregular changes of the orbital
elements, and therefore its presence can be overlooked by studying
variations of the elements. We confirm explicitly the strong and sensitive
dependence of the dynamical behaviour on the companion masses.
Key words: celestial mechanics, stellar dynamics  methods: numerical, Nbody simulations  planetary systems  stars: individual ( Andromedae)
Although the other orbital parameters are supposed to be determined without significant ambiguity (it is in fact a very rare situation), we have to explore a huge volume of parameter space in order to determine possible dynamical states and different zones of stability of planetary systems. Relatively large masses, highly eccentric orbits and small semimajor axes may cause strong gravitational interactions between planets, what may end up in chaotic states and unpredictability. Consequently, the factors may lead in some cases to the disintegration of the planetary system. This is why an efficient method, which can help to identify and localize such states in the parameter space is highly desirable. This method should provide additional dynamical constraints to the conditions under which the system can survive on an evolutionary time scale. Among such tools for investigation of global dynamics, the classical method is the Lyapunov Characteristic Number (LCN hereafter). It is widely applied and well understood (e.g., Benettin et al. 1976; Benettin et al. 1980; Wisdom 1983; Froeschle 1984; Sussman & Wisdom 1988; Mikkola & Innanen 1995; Froeschlé & Lega 2000; Tancredi et al. 2001). In this paper we consider an alternative technique invented by Cincotta et al. (Cincotta & Simó 1999, 2000; Cincotta & Giordano 2001; Cincotta & Núnez 2000). This method seems to have better computational properties than the LCN, and can be very useful for efficient discrimination of regions in parameter space corresponding to quasiperiodic and chaotic motions. It will be called MEGNO hereafter (the acronym of Mean Exponential Growth of Nearby Orbits). We describe the theory of this method in Sect. 2.
The main purpose of this work is to apply the MEGNO method to the study of the dynamics of planetary systems. As the object of our investigations we choose the And planetary system. In Sect. 3 we present two sets of observational parameters (Lick and AFOE) related to this system that we use in this paper. In Sect. 4 we formulate the equations of motion and discuss details of their numerical solutions. Until now the MEGNO was applied only to systems with two or three degrees of freedom. The dynamics is much more complicated in higher dimensional systems. Therefore a careful investigation of the MEGNO properties is required when the method is applied to such systems. In Sect. 5 we present preliminary numerical tests performed on a fictitious planetary system, and also on the Solar system. Section 6 contains results for the nominal And orbital elements (Stepinski et al. 2000) as well as tests that we performed on the MEGNO indicator. Section 7 presents a look at the global dynamics of the And system. We explore to a certain extent its orbital parameter space, in order to locate possible regular motions of the system, as well as to investigate the properties of the MEGNO technique.
Considering a planetary system as N point masses in gravitational
interactions, we can look for its regular or irregular states by computing the
LCN. This quantity is known to be a measure of the hyperbolicity degree (the
divergence rate of nearby orbits) in the dynamical system, described here as a
set of ordinary differential equations:
The MEGNO method introduced by Cincotta and Simó belongs to this group of methods, called fast indicators. Detailed inspection of its properties (see CSGN) inspired us to expect that it can be very useful for studying the dynamics of planetary systems. The MEGNO provides more informations than other indicators about the dynamics of the system (see CSGN), and its properties are justified by solid analytical investigations and numerical tests. It originates from the concept of the symmetric conditional entropy of nearby orbits (Cincotta & Simó 1999), and an heuristic observation that close regular trajectories are strongly correlated with each other, but in chaotic domains of the phase space the nearby orbits loose this correlation very quickly.
To define the MEGNO indicator following Cincotta & Simó (2000), let us notice
that the definition of the LCN can be rewritten in the following way
This property can also be used to estimate efficiently the LCN, which can be recovered by the linear least square fit to . This procedure uses dynamical information collected over the whole time of calculations (see CSGN for details). It requires also much shorter motion times than the classical methods of neighboring trajectories (twoparticle method) or the method of variational equations. The LCN derived on the basis of the MEGNO is very sensitive to the fine structure of phase space. In this work we do not make explicit use of this feature, but it is certainly a fundamental one, because it can be used, for example, to find resonances and their widths. For details and discussion see CSGN.
Comparing the convergence rates of Y(t) and the LCN, one obtains the following
estimates when
:
From the computational point of view the MEGNO is relatively easy to implement. We have to write
down the equations of motion together with the variational equations, and
then Y(t) and
can be found by solving this system of equations,
complemented by the following two supplementary equations:
The recent paper by Tancredi et al. (2001) shows that the variational method of computing the LCN provides more reasonable results than computationally less expensive twoparticle method (Benettin et al. 1976). Thus an estimation of the MEGNO requires almost the same computational efforts per unit of dynamical time as the LCN but the required motion time is much shorter for the MEGNO.
By evaluating the LCN, we finally get only one number and after long computations all the dynamical information (which is in fact hidden in the evolution of variationals ), is neglected. Also the discrimination of regular and chaotic orbits with the LCN is not always obvious. It may happen that in a case of chaotic behaviour we get a clear sign of it. But it is very hard or even impossible to prove that, when one gets the LCN converging typically for regular motion, its limit is really nonzero (Murray & Dermott 2000). A striking example is provided by the computation of the LCN for Pluto and the outer Solar system (Sussman & Wisdom 1988). After indicating over more that 10^{8} yr a quasiregular evolution, the LCN of this system actually takes a nonzero value, then giving a divergence rate of the order of .
Compared to the LCN, the MEGNO offers much more detailed identification of the regular and chaotic orbits in a system. Because of the "amplification" of chaotic effects in the evolution of the variationals (see the definition, Eq. (6)) the MEGNO allows to detect signs of chaotic behaviour much earlier than the LCN. Thus it helps to resolve the stability question more efficiently  we shall demonstrate this later in this paper. This efficiency is highly appreciated when one deals with the global dynamics and requires to examine large sets of initial conditions (of order of 10^{3}10^{4}).
We refer to the paper by Stepinski et al. (2000), SMB paper hereafter, where we
have found the initial conditions of the And. By independent
analysis of the observations published by Butler et al. (1999), SMB found two sets
of planetary orbital parameters (see Table 3 in SMB and Table 1
here),
Planet  Observatory  mass [ ]  a [AU]  P [d]  e  [deg]  [deg]  M [deg] 
B  Lick  0.68811  0.0592  4.617  0.0430  0.0  0.58  266.48 
C  Lick  1.8894  0.8282  242.0  0.3478  0.0  248.21  123.13 
D  Lick  3.9087  2.5334  1269.0  0.2906  0.0  242.99  354.78 
B  AFOE  0.736  0.0592  4.617  0.042  0.0  63.827  203.515 
C  AFOE  1.905  0.8326  234.5  0.230  0.0  219.214  151.605 
D  AFOE  5.150  2.7865  1481.2  0.426  0.0  245.627  348.759 
It is claimed (Stepinski et al. 2000; Jiang & Ip 2001) that the innermost planet, being the least massive and very close to the star, does not affect in the first approximation the longterm dynamical evolution of the system. By neglecting this planet it is possible to perform longterm integrations with much larger time step ( days) than it is permitted when the all the planets are taken into account ( days). Thus the dynamical timescale of the system increases by 2 orders of magnitude.
For that reasons we focus mostly on twoplanet systems containing the central mass of the star ( ) and two companion point masses labeled hereafter by C and D (the innermost planet is called B). However in the light of our results (presented below), the presence of the innermost planet may change qualitatively the dynamical behaviour of the system in the strict sense of distinction between quasiperiodic and chaotic motion.
Our study is not complete in the sense that we do not explore all the available parameter space, and we do not take into account formal errors of the orbital parameters and other factors, which may reflect the dynamical behaviour of the And (such as relativistic precession, quadrupole distortion of the system bodies due to their mutual gravity and intrinsic spin, tidal friction). Our main objective is to present the MEGNO technique in a practical way, and to compare the results with those given by the LCN method. We think that application of the MEGNO will greatly extend conclusions obtained by investigating the dynamics of planetary systems by numerical Nbody integrations and by studying variations of orbital elements. The latter approach is used as the general basis of many papers devoted to investigating the dynamical stability of the And, see e.g. papers of Laughlin & Adams (1999), Stepinski et al. (2000), Rivera & Lissauer (2000a,b), Jiang & Ip (2001), Barnes & Quinn (2001), Ito & Miyama (2001), Chiang et al. (2001).
In order to compute the MEGNO one has to solve the system of first order
differential equations built from the equations of motion and the variational
Eq. (3). To avoid explicit computations of the equations of
motion of the central mass, the equations are transformed to the barycentric
reference frame, defined through
The preservation of the integrals of the total energy and the angular momentum, as well as their numerical shifts, have been controlled within all the numerical integrations. The integrals are preserved with a relative accuracy at the level of 10^{8}  10^{12} over the timescale of the order of 10^{4}10^{5} orbital periods of the outermost body, depending on the given test. Integrations of the 3planet Andr systems lead to quite significant degradation of the relative accuracy (i.e. of the order of 10^{9}), probably due to a large number of equations and a large range of the orbital periods (from 4 to 1300 days).
Estimation of a proper period of time required for computing the MEGNO is crucial for both the computational efficiency and the classification of the dynamical state of the studied system. According to CSGN, this time should be of the order of , where is the characteristic period of the system. In the case of the And the period of the outermost planet is equal to about 1300d, and so the integration time should be of the order of yr. We found that the upper limit was usually overestimated in the case of chaotic evolution, but we have found some border cases for which the period of time is not long enough to get a definitive conclusion. This concerns the situation when the motion evolves on the border between regular and chaotic states. But this is not a very surprising fact as it is a common problem for every known numerical indicator of dynamical state.
Because the variational (see the definition of the MEGNO) grows at an exponential rate in the case of chaotic motion, it is necessary to renormalize after a period of time, which we call here the renormalization step .
The normalization of follows from the definition, Eqs. (6), (5) and the fact that depends linearly on . The total integration time T consists of renormalization intervals .
We mentioned already that the MEGNO technique has been applied up to now to relatively simple 2D and 3D Hamiltonian models (Cincotta & Simó 2000; Cincotta & Núnez 2000). Thus from the numerical standpoint, it is important to understand the behaviour of the MEGNO when applying it to higherdimensional systems, as for instance to the 4body model describing the And. This is why, at first, we studied many simplified, and in some sense artificial models. One of these tests is presented in Fig. 1.
Planet  mass [ ]  a [AU]  e  i [deg]  [deg]  [deg]  L [deg] 
Earth  0.003146  1.00000011  0.01671022  0.00005  11.26064  102.94719  100.46435 
Jupiter  1.0  5.20336301  0.04839266  1.30530  100.55615  14.75385  34.40438 
Saturn  0.299410  9.53707032  0.05415060  2.48446  113.71504  92.43194  49.94432 
Pluto  39.48168677  0.24880766  17.14175  110.30347  224.06676  238.92881 
Note 1: Reference: Explanatory Supplement to the Astronomical Almanac, 1992,
K. P.
Seidelmann,
Ed., p. 316 (Table 5.8.1), University Science Books, Mill Valley,
California.
Figure 1: The MEGNO calculated for the Lick Keplerian elements of the And, but for very small masses of the "planets" . a) The upper panel shows the time evolution of Y(t), b) the mean over the motion time of 10^{4} periods of planet D is marked by dashed line, c) the LCN computed by the direct variational method, and its estimation by the evolution law of the MEGNO. The integrator accuracy is set to , the renormalization step is days, and the number of steps is .  
Open with DEXTER 
For a similar test on two planetary systems with three planets each, we used
initial conditions describing orbits of some bodies of the Solar system. The
data are given in Table 2.
The first computation concerns Earth,
Jupiter and Saturn, and the results are shown in Fig. 2.
Figure 2: The MEGNO calculated for the 3planet system of Earth (E), Jupiter (J) and Saturn (S). The two upper panels show the time evolution of Y(t) computed with two different choices of the initial tangent variation . The two middle panels show the resulting plots of the mean . The two bottom panels illustrate the evolution of the eccentricities and inclinations of the orbits. Integration was performed over yr. Relative errors of the total energy the angular momentum do not exceed 10^{12} in the runs.  
Open with DEXTER 
A similar experiment on a part of the outer Solar system (Jupiter, Saturn and Pluto) leads however to serious difficulties of interpretation. In fact we choose it as a challenging problem for the method, knowing that Pluto considered as a member of the outer Solar system exhibits chaotic orbital behaviour, which was confirmed after very careful, longterm integrations (Sussman & Wisdom 1988). In the JupiterSaturnPluto system, the roughlyestimated integration time should be of the order of  yr. In order to have a clear answer we continued the integration up to 10^{7} yr. The MEGNO was computed on two different processor architectures  a DEC/Alpha machine with 16 digits double precision accuracy, and a Linux/Intel machine with the same double precision accuracy but the Intel's coprocessor works internally on 19 digits, what decreases the roundoff errors.
Our results are shown in Fig. 3.
Figure 3: The MEGNO calculated for the 3planet system of Jupiter, Saturn and Pluto. Panels ac) show the time evolution of the MEGNO computed with the DEC/Alpha and the Intel machines with integration accuracy yr a), on the Alpha with yr b), and on the Intel with yr, (thin line) as well as (thick line), panel c). Relative errors of the energy and the angular momentum are of the order of 10^{9}and 10^{10}, depending on the test. Panel d) shows the maximal Lyapunov exponent computed over time interval of 1 Gyr.  
Open with DEXTER 
Figure 4: The MEGNO indicator computed for the nominal Lick data set, when 2 [CD] and 3 [BCD] planets are taken into account, respectively. From the top: a) the time evolution of the MEGNO Y(t) for the system including planets C and D, b) Y(t) for the system of 3 planets, c) the mean value of the MEGNO , d) variations of the eccentricity of planet C in the two models, e) the same quantity for the outer planet D.  
Open with DEXTER 
Figure 5: The MEGNO computed for the nominal AFOE data set, when 2 [CD] and 3 [BCD] planets are taken into account, respectively. From the top: a) the time evolution of the MEGNO Y(t) for the system including planets C and D, b) the mean value of the MEGNO for the 2planet system, with different choices of the initial tangent vector (at random), c) Y(t) for the system of 3 planets (dotted curve) with the graph of (thick, smooth curve)  the factor 2 was added due to the averaging of Y (see papers by Cincotta and coauthors for details), d) the MEGNO computed with the data of the c) case, e) the LCN computed by the direct variational method, and also estimated by the MEGNO evolution law .  
Open with DEXTER 
Let us note that the system of JupiterSaturnPluto considered here is in some sense artificial thus the conclusions concerning its behaviour cannot be strictly related to the whole outer Solar system.
Unfortunately in such unclear situations the theory of stability of the Nbody system cannot help. The analytical estimations of the orbital parameters and masses that provide that the bodies follow bounded trajectories, are very restrictive on the assumptions. Recently Niederman (1996) proved theorems extending the Nekhoroshev stability theory. When they are applied to the SunJupiterSaturn system, they ensure its time of stability comparable with the age of the Solar system (i.e., 6 Gyr) whether the ratio of the mass of the heavier planet over mass of the central body is lower than 10^{13}. This is very far from the actual ratio of 10^{3} between Jupiter and the Sun. It is even worse in the case of new planetary systems, because the theories concerning planetary stability are applicable to systems close to integrable ones. This assumption is no longer valid when the planetary masses, eccentricities and inclinations of orbits are large.
In the following tests of the MEGNO, for the nominal sets of initial conditions, we used the Lick and the AFOE orbital parameters (see Table 1). The computation of the MEGNO over yr (corresponding to about 10^{4} revolutions of planet D) is presented in Fig. 4 for the Lick data, and Fig. 5 for the AFOE parameters, respectively. The indicator was computed for a system firstly with the two outer planets, and secondly with all the three companions included, in order to check whether planet B influences the dynamical regime of the system. For the Lick data in both cases we obtained an indication of the quasiperiodic regime, as the mean value of Y(t) saturates very fast to 2with characteristic oscillations of decreasing amplitude. The oscillations of may be interpreted as an effect of close encounters of the system orbit (as a whole) with unstable manifolds of invariant objects existing in the phase space. The qualitative behaviour of the system does not change in both cases. The presence of planet B induces some small changes in periods of orbital elements  it is illustrated in panels d) and e) of Fig. 4.
The case of the AFOE data is partly different. The 2planet system has the MEGNO signature of the quasiperiodic regime (Figs. 5a,b). We tested it by repeating the calculations with the variational chosen at random, and always found behaviour characteristic for regular orbits.
For the 3planet model we found clear signs of chaotic behaviour. In this
example the system remains close to the quasiperiodic regime during less than
10^{4} yr, as
goes very close to 2. However, after this time a large
jump of
occurs. This effect was verified by integration with
different accuracies, the renormalization step size, and the total time of
integration. The LCN has been simultaneously evaluated by the variational
method in the same program runs. The results are shown in
Fig. 5e  both the LCN and its estimation by the evolution
law
 are in a good agreement since they converge to the same
value. The estimate of the LCN by the MEGNO converges faster to zero than the
estimation provided by the direct computation.
Figure 6: Effect of sensitivity to initial conditions for the system related to the AFOE elements. The plots in the left column show the data computed with the initial conditions differing from the 14th digit in the heliocentric xposition of planet C  the dotted graphs are for the unchanged data. After yr the eccentricities as well as other elements (not shown in the plots) differ by values of the first order of magnitude. The plots in the right column show absolute differences between the eccentricities obtained by the same way when using the Lick parameters (straight lines centered at 0) and the AFOE data.  
Open with DEXTER 
Figure 7: a) The MEGNO indicator computed for the nominal geometrical Lick elements ( , d, and ). The levels of the nominal (minimal) mass of planet C and the higher value are indicated by the dashed lines. b) the MEGNO computed with a better accuracy, namely , the same renormalization step , and (1000 points were examined). c) Zoom of panel b). d) Estimation of the LCN by a linear fit to the data taken into account in panel b).  
Open with DEXTER 
Figure 8: The MEGNO computed for the nominal, geometrical Lick elements, with another mass of planet C, namely (it is marked in Fig. 7a). Panel a) shows four different plots of the MEGNO obtained with different choices of the initial tangent vector . Panel b) is a magnification of plot a). Panel d) shows time variation of Y(t) for one of the seemingly converging plots shown in panel b). Panel c) is for the diverging case. Panels e) and f) illustrate computations of over a longer time span than seen on a), as well as an independent estimation of the LCN in comparison with predictions obtained by the MEGNO (curve labeled by ).  
Open with DEXTER 
Moreover the time evolution of the AFOE orbital elements appears to be quasiperiodic, even over a long period of time (see for example SMB). In order to explain this paradox we computed the following test. Using the CowellStormer integrator (as an independent integration scheme) by Varadi (1996), we followed the evolution of the system defined by two sets of very similar initial conditions. One of them was defined by the nominal AFOE elements. As the second set we get the nominal data modified by a change of the x heliocentric coordinate of planet C in the 14th digit (specifically from 5 to 6). We run two integrations over the time span of 10^{5} yr with time step 0.004d, following the code author's recommendation to use the integration time step of , where is the shortest orbital period. That choice provided a local error of the order of 10^{19}. Relative errors of the energy and angular momentum did not exceed 10^{12} in the whole run. Results of the comparison are given in Fig. 6. After less than yr the eccentricities as well as other elements (not shown) differ by values of the first order of magnitude  it is illustrated by the plots. In a similar experiment performed on the Lick data we did not find any notable change in the orbital elements. Estimation of the Lyapunov time for the AFOE system, which is very short, of the order of 10^{3} yr (see Fig. 5d), is consistent with other estimations quoted by Laughlin & Adams (1999), although these authors deal with the exponential divergence of eccentricities instead of the strictly defined LCN.
Our results are in agreement with Laughlin & Adams (1999) and Stepinski et al. (2000), who found that the innermost planet produces an effect rather destabilizing than stabilizing the system. We would like to emphasize that as soon as one considers the 2planet configurations, one describes only necessary but not sufficient constraints on the system stability.
As mentioned already, the method of radial velocity measurements does not
allow in general to determine real values of the planetary masses. The second
main indeterminacy follows from the unknown relative inclinations between
planetary orbits. In the case of a 2planet system the relative inclination of
orbits
is related to the longitudes of nodes
by the following formulae (used also
by SMB):
In a first simple examination of the space of the undetermined orbital
parameters, we may assume that one of the orbits has ,
i.e.
(so, the corresponding planet has the nominal minimal mass), but
the inclination of the second planet can be varied, and its mass as well. In
this test the lines of nodes of planets C and D are parallel and perpendicular
to the line of sight. We performed such an experiment on the basis of the Lick
data and its results are illustrated in Fig. 7.
We varied the
mass of the planet C from an artificial limit
up to the
value of
,
(so the mass hierarchy of the planets C and D can
be switched in some cases), that corresponds to relative inclinations
in the range

.
The first run we
did with days and
steps with integrator
accuracy
.
That gave us
yr. The resulting graph of
as a function of mass
is shown in panel (a) of Fig. 7. It is easy
to see that the minimal mass of C, namely
,
is located within a
wide range of masses for which the MEGNO converges very close to the value of
2, and so implies regular motion. In order to check the resolution ability of
the technique, and its application on a borderline between chaotic and regular
regimes, we examined in detail its time behaviour for the value
(see a narrow valley in panel 7a). We run a few computations of
with different
values of
.
In some cases we obtained apparently very good and
fast convergence, but one of the runs
ended up with a value about 2.5,
indicating no convergence at all (Fig. 8a).
A closer look at
the "converging'' plots (Fig. 8b) shows that in fact they are
not smooth as one would expect in the case of regular motion. In fact the time
evolution of Y(t)  panels c) and d)  also shows this very clearly. The
calculation extended over a longer period of time (panel f) shows the
indication of a chaotic regime, which can be also verified by the LCN computed
by the variational method  panel e). Unfortunately, the integration time,
which is required to get the true answer about the chaotic nature of this
system, is simply
longer than we originally used. The temporal plateau seen in panel b) correspond to almost quasiperiodic regime, i.e., the
system orbit is evolving on a border of regular and irregular states. It is a
good example showing that the MEGNO (as other fast indicators) may fail in
some borderline cases, and one has to use it carefully. It illustrates also a
possible dependence of the MEGNO on the choice of
:
we obtain an
appearance of different kinds of
evolution when the time span is too
short.
Figure 9: Examples of typical behaviour of the MEGNO and its comparisons to the LCN (computed in the same program runs by the variational method). Keplerian elements are defined by the Lick data. The choices of the planet C mass are coherent with Fig. 7. The right panels show the estimation of the LCN both by the law (continuous lines) and by the direct variational method (dotted lines).  
Open with DEXTER 
Figure 10: The MEGNO evaluated for the nominal geometrical Lick elements (upper plots) and AFOE data (lower plots), as a function of the inclination i. Planets C and D are included in the test. Plots on the right hand side show the estimation of the LCN based on the linear fit .  
Open with DEXTER 
The example that we examined above, brings in a question the validity of the scan, so we have calculated it again for the whole range of masses between 0 and 10 with a better accuracy ( ) and doubled steps. In that case the final plot containing 1000 points is shown in Fig. 7b. Consequently, we do not notice any qualitative differences between the two calculations; however in the range of between and , the MEGNO changed its quantitative character.
A closeup view of Fig. 7c shows an excellent convergence of
the MEGNO in the mass ranges corresponding to quasiperiodicity. As we noticed
already, it is easy to obtain an estimate of the LCN on the MEGNO basis, by
using the law
.
Figure 11: Behaviour of the MEGNO for the initial conditions discussed in the SMB paper. On the left hand side: (i.e. the nominal Lick masses are doubled), , relative inclination . On the right hand side: ( ). See the SMB paper (Fig. 7) for a comparison, and the text for details. Graphs in the two last rows show eccentricities and absolute inclinations of planets C and D for the both initial conditions. Let us observe the long periods in the two eccentricity fluctuations in the 6th panel computed with high relative inclination ; it could be caused by the Kozai cycles (Kozai 1962). ( , d, and ).  
Open with DEXTER 
This experiment shows the possibility of changing the mass hierarchy in the And. Fixing at its minimal value, we can select up to almost 2 times of , and the system still remains in the quasiperiodic regime. This conclusion is not in contradiction with the stability limits found by SMB (confirmed by us, see bellow) and by Ito & Miyama (2001), because the authors adopted equal inclinations, , for both outer planets and therefore did not switch their mass hierarchy. However, let us note, that if one simply interchanges masses and leaving all other initial orbital parameters unchanged, the Lick 2planet And system becomes unstable unless additional perturbing forces are taken into account (KiselevaEggleton & Bois 2001). Their results confirm that not only the mass hierarchy (or mass ratio of and ) is important for the dynamical stability, but also the real values of these masses.
Finally, we examine the time behaviour of the MEGNO for different masses , which correspond to different dynamical regimes, taking as a reference the data shown in Fig. 7. The results are illustrated in Fig. 9.
We selected (in the middle of the second pick in Fig. 7a), two values in the wide regular zones between and , and a large mass located in another chaotic region. The corresponding time changes of the MEGNO are shown in the left panels of Fig. 9, the right panels illustrate the estimation of both by the direct variational method (plots labeled by the LCN) and from . A very good agreement of the two methods of the LCN derivations is obvious. In the case of chaotic regime, the MEGNO shows signs of chaos by two orders of magnitude faster than a direct estimation of the LCN obtained by the variational method. Also the plots corresponding to the quasiperiodic regime show rapid, regular convergence to values very close to 2. Further, the convergence of provided by the MEGNO is faster than in the case of the LCN  it is the effect predicted by the theory.
Another onedimensional test of the parameter space available for the And system is illustrated in Fig. 10. We assumed that the longitudes of nodes are equal, the system is coplanar, and we simulated the effect of changing its observationally undetermined inclination i. The two upper plots correspond to the Lick data, and the lower plots are for the AFOE data. In the right panels, the LCN estimation is presented. The Lick parameters are much less sensitive to the inclination factor (and to much larger masses than their minimal values). Even in the range of to related to and , the system still may be found in a quasiperiodic regime. The AFOE data do not permit for such large masses. In both cases the minimal masses are located in the wide regular range.
In the next experiment, we simulated the change of the absolute inclination i, as well as the relative inclination of orbits defined through the relation involving the difference in the longitudes of the line of nodes (Eq. (12)). Let us notice that in these simulations, involving mutually inclined planets, we assumed that both the companions move on prograde orbits, thus we do not account for the general dependence of on the viewing geometry, as detailed in Chiang et al. (2001).
At first, in order to have a quantitative reference, we tested the same initial conditions as examined in the SMB paper (see its Fig. 7): SMB found that the configuration with , corresponding to relative inclination , is regular while that with ( ) leads to irregular changes of the orbital elements and a chaotic state of the system. Our computations are illustrated in Fig. 11 which shows the MEGNO variations, as well as the eccentricities and inclinations of the planets. In the both examples, the MEGNO gives an answer which is coherent with the results of 10^{6} yr integration. The irregular behaviour is detected very quickly, only after yr.
SMB examined the stability of the And system for the eccentricity ranges of planet C induced by different initial values of (and so, different masses and ) as well as by different relative inclinations . We have performed the same test using the MEGNO. The results are illustrated in a 3D plot (see Fig. 12), and its 2D version focusing on the areas of quasiperiodic states is presented in Fig. 13. Comparing the results with those obtained by SMB, we may conclude that there is a good agreement: for instance there are exist very few initial conditions related to regular behaviour for , which was determined as a border of stability by SMB for the Lick parameters. But the MEGNO permits to determine the ranges of parameters corresponding to regular states more precisely and in much shorter integration times  of the order of yr for one initial condition compared to 10^{6} yr integrations by SMB.
Our detailed scan over the undefined parameters of the
And entirely
confirms the results of SMB. In most cases the quasiperiodic behaviour of the
system provide initial conditions such that the lines of nodes are nearly
parallel, i.e. for the relative inclinations rather small. This is clearly seen
in Figs. 1214.
Figures 12, 13 present the results for the Lick data,
and Fig. 14 corresponds to the same experiments performed on the
AFOE elements.
Figure 12: The MEGNO calculated for the Lick elements of planets C and D in threedimensional form. The iabscissa is the inclination of the line of sight, so that both masses are . The second abscissa is the relative inclination of the orbits through the formula 12. The flat areas (without any peaks) indicate quasiperiodic zones and the others chaotic ones. Resolution of the grid is points. ( , d, and ).  
Open with DEXTER 
Although the coincidence of the conclusions in most cases is striking, the comparison of results given by the MEGNO with those obtained by looking for large and irregular changes of the Keplerian elements is a subtle task. Chaotic states do not necessarily mean that the system behaves in some erratic and irregular way: see the example of the 3planet AFOE system described earlier in this paper. In this sense the MEGNO rejects initial conditions, which induce a seemingly regular evolution of the orbital elements, although in a strictly dynamical meaning they generate an unpredictable behaviour. The MEGNO gives then strong constraints on the initial conditions. This may explain the difference between the stability limits of for the AFOE data found by us and by SMB. One may argue indeed that it is difficult to extrapolate the conclusions based on the 10^{6} yr integrations, when taking into account the orbital element changes, onto the evolutionary time scale. Therefore, the MEGNO can be an efficient indicator in such unclear situations.
Finally we checked whether an error of the pericenter longitude has any effect on the distinction between the chaotic and regular evolution of the 2planet system, defined by the nominal Lick parameters. From the results illustrated by Fig. 15, we may conclude that errors of the order of in mostly would not change the qualitative state of the system. Moreover there are three zones of instability represented by peaks in . They correspond to secular resonances in the pericenter motion. For detailed studies of the effects of the secular resonance on the stability of the And, see recent papers by Ito & Miyama (2001) and Chiang et al. (2001). The instabilities are much more wide when one takes larger masses; we did a test selecting for both of the planets, following Ito & Miyama (2001). The results illustrated in Fig. 16 are in a good agreement with their conclusions. The secular resonance line ( ) lies in relatively narrow zone of quasiperiodic motion. We can also see other very narrow stable zones. In these experiments we also found very good convergence of the MEGNO: in most of the points of the grid of the parameters, which were examined, the differences from the value of 2 are less than 0.005.
The MEGNO method seems to be promising and an efficient tool for studying the longtime stability and global dynamics of new planetary systems. This method belongs to the class of fast indicators, which permit efficiently distinguish between regular and chaotic evolutions of a dynamical system. Its main features, in the context of planetary dynamics, may be summarized as follows:
Figure 13: The MEGNO presented in Fig. 12 in the form of a contour plot. Values such as are marked with filled black points. Open circles indicate , and dots . Two large crosses coincide with the initial data used in Fig. 11. The value of is the stability limit quoted in Stepinski et al. (2000).  
Open with DEXTER 
Figure 14: The MEGNO plotted in the same way as in Fig. 12, but for the AFOE parameters. Filled black points mean , open circles , and dots .  
Open with DEXTER 
Figure 15: The MEGNO computed for the nominal Lick data for the two planets C and D, when taking into account a possible error of the periastron longitude. The nominal data are pointed out on the left panel by a dashdotted line. On the right panel, the same data are represented in the form of a contour plot: filled black points indicate , open circles , and dots . The nominal arguments of periastron are marked by ( ) on the left panel and by a cross on the right panel. The exact secular resonance is marked by a thick line on the right panel.  
Open with DEXTER 
Figure 16: The MEGNO computed for the geometrical Lick elements for the two planets C and D, when taking into account a possible error of the periastron longitude and assuming . The nominal data are marked on the left panel by a dashdotted line. On the right panel, the same data are illustrated in the form of a contour plot: filled black points indicate , open circles , and dots . The nominal arguments of periastron are marked by ( ) on the left panel and by a cross on the right panel. The exact secular resonance is marked by a thick line on the right panel.  
Open with DEXTER 
In this work, we did not make use of important features of the technique, as for example, its ability to resolve the fine structure of the phase space. It is provided by the LCN estimations derived from the MEGNO. These estimations are a key for localizing hyperbolic invariant objects (e.g. unstable periodic solutions, unstable tori), resonances in the phase space and their widths.
By using the method we may confirm and precise the conclusions of the SMB paper. For the Lick data and the AFOE data, we also found that the system with the two outer planets remains most probably in quasiperiodic regime under the assumption that the lines of nodes are initially closely aligned.
Moreover the technique seems to be much more efficient in bringing out the same qualitative and quantitative conclusions that those based on studying variations of the orbital elements.
From our calculations presented in this paper, one can derive the following general conclusions, which seem to be often overlooked in other studies of the And planetary system:
Acknowledgements
We would like to thank Pablo Cincotta for introducing us to the method, for immediate access to his new papers, as well as for discussions and numerous explanations concerning the technique.
We would like to thank the referee, Dr. G. Laughlin, for his comments, which have improved the paper. The postdoctoral position of K. G. is supported by a grant from the John Templeton Foundation (Agreement ID#938COS272).