A&A 398, 315-325 (2003)
DOI: 10.1051/0004-6361:20021602

A dynamical analysis of the HD 37124 planetary system

K. Gozdziewski

Torun Centre for Astronomy, N. Copernicus University, Gagarin Str. 11, 87-100 Torun, Poland

Received 29 July 2002 / Accepted 7 October 2002

Abstract
In this paper we study the stability of the HD 37124 planetary system. Using initial conditions found by Butler et al. (2002) we estimate the dynamical bounds on orbital parameters that provide stable (quasi-periodic) motions of the system. The stability analysis has been performed with the help of the MEGNO technique. This method, invented by Cincotta & Simó (2000), makes it possible to distinguish efficiently between chaotic and regular dynamics of a dynamical system. The MEGNO analysis helps us to confirm the result of Butler et al. (2002) who found that the critical factor for system stability is the eccentricity of the outer planet. This parameter is not well constrained by the current set of observations. For coplanar configurations, the limiting value of the eccentricity, providing quasi-periodic motion of the system, is roughly equal to 0.55. This value is typically slightly smaller when the system inclination increases (and companions' masses grow) but for $i \simeq 25^\circ$ it can be as low as $\simeq 0.25$. The system is located in a wide stable zone in the plane of eccentricities. The dynamics are insensitive to the initial phases of the planets. If the eccentricity of the outer planet is close to the current fit value $\simeq$0.4, then the system is regular over wide ranges of the relative inclination of the planets. We analyze whether a telluric planet can survive orbitally in the habitable zone of the system. This zone is covered by the investigated region of semi-major axes between 0.6 AU and 2.4 AU. The dynamics in this zone strongly depend on the eccentricity of the outer companion. For moderately low values of this parameter, $e_{{\mbox{\scriptsize c}}} \simeq 0.1$-0.2, the habitable zone is a dynamical analogue of the Asteroid belt in the Solar system, and the habitable zone of the 47 UMa system. Moreover, in this instance, relatively wide stable areas exist. For increasing eccentricity of the outer planet the stable zones shrink rapidly, and if the parameter is larger than 0.4 they practically disappear. We describe an alternative integration of the MEGNO indicator in the planetary N-body problem, based on the symplectic leapfrog scheme and the tangent map by Mikkola & Innanen (1999).

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


1 Introduction

The search for extrasolar planetary systems conducted in the past years by the California & Carnegie Planet Search Team has recently revealed new multi-planetary systems (Butler et al. 2002). Our work is devoted to the system orbiting the star HD 37124. The variability of its radial velocity was detected at the beginning of the survey initiated by Mayor and Queloz in 1994 (Udry 2002). A first planet in $\simeq$155 day orbit was found in the frame of the Lick survey (Vogt et al. 2000). The second planet in $\simeq$2000 day orbit has been reported by Butler et al. (2002). The discovery of the second planet has been also announced by the Geneva Planet Search Team (Udry 2002).

The observational time span of HD 37124 is short with respect to the orbital period of the outer companion. During the years 1996-2002 only 30 observations have been collected. The orbital fit to the radial velocity curve of HD 37124 is not well constrained with these data (Butler et al. 2002). A double Keplerian solution found by these authors allows for the eccentricity of the outer planet to vary from 0.3 to 0.8, giving residuals consistent with the estimated velocity uncertainty 4.8 m  $\mbox{s}^{-1}$. In this work we conduct a preliminary study of the orbital stability of the HD 37124 system using the 2-Keplerian fit published in (Butler et al. 2002). The data are reproduced in Table 1. We consider them as formal osculating elements at the epoch of the periastron passage of the outer planet. Because the orbital fit is uncertain, we try to recover some qualitative dynamical properties of the HD 37124 planetary system, and to set dynamical constraints on undetermined or badly fixed orbital parameters.

The recent discoveries of the second companion of HD 37124, the system of 47 UMa (Fischer et al. 2002b), 55 Cnc (Marcy et al. 2002), HD 12661, HD 38529 (Fischer et al. 2002a), and (possibly) HD 160169 (Jones et al. 2002), revealed exosystems that can be considered more or less analogues of the Solar system. In these systems, except 47 UMa, which is by far the closest dynamical relative of the Solar system, the inner planets' semi-major axes are much smaller than 1 AU, and the outer planet in an orbit $\simeq$(3,5) AU plays the dynamical role of Jupiter. A common feature of these systems is that they do not possess giant planets at distances $\simeq$1 AU, which roughly coincide with habitable zones of Sun-like stars (Franck et al. 2000). This leads to a natural question: can an Earth-like planet form and survive in an orbit confined to such zone? A positive answer to this question is important for selecting targets of terrestrial planets' search missions. In this context, the system of 47 UMa has been studied by many authors (Jones et al. 2001; Laughlin et al. 2002; Thébault et al. 2002; Noble et al. 2002; Gozdziewski 2002). A common conclusion of these works is that stable orbits are possible in the habitable zone of 47 UMa, but the creation of terrestrial planets in the dynamical environment created by the two planets is rather unlikely. It does not seem to be a rule in the case of new planetary systems. Butler et al. (2002) announced the results of their investigations of 55 Cnc. They found a wide and stable habitable zone in this system. In the present paper we search for a possibility of stable orbits in the habitable zone of HD 37124.

The orbital parameters space of this exosystem is explored with the help of the Mean Exponential Growth factor of Nearby Orbits (MEGNO) (Cincotta & Simó 2000). This technique seems to be well designed for studying planetary dynamics in the framework of the N-body problem. The basic idea of our approach relies on a rapid determination whether a given initial condition leads to a quasi-periodic or irregular (chaotic) motion of the planetary system. We applied it already, studying the global dynamics of $\upsilon$ Andr (Gozdziewski et al. 2001), HD 82943 (Gozdziewski & Maciejewski 2001), GJ 876 (Gozdziewski et al. 2002), and, recently, to the 47 UMa system (Gozdziewski 2002). These papers contain an explanation of our motivations and technical aspects of the calculations. Summarizing the idea of the calculations briefly, MEGNO $\langle Y\rangle$ is calculated during numerical integration of the equations of motion and the respective variational equations. As was shown by Cincotta & Simó (2000), if $\langle Y\rangle(t)$converges to 2, and after a chosen integration time $t_{{\rm max}}$ we have $\vert\langle Y\rangle(t_{{\rm max}})-2\vert<\epsilon$, (where $\epsilon$ is a small quantity, of the order of 10-2-10-3), then the tested initial condition leads to a quasi-periodic, stable solution of the system; otherwise it corresponds to chaotic behavior, with a non-zero maximal Lyapunov exponent. MEGNO is related to the maximal Lyapunov exponent $\lambda_{{\mbox{\scriptsize max}}}$ through the asymptotic, linear relation $ \langle Y\rangle(t) \simeq a t + b,$where $a = \lambda_{{\mbox{\scriptsize max}}}/2$, $b \simeq 0$ for chaotic motion and a = 0, $b \simeq 2$ for a quasi-periodic solution. This relation makes it possible to estimate efficiently the Lyapunov exponent by the linear fit to $\langle Y\rangle(t)$ (Cincotta & Simó 2000). The typical time span ( $t_{{\rm max}} \simeq 10^4$ orbital periods of the outermost planet), required for the calculation of MEGNO is 10-102 times shorter than the time required to estimate the maximal Lyapunov exponent by any classical approach. MEGNO provides computational efficiency, which is crucial in the examination of large sets of initial conditions. The MEGNO concept presented by Cincotta & Simó (2000) has been generalized for the case of multidimensional systems. It is described in a recent paper by Cincotta et al. (2002).

In this work we basically follow Gozdziewski (2002). MEGNO is used for constructing one- and two-dimensional sections of the orbital parameters space in the neighborhood of the initial condition. Such maps make it possible to find dynamical constraints on the orbital parameters that cannot be determined by the fits, for instance, the inclination of the system, the masses of the planets, and the relative inclination of orbits. We are not interested in studying the orbital stability with long-term integrations. Instead, we explore the space of initial conditions in order to find regions in the phase space, which are filled up with quasi-periodic, stable evolutions. This allows us to specify whether chaotic motions are correlated with changes of orbital elements that appear in a relatively short period of MEGNO integrations, of the order of $N_{{\mbox{\scriptsize P}}}=10^4$ orbital periods of the outer companion, $P_{{\mbox{\scriptsize c}}}$. Such short time-scale, $T_{{\mbox{\scriptsize M}}} = N_{{\mbox{\scriptsize P}}} P_{{\mbox{\scriptsize c}}}$, is typically required to derive reliable estimates of MEGNO, and guarantees computational efficiency of the tool. The analysis of the short-term orbital dynamics is also very helpfull for interpreting the MEGNO results.

Most of the MEGNO integrations in this paper are 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.

2 MEGNO signature of the initial condition

First we examined the stability of the 2-Keplerian fit, found by Butler et al. (2002). The results of this test are illustrated in Figs. 1a,c, where a number of MEGNO graphs, derived for randomly selected initial variations[*], are plotted as functions of time. We observe that for the integration time span $\simeq T_{{\mbox{\scriptsize M}}}$, the character of MEGNO convergence is not always quite clear. In particular, this concerns MEGNO variations shown in two plots, labeled by (1) and (2), in Fig. 1c.

Integrations continued over $10 T_{{\mbox{\scriptsize M}}}$ do not give clear answers for all choices of the initial variations, either. Actually, the two $\langle Y\rangle$-graphs, labeled with (1) and (2) in Fig. 1c, seem to be slowly divergent. We have observed this kind of MEGNO behavior in many runs. Note, that oscillations of $\langle Y\rangle$ are correlated with variations of the eccentricities (compare Figs. 1a,c and Fig. 2). In all examined runs the energy integral grows linearly with time over a few orders of magnitude (from 10-13 to 10-9) in the time span of integrations. We already noticed the effect of divergent evolution $\langle Y\rangle$ in some models (Gozdziewski et al. 2001). It can be explained by an accumulation of numerical errors. A wrong classification of orbits is rather unlikely with the mean value of MEGNO being in fact very close to the value expected for a quasi-periodic orbit. Moreover, to verify the results in such unclear cases we propose to use an alternative method of MEGNO integration, which is based on the mixed-variable, symplectic algorithm by Wisdom & Holman (1991). It is described in the next section.


 

 
Table 1: Orbital parameters of the HD 37124 system (Butler et al. 2002). Mass of the central star is $0.91~M_{\odot}$.
Orbital parameter Planet b Planet c
$\mbox{m}_2 \sin i$ [M $_{{\mbox{\scriptsize J}}}$] .............. 0.86 1.01
a [AU] ............................. 0.54 2.95
P [d] ............................... 153 (1) 1942 (400)
e ...................................... 0.10 (0.06) 0.4 (0.2)
$\omega$ [deg] ......................... 97 (40) 256 (120)
$T_{{\mbox{\scriptsize p}}}$ [JD-2450000] ....... 1227 (300) 1828 (400)
$M(\tau_{{\mbox{\scriptsize c}}})$ [deg] ............... 334.12 0.00



  \begin{figure} \par\includegraphics[width=16cm,clip]{h3884f1.eps}\end{figure} Figure 1: MEGNO signatures of the HD 37124 system. The plots on the left are for the Bulirsh-Stoer integrator, the plots on the right are derived with the help of the tangent map and Yoshida's sixth-order symplectic integrator. Panels a) and b) show typical temporal variations of Y(t). Panels c) and d) illustrate mean value of MEGNO $\langle Y\rangle$, obtained from a number of program runs (with the initial variations chosen at random). Thick vertical lines mark 104 periods of the outer companion. The time step of the symplectic integrator is 1/20 of the period of the innermost planet.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=7.4cm,clip]{h3884f2.eps}\end{figure} Figure 2: Time evolution of eccentricities for the system defined in Table 1.
Open with DEXTER

3 Symplectic calculation of MEGNO

Nowadays, the W-H algorithm is widely used for solving the equations of motion for the planetary N-body problem. The W-H mapping has the remarkable property of conserving the total energy integral, and is more efficient than other integration algorithms. For our purposes, it is important that the variational equations of the N-body problem can be solved within the same, symplectic scheme (Mikkola & Innanen 1999). These authors have shown that the leapfrog algorithm can be differentiated analytically. The tangent map obtained in this way can be integrated with the W-H algorithm, simultaneously with the equations of motion, instead of solving the variational system directly. This method was designed for computing the maximal Lyapunov exponent (Mikkola & Innanen 1995). It is straightforward to realize that it can also be applied to computations of MEGNO.

Let us assume that the initial condition is given at time t0: ${\bf x}(t_0)$ are the Cartesian coordinates of planets, and ${\bf v}(t_0)$ their velocities. The vector of initial variations $\vec{\delta}$ is combined from variations $\delta {\bf x}(t_0), \delta {\bf v}(t_0)$. After the leapfrog step h we obtain the orbital solution ${\bf x}(t_0+h), {\bf v}(t_0+h)$, and the variations $\delta {\bf x}(t_0+h), \delta {\bf v}(t_0+h)$ of the coordinates and velocities, respectively. A numerical procedure computing the right hand sides of the variational equations provides the time derivatives of the variations, $\dot {\vec{\delta}}(t_0+h)$. Let us assume that N integration steps are performed in the interval of time [0,t], at moments tk= t0 + k h, $k=0,\ldots,N$.

The temporal value of MEGNO, written in the integral form (Cincotta & Simó 2000),

\begin{eqnarray*}Y(t) = \frac{2}{t} J(t), \quad J(t)= \int_0^t \frac{\dot \delta}{\delta} s {\rm d}s, \end{eqnarray*}


can be approximated by the extended trapezoidal rule (Press et al. 1992) as follows:

\begin{eqnarray*}J(t) = \tilde J(t) + O(t^3 f''/N^2), \end{eqnarray*}


where

\begin{eqnarray*}&& \tilde J(t) = h \left( \frac{1}{2} f_0 + f_1 + \ldots + \fr... ...s, \quad \delta = \vert\vec{\delta}\vert, \quad f_k= f(t_k), \end{eqnarray*}


and f'' is evaluated at $t_{{\mbox{\scriptsize c}}} \in [0,t]$. Then

\begin{eqnarray*}Y(t) = \tilde Y(t) + O(h^2 f''), \end{eqnarray*}


where $\tilde Y(t)$ is the numerical approximation of Y(t).

Finally, MEGNO is evaluated as the mean over Y(t):

\begin{eqnarray*}\langle Y\rangle(t) = \frac{1}{t} \int_{0}^t Y(t) = \frac{1}{t} \sum_{k=0}^N \tilde Y(t_k) + O(h f''). \end{eqnarray*}


Because $\langle Y\rangle(t)$ is a quantity of the first order of magnitude, this rather crude approximation is sufficient in practice. To test this simple algorithm we applied the second order MVS integrator of Wisdom & Holman (1991), and the higher-order schemes of Yoshida (1990). These later algorithms are combined from the second order leapfrog, and their application is straightforward. The equations of motion have been integrated in the Jacobi coordinates.

Figures 1b,d show a number of MEGNO signatures of the nominal system (see Table 1), derived with the 6th order integrator of Yoshida. The time step has been set to $P_{{\mbox{\scriptsize b}}}/20$. ( $P_{{\mbox{\scriptsize b}}}$ is the orbital period of planet b). In this case we obtained a perfect convergence of $\langle Y\rangle$ in every run. The total energy errors are smaller by about two orders of magnitude than in the case of ODEX. We did not find problematic cases of MEGNO divergence which are described in the previous section. This result confirms that the orbital fit is related to a regular, stable system.

The algorithm has been tested through another set of simulations. For the tests we selected the planetary system described in Sect. 5: it contains two planets (their initial parameters are the same as those given in Table 1, except that $e_{{\mbox{\scriptsize c}}}=0.1$), and an Earth-like planet. The scans are performed over the initial semi-major axis of the small planet, and the other initial orbital parameters are the same as those selected for the experiments described in Sect. 5.

Figure 3 illustrates a comparison of MEGNO scans obtained with the symplectic integrators of the second-, the fourth-, and the sixth-order, and with the Bulirsh-Stoer integrator (ODEX). The time step has been set to $P_{{\mbox{\scriptsize b}}}/20$ for the sixth-order scheme, and $P_{{\mbox{\scriptsize b}}}/30$ for the second- and fourth-order schemes. The results of these calculations coincide very well. Qualitative differences between the MEGNO estimations, which are present at a few isolated $a_{{\mbox{\scriptsize0}}}$ points, can be caused by the relatively short integration time, equal to $T_{{\mbox{\scriptsize M}}}$. The observed differences are also related to the fact that the initial tangent vector has been selected at random, and MEGNO convergence depends, in general, on this selection (Gozdziewski et al. 2001).

Other simulations have shown that the results of MEGNO calculations obtained with ODEX and with the leapfrog integrators can differ substantially if the Earth-like bodies are evolving into highly eccentric orbits, with e0 about 0.6-0.8. It is well known that the symplectic integrators applied with a fixed-step size are not well suited for systems permitting highly eccentric orbits and close encounters between the bodies. In such an instance, the errors in the obtained solution are substantial and a reliable approximation of MEGNO is not possible.


  \begin{figure} \par\includegraphics[width=17cm,clip]{h3884f3.eps}\end{figure} Figure 3: Comparison of MEGNO scans obtained with the help of the ODEX integrator (lines) and with the help of symplectic algorithms of Yoshida: the second-order (filled circles), the fourth-order (open circles), and the sixth-order (triangles). The dynamical model describes the motion of a planetary system with two giant planets and an Earth-like object. The initial orbital parameters of the planets are the same as given in Table 1, with the eccentricity of the outer planet lowered to $e_{{\mbox{\scriptsize c}}}=0.1$. The initial value of semi-major axis $a_{{\mbox{\scriptsize0}}}$ of the small planet ( $m_{{\mbox{\scriptsize EB}}}=0.001 \mbox{M}_{{\mbox{\scriptsize J}}}$) is marked on the x-axis, its initial eccentricity $e_{{\mbox{\scriptsize 0}}}=0$. See Sect. 5 for more details. The resolution of the scan is 0.0015 AU.
Open with DEXTER

This is the main reason that the tangent map approach requires a more detailed analysis than reported here, in order to be applied "productively''. (Note that also for this reason the further simulations in this paper have used the ODEX-based integrations.) However, we foresee that the tangent map method can provide many advantages over a general purpose integrator. The symplectic algorithms of integration help to avoid artificial changes of the regime of motion in borderline cases, which are caused by the errors in the energy. Thus their application should help to classify orbits more precisely, with generally less computational effort than that required by the general purpose integrators. Calculations of MEGNO, which utilize such integrators, become very CPU-expensive if the number of bodies grows. In this context the symplectic integrators can be rather easily parallelized (Saha et al. 1997). This serves the purpose of running the MEGNO code in computer-cluster environments, and increases the computational efficiency. The problem of numerical difficulties caused by the presence of highly eccentric orbits can be avoided through an application of regularized algorithms of integration. A very promising way of dealing with this issue seems to be the logarithmic Hamiltonian approach (Mikkola & Wiegert 2002; Mikkola et al. 2002). This variant of MEGNO integrations is under development (Gozdziewski 2002).

4 Dynamical constraints on the orbital parameters

To estimate the dynamical constraints on the orbital parameters, we computed a number of stability maps in wide ranges around the initial condition. (The position of the initial data is marked in the 2D maps by the intersection of two lines). Assuming the minimum masses we computed the $(M_{{\mbox{\scriptsize b}}},M_{{\mbox{\scriptsize c}}})$-, and $(\omega_{{\mbox{\scriptsize b}}},\omega_{{\mbox{\scriptsize c}}})$- MEGNO maps. The latter is shown in Fig. 4. We did not find any unstable zone in this map. A very similar picture is obtained in the $(M_{{\mbox{\scriptsize b}}},M_{{\mbox{\scriptsize c}}})$-map. With these data, we can conclude that the system's behavior is insensitive to changes of the initial phases of the planets.

The initial condition is located in a wide stable zone which is seen in the $(e_{{\mbox{\scriptsize b}}},e_{{\mbox{\scriptsize c}}})$-plane, see Fig. 5. We already noticed that the eccentricity of the outer planet is not well constrained by the radial velocity observations (Butler et al. 2002). The MEGNO map helps us to determine the dynamical constraints on this parameter. The border line between regular and chaotic zones is generally complex, but we can conclude that if the eccentricity of the inner planed is close to the fitted value, then the dynamics preclude regular systems with $e_{{\mbox{\scriptsize c}}} > 0.55$. A similar estimate, $e_{{\mbox{\scriptsize c}}} > 0.65$, derived on the basis of direct integrations, is reported by Butler et al. (2002). These authors found that coplanar systems are destabilized if $e_{{\mbox{\scriptsize c}}}$ is grater than this critical value.


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{h3884f4.eps}\end{figure} Figure 4: Stability map of the HD 37124 system when both the arguments of periastron of the planets are changed. The initial conditions are defined in Table 1, and considered as osculating elements at the epoch of the periastron passage of the outer planet. The position of the initial condition is marked by the intersection of two lines. The resolution of the plot is $6^{\circ }\times 6^{\circ }$. Note that the overall range of $\langle Y\rangle$ is close to 2, and small fluctuations about this value indicate globally stable system.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{h3884f5.eps}\end{figure} Figure 5: Stability map of the HD 37124 system in the $(e_{{\mbox{\scriptsize b}}},e_{{\mbox{\scriptsize c}}})$-plane. The data grid has a resolution of $50 \times 50$ points. Dark shaded areas correspond to chaotic zones, and light-gray areas, centered about 2, represent quasi-periodic motions of the system.
Open with DEXTER

We examined the dependence of this dynamical constraint on the inclination, assuming coplanar systems. The results are shown in Fig. 6. MEGNO was scanned in the $(a_{{\mbox{\scriptsize c}}},e_{{\mbox{\scriptsize c}}})$-plane. This plane was selected because the orbital period (and the semi-major axis) of the outer companion is not determined precisely by the orbital fit, either. For the examined inclinations: $90^{\circ}$, $60^{\circ}$, $30^{\circ}$, and $15^{\circ}$, the system is located in an extended stable zone. If the inclination decreases (and the masses grow), the border of stability is shifted toward low eccentricities. This shift is substantial for large mass-factors. For $f\simeq 4$ the limiting value of $e_{{\mbox{\scriptsize c}}}$ is typically about 0.5. Moreover, in all cases the critical eccentricity can be as low as $\simeq$0.25-0.3 in a few narrow ranges of $a_{{\mbox{\scriptsize c}}}$, specifically in a region close to the nominal value of the semi-major axis. These values of $e_{{\mbox{\scriptsize c}}}$ are substantially smaller than the fit value, $e_{{\mbox{\scriptsize c}}}=0.4$. Similar estimates of the critical eccentricity can be derived from the $(M_{{\mbox{\scriptsize c}}},e_{{\mbox{\scriptsize c}}})$-maps, shown in Fig. 7. These maps have been calculated for inclinations $90^{\circ}$, and $30^{\circ}$. In the case of the minimum masses, the dynamics are regular and practically insensitive to the initial position of the outer companion, if the eccentricity is less than $\simeq$0.55-0.6. This estimation agrees with the value obtained in the previous simulation. Let us note that for $i=30^{\circ}$ the critical eccentricity is reduced to $\simeq$0.4 for $M_{{\mbox{\scriptsize c}}} \simeq 0^{\circ}$.


  \begin{figure} \par\includegraphics[width=7.4cm,clip]{h3884f6a.eps}\hspace*{5mm}... ...c.eps}\hspace*{5mm} \includegraphics[width=7.4cm,clip]{h3884f6d.eps}\end{figure} Figure 6: Stability maps of the HD 37124 system for varied inclinations: a) $i=90^{\circ },f=1$; b) $i=60^{\circ },f=1.15$; c) $i=30^{\circ },f=2$; d)  $i=15^{\circ },f=3.86$. The resolution of the plots is $50 \times 50$ points.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=7.4cm,clip]{h3884f7a.eps}\hspace*{5mm} \includegraphics[width=7.4cm,clip]{h3884f7b.eps}\end{figure} Figure 7: Stability maps of the HD 37124 system for varied eccentricity and the mean anomaly of the outer companion: a) for the minimal masses, b) for $i=30^{\circ },f=2$. The resolution of the maps is $6^{\circ } \times 0.01$.
Open with DEXTER

The relation between the relative inclination, the inclinations of orbits to the line of sight, and the longitudes of their lines of nodes is given, for a 2-planet system, through $\cos (i_{{\rm rel}}) = \cos(i_{{\mbox{\scriptsize b}}}) \cos(i_{{\mbox{\scripts... ..._{{\mbox{\scriptsize b}}}) \sin(i_{{\mbox{\scriptsize c}}}) \cos \Delta \Omega,$where $\Delta\Omega=\Omega_{{\mbox{\scriptsize c}}}-\Omega_{{\mbox{\scriptsize b}}}$, $i_{{\mbox{\scriptsize b}}}$, $i_{{\mbox{\scriptsize c}}}$ are the inclinations of the planets, and $\Omega_{{\mbox{\scriptsize b}}}$, $\Omega_{{\mbox{\scriptsize c}}}$ are the longitudes of their lines of nodes. In this test $i=i_{{\mbox{\scriptsize b}}}=i_{{\mbox{\scriptsize c}}}$, and changes of i alter the planetary masses according to the mass factor $f=\sin^{-1}i$. A sequence of MEGNO maps in the $(i,\Delta\Omega)$-plane, calculated for $e_{{\mbox{\scriptsize c}}}=0.4$ (the value of the orbital fit), $e_{{\mbox{\scriptsize c}}}=0.5$, $e_{{\mbox{\scriptsize c}}}=0.55$, and $e_{{\mbox{\scriptsize c}}}=0.6$, are shown in Fig. 8. Maps corresponding to the bounding eccentricities are accompanied by their 3D-versions (Fig. 9) to illustrate that the borders between regular and chaotic zones are very sharp.

If the eccentricity of the outer planet is lower than $\simeq$0.55, then the stability is preserved in extended ranges of the relative inclination and masses of the companions. Instabilities are confined to a few narrow areas. We did not observe such extended zones of stability in other exosystems analyzed up to now. If the eccentricity of the outer planet becomes large, about 0.6, the distribution of stable zones is reversed - the stability is preserved only in a few relatively narrow areas. In this case the system becomes very active dynamically. Figures 8a-d illustrate clearly how the unstable zones expand with increasing $e_{{\mbox{\scriptsize c}}}$. The last panel, for $e_{{\mbox{\scriptsize c}}}=0.6$, makes it possible to reach an interesting conclusion: although the system is unstable for $i=90^{\circ}$, a small change of the relative inclination, or an increase of masses can stabilize its motion. Finally, we examined the dependence of the stability on the masses of the planets, and on the eccentricity of the outer companion, in coplanar systems. The results of this test are shown in Fig. 10. These data help us conclude that the critical value of $e_{{\mbox{\scriptsize c}}}$ strongly depends on the masses. This limit can be as low as 0.25 for $i \simeq 25^\circ$. Note that the map is closely related to scans shown in Fig. 6.


  \begin{figure} \par\includegraphics[width=7.6cm,clip]{h3884f8a.eps}\hspace*{5mm}... ...c.eps}\hspace*{5mm} \includegraphics[width=7.6cm,clip]{h3884f8d.eps}\end{figure} Figure 8: Stability of the HD 37124 system when the inclination i and the relative inclination are varied, through changes of $\Delta\Omega=\Omega_{{\mbox{\scriptsize c}}}-\Omega_{{\mbox{\scriptsize b}}}$. Plots  a), b), c), and d) are for $e_{{\mbox{\scriptsize c}}}=0.4,\; 0.5,\; 0.55,\; 0.6$, respectively. The resolution of the maps is $6^{\circ } \times 1^{\circ }$.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{h3884f9a.eps}\hspace*{2mm} \includegraphics[width=8.8cm,clip]{h3884f9b.eps}\end{figure} Figure 9: Three-dimensional stability maps for the initial condition, when the system inclination i, and $\Delta \Omega $ are varied. The left panel is for the fit value, $e_{{\mbox{\scriptsize c}}}=0.4$, and the right panel is for $e_{{\mbox{\scriptsize c}}}=0.6$. The resolution of the maps is $6^{\circ } \times 1^{\circ }$.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{h3884f10.eps}\end{figure} Figure 10: Stability of the HD 37124 system as a function of the system inclination i, and the eccentricity of the outer companion. The resolution of the map is $1^{\circ } \times 0.02$.
Open with DEXTER

5 Earth-like bodies in the habitable zone of HD 37124

The habitable zone (HZ) in a main-sequence star planetary system can be defined through the requirement of moderate planetary surface temperatures, which imply the liquid phase of water (Lissauer 1999). A more restrictive and general definition of habitability, based on the long-term possibility of photosynthetic life on extrasolar planets, is considered by Franck et al. (2000). These authors developed a model of the habitable zone of a main-sequence star, taking into account astrophysical, climatological, bio-geochemical and geodynamical processes providing a possibility of photosynthesis. The exact solution for this model is complex. For our purposes it is enough to interpolate rough boundaries of the HZ using the known solutions for stellar-mass cases 0.8, 1.0 $M_{\odot}$ (Franck et al. 2000). For $\simeq$ $0.91~M_{\odot}$ of HD 37124, the inner boundary of the HZ  $\simeq 0.7$ AU, and its outer boundary does not expand beyond $\simeq 1.5$ AU. This zone is preserved over a few Gyr of time.

The HZ of HD 37124 is free from giant planets, thus a natural question, which we have already stated, arises: can an Earth-like body (EB) form and survive in a stable orbit within the zone? Here we assume that the EB has already been created, and we investigate the dynamical structure of the HZ. Essentially, we apply the same method that has been used in Gozdziewski (2002). It is based on MEGNO analysis of the orbital dynamics of a test planet which is initially placed in the HZ.


  \begin{figure} \par\includegraphics[width=16.5cm,clip]{h3884f11.eps}\end{figure} Figure 11: Stability of an Earth-like body in the habitable zone of the HD 37124 system when the initial semi-major axis of the terrestrial planet is varied. The resolution of the scan is 0.0015 AU. The upper plot is for MEGNO, the lower graph illustrates maximal values of the eccentricity of the test planet that are reached during the integration time (at most 104 orbital periods of the outer planet). Small triangles about 5 mark the nominal positions of mean motion resonances with the outer planet, and they are labeled. Labels (a), (b), (c), and (d), accompanying large triangles, mark positions of the resonances examined in more detail in Fig. 12.
Open with DEXTER


  \begin{figure} \includegraphics[width=16cm,clip]{h3884f12.eps}\end{figure} Figure 12: Examples of resonant motion of the EB. a) The critical argument of the secular resonance $\nu _5$ is plotted as a function of time. The initial semi major-axis a0=0.86 AU is marked in the previous figure. b) Time evolution of the longitudes of periastron $\varpi_{{\mbox{\scriptsize b}}}$, $\varpi _0$indicate the presence of the secular resonance $\nu _5$ in the neighborhood of a0=1.12 AU; thick line is for the planet, dots are for the EB. c) The critical arguments of the secular resonances $\nu _5$, $\nu _6$ are plotted as a function of time; these resonances explain the eccentricity jump at $a_0\simeq 1.577$ AU, and an ejection of the EB from the system after 20 000 yr only. d) The critical argument $\theta =2\lambda _0-3\lambda _c+\varpi _0$ of the mean motion resonance 3:2 with the outer planet, indicating stable orbital evolution of the test body.
Open with DEXTER

The resonance overlapping criterion for the circular restricted three body problem predicts large scale chaos for small bodies, with initially low-eccentric orbits, for semi-major axes $a_0 \in [a_{{\mbox{\scriptsize P}}}-\Delta a,a_{{\mbox{\scriptsize P}}}+\Delta a]$, where $ \Delta a \simeq 1.5 a_{{\mbox{\scriptsize P}}} \mu^{2/7},$ $a_{{\mbox{\scriptsize P}}}$ is the semi-major axis of the perturber, and $\mu$ is the mass ratio of the primary and the star (Holman & Murray 1996). Applying this criterion, we found that $\Delta a_{{\mbox{\scriptsize b}}} \simeq 0.1$ AU, and $\Delta a_{{\mbox{\scriptsize c}}} \simeq 0.6$ AU. This criterion permits stable motions in the zone between $\simeq$0.65 AU and 2.35 AU. This region overlaps the HZ, and we restrict our interest to this area.


  \begin{figure} \par\includegraphics[width=6.5cm,clip]{h3884f13a.eps}\hspace*{5mm... ....eps}\hspace*{5mm} \includegraphics[width=6.5cm,clip]{h3884f13d.eps}\end{figure} Figure 13: The stability of an Earth-like body in the habitable zone of the system when the eccentricity of the outer companion is varied. At the top: the left panel is for the minimal masses, the right panel is for $i=30^{\circ },f=2$. At the bottom: maximal values of the eccentricity of the EB reached during the integration time, at most 104 periods of the outer planet. The resolution of the scans is $0.015~\mbox{AU} \times 0.05$. The irregular patterns visible in the bottom maps, in the chaotic areas, are caused by short integration time due to the limit of 5, which has been set to MEGNO.
Open with DEXTER

In the numerical experiments, we assumed initial conditions of the planets that provide quasi-periodic evolution. In this instance, dynamics of the EB are unaffected by (possibly) chaotic motions of the giant planets. The system has been assumed to be edge-on and coplanar. The EB mass has been set to 0.001  $M_{\mbox{\scriptsize J}}$, its eccentricity at the initial time to e0=0, the longitude of periapse $\varpi_0=0^{\circ}$, and the mean anomaly $M_0=0^{\circ}$. Variations of MEGNO and orbital elements were calculated over $\simeq$ $6 \times 10^4$ yr, i.e., over a time span corresponding to 104 orbital periods of the outer planet.

The results of a scan over the initial semi-major axis of the EB are shown in Fig. 11. MEGNO is marked with thin lines in the upper plot. The lower plot in this figure illustrates the maximum values of the eccentricity e0 reached by the EB during the integration time. (The integrations were stopped if MEGNO was larger than 5.) Note that in this experiment $e_{{\mbox{\scriptsize c}}}=0.1$; it will be clearer later on why it was reduced to such a small value in the simulation of the HZ.

With the help of this simulation we can conclude that the region between the planets recalls the Asteroid belt in the Solar system. In the HD 37124 system the outer planet is acting as Jupiter. The dynamical structure of the HZ is also very similar to the HZ in the 47 UMa system (Gozdziewski 2002). The global instabilities, predicted by the resonance overlap criterion, are clearly present in the MEGNO scan. A large number of MEGNO spikes in the investigated zone can be identified with resonances of the EB with the planets. The nominal positions of mean motion resonances, (i.e., unaffected by the precession of the orbits), with the outer planet are marked by small triangles, and they are labeled. It can be seen that these positions are substantially shifted from the actual locations of the resonances. The strongest mean motion resonances are 4:1, 3:1, 5:2, 2:1. They cause a rapid increase of the eccentricity of the EB. A few MEGNO spikes and relatively wide plateau can be identified with positions of secular resonances. This identification is illustrated and labeled in Fig. 12. A broad instability, centered around 0.86 AU, (this value is marked by a triangle in Fig. 11, and labeled (a)), is caused by the analogue of the $\nu _5$ secular resonance in the Solar system. The critical argument of this resonance is plotted as a function of time in Fig. 12a. Because the maximal eccentricity of the EB remains moderately low in this area, in fact the resonance can stabilize the motion. This needs further explanation. Panel (b) in Fig. 12 is for $a_{{\mbox{\scriptsize0}}}=1.12$ AU. Also in this case the argument of periastron of the inner planet and the argument of perihelion of the EB are clearly correlated. A relatively narrow instability marked by (b) in Fig. 11 is possibly an effect of interaction of this resonance with the 13:3 mean motion resonance. A rapid jump of eccentricity of the EB, centered at $a_{{\mbox{\scriptsize0}}}=1.577$, is correlated with MEGNO spike too. In this case two secular resonances, $\nu _5$ and $\nu _6$, act simultaneously, and this is illustrated in Fig. 12c. These resonances are also present in a wide unstable zone, centered at 2.1 AU, and cause rapid ejection of the EB from the system, on the time-scale $\simeq$5 $\times 10^4$ yr. But the orbital resonances do not always have a catastrophic influence on the motion of the EB. The last example, shown in Fig. 12d, illustrates the dynamical protection of the EB by the 3:2 mean motion resonance. The plot of MEGNO helps us to precisely determine its width, as well as it is revealing some details of its structure.

The MEGNO scan makes it possible to find a relatively extended zone of stable motions, confined to $a_{{\mbox{\scriptsize0}}} \simeq (1,1.6)$ AU. This area is divided by only a few narrow regions of mean motion resonances. In this zone the eccentricity of the EB, which was initially set to zero, actually is forced to $\simeq$0.1. Note that a similar kind of short-term dynamics is present in the 47 UMa system (Gozdziewski 2002). The influence of the inner planet on the motion of the EB is weak, except for the regions of the secular resonances. The strongest instabilities are caused by the outer planet.

The results could be encouraging for a search of terrestrial planets in the HD 37124 system, but one should remember that the fitted eccentricity of the outer planet is much larger than 0.1, which has been selected in the simulation. We computed the MEGNO map in the $(a_{{\mbox{\scriptsize0}}},e_{{\mbox{\scriptsize c}}})$-plane to examine the influence of this parameter on the dynamics of the EB. The results are shown in Fig. 13, which is the 2-dimensional analogue of Fig. 11. The panels at the top show MEGNO for minimum masses of the planets, and for a system inclined by $30^{\circ}$ (f=2), respectively. This test confirms that the eccentricity of the outer planet is the critical factor for stability of the EB. The width of unstable zones related to the mean motion resonances develops rapidly with growing eccentricity of the outer giant planet. If  $e_{{\mbox{\scriptsize c}}}$ is about 0.4 then the stable zones practically disappear. The plots at the bottom of Fig. 13 illustrate a fine correlation of the maximal eccentricity of the EB, attained during the integration time, with the MEGNO signatures.

The conclusions derived with the help of MEGNO scans can be further refined by direct long-term integrations of systems with the eccentricity of the outer planet set relatively close to the fitted value. In this case the formally chaotic motions of the EB can be bounded over very long periods of time. Possibly, outside the strongest resonances the EB can survive. However, here we meet the problem of forming the EB in the unstable dynamical environment which is present in the HD 37124 system. A similar question has been already addressed in the instance of 47 UMa. Taking into account the dynamical analogies between the HZ's of these two planetary systems, we can claim that also in the case of HD 37124 the process of the EB formation has been hindered or even made impossible if the planets formed quickly (Kortenkamp & Wetherill 2000; Thébault et al. 2002; Laughlin et al. 2002). Possibly, this problem requires a separate, detailed study. The HD 37124 system seems to be a good candidate for investigating the creation of terrestrial planets.

6 Conclusions

In this paper we analyze the orbital dynamics of the HD 37124 planetary system, and we search for stable orbits in the habitable zone of this system. The stability of the planetary system is considered in the strict sense of quasi-periodic motions. Our analysis follows mainly the previous work devoted to the 47 UMa system (Gozdziewski 2002).

The dynamical analysis of the initial condition, found by the 2-planet Keplerian fit to the radial velocity observations of HD 37124 (Butler et al. 2002), makes it possible to determine some dynamical characteristics of this system. The crucial factor for the system stability is the eccentricity of the outer planet. This parameter is not well constrained by the current set of observations. These data permit the eccentricity for the outer planet to vary from 0.3 to 0.8 (Butler et al. 2002). The MEGNO analysis makes it possible to determine the upper limit of the eccentricity. For the minimum masses, this limit is about 0.55. If the system inclination is decreasing, then the critical value can be reduced, even to $\simeq$0.25 for $i \simeq 25^\circ$. A precise value of this limit depends on many parameters, e.g., the eccentricity of the inner planet, the semi-major axis of the outer planet, and the relative inclination of the companions. Varying the system inclination and the semi-major axis of the outer planet within reasonable ranges, we found a few narrow unstable zones in the $(a_{{\mbox{\scriptsize c}}},e_{{\mbox{\scriptsize c}}})$-plane, which extend toward $e_{{\mbox{\scriptsize c}}} \simeq 0.25$. On the other hand, small changes of the relative inclination of the planets can stabilize the system, even for high eccentricities of the outer planet, and large mass factors.

If the eccentricity of the outer planet is smaller than $\simeq 0.55$ then the dynamics are regular within wide ranges of the relative inclination and companions' masses. In this case the regime of motion is basically not sensitive to the initial phases of the planets. For the nominal system no instabilities are found within the whole range of the initial mean anomalies and the arguments of periastron. These two characteristics can be explained by a relatively large separation of the companions and a lack of strong resonances in their motion.

The habitable zone of HD 37124 seems to be a close analogue of the Asteroid belt, and of the habitable zone of 47 UMa. If the eccentricity of the outer planet is relatively small, about 0.1-0.2, then chaotic zones are confined to narrow ranges of the initial semi-major axes of an Earth-like planet. The main sources of unstable motions of the EB can be identified with the low-order mean motion resonances with the outer planet, and with the secular resonances $\nu_5,\nu_6$ (analogues of the secular resonances in the Solar system). If the eccentricity grows, then the regular zones shrink rapidly. For $e_{{\mbox{\scriptsize c}}}\simeq 0.4$, comparable to the value of the orbital fit, they disappear almost completely. This result leads to a strong dynamical difficulty in maintaining a stable orbit in the region between the planets. This claim can be further refined by long-term numerical integrations of test bodies. Such integrations could help to verify whether the chaotic motions are bounded. A few open questions remain. If the eccentricity of the outer planet is about 0.1-0.2, could a terrestrial planet form in a zone around $\simeq$1 AU? If yes, the HZ would be stable in a quite broad area, and the HD 37124 exosystem might be a target for a search of such planet.

This work gives one more example of the application of the MEGNO technique for studying planetary dynamics. To increase the efficiency and precision of the MEGNO computations, we propose to use the symplectic integration schemes and the tangent map (Mikkola & Innanen 1999) for solving the variational equations of motion. Preliminary tests of the symplectic algorithm, performed in this paper, give encouraging results. The use of the logarithmic Hamiltonian regularization (Mikkola et al. 2002) is a promising way of avoiding computational difficulties when the eccentricities become large.

Acknowledgements
I am very grateful to Dr. Claudia Giordano and Dr. Pablo Cincotta for a critical review of the manuscript and vital remarks that helped me to improve its contents. I am indebted to Zbroja for correcting the manuscript. I would like to thank Dr. Jet Katgert, the Editor, for extensive language revision of the final version of the manuscript. The numerical calculations in this paper have been performed on the HYDRA computer-cluster supported by the Polish Committee of Scientific Research, Grant No. 5PO3D 006 20. This work has been supported by the Polish Committee of Scientific Research, Grant No. 2P03D 001 22, and by the N. Copernicus University Research Grant No 362-A.

References

 

Copyright ESO 2003