A&A 378, 569-586 (2001)
DOI: 10.1051/0004-6361:20011189

Global dynamics of planetary systems with the MEGNO criterion

K. Gozdziewski1,2 - E. Bois1 - A. J. Maciejewski3 - L. Kiseleva-Eggleton4,5


1 - Observatoire de Bordeaux, UMR/CNRS/INSU 5804, BP 89, 33270 Floirac, France
2 - Torun Centre for Astronomy, N. Copernicus University, Gagarina 11, 87-100 Torun, Poland
3 - Institute of Astronomy, Zielona Góra University, Lubuska 2, 65-265 Zielona Góra, Poland
4 - IGPP, LLNL, L-413, 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 N-body 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 $\upsilon $ And planetary system, which consists of a star of $1.3~M_{\odot}$ and three Jupiter-size 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 104 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 103-104 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, N-body simulations - planetary systems - stars: individual ($\upsilon $ Andromedae)


1 Introduction

One of the main difficulties in dynamical studies of extrasolar planetary systems follows from the undetermined orbital parameters. The most successful observational technique - radial velocity measurements - allows us to estimate the planetary companion masses up to the line-of-sight inclination factor $\sin(i)$ (Cumming et al. 1999), and thus it gives only the lower mass limit $M \sin(i)$. Let us note that traditional Keplerian fits to radial velocity variations of the central star, used for determination of planetary orbits cannot, strictly speaking, be applied to systems with multiple planets. This problem was recently addressed by Laughlin & Chambers (2001). Using Keplerian elements as a representation of the planetary orbits, we are also not able to determine from observations the orientation of the lines of nodes and, as a consequence, the relative inclinations of orbits in systems with more than one planet.

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 semi-major 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 $\upsilon $ 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 $\upsilon $ 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 $\upsilon $ 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.

2 Overview of the method

A detailed description including examples of the application of the MEGNO method may be found in series of papers (CSGN hereafter) (Cincotta & Simó 2000; Cincotta & Giordano 2001; Cincotta & Núnez 2000). Here we give a simple overview of the method in order to make the paper self-contained.

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:

 \begin{displaymath} \frac{{\rm d}\phantom{t} }{{\rm d}t}\vec{x}=\vec{v}(\vec{x}), \qquad \vec{x}\in{\rm I\!R}^n. \end{displaymath} (1)

For a solution $\vec{\phi}=\vec{\phi}(t)$ of (1) the LCN is defined as (Benettin et al. 1980)

 \begin{displaymath} \lambda:= \lim_{t\rightarrow\infty}\lambda_1(t), \qquad \m... ...ac{\Vert \vec{\delta}(t) \Vert}{\Vert\vec{\delta}(t_0)\Vert}, \end{displaymath} (2)

where $\vec{\delta}(t)$ satisfies the variational equations

 \begin{displaymath} \frac{{\rm d}\phantom{t} }{{\rm d}t}\vec{\delta}= \vec{A}(t... ...) := \frac{\partial \vec{v}}{\partial \vec{x}}(\vec{\phi}(t)). \end{displaymath} (3)

Since the the LCN can be calculated analytically only in very special cases, in practice it is computed by numerical methods. These calculations usually require quite a long time due to a slow convergence of $\lambda_1(t)$ to $\lambda$, which is often much longer than the lifetime of the studied model (Cincotta & Simó 2000). Consequently, different new approaches allowing to explore the phase space on relatively short time-scales have been proposed, e.g. an elaborate method of frequency analysis (Dumas & Laskar 1993; Laskar 1993), the fast Lyapunov indicators (Froeschle et al. 1997; Froeschlé & Lega 2000; Froeschlé & Lega 2000), dynamical spectra, helicity and twist angles Voglis et al. (1999), just to mention a few.

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

 \begin{displaymath} \lambda:= \lim_{t\rightarrow\infty}\frac{1}{t}\int_0^t \fr... ...{\delta(s)}\,{\rm d}s= \Bigl<\frac{\dot\delta}{\delta}\Bigr>, \end{displaymath} (4)

where

 \begin{displaymath} \delta = \Vert\vec{\delta}\Vert, \qquad \dot\delta=\frac{{\... ...ec{\delta}} \cdot\vec{\delta}}{\Vert\vec{\delta}\Vert^2}\cdot \end{displaymath} (5)

The MEGNO is defined by

 \begin{displaymath} Y(t):= \frac{2}{t}\int_0^t \frac{\dot\delta(s)}{\delta(s)}s\,{\rm d}s, \end{displaymath} (6)

and it is used together with its mean value $\langle Y\rangle(t)$:

 \begin{displaymath} \langle Y\rangle(t):= \frac{1}{t}\int_0^t Y(s)\,{\rm d}s. \end{displaymath} (7)

Then, as it was shown by CSGN, we have:
a)
If $\vec{\phi}(t)$ is a regular solution with a linear divergence of nearby orbits then $\lim_{t\rightarrow\infty}\langle Y\rangle(t)=2$.
b)
If $\vec{\phi}(t)$ is chaotic then $\langle Y\rangle(t)\approx (\lambda/2)t$ as $t\rightarrow \infty$, where $\lambda$ is the LCN of the solution $\vec{\phi}(t)$.
Moreover, when $\langle Y\rangle(t)$ tends towards a fixed value different from 2 then it indicates that close trajectories diverge according to a certain power law. In fact the statements (a) and (b) can be merged; this is a very important feature of the MEGNO, distinguishing it from the LCN. The time evolution of $\langle Y\rangle$may be written in an uniform way for any kind of motion. The asymptotic behaviour of $\langle Y\rangle$ is given as

 \begin{displaymath}\langle Y\rangle(t) \simeq a t + d, \end{displaymath} (8)

where a=0, $d \simeq 2$ for stable quasiperiodic motion while $a=\lambda/2$, $d\simeq 0$ for irregular and chaotic motion.

This property can also be used to estimate efficiently the LCN, which can be recovered by the linear least square fit to $\langle Y\rangle(t)$. 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 (two-particle 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 $t\rightarrow \infty$:

\begin{displaymath}\frac{1}{t}Y(t) \simeq \frac{2}{t} \end{displaymath}

for a regular solution, and

\begin{displaymath}\frac{1}{t}Y(t) \simeq \lambda \end{displaymath}

for a chaotic solution. This means that in the case of regular motion Y(t)/t converges faster to 0 than $\lambda_1(t)$ (which converges as $\ln t/t$). Because $\langle Y\rangle(t)$ grows at the rate $\lambda/2$, we use this property for an alternative estimation of the LCN to the classical $\log\lambda(t)$vs. $\log t$ graphs. When the both quantities will be required, we will mark the Lyapunov Characteristic Number by $\lambda_{\max}$ on the plot's y-axis.

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 $\langle Y\rangle(t)$ can be found by solving this system of equations, complemented by the following two supplementary equations:

 \begin{displaymath}\frac{{\rm d}\phantom{t} }{{\rm d}t}y = \frac{\dot{\vec{\delt... ...ad \frac{{\rm d}\phantom{t} }{{\rm d}t}w = 2 \frac{y}{t}\cdot \end{displaymath} (9)

Then Y(t)=2y(t)/t and $\langle Y\rangle(t)=w(t)/t$.

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 two-particle 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 $\vec{\delta}$), 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 non-zero (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 108 yr a quasi-regular evolution, the LCN of this system actually takes a non-zero value, then giving a divergence rate of the order of $10^{-7.3}~\mbox{yr}^{-1}$.

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 103-104).

3 An object to study

For our dynamical studies, we choose the recently discovered $\upsilon $ And system (Butler et al. 1999). It attracts the attention of many researchers as the first confirmed extrasolar multi-planetary system consisting of three Jupiter-mass companions. The orbits of two outer bodies have relatively high eccentricities ( $e \simeq 0.2$-0.4), and the least massive innermost body is orbiting the central star very closely ( $P \simeq 4.6$ d) in an almost circular orbit. Thus it provides an excellent opportunity to study the planetary system stability. In spite of extensive computational efforts devoted to the study of its dynamical behaviour, it is difficult to determine whether the system is dynamically stable on the time span of the lifetime of the star, which is estimated to be 2-3 Gyr old. Some authors, e.g. Rivera & Lissauer (2000b) claim that there are many observationally allowed regions of parameter space in which the system can indeed survive on a stellar evolutionary timescale. There is still a fundamental question: which factors would allow for long-term stability of the system, in spite of strong mutual gravitational interaction between the two outer planets?

We refer to the paper by Stepinski et al. (2000), SMB paper hereafter, where we have found the initial conditions of the $\upsilon $ 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),

 

 
Table 1: Parameters and Keplerian geometrical elements of the $\upsilon $ Andr (1.3 $M_{\odot }$).
Planet Observatory mass [ $M_{{\rm J}}$] a [AU] P [d] e $\Omega$ [deg] $\omega$ [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


based on Lick and on AFOE observations. In some of our tests, a tiny relative inclination ${\simeq}0.001^{{\rm o}}$ is added in order to allow the system explore the third spatial dimension. There are other known estimations of the initial conditions of the system (see Chiang et al. 2001 for the latest estimates), but their inspection and comparison is out of the scope of the present paper.

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 long-term dynamical evolution of the system. By neglecting this planet it is possible to perform long-term integrations with much larger time step ( ${\simeq}10$ days) than it is permitted when the all the planets are taken into account ( ${\simeq}0.2$ days). Thus the dynamical timescale of the system increases by 2 orders of magnitude.

For that reasons we focus mostly on two-planet systems containing the central mass of the star ( $1.3~M_{\odot}$) 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 $\upsilon $ 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 N-body 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 $\upsilon $ 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).

4 Equations of motion

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

 \begin{displaymath}M_* \vec{r}_0 + \sum_{p=1}^{N} m_p \vec{r}_p = \vec{0}, \end{displaymath} (10)

where M* denotes the mass of the central body and mp, $p=1,\ldots,N$are masses of N planets. Setting

\begin{displaymath}\mu_p = \frac{m_p}{M_*}, \quad p=1,\ldots, N, \end{displaymath}

we obtain the position of the star M* in the barycenter reference frame:

\begin{displaymath}\vec{r}_0 = - \sum_{p=1}^{N} \mu_p \vec{r}_p, \end{displaymath}

and the vector $\vec{R}_i$ of the relative position of a planet Pwith respect to the star in that frame is

 \begin{displaymath} \vec{R}_{i} = \vec{r}_i - \vec{r}_0 = \sum_{p=1}^{N} (\delta_{pi}+\mu_p) \vec{r}_p. \end{displaymath} (11)

The Hamiltonian of the system is given by

\begin{eqnarray*}H = \frac{1}{2} M_* \vec{v}_0^2 - k^2 \sum_{p=1}^N \frac{M_... ...p,i\neq p}^N \frac{m_i m_p}{\vert\vec{r}_p-\vec{r}_i\vert}\cdot \end{eqnarray*}


The equations of motion of a planet i in the barycenter frame are then as follows

\begin{eqnarray*}\frac{{\rm d}\,{\vec r}_i}{{\rm d}\,t} &=& \vec{v}_i, \\ \fra... ...frac{\vec{r}_{i}-\vec{r}_p}{\vert\vec{r}_{i}-\vec{r}_p\vert^3}, \end{eqnarray*}


where $\vec{R}_i$ is given by (11) and $ i=1,\ldots,N$. These equations should be completed with the first order variational equations written in the general form of Eq. (3). The system of units is derived from the assumption k2=1: the length unit is [1 AU], the mass unit is [1 $M_{\odot }$] and time is expressed in days [1 d], with 1 Julian Year (denoted here by $\mbox{[yr]}$) equal to 365.25d. In these units the solar mass, derived from the third Keplerian law, is

\begin{displaymath}M_{\odot} = \frac{4 \pi^2}{\mbox{[yr]}^2} \mbox{[AU]}^3 \equiv 0.00029592\frac{\mbox{[AU]}^3 }{\mbox{[d]}^2 }\cdot \end{displaymath}

Because the geometrical orbital elements are usually given in the heliocentric reference frame, we need to change the heliocentric coordinates to the barycentric frame at the beginning of integration. Also the transformation of the results of computations to heliocentric orbital elements is required. This can be effectively handled by the following simple algorithm. We introduce 3Nvectors $\vec{R}=[\vec{R}_1^T,\ldots,\vec{R}_N^T]^T$ and $\vec{r}=[\vec{r}_1^T,\ldots,\vec{r}_N^T]^T$; then we rewrite Eq. (11) in the matrix form


\begin{displaymath}\vec{R} = \vec{A} \vec{r}, \quad \vec{r} = \vec{A}^{-1} \vec{R}, \end{displaymath}

where

\begin{displaymath}\vec{A}=\vec{E}_{3N} + \vec{M}, \end{displaymath}

$\vec{E}_{3N}$ is the $3N\times 3N$ identity matrix, and

\begin{displaymath}\vec{M}=\left[\begin{array}{ccc} \mu_1\vec{E}_3&\ldots& \mu_... ...mu_1\vec{E}_3&\ldots& \mu_N\vec{E}_3 \end{array}\right]\cdot \end{displaymath}

It is easy to demonstrate, the inverse of $\vec{A}$ is given by the formula:

\begin{displaymath}\vec{A}^{-1} = \vec{E}- \frac{1}{\det \vec{A}}\vec{M}, \end{displaymath}

where

\begin{displaymath}\det\vec{A} = 1 + \sum_{p=1}^N \mu_p. \end{displaymath}

Specifically, in the case of the 3-body problem (i.e. a star with 2 planets) we get

\begin{eqnarray*}\vec{r}_1 &=& \frac{1}{1+\mu_1+\mu_2} \left[ (1+\mu_2) \vec{R... ...mu_2} \left[ -\mu_1 \vec{R}_1 + (1 + \mu_1) \vec{R}_2 \right]. \end{eqnarray*}


For solving the equations of motion mixed with the variational equations, we used the code ODEX by Hairer & Wanner (1995), which incorporates the Bulirsh-Stoer-Gragg method. The absolute and relative error tolerances we set in the range of 10-14-10-16, depending on the machine architecture. After the reduction to the barycentric frame for N-planet system, the number of equations is 12N+2. For instance in the case of three planets we deal with 38 equations. The initial variation $\vec{\delta}_0$ was chosen as follows: with a random number generator we selected 6N positive components of $\vec{\delta}_0$, and a second sequence of 6N random numbers decided on the sign changes of the components. Finally the vector was normalized. At the beginning of calculations, utilizing the machine clock, we shifted the generator sequence, in order to have different $\vec{\delta}_0$ components in every program run. This is important because in systems having more than 3 degrees of freedom, the evolution of the MEGNO can depend on the choice of the initial tangent variational. This will be illustrated and discussed later.

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 104-105 orbital periods of the outermost body, depending on the given test. Integrations of the 3-planet $\upsilon $ 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).

5 Preliminary tests

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 $10^3T_{{\rm D}}-10^4T_{{\rm D}}$, where $T_{{\rm D}}$ is the characteristic period of the system. In the case of the $\upsilon $ And the period of the outermost planet is equal to about 1300d, and so the integration time should be of the order of $3.5 \times 10^3-3.5\times 10^4$ 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 $\vec{\delta}$ (see the definition of the MEGNO) grows at an exponential rate in the case of chaotic motion, it is necessary to renormalize $\vec{\delta}$ after a period of time, which we call here the renormalization step $\tau $.

The normalization of $\vec{\delta}$ follows from the definition, Eqs. (6), (5) and the fact that $\dot{\vec\delta}$ depends linearly on $\vec{\delta}$. The total integration time T consists of $N_{{\rm s}}$renormalization intervals $\tau $.

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 higher-dimensional systems, as for instance to the 4-body model describing the $\upsilon $ 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.


 

 
Table 2: Planetary mean orbits (J2000) (epoch = J2000 = 2000 January 1.5).
Planet mass [ $M_{{\rm J}}$] a [AU] e i [deg] $\Omega$ [deg] $\omega$ [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 $7.9006\times 10^{-06}$ 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.



  \begin{figure} \par\includegraphics[width=6.6cm,clip]{h2931f1.eps}\end{figure} Figure 1: The MEGNO calculated for the Lick Keplerian elements of the $\upsilon $ And, but for very small masses of the "planets" ${\simeq}10^{-5}~{M}_{{\rm J}}$. a) The upper panel shows the time evolution of Y(t), b) the mean $\langle Y\rangle(t)$ over the motion time of 104 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 $2\langle Y\rangle/t$ of the MEGNO. The integrator accuracy is set to $\epsilon _{{\rm rel}}=10^{-14}$, the renormalization step is $\tau =2\pi $ days, and the number of steps is $N_{{\rm s}}=5\times 10^6$.
Open with DEXTER

We followed the evolution of the 2-planet system by taking the Lick geometrical elements and lowering the planetary companion masses to very small values ${\simeq}10^{-5}~{M}_{{\rm J}}$. In such a model, one expects quasi-Keplerian orbits of the companions since the mutual perturbations are almost negligible. In this case indeed, the MEGNO time evolution (see Fig. 1a), gives the signature of a quasiperiodic regime of the system (see Fig. 1a), with the mean value $\langle Y\rangle$ converging very fast to the value of 2. Let us note that we performed the integration over motion time much longer than theoretically required ( ${\simeq}5 \times 10^4$ of the characteristic periods), in order to examine the character of convergence of the MEGNO.

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.

  \begin{figure} \par\includegraphics[width=10cm,clip]{h2931f2.eps}\end{figure} Figure 2: The MEGNO calculated for the 3-planet 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 $\vec{\delta}_0$. The two middle panels show the resulting plots of the mean $\langle Y\rangle$. The two bottom panels illustrate the evolution of the eccentricities and inclinations of the orbits. Integration was performed over ${\simeq }7\times 10^5$ yr. Relative errors of the total energy the angular momentum do not exceed 10-12 in the runs.
Open with DEXTER

Because the period of the outermost body (Saturn) is 29.46 yr, by the "rule of thumb'' choice the integration time should be of the order of $3\times 10^4$- $3\times 10^{5}$ years. The behaviour of the MEGNO in this case is characteristic for quasiperiodic regular orbits. Figure 2 illustrates two from a few tests that we computed with this model, choosing the initial tangent vector $\vec{\delta}_0$ at random. The bottom panels demonstrate that, depending on the choice of $\vec{\delta}_0$, the mean $\langle Y\rangle(t)$ can approach the limit 2 from above or from below. This effect was not observed in 2D and 3D-systems, where these two different time evolution patterns were used to distinguish whether the investigated solution is close to an unstable or a stable periodic solution.

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, long-term integrations (Sussman & Wisdom 1988). In the Jupiter-Saturn-Pluto system, the roughly-estimated integration time should be of the order of $2.5\times 10^5$- $2.5\times 10^6$ yr. In order to have a clear answer we continued the integration up to 107 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 round-off errors.

Our results are shown in Fig. 3.

  \begin{figure} \par\includegraphics[width=10cm,clip]{h2931f3.eps}\end{figure} Figure 3: The MEGNO calculated for the 3-planet system of Jupiter, Saturn and Pluto. Panels  a-c) show the time evolution of the MEGNO computed with the DEC/Alpha and the Intel machines with integration accuracy $\epsilon _{{\rm rel}}=10^{-14},\tau =10$ yr a), on the Alpha with $\epsilon _{{\rm rel}}=10^{-14}, \tau =2$ yr b), and on the Intel with $\tau =2$ yr, $\epsilon _{{\rm rel}}=10^{-16}$ (thin line) as well as $\epsilon _{{\rm rel}}=10^{-14}$ (thick line), panel  c). Relative errors of the energy and the angular momentum are of the order of 10-9and 10-10, depending on the test. Panel  d) shows the maximal Lyapunov exponent computed over time interval of 1 Gyr.
Open with DEXTER

Tests on the Alpha and Intel machines (shown in panel a), continued with $\tau=10~\mbox{yr}$ and $\epsilon _{{\rm rel}}=10^{-14}$ give initially fast convergence of the indicator to 2 but after $2\times 10^6$ yr the MEGNO seems to be slowly divergent from this value. Suspecting that the instability appears because of a too long renormalization step $\tau $ and accumulation of numerical errors in the variationals, we repeated the calculation with $\tau =2$ yr on the Alpha machine (panel b), and the Intel (panel c), which indeed gave much better convergence. In these cases we obtained a behaviour characteristic for quasi-regular orbits. Differences between the graphs are due to different $\vec{\delta}_0$. In the cases when we observe a slow divergence of the MEGNO, we cannot definitely conclude that the motion is regular. In all the calculations we get the LCN converging to zero, but this cannot be used as a definitive argument - the integration time 10 Myr might be too short to get the conclusion (as it is suggested by the investigations of the whole outer Solar system).
  \begin{figure} \par\includegraphics[width=6.37cm,clip]{h2931f4.eps}\end{figure} 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 $\langle Y\rangle(t)$, d) variations of the eccentricity of planet C in the two models, e) the same quantity for the outer planet D.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=6.5cm,clip]{h2931f5.eps}\end{figure} 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 $\langle Y\rangle(t)$ for the 2-planet system, with different choices of the initial tangent vector $\vec{\delta}_0$ (at random), c) Y(t) for the system of 3 planets (dotted curve) with the graph of $2\langle Y\rangle$ (thick, smooth curve) - the factor 2 was added due to the averaging of Y (see papers by Cincotta and co-authors 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 $2\langle Y\rangle/t$.
Open with DEXTER

For this reason we computed the LCN (using an independent program) over much longer period of time 1 Gyr (panel d) - in this test we found that the LCN is converging to zero too, what indicates that the motion is regular, indeed. Anyway, this example shows that the accumulation of numerical errors (due to a too large renormalization interval or the integrator errors) may result in an output that can be wrongly interpreted, for example as the presence of weak chaos (see the thick graph in panel b). This problem is however independent on what kind of indicator we use.

Let us note that the system of Jupiter-Saturn-Pluto 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 N-body 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 Sun-Jupiter-Saturn 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.

6 Test of the nominal parameters of the $\upsilon $ And

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 $43\,000$ yr (corresponding to about 104 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 $\langle Y\rangle(t)$ 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 2-planet system has the MEGNO signature of the quasiperiodic regime (Figs. 5a,b). We tested it by repeating the calculations with the variational $\vec{\delta}_0$ chosen at random, and always found behaviour characteristic for regular orbits.

For the 3-planet model we found clear signs of chaotic behaviour. In this example the system remains close to the quasiperiodic regime during less than 104 yr, as $\langle Y\rangle(t)$ goes very close to 2. However, after this time a large jump of $\langle Y\rangle(t)$ 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 $2\langle Y\rangle(t)/t$ - 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.

  \begin{figure} \par\includegraphics[width=10.6cm,clip]{h2931f6.eps}\end{figure} 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 $3\times 10^4$ 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


  \begin{figure} \par\includegraphics[width=10cm,clip]{h2931f7.eps}\end{figure} Figure 7: a) The MEGNO indicator $\langle Y\rangle(T)$ computed for the nominal geometrical Lick elements ( $\epsilon _{{\rm rel}}=10^{-14}$, $\tau =4\pi $d, and $N_{{\rm s}}=10^6$). The levels of the nominal (minimal) mass of planet C and the higher value ${{M}}_{{\rm C}}=8.52~{{M}}_{{\rm J}}$ are indicated by the dashed lines. b) the MEGNO computed with a better accuracy, namely $\epsilon _{{\rm rel}}=10^{-16}$, the same renormalization step $\tau $, and $N=2\times 10^6$ (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


  \begin{figure} \par\includegraphics[width=10cm,clip]{h2931f8.eps}\end{figure} Figure 8: The MEGNO $\langle Y\rangle$ computed for the nominal, geometrical Lick elements, with another mass of planet C, namely $8.52~ {{M}}_{{\rm J}}$ (it is marked in Fig. 7a). Panel  a) shows four different plots of the MEGNO obtained with different choices of the initial tangent vector $\vec{\delta}_0$. 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 $\langle Y\rangle(t)$ 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 $2\langle Y\rangle/t$).
Open with DEXTER

Also the estimation of the LCN, obtained by measuring the slope of $\langle Y\rangle(t)$, gives the consistent value of $\lambda_{\max} \simeq 8.6\times 10^{-4}~\mbox{yr}^{-1}$. These relationships between the LCN and the MEGNO are supported by the theory.

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 Cowell-Stormer 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 105 yr with time step 0.004d, following the code author's recommendation to use the integration time step of $T_{{\rm S}}/1000$, where $T_{{\rm S}}$ 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 $4\times 10^4$ 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 103 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 2-planet configurations, one describes only necessary but not sufficient constraints on the system stability.

7 A look at the global dynamics

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 2-planet system the relative inclination of orbits $i_{\rm r}$ is related to the longitudes of nodes $\Omega_{{\rm C}},\Omega_{{\rm D}}$ by the following formulae (used also by SMB):

 \begin{displaymath}\cos (i_{\rm r}) = \cos(i_{{\rm C}}) \!\cos(i_{{\rm D}}) + \... ...\sin(i_{{\rm D}}) \!\cos (\Omega_{{\rm D}}-\Omega_{{\rm C}}), \end{displaymath} (12)

where $i_{{\rm C}}, i_{{\rm D}}$ are the inclinations of the planetary orbits to the line of sight. In this section we apply the MEGNO for determining those parameter ranges, which would correspond to regular and to chaotic motions of the $\upsilon $ And. As a source of the initial conditions we use again the SMB paper.

In a first simple examination of the space of the undetermined orbital parameters, we may assume that one of the orbits has $\sin i=1$, i.e. $i=90^{{\rm o}}$ (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 $0~M_{{\rm J}}$ up to the value of $10~M_{{\rm J}}$, (so the mass hierarchy of the planets C and D can be switched in some cases), that corresponds to relative inclinations $i_{\rm r}$ in the range ${\simeq}0^{{\rm o}}$- $80^{{\rm o}}$. The first run we did with $\tau =4\pi $ days and $N_{{\rm s}}=10^6$ steps with integrator accuracy $\epsilon _{{\rm rel}}=10^{-14}$. That gave us $T\simeq 3.4\times 10^4$ yr. The resulting graph of $\langle Y\rangle(T)$ as a function of mass $M_{{\rm C}}$ is shown in panel (a) of Fig. 7. It is easy to see that the minimal mass of C, namely $M_{{\rm C}}^{{\min}}\approx 1.9~ M_{{\rm J}}$, 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 ${{M}}_{{\rm C}}=8.52~{{M}}_{{\rm J}}$ (see a narrow valley in panel 7a). We run a few computations of $\langle Y\rangle(t)$ with different values of $\vec{\delta}_0$. In some cases we obtained apparently very good and fast convergence, but one of the runs $\langle Y\rangle(t)$ 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 $20\%$ 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 $\vec{\delta}_0$: we obtain an appearance of different kinds of $\langle Y\rangle(t)$ evolution when the time span is too short.

  \begin{figure} \par\includegraphics[width=9cm,clip]{h2931f9.eps}\end{figure} 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 $2\langle Y\rangle/t$ (continuous lines) and by the direct variational method (dotted lines).
Open with DEXTER


  \begin{figure} \par\includegraphics[width=9cm,clip]{h2931f10.eps}\end{figure} 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 $\langle Y\rangle(t)=(\lambda/2) t + d$.
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 $M_{{\rm C}}$ between 0 and 10  $M_{{\rm J}}$ with a better accuracy ( $\epsilon _{{\rm rel}}=10^{-16}$) and doubled $N_{{\rm s}}=2\times 10^6$ 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 $M_{{\rm C}}$ between $8~M_{{\rm J}}$ and $9~M_{{\rm J}}$, the MEGNO changed its quantitative character.

A close-up 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 $2\langle Y\rangle/T$.

  \begin{figure} \par\includegraphics[width=11cm,clip]{h2931f11.eps}\end{figure} Figure 11: Behaviour of the MEGNO for the initial conditions discussed in the SMB paper. On the left hand side: $\sin(i_{{\rm C}})=\sin(i_{{\rm D}})=0.5$ (i.e. the nominal Lick masses are doubled), $\Omega _{{\rm D}}-\Omega _{{\rm C}}=15^{{\rm o}}$, relative inclination $i_{\rm r}=7.5^{{\rm o}}$. On the right hand side: $\Omega _{{\rm D}}-\Omega _{{\rm C}}=90^{{\rm o}}$( $i_{\rm r}=41{\hbox{$.\!\!^\circ$ }}4$). 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 $i_{\rm r}$; it could be caused by the Kozai cycles (Kozai 1962). ( $\epsilon _{{\rm rel}}=10^{-14}$, $\tau =4\pi $d, and $N_{{\rm s}}=2\times 10^6$).
Open with DEXTER

The similarity of the plots of $\langle Y\rangle(T)$(Fig. 7b) and LCN (Fig. 7d) as functions of $M_{{\rm C}}$ is obvious.

This experiment shows the possibility of changing the mass hierarchy in the $\upsilon $ And. Fixing $M_{{\rm D}}$ at its minimal value, we can select $M_{{\rm C}}$ up to almost 2 times of $M_{{\rm D}}$, 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, $i_{{\rm C}} = i_{{\rm D}}$, for both outer planets and therefore did not switch their mass hierarchy. However, let us note, that if one simply interchanges masses $M_{{\rm C}}$ and $M_{{\rm D}}$ leaving all other initial orbital parameters unchanged, the Lick 2-planet $\upsilon $ And system becomes unstable unless additional perturbing forces are taken into account (Kiseleva-Eggleton & Bois 2001). Their results confirm that not only the mass hierarchy (or mass ratio of $M_{{\rm C}}$ and $M_{{\rm D}}$) 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 $M_{{\rm C}}$, 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 $M_{{\rm C}}=2.74~M_{{\rm J}}$ (in the middle of the second pick in Fig. 7a), two values in the wide regular zones between $M_{{\rm C}}=3~M_{{\rm J}}$ and $M_{{\rm C}}=7~M_{{\rm J}}$, and a large mass $M_{{\rm C}}=9~M_{{\rm J}}$ 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 $\lambda_{\max}$ both by the direct variational method (plots labeled by the LCN) and from $2\langle Y\rangle(t)/t$. 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 $\lambda_{\max}$ provided by the MEGNO is faster than in the case of the LCN - it is the effect predicted by the theory.

Another one-dimensional test of the parameter space available for the $\upsilon $ 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 $25^{{\rm o}}$ to $30^{{\rm o}}$ related to $M_{{\rm C}}\simeq4.5~M_{{\rm J}}$ and $M_{{\rm D}}\simeq9~M_{{\rm J}}$, 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 $\omega$ 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 $\sin(i_{{\rm C}})=\sin(i_{{\rm D}})=0.5,\ \Omega_{{\rm D}}-\Omega_{{\rm C}}=15^{{\rm o}}$, corresponding to relative inclination $i_{\rm r}=7.5^{{\rm o}}$, is regular while that with $\Omega _{{\rm D}}-\Omega _{{\rm C}}=90^{{\rm o}}$( $i_{\rm r}=41{\hbox{$.\!\!^\circ$ }}4$) 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 106 yr integration. The irregular behaviour is detected very quickly, only after $2.5\times 10^4$ yr.

SMB examined the stability of the $\upsilon $ And system for the eccentricity ranges of planet C induced by different initial values of $\sin i$ (and so, different masses $M_{{\rm C}}$ and $M_{{\rm D}}$) as well as by different relative inclinations $i_{\rm r}$. 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 $\sin i<0.23$, 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 $3\times 10^4$ yr for one initial condition compared to 106 yr integrations by SMB.

Our detailed scan over the undefined parameters of the $\upsilon $ 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. 12-14. Figures 12, 13 present the results for the Lick data, and Fig. 14 corresponds to the same experiments performed on the AFOE elements.


  \begin{figure} \par\includegraphics[width=7.5cm,clip]{h2931f12.eps}\end{figure} Figure 12: The MEGNO calculated for the Lick elements of planets C and D in three-dimensional form. The i-abscissa is the inclination of the line of sight, so that both masses are $\propto1/\sin(i)$. 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 $50\times 50$ points. ( $\epsilon _{{\rm rel}}=10^{-14}$, $\tau =2\pi $d, and $N_{{\rm s}}=10^6$).
Open with DEXTER

The system remains in quasiperiodic regime if the relative inclinations of the orbits are bounded by raw limits of ${\simeq}60^{{\rm o}}$ for $\sin i$ close to 1 and $20^{{\rm o}}$- $30^{{\rm o}}$ for $\sin i \simeq 1/2$ for the Lick orbital parameters, and ${\simeq}30^{{\rm o}}$ and ${\simeq}20^{{\rm o}}$ for the AFOE parameters, respectively. For the Lick parameters, the MEGNO gives almost the same stability limit $\sin i \simeq 0.23$, as it was found by SMB, but for the AFOE data we obtain more restrictive stability limit, namely $\sin i \simeq 0.65$, instead of ${\simeq}0.33$, determined by the authors.

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 3-planet 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 $\sin i$ 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 106 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 $\omega$ has any effect on the distinction between the chaotic and regular evolution of the 2-planet system, defined by the nominal Lick parameters. From the results illustrated by Fig. 15, we may conclude that errors of the order of $20^{{\rm o}}$ in $\omega$ mostly would not change the qualitative state of the system. Moreover there are three zones of instability represented by peaks in $\langle Y\rangle$. They correspond to secular resonances in the pericenter motion. For detailed studies of the effects of the secular resonance on the stability of the $\upsilon $ 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 $\sin(i)=0.4$ 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 ( $\omega _{{\rm C}}=\omega _{{\rm D}}$) 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 $1^{{\rm o}}\times1^{{\rm o}}$ grid of the parameters, which were examined, the differences from the value of 2 are less than 0.005.

8 Conclusions

The MEGNO method seems to be promising and an efficient tool for studying the long-time 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:

Considering the application of the method to the $\upsilon $ And case, we claim that our work demonstrates the usefulness of the MEGNO in planetary dynamics. The experiments that we have performed on the 3-body and 4-body models, i.e. on multidimensional dynamical systems, confirm numerically its useful features predicted by the theory.
  \begin{figure} \par\includegraphics[width=9cm,clip]{h2931f14.eps}\end{figure} Figure 14: The MEGNO plotted in the same way as in Fig. 12, but for the AFOE parameters. Filled black points mean $\vert(\langle Y\rangle-2\vert \leq 0.001$, open circles $\vert\langle Y\rangle-2\vert \leq 0.005$, and dots $\vert\langle Y\rangle-2\vert<0.01$.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=9cm,clip]{h2931f15.eps}\end{figure} 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 dash-dotted line. On the right panel, the same data are represented in the form of a contour plot: filled black points indicate $\vert\langle Y\rangle-2\vert\leq 0.001$, open circles $\vert\langle Y\rangle-2\vert \leq 0.005$, and dots $\vert\langle Y\rangle-2\vert<0.01$. The nominal arguments of periastron are marked by ( $\omega _{{\rm C}}^0,\omega _{{\rm D}}^0$) on the left panel and by a cross on the right panel. The exact secular resonance $\omega _{{\rm C}}=\omega _{{\rm D}}$ is marked by a thick line on the right panel.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=9cm,clip]{h2931f16.eps}\end{figure} 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 $\sin(i)=0.4$. The nominal data are marked on the left panel by a dash-dotted line. On the right panel, the same data are illustrated in the form of a contour plot: filled black points indicate $\vert\langle Y\rangle-2\vert\leq 0.001$, open circles $\vert\langle Y\rangle-2\vert \leq 0.005$, and dots $\vert\langle Y\rangle-2\vert<0.01$. The nominal arguments of periastron are marked by ( $\omega _{{\rm C}}^0,\omega _{{\rm D}}^0$) on the left panel and by a cross on the right panel. The exact secular resonance $\omega _{{\rm C}}=\omega _{{\rm D}}$ is marked by a thick line on the right panel.
Open with DEXTER

However, as any numerical method, the MEGNO needs some tuning and experience in using it.

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 $\upsilon $ 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#938-COS272).

References

 

Copyright ESO 2001