A&A 398, 1151-1161 (2003)
DOI: 10.1051/0004-6361:20021713

Stability of the HD 12661 Planetary System

K. Gozdziewski

Torun Centre for Astronomy, N. Copernicus University, Poland

Received 26 August 2002 / Accepted 14 November 2002

Abstract
In this paper we study the stability of the HD 12661 planetary system in the framework of the N-body problem. Using the initial conditions found and announced by the California & Carnegie Planet Search Team, (http://exoplanets.org/almanacframe.html), we estimate the dynamical limits on orbital parameters that provide stable (quasi-periodic) motions of the system. We investigate the orbital stability by combining the MEGNO indicator analysis with short-term integrations of the orbital dynamics. The MEGNO technique, invented by Cincotta & Simó (2000), makes it possible to distinguish efficiently between chaotic and regular dynamics of a conservative dynamical system. The orbital evolution, derived simultaneously with MEGNO, helps to identify sources of instability. The nominal initial condition leads to a chaotic solution, with a Lyapunov time $\simeq $1300 yr. In spite of this, the system motion seems to be bounded. This was examined directly, by long-term, 1 Gyr integrations. During this time, the eccentricities vary in the range (0.1, 0.4). The system is locked in apsidal resonance with the critical argument librating about $180^{\circ }$, with a full amplitude varying typically between $40^{\circ}$ and $180^{\circ }$. Using MEGNO, we found that the HD 12661 system evolves on a border of the 11:2 mean motion resonance. This resonance is stable and results in a quasi-periodic evolution of the system. From the viewpoint of global dynamics, the crucial factor for system stability is the presence of the apsidal resonance. We detected this resonance in a wide neighborhood of the initial condition in the space of orbital parameters of the system, and in wide ranges of relative inclination and masses of the planets. The center of libration can be $180^{\circ }$ (as in the nominal system) or $0^{\circ }$. The regime depends on the initial values of the apsidal longitudes. Statistically, the system prefers almost exclusively one of these two resonance regimes. The HD 12661 system gives a very evident example of the dynamical role of secular resonances, and their influence on the stability of exosystems containing Jupiter-like planets. Data derived by numerical experiments are compared with the results of Laplace-Lagrange secular theory. The analytical theory gives a crude approximation of the secular dynamics, because the eccentricities and masses of the planets are large, and the nominal system is near the 11:2 mean motion resonance.

Key words: celestial mechanics, stellar dynamics - methods: numerical, N-body simulations - planetary systems - stars: individual: HD 12661


1 Introduction

The search for extrasolar planetary systems conducted in the past several years by the California & Carnegie Planet Search Team has recently revealed a number of new multi-planetary systems (Butler et al. 2002). This work is devoted to a preliminary analysis of the orbital stability of the HD 12661 system, and its dependence on the orbital parameters that are unconstrained by radial velocity observations. In our numerical experiments, we use the 2-Keplerian fit published by the California & Carnegie Planet Search Team on their web site[*] and given in Table 1. Formally, we consider the orbital parameters as osculating elements at the epoch of the periastron passage of the outer planet. Because the orbital solution is uncertain, we try to recover some qualitative dynamical properties of the HD 12661 planetary system.


 

 
Table 1: Orbital parameters of the HD 12661 system, adopted from data published in http://exoplanets.org/almanacframe.html (data revised on 8 July 2002). The mass of the central star is equal to $1.07~{M}_{\odot}$.
Orbital parameter Planet b Planet c
${m}_2 \sin i$ [M$_{{\mbox{\scriptsize J}}}$] ...................... 2.3 1.56
a [AU] .................................. 0.82 2.56
P [d] ..................................... 263.3 1444.5
e ............................................ 0.35 0.20
$\omega$ [deg] ................................ 292.6 147.0
$T_{{\mbox{\scriptsize p}}}$ [JD-2 450 000] .......... 9943.7 9673.9
$M(T_{{\mbox{\scriptsize c}}})$ [deg] ....................... -8.89 0.00


To explore the orbital parameters' space of the HD 12661 system, we use the Mean Exponential Growth factor of Nearby Orbits (MEGNO), the fast indicator invented by Cincotta & Simó (2000). This technique seems to be well designed for studying the planetary dynamics in the framework of the N-body problem. The basic idea of our approach relies on rapid determination whether an investigated initial condition leads to quasi-periodic or irregular (chaotic) motion of the planetary system. In Gozdziewski et al. (2001) we explain our motivations and technical aspects of the calculations.

In this work, we generally follow Gozdziewski (2002). MEGNO is utilized to construct one- and two-dimensional sections of the orbital parameters' space in the neighborhood of the nominal initial condition. Such maps make it possible to determine dynamical bounds on the orbital parameters that are unconstrained by the fits, for instance, the inclination of the system, the companions' masses, and the relative inclination of orbits. Basically, we are not interested in studying the orbital stability with long-term integrations. Instead, we search for regular motions that constitute a skeleton of the phase space filled up with quasi-periodic, stable evolutions. Simultaneously, we monitor the evolution of orbital elements. The timescale of the integrations is about $N_{{\mbox{\scriptsize P}}}=10^4$orbital periods of the outer companion. This time is long enough to derive reliable estimates of both MEGNO as well as the Lyapunov Characteristic Number (LCN, $\lambda{{\mbox{\scriptsize max}}}$) and to obtain valuable data that help to identify the main sources of instabilities. The computational efficiency of the method is a crucial factor that makes it possible to examine large sets of initial conditions.

Most of the MEGNO integrations in this paper have been driven by a Bulirsh-Stoer integrator. We used the ODEX code (Hairer & Wanner 1995). The relative and absolute accuracies of the integrator have been set to 10-14 and $5\times 10^{-16}$, respectively. The position of the nominal condition is marked in contour plots by the intersection of two thin lines. The stable areas are marked in MEGNO maps by light-gray (with value $\simeq $2), and chaotic zones are dark-shaded. (Colour versions of these maps are published in the electronic edition of the Journal.)

2 Analysis of the initial condition

Using the initial condition given in Table 1 to examine the evolution of the HD 12661 planetary system, the system immediately reveals its inquiring nature. The results of a 300 Myr integration, performed with the help of the ODEX integrator, are shown in Figs. 1a-f. Apparently, the system behaves chaotically, and this can be clearly seen on graphs representing the time-evolution of the semi-major axes and eccentricities of the companions. This observation has a strong theoretical confirmation. Using the method announced in (Gozdziewski 2003), we calculated the orbital evolution of the system over 10 Myr, and we found simultaneously its MEGNO and LCN signatures. The results are illustrated in Figs. 1g, h. MEGNO grows linearly with time, at the rate $3.8 \times 10^{-4}~\mbox{yr}^{-1}$, obtained by the linear fit to the theoretical asymptotic behaviour of MEGNO, described by $\langle Y\rangle(t)= (\lambda_{\mbox{\scriptsize max}}/2) t + d$ (Cincotta & Simó 2000). This fit provides a very good estimate of LCN - using the direct variational algorithm, we obtained $\lambda_{{\mbox{\scriptsize max}}} = 7.5 \times 10^{-4}~\mbox{yr}^{-1}$ (see Fig. 1h). Thus, the Lyapunov time scale of the HD 12661 system is very short, $\simeq $1300 yr. This means that, integrating the equations of motion, one obtains the orbital solution with 100% undeterminancy after $\simeq $40 000 yr, assuming that calculations are performed with a typical floating precision of 14-15 significant digits. In spite of the strong chaotic behaviour, the system motion seems to be bounded. A dynamical mechanism, which certainly helps to maintain the stability, is the secular apsidal resonance (SAR). This is clearly illustrated in Fig. 1e. The apsidal lines remain antialigned, and the critical argument of the SAR, $\theta_3=\varpi_{\mbox{\scriptsize b}}-\varpi_{\mbox{\scriptsize c}}$, librates with varying full amplitude $\simeq $40- $120^{\circ}$ about $180^{\circ }$. For comparison, Fig. 1f shows the evolution of $\theta _3$ when the initial mean anomaly of the inner planet is increased by $20^{\circ }$. In this case, we obtain perfectly regular oscillations of the critical argument. (This choice will be clarified later.) The long term evolution of the HD 12661 system has been further examined through a number of numerical integrations over 1 Gyr. The results of these experiments are illustrated in Figs. 1i, j, showing variations of the critical argument of the SAR. Panel 1i is for the Bulirsh-Stoer integrator and panel 1j is for the RMVS integrator (Levison & Duncan 1994), used with a time-step of 8 d. Actually, we found that $\theta _3$ can circulate, and the chaotic nature of the nominal HD 12661 system is very evident.

Possibly, in the case of the nominal HD 12661 system, we have convincing evidence of the stability in the astronomical sense (Lissauer 1999). This notion of stability implies that the system will remain bound, without ejections or collisions, for a possibly long but finite period of time, and this result is robust against sufficiently small perturbations. Although the motion of the HD 12661 system is chaotic, its planets avoid close encounters due to the protective dynamical mechanism of the secular resonance, and the motion of planets can be bounded over a very long period of time.


  \begin{figure} \par\includegraphics[width=15.5cm,clip]{h3945f1.eps}\end{figure} Figure 1: Orbital evolution of the HD 12661 system for the nominal initial condition given in Table 1. Panels a), b) show changes of the semi-major axes. Panels c), d) are for the eccentricities. Panel e) shows the critical argument of the apsidal resonance $\theta_3=\varpi_{\mbox{\scriptsize b}}-\varpi_{\mbox{\scriptsize c}}$. Panel f) shows the evolution of the critical argument when the mean anomaly of the outer planet is increased by $20^{\circ }$. Panel g) is for the MEGNO $\langle Y\rangle$. Panel  h) is for the Lyapunov Characteristic Number (LCN) calculated with the help of the direct variational algorithm. Panels i) and j) illustrate variations of $\theta _3$ over 1 Gyr derived with the Bulirsh-Stoer integrator, and with the RMVS integrator by Levison & Duncan (1994) from their SWIFT package, respectively. Time is in Myr's.
Open with DEXTER

To examine where possible stable systems are located in the orbital parameters' space, we computed a number of MEGNO sections of this space. A one-dimensional scan of MEGNO, over the semi-major axis of the inner planet, computed for the fixed, nominal value $a_{{\mbox{\scriptsize c}}}$, is shown in Fig. 2. We plotted both MEGNO (the upper graph), as well as the maximal values of the eccentricities attained during the integration time (equal to 104 periods of the outer companion). The later values are shown in the two lower graphs. The MEGNO plot is accompanied by small triangles that indicate the nominal positions of mean motion resonances (MMR). These resonances are labeled by i:j, where $j n_{{\mbox{\scriptsize b}}} - i n_{{\mbox{\scriptsize c}}} \simeq 0$, and $n_{{\mbox{\scriptsize b,c}}}$ denote the mean motions of the planets, respectively. This experiment reveals that the nominal initial condition is located in a chaotic zone, which is very close to the 11:2 MMR. The resonance can be characterized by librations of the critical argument 211360, with a full amplitude of $\simeq $ $180^{\circ }$ (see Fig. 3a). It is accompanied by the SAR: its critical argument librates with a full amplitude of $\simeq $ $80^{\circ}$ (Fig. 3a'). In this instance, MEGNO is close to 2, indicating a stable, quasi-periodic system. In the vicinity of this stable MMR we see two symmetric peaks of MEGNO, and the nominal system is located at one of these peaks[*]. They can be identified with the components of the 11:2 MMR separatrix.

  \begin{figure} \par\includegraphics[width=16.7cm,clip]{h3945f2.eps}\end{figure} Figure 2: A scan of MEGNO for the nominal initial condition of the HD 12661 system, over the semi-major axis of the inner planet (in AU). The upper graph is for MEGNO. Two graphs at the bottom (in the gray shaded area) are for the maximal values of the eccentricity of the inner planet (thick line) and the outer planet (thin line), attained during the integration time equal to 104 periods of the outer planet. The resolution of the plots is $8 \times 10^{-5}$ AU. Triangles at the level of 3.5 indicate the nominal positions of mean motion resonances, and they are labeled. The nominal value of $a_{{\mbox{\scriptsize b}}}$ is marked with a thin vertical line at 0.82 AU. Regions of mean motion resonances analyzed in Fig. 3 are shaded in light-gray.
Open with DEXTER

The nominal HD 12661 system lies on the border of a regular regime of motion. This fact can provide another argument for the long-term, bounded behaviour of the system, although without its confirmation by direct integrations, for now it should be treated as speculative. Because the motion is close to a separatrix (or takes place within it), it can be influenced by the stickiness phenomenon. This effect was studied, e.g., by Murison et al. (1994), and recently, by Tsiganis et al. (2000, 2002), for the models of asteroid motion in the Solar system. The former authors found that the long-term evolution of chaotic orbits initiated in the vicinity of some high-order mean motion resonances with Jupiter is possible if there exist no resonant periodic orbits in the framework of the elliptic-three body problem. They showed that, in more complicated dynamical models, the long-term evolution of chaotic orbits may also be governed by secular resonances. The "stable-chaotic'' orbits found for some asteroids behave in a way characteristic for trajectories being "sticky'' around an invariant manifold of lower-than-full dimensionality in phase space. The dynamical effects observed in the nominal HD 12661 system is reminiscent of the models analyzed by Tsiganis et al. (2000, 2002). By analogy, the hypothesis of the "stable-chaotic'' motion of the HD 12661 system can be supported by the presence of the stickiness effect, which takes place in addition to the secular resonance.

In the scan shown in Fig. 2, we can identify other MMR's which are characterized by structures similar to that found in the case of the 11:2 MMR, i.e., a stable zone bounded by two symmetric spikes, indicating the separatrix regions. This is clearly seen for the 8:1, 7:1, 6:1, 13:2, and 9:2 MMR. (Note that qualitatively similar structures, indicating the presence of resonances, are observed in low-dimensional Hamiltonian systems (Cincotta & Simó 2000).) The identification is illustrated by a few examples shown in Fig. 3. This figure illustrates critical arguments of the MMR's. In all these cases the apsidal resonance is also present. The stabilizing influence of these resonances can be clearly recognized in the evolution of the eccentricities (see the two graphs at the bottom of Fig. 2). The maximal eccentricity of the outer, smaller planet, is significantly lowered in the stable zones of the MMR's. The MEGNO scan reveals a very interesting structure of the 5:1 MMR. In this instance, we found an unstable zone in the middle of the resonance, together with surrounding it, two unstable "wings''. The SAR is also present. The scan makes it possible to estimate the width of the mean motion resonance as $\simeq $0.025 AU. In Fig. 2, a large number of other spikes of MEGNO can be identified with high order resonances.


  \begin{figure} \par\includegraphics[width=16cm,clip]{h3945f3.eps}\end{figure} Figure 3: Evolution of the critical arguments of mean motion resonances (the left panels) and the critical argument of the apsidal resonance (the right panels), calculated for a few initial conditions localized within the resonances identified in Fig. 2. Panels a), a') are for the stable 11:2 MMR. Panels b), b') and c), c') are for the middle unstable and the left stable region of the 5:1 MMR, respectively. Note the chaotic changes of $\theta _3$ in panel b'). Panels d), d'), and e), e') are for the stable center and for the unstable border of the 6:1 MMR, respectively.
Open with DEXTER

3 Stability as a function of eccentricities and semi-major axes

In the next experiment, we calculated 2D-MEGNO maps in the plane of semi-major axes. The result is shown in the left panel of Fig. 4. In this plot we can easily recognize MMR's already identified in the one-dimensional scan. The middle panel of Fig. 4 helps us to see the SAR and to estimate a semi-amplitude of its critical argument $\theta _3$. It has been derived by the following simple method. In the program computing MEGNO, we evaluated maximal value $\theta_3^{{\mbox{\scriptsize max}}}$ of $\theta _3$ after every step of renormalization of the variational equations. (The time step was equal to the orbital period of the outer planet.) The maximal value of the critical argument was taken relative to the center of libration $0^{{\rm o}}$, or $180^{\circ }$. To avoid the effects of a possible transition into the apsidal resonance, the determination of $\theta_3^{{\mbox{\scriptsize max}}}$ was started after the first half of the integration period. Finally, if $\theta_3^{{\mbox{\scriptsize max}}}<180^{\circ}$, then we treated this value as a semi-amplitude of apsidal librations. Although the period of integrations is relatively short, such an approximation of the semi-amplitude gives much insight into the global dynamics of the system. One should be aware that the apsidal resonance, in the sense considered here, can appear in chaotic zones temporarily (for an illustration see Fig. 3b'). The result is also false if the system is trapped into the secular resonance after a time which is longer than the integration period. Moreover, a combination of MEGNO- with $\theta_3^{{\mbox{\scriptsize max}}}$-maps still gives us much additional, valuable information.

With this approach, we found that the signature of the apsidal resonance covers a very wide neighborhood of the nominal initial condition in the $(a_{{\mbox{\scriptsize b}}},a_{{\mbox{\scriptsize c}}})$-plane and extends over 0.2 AU in $a_{{\mbox{\scriptsize b}}}$and $\simeq 0.8$ AU in $a_{{\mbox{\scriptsize c}}}$, respectively.


  \begin{figure} \par\includegraphics[width=5.6cm]{h3945f4a.eps}\hspace*{4mm} \in... ...f4b.eps}\hspace*{4mm} \includegraphics[width=5.6cm]{h3945f4c.eps} \end{figure} Figure 4: Stability of the HD 12661 system in the semi-major axes plane. The left panel is for MEGNO. The middle panel is for the semi-amplitude of apsidal librations, $\theta_3^{{\mbox{\scriptsize max}}}$, relative to the libration center $180^{\circ }$. The right panel shows the maximal eccentricity of the outer planet. The data grid has the resolution $200 \times 50$ points.
Open with DEXTER

The stability map was also computed in the $(e_{{\mbox{\scriptsize b}}},e_{{\mbox{\scriptsize c}}})$-plane (see Fig. 5). The left panel is for MEGNO, and the middle panel helps us to localize a region of the SAR. This simulation makes it possible to determine the border of a global instability of the system. This result is illustrated in the right panel of Fig. 5. In the zone of apsidal resonance, the eccentricity of the outer planet can be as large as $\simeq $0.4 for small values of $e_{{\mbox{\scriptsize b}}}$, and it cannot be larger than $\simeq $0.35 for $e_{{\mbox{\scriptsize b}}} \simeq 0.5$. After this limit, global instability occurs, leading to collisions between the companions or ejections of a planet from the system. These limits are more extended than borders determining quasi-periodic systems. Again, a fine correlation is evident among the semi-amplitude of the apsidal resonance, the maximal eccentricity of the outer planet, and structures seen in the MEGNO map. For example, let us note an island of regular motions in the MEGNO plot around ( $e_{{\mbox{\scriptsize b}}}\simeq 0.45,e_{{\mbox{\scriptsize c}}}\simeq 0.28$), and its mirror images in the $\theta_3^{{\mbox{\scriptsize max}}}$-, and $e_{{\mbox{\scriptsize c}}}^{{\mbox{\scriptsize max}}}$-plots. This is even better seen in the $(a_{{\mbox{\scriptsize c}}},e_{{\mbox{\scriptsize c}}})$-plane, shown in Fig. 6. In this case, the characteristic triangular -patterns, which are present for both $\theta_3^{{\mbox{\scriptsize max}}}$ and $e_{\rm c}^{{\mbox{\scriptsize max}}}$, have their counterparts in the plot of MEGNO.


  \begin{figure} \par\includegraphics[width=5.6cm]{h3945f5a.eps}\hspace*{4mm} \in... ...b.eps}\hspace*{4mm} \includegraphics[width=5.6cm]{h3945f5c.eps}\par\end{figure} Figure 5: Stability of the HD 12661 system in the plane of the eccentricities. The left panel if for MEGNO. The middle panel gives an estimation of the semi-amplitude of apsidal librations, $\theta_3^{{\mbox{\scriptsize max}}}$, about $180^{\circ }$. The right panel shows the maximal eccentricity of the outer planet. The data grid has the resolution $100 \times 100$ points.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=5.6cm]{h3945f6a.eps}\hspace*{4mm} \in... ...b.eps}\hspace*{4mm} \includegraphics[width=5.6cm]{h3945f6c.eps}\par\end{figure} Figure 6: Stability of the HD 12661 system in the $(a_{{\mbox{\scriptsize c}}},e_{{\mbox{\scriptsize c}}})$-plane. The left panel is for MEGNO. The middle panel gives an estimation of the semi-amplitude of apsidal librations, $\theta_3^{{\mbox{\scriptsize max}}}$, about the line of apsides antialignment. The right panel shows the maximal eccentricity of the outer planet. The data grid has the resolution $0.006\mbox{AU} \times 0.005$.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=5.75cm]{h3945f7a.eps}\hspace*{3mm} \i... ...7b.eps}\hspace*{3mm} \includegraphics[width=5.75cm]{h3945f7c.eps} \end{figure} Figure 7: Stability map of the HD 12661 system in the plane of initial mean anomalies. The left panel is for MEGNO, the middle panel shows the semi-amplitude of the librations of $\theta _3$ about the line of antialignment of apsides, the right panel shows the maximal eccentricity of the outer planet. The data grid has the resolution $90\times 90$ points.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=5.75cm]{h3945f8a.eps}\hspace*{3mm} \i... ....eps}\hspace*{3mm} \includegraphics[width=5.75cm]{h3945f8f.eps}\par\end{figure} Figure 8: Stability maps and zones of the apsidal resonance in the ( $\omega_{{\mbox{\scriptsize b}}},M_{{\mbox{\scriptsize b}}})$-plane (top), and in the ( $\omega_{{\mbox{\scriptsize c}}},M_{{\mbox{\scriptsize c}}})$-plane (bottom). The left column is for MEGNO. The middle columns give estimates of the semi-amplitude $\theta_3^{{\mbox{\scriptsize max}}}$ for the case of apsides aligned in the resonance. The right column is for the semi-amplitude of the librations when the apsides are antialigned in the exact resonance. The resolution of the plots is $90\times 90$ points.
Open with DEXTER

4 Stability vs. the initial phase of planets

An inspection of stability maps calculated for different combinations of orbital parameters, which are related to the initial orbital phase of the planets, makes it possible to determine one more important characteristic of the HD 12661 system. A map in the plane of the mean anomalies is shown in Fig. 7. It uncovers an extremely complicated pattern of stable and unstable zones. The nominal initial condition lies very close to a relatively extended stable region. The signature of the apsidal resonance covers the whole scan. The semi-amplitude of librations is $\simeq $ $ 60^{\circ}$, with much lower values in areas correlated with stability zones that are present in the MEGNO map. Also, the scan of the maximal eccentricity of the outer planet is clearly correlated with this map.

This experiment shows that the HD 12661 system's dynamics are strongly dependent on the orbital phases of the planets. A very small change of the mean anomalies can lead to both stable and unstable motions. (For an illustration, we changed the mean anomaly of the inner planet by $+20^{\circ}$, with respect to the nominal data, and the resulting, quasi-periodic evolution of $\theta _3$ is shown in Fig. 1f.) Because the fit errors of the mean anomalies are typically substantial, a conclusion whether the system behaviour is quasi-periodic or chaotic, actually seems to be very difficult to reach. Additionally, a small shift of these initial parameters results in essential deformations of stability maps in the planes of other orbital parameters. This is a consequence of the sophisticated, multi-dimensional structure of the phase space in which the HD 12661 system evolves.

Another set of simulations, performed in the planes of $(\omega_{{\mbox{\scriptsize b,c}}},M_{{\mbox{\scriptsize b,c}}})$, is illustrated in Fig. 8 and helps us to determine a dependence of the width of the apsidal resonance on the initial longitudes of the planets. This test reveals that the librations of $\theta _3$ can take place not only about $180^{\circ }$(semi-major axes are antialigned, as in the nominal system), but also about $0^{\circ }$ (semi-major axes are aligned). The apsidal resonance covers almost all the scanned area. It can be seen more clearly in the one-dimensional scan over $\omega_{{\mbox{\scriptsize b}}}$, calculated for the nominal $\omega_{{\mbox{\scriptsize c}}}$, and shown in Fig. 9. With the help of this plot we discover that the zone of circulation is limited to narrow areas around $\omega_{{\mbox{\scriptsize b}}} \simeq 60^{\circ}$, and $\simeq 240^{\circ}$. The semi-amplitude of the librations can be as low as $ 60^{\circ}$. The MEGNO scan, accompanying the $\theta_3^{{\mbox{\scriptsize max}}}$-graphs, helps us to predict that configurations corresponding to an alignment of apsidal lines would be mostly quasi-periodic and stable.


  \begin{figure} \includegraphics[width=8.8cm,clip]{h3945f9.eps}\end{figure} Figure 9: An estimation of the semi-amplitude of the librations of $\theta _3$. The scan is over $\varpi_{{\mbox{\scriptsize b}}}$, with $\varpi_{{\mbox{\scriptsize c}}}$ fixed at the nominal value. The nominal  $\omega_{{\mbox{\scriptsize b}}}$ is marked by a thin vertical line. The upper graphs correspond to libration centers $0^{\circ }$ (thin line) and $180^{\circ }$ (thick line), respectively. (Note, that to make the figure more clear, the plot for libration center $0^{\circ }$ is artificially shifted by $\simeq $ $+10^{\circ }$.) The graph at the bottom is for MEGNO. The resolution of the scan is $0.5^{\circ }$.
Open with DEXTER

5 Comparison with secular theory

In a recent paper, Malhotra (2002) analyzed a dynamical mechanism for establishing apsidal resonance in a planetary system with two planets initially in circular orbits. The author found that a perturbation (caused by the ejection of a planet from the system, or by torques from a protoplanetary disk), that changes the eccentricity of one planet's orbit, causes the other planet's orbit to be eccentric as well. Such a mechanism can explain the libration of the relative longitudes of apsides for a wide range of initial conditions. The HD 12661 system seems to be a good candidate for applying the model. The paper is based on the classical Laplace-Lagrange secular theory (Murray & Dermott 2000). Assuming that the mean motion resonances are absent, and the eccentricities and inclinations are small, we can describe the secular variation of eccentricities and apsidal longitudes by a set of first-order linear differential equations, in terms of the quantities

\begin{displaymath}h = e \sin \varpi, \quad k = e \cos \varpi, \end{displaymath}

called the eccentricity vectors. The general solution for two planets is a superposition of two eigenmodes (for each planet, j=1,2)

\begin{eqnarray*}h_j(t) &=& e_{j1} \sin( g_1 t + \beta_1 ) + e_{j2} \sin( g_2 t ... ... e_{j1} \cos( g_1 t + \beta_1 ) + e_{j2} \cos( g_2 t + \beta_2). \end{eqnarray*}


The amplitudes, eij, and fundamental frequencies, gj, of the solution components are functions of the ratio of semi-major axes, $\alpha=a_1/a_2$, and masses of the companions. The phases, $\beta_{j}$, are determined by the initial conditions (eccentricities and apsidal longitudes). Using this theory, Laughlin et al. (2002) analyzed the secular resonance in the 47 UMa system. They defined a quantity S

\begin{displaymath}S = \left\vert \frac{e_{11} e_{22} + e_{12} e_{21}} {e_{11} e_{21} + e_{12} e_{22}} \right\vert, \end{displaymath}

which determines the behaviour of $\Delta \varpi = \varpi_{{\mbox{\scriptsize c}}}-\varpi_{{\mbox{\scriptsize b}}}$. If S<1 then $\Delta \varpi$ librates, and if S>1 then $\Delta \varpi$ circulates. In the 47 UMa system the eccentricities are moderately low, and the system is rather far from strong resonances. These factors helped the authors to obtain remarkably good qualitative agreement of the secular theory with the results of numerical integrations.
  \begin{figure} \par\includegraphics[width=8.8cm,clip]{h3945f10.eps}\par\end{figure} Figure 10: Time evolution of the eccentricities as predicted by the secular theory of Laplace-Lagrange (thin lines) and obtained by the numerical integration of the equations of motion (thick lines), for the nominal HD 12661 system.
Open with DEXTER

We performed the same test for the HD 12661 system. One should be aware that in this case the validity of Laplace-Lagrange theory is very problematic, because the eccentricities are large and the nominal system is near the 11:2 mean motion resonance. In fact, when applied to the nominal HD 12661 system, the secular theory gives a very crude approximation of amplitudes and modulation frequencies of the eccentricities (see Fig. 10). To illustrate a comparison of the predictions of the secular theory with our previous estimations of the zones of the SAR, we show in Fig. 11 scans of the quantity S. Formally, the initial condition is well located within areas related to the apsidal resonance. The scan, presented in Fig. 9, is surprisingly well reproduced in the middle panel of Fig. 11. Both the width and the location of the secular resonance coincide very well. However, the observed coincidence of the results seems to be rather a matter of chance than a real value of the secular theory. Possibly, a more elaborate secular theory, incorporating higher terms in the masses and eccentricities, and terms related to the 11:2 mean motion resonance, will give a better analytical approximation of the secular dynamics.

In the end, we return to the mechanism of establishing the apsidal resonance investigated by Malhotra (2002). The author found that the amplitude of secular variation of $\Delta \varpi$ provides an observational constraint on the dynamical history of planetary systems. Because in the HD 12661 system we observe rather large full amplitudes of librations $\simeq $ $140^{\circ}$- $180^{\circ }$ and relatively large values of $e_{{\mbox{\scriptsize b}}} \simeq 0.4$, the impulse mechanism could explain the eccentric orbits and the presence of the apsidal resonance in the system. Another theory, which could explain these features of the HD 12661 system is presented in a recent paper by Chiang & Murray (2002).


  \begin{figure} \par\includegraphics[width=5.75cm]{h3945f11a.eps}\hspace*{3mm} \... ...eps}\hspace*{3mm} \includegraphics[width=5.75cm]{h3945f11c.eps}\par\end{figure} Figure 11: The quantity S calculated for the nominal HD 12661 system in the $(e_{{\mbox{\scriptsize b}}},e_{{\mbox{\scriptsize c}}})$-, $(\varpi_{{\mbox{\scriptsize b}}},\varpi_{{\mbox{\scriptsize c}}})$-, and $(e_{{\mbox{\scriptsize c}}},\varpi_{{\mbox{\scriptsize c}}})$-planes. Values S<1 indicate the presence of the apsidal resonance, as predicted by the secular theory of Laplace-Lagrange. The data grids have the resolution $90\times 90$ points.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=5.75cm]{h3945f12a.eps}\hspace*{3mm} \... ...12b.eps}\hspace*{3mm} \includegraphics[width=5.75cm]{h3945f12c.eps}\end{figure} Figure 12: Stability map of the HD 12661 system when the relative inclination of the planets and the absolute inclination i are varied. The left panel is for MEGNO. The middle panel shows the semi-amplitude of the librations of $\theta _3$ about $180^{\circ }$. The right panel shows the maximal eccentricity of the outer planet. The data grid has the resolution $4^{\circ }\times 1^{\circ }$.
Open with DEXTER

6 Dependence of the stability on companion masses

The last series of numerical tests give some insight into the dependence of system stability on planetary masses and the relative inclination of the orbits. These parameters are unconstrained by the radial velocity measurements. Not knowing even approximate values of these parameters, we can only speculate on the system dynamics when they are changed within allowable limits. We assumed in the tests that the inclinations of the orbits are equal $i=i_{{\mbox{\scriptsize b}}}=i_{{\mbox{\scriptsize c}}}$, and the relative inclination of the orbits is primarily affected by the difference of nodal longitudes, $\Delta \Omega =\Omega_{{\mbox{\scriptsize c}}}-\Omega_{{\mbox{\scriptsize b}}}$, according to the formula:

\begin{displaymath}\cos (i_{\rm rel}) = \cos(i_{{\mbox{\scriptsize b}}}) \cos(i_... ...ize b}}}) \sin(i_{{\mbox{\scriptsize c}}}) \cos \Delta \Omega. \end{displaymath}

We also assumed that both orbits are prograde. By changing i, we alter the masses because the minimum masses, $m_{{\mbox{\scriptsize b,c}}} \sin i$, corresponding to edge-on orbits, must be preserved.

The results are shown in Fig. 12. It is not surprising to find the nominal system in a chaotic zone in the MEGNO scan. Moreover, within this zone we detect an extended region of apsidal resonance. As in previously analyzed scans, we can recognize a clear correlation of the $\theta_3^{{\mbox{\scriptsize max}}}$-, $e_{{\mbox{\scriptsize c}}}^{{\mbox{\scriptsize max}}}$-, and MEGNO-plots in extended ranges of inclinations. In the plot illustrating the maximal eccentricity of the outer planet, we find a striking evidence that the secular resonance is a crucial factor in system stability: without its presence, the system can be strongly unstable. In chaotic zones related to high relative inclinations and i confined to regions around $50^{\circ}$, the system disrupts on a timescale of a few thousands of years.

Finally, we examined the stability and the width of the apsidal resonance in ( $i,\omega_{{\mbox{\scriptsize b}}}$)-, and ( $i,\omega_{{\mbox{\scriptsize c}}}$)-planes. The results of this test are shown in Fig. 13. A clear correlation of MEGNO maps with the maps showing $\theta_3^{{\mbox{\scriptsize max}}}$ and $e_{{\mbox{\scriptsize c}}}^{{\mbox{\scriptsize max}}}$ is evident. The apsidal resonance is present in very extended ranges of the system inclination. Librations of $\theta _3$ can take place with the semi-major axes antialigned or aligned. The nominal system seems to be deeply locked in the apsidal resonance, with $\theta _3$librating about $180^{\circ }$ for a wide range of companion masses.

  \begin{figure} \par\includegraphics[width=5.6cm]{h3945f13a.eps}\hspace*{3mm} \in... ....eps}\hspace*{3mm} \includegraphics[width=5.6cm]{h3945f13f.eps}\par\end{figure} Figure 13: Stability maps and zones of the apsidal resonance in the ( $i,\omega_{{\mbox{\scriptsize b}}})$-plane (top), and in the ( $i,\omega_{{\mbox{\scriptsize c}}})$-plane (bottom). The left column is for MEGNO. The plots in the middle column give estimates of the semi-amplitude of the librations, $\theta_3^{{\mbox{\scriptsize max}}}$, for apsidal lines aligned in the exact resonance. The right column is for the semi-amplitude of apsidal librations when the apsides are antialigned in the resonance. Thin lines mark the initial values of the apsidal longitudes in the nominal system. The resolution of the plots is $90 \times 120$ data points.
Open with DEXTER

7 Conclusions

The dynamical analysis of the HD 12661 system provides some interesting conclusions and raises a number of dynamical questions. The system evolution, determined by the 2-Keplerian fit to the radial velocity observations, is essentially chaotic, with a very short Lyapunov time of $\simeq $1300 yr. In spite of this strong chaotic behaviour, which can be easily recognized in the evolution of the orbital elements, the system motion seems to be bounded. This conclusion is supported by the 1 Gyr integrations of the equations of motion. The stability of the system, in the astronomical sense, is likely provided by the dynamical mechanism of the secular apsidal resonance. It acts against the chaotic diffusion, and protects the planets from an encounter. The secular resonance can be detected in wide ranges of orbital elements. The current orbital fit implies that the apsidal lines remain quite closely antialigned. By a change of initial apsidal longitudes this configuration may be switched to librations of $\theta _3$about $0^{\circ }$. Statistically, the possibility of $\Delta \varpi$circulating is much less likely than librating about $0^{\circ }$ or $180^{\circ }$.

Applied to the system, the classical Laplace-Lagrange secular theory gives only a crude description of planetary orbits calculated by numerical integrations. Because the eccentricities of the HD 12661 planets are large and the system is near the 11:2 mean motion resonance, this theory seems to be not well suited for this case. Surprisingly, theoretical predictions of the secular resonance's width and its localization in the parameters space agree qualitatively very well with the results obtained by the numerical integrations. Possibly, this agreement can be confirmed by a more elaborate version of the secular theory, accounting for higher orders of masses and eccentricities, and terms related to the 11:2 mean motion resonance.

Recently, Malhotra (2002) has shown that in a two planet-system, an impulse perturbation on the eccentricity of the outer planet, caused, for example, by planet-planet scattering, can excite the eccentricity of the inner planet and will result in apsidal libration with $\simeq $0.5 probability. In the opposite limit of a slow adiabatic perturbation (caused for example by torques from the protoplanetary disk) that increases one planet's eccentricity on a timescale much longer than the secular timescale $\simeq $ |g1-g2|-1, the apsidal resonance would occur with probability nearly equal to 1. In this case, the resulting libration amplitude will be small for nearly circular initial orbits. According to the theory, the apsidal libration amplitude provides an observational diagnostic: large amplitudes favour the impulse mechanism, whereas small amplitudes favour the adiabatic mechanism for establishing the secular resonance. In the HD 12661 system, we observe large values of $e_{{\mbox{\scriptsize b}}} \simeq 0.4$, and large full libration amplitudes of $\vert\Delta \omega\vert\simeq 140^{\circ}$- $180^{\circ }$. These factors would favour the impulse mechanism.

Using the combined MEGNO and short-term analysis of orbital elements, we found that the system evolves on the border of the 11:2 mean motion resonance. Within this resonance, the apsidal resonance acts simultaneously, providing strictly stable, quasi-periodic configurations. A question of whether the system may be captured into the 11:2 mean motion resonance, and which dynamical mechanism can lead to this capture, needs a further study. We can speculate that precise fits to the radial velocity observations, taking into account the mutual gravitational interaction between the planets (Laughlin & Chambers 2001), may provide systems that actually are locked in this resonance.

  \begin{figure} \par\includegraphics[width=5.6cm]{h3945f14a.eps}\hspace*{3mm} \in... ...5f14b.eps}\hspace*{3mm} \includegraphics[width=5.6cm]{h3945f14c.eps}\end{figure} Figure 14: Stability of the HD 12661 system for the revised initial condition by Fischer et al. (2003). The left panel is for the $(a_{{\mbox{\scriptsize b}}},a_{{\mbox{\scriptsize c}}})$-plane. The middle panel is for the $(e_{{\mbox{\scriptsize b}}},e_{{\mbox{\scriptsize c}}})$-plane. The right panel is for the $(a_{{\mbox{\scriptsize c}}},e_{{\mbox{\scriptsize c}}})$-plane. Compare with Figs. 4, 5, and 6, respectively.
Open with DEXTER

The regime of the system motion appears to be extremely sensitive to small changes in orbital elements. Thus, the determination of its real dynamical state seems to be very difficult, unless we have a very precise fit to the radial velocity data. In this case, the MEGNO stability maps change rapidly with small changes in the initial conditions. The HD 12661 system seems to give an example of a planetary configuration for which the formal distinction into chaotic and quasi-periodic motion does not seem to be satisfactory for deciding on the long-term stability in the astronomical sense. In such cases, we propose to accompany MEGNO with simultaneous analysis of the orbital elements, which are obtained as "a side effect'' of the MEGNO calculations. This approach makes it possible to identify the main sources of instabilities, and to determine dynamical factors that can stabilize the motion. We stress that the distinct behaviours of the system, detected with MEGNO, are definitely recognizable in the evolution of orbital elements. For example, the MEGNO stability maps are very well correlated with maps of the maximal eccentricities attained by the planets, in spite of a relatively short period of integrations.

The secular apsidal resonance has been found as an important factor in stabilizing the motion in a few newly discovered exosystems, e.g., $\upsilon$ Andr (Chiang et al. 2001; Chiang & Murray 2002), Gliese 876 (Lee & Peale 2002), possibly in HD 82943 (Gozdziewski & Maciejewski 2001), and in the 47 UMa system (Laughlin et al. 2002; Gozdziewski 2002). The case of HD 12661 seems to give one of the most evident examples of this resonance, and its protective role for the dynamics of exosystems containing giant planets. We hope that when the orbital fit for the HD 12661 system is more precise, the inquiring dynamics of this system can be investigated more deeply.

Recently, Fischer et al. (2003) published a revised fit to the HD 12661 radial velocity curve. The new solution does not lead to qualitative changes of the system dynamics. For a reference, we show scans of MEGNO in the plane of $(a_{{\mbox{\scriptsize b}}},a_{{\mbox{\scriptsize c}}})$, $(e_{{\mbox{\scriptsize b}}},e_{{\mbox{\scriptsize c}}})$, and $(a_{{\mbox{\scriptsize c}}},e_{{\mbox{\scriptsize c}}})$(see Fig. 14). The orbital parameters are as follows: $i=90^{\circ}$, $m_{{\mbox{\scriptsize b}}}=2.3 M_{{\mbox{\scriptsize J}}}$, $a_{{\mbox{\scriptsize b}}}=0.83$ AU, $e_{{\mbox{\scriptsize b}}}=0.35$, $\omega_{{\mbox{\scriptsize b}}}=293.1^{\circ}$, $M_{{\mbox{\scriptsize b}}}=75.5^{\circ}$, $m_{{\mbox{\scriptsize c}}}=1.57 M_{{\mbox{\scriptsize J}}}$, $a_{{\mbox{\scriptsize c}}}=2.56$ AU, $e_{{\mbox{\scriptsize c}}}=0.2$, $\omega_{{\mbox{\scriptsize c}}}=162.4^{\circ}$, $M_{{\mbox{\scriptsize c}}}=0^{\circ}$ (see Table 5 in the paper cited).

Note added in proofs: Recently, we found a new, better orbital fit to the radial velocity curve published in the preprint by Fischer et al. (2003). This fit has $\sqrt{\chi^2_{\nu}}(12) \simeq 1.37$and RMS $\simeq $ 7.9 m-1. A preliminary dynamical analysis of this solution reveals that the HD 12661 system can be locked in a stable 1:6 mean motion resonance (Gozdziewski & Maciejewski 2003, in preparation).

Acknowledgements
I am indebted to Dr. Eugenio J. Rivera for a detailed review of the manuscript and vital remarks that improved the contents. I am very grateful to Zbroja for corrections of the manuscript. Calculations in this paper have been performed on the HYDRA computer-cluster, supported by the Polish Committee for Scientific Research, Grants No. 5P03D 006 20 and No. 2P03D 001 22. This work is supported by the Polish Committee for Scientific Research, Grant No. 2P03D 001 22, and by the N. Copernicus University Research Grant No. 362-A.

References

 

Copyright ESO 2003