© ESO, 2011

1. Introduction

More than 400 extrasolar planets have been detected, among which there are at least 40 multi-planetary systems. The majority of the detections were achieved using the radial velocity technique, which is incapable of measuring the mutual inclination of the orbits. Because of this lack of knowledge, as well as a general belief that the orbits in planetary systems are coplanar like those in our Solar System, only a few studies of the dynamics of extrasolar systems have considered the three-dimensional case. Some analytical works have highlighted the possibility that stable highly non-coplanar planetary systems may exist (e.g. Libert & Henrard 2007; Libert & Tsiganis 2009a), while studies of the formation of these systems remain in progress (e.g. Levison et al. 1998; Thommes & Lissauer 2003; Chatterjee et al. 2008; Libert & Tsiganis 2009b).

Terrestrial planets have not – until now – been easily detected because of the limitations of the applied methods. Hence it is quite important to perform numerical studies to find out whether these small planets can be in stable configurations with the giant planets observed so far. In the same way, many studies of the restricted problem have investigated the stability of additional planets in known single and multi-planetary systems (e.g. Menou & Tabachnik 2003; Sándor et al. 2007) under the assumption of coplanar orbits. In Funk et al. (2008, 2010), a possible inclination of the fictitious test-planet was considered in a parametric study of the dynamical stability of potential additional terrestrial planets inside the habitable zone¹ of nearby (within 30 pc) extrasolar planetary systems². In this present paper, we highlight the influence of the eccentricity and the mass of the gas giant on the evolution of fictitious, inclined, earth-mass companions that are initially on circular orbits. These results are analyzed thoroughly in the present contribution.

On the one hand, we consider gas giants on circular orbits and show that our results are in agreement with the well-known results of the Kozai mechanism (Kozai 1962) as expected. On the other hand, we study the influence of increasing the eccentricity of the gas giant on the Kozai mechanism. In particular, we pay special attention to a stability region close to 35° of inclination of the fictitious planet, which – to our knowledge – has not yet been reported in the literature. Since we used the dimensionless form of the equations of motion, by scaling the results of our parametric study, we can make it valid for many real extrasolar planetary systems.

The paper is organized as follows. In Sect. 2, we describe our parametric study of the restricted three-body problem (here the fictitious planet is regarded as massless). The results of our simulations are detailed in Sect. 3.1 for circular gas giants and in Sect. 3.2 for eccentric gas giants. An application to extrasolar planetary systems is discussed in Sect. 4. Finally, our results are summarized in Sect. 5.

2. The dynamical model and the methods

To study the dynamical stability, we applied the model of the spatial elliptic restricted three-body problem. We integrated the dimensionless equations of motion. An obvious advantage of using these equations is that the results are independent of the exact value of the semi-major axis of the primary, i.e. the gas giant in single extrasolar planetary systems. The unit of length was chosen such that the separation of the star and the gas giant (the primaries) is unity, i.e. the semi-major axis of the gas giant a₁ = 1 in all computations. We defined the unit of mass to be the sum m₀ + m₁ = 1, where m₀ and m₁ are the masses of the star and the gas giant, respectively. We choose the unit of time such that k²(m₀ + m₁) = 1, where k is the gravitational constant. The orbital plane of the primaries was used as a reference plane, in which the line connecting the primaries at t = 0 defines a reference x-axis. The longitude of the ascending node Ω is the angle between the line of nodes and the x axis. The argument of pericenter ω is measured between the line of nodes and the radius vector of the pericenter. The values of Ω = 0, ω = 0, and the mean anomaly M = 0 places the test particle on the x axis at a distance of a(1 − e) from the astrocenter, where a is the semi-major axis and e is the eccentricity of the test particle. The mass parameter of the system is μ = m₁/(m₀ + m₁).

As a dynamical model, we used the aforementioned elliptic restricted three-body problem, hence the results are valid for massless test-planets that move in the gravitational field of a star and a gas giant. Moreover, we tested whether the results are still valid for low-mass planets (m = 1,...,10   M _⊕ ), where M _⊕  is the mass of the Earth, and applied them to the extrasolar system HD 154345. In Table 1, we summarize the initial conditions (semi-major axis (a), eccentricity (e), inclination (i), argument of the pericenter (ω), longitude of the ascending node (Ω), and the mean anomaly (M)) for the gas giant and the test-planet. The studies were performed for 10⁵ periods of the primary.

Three mass ratios were considered μ = 0.0005,0.001,0.003, which are among others equivalent to a giant planet of respectively 0.5,1, and 3 Jupiter mass orbiting a solar mass star. For the case of μ = 0.001, we performed calculations for eccentricities of the gas giant (e₁) up to 0.5 and for μ = 0.0005 even up to 0.9, hence we obtained a quite good estimate of the stable region even for high eccentric orbits.

All calculations were performed with the Lie integration method (see Hanslmeier & Dvorak 1984; Lichtenegger 1984), which can deal with high eccentricities and close encounters between the bodies. In addition we performed some test-calculations with the Bulirsch-Stoer integration method in combination with the LCI (see Fig. 3), an efficient variable-time step algorithm. The most important feature of the Bulirsch-Stoer algorithm for N-body simulations is that it is capable of keeping an upper bound on the local errors introduced by taking finite time steps when adaptively reducing the step size as interactions between the particles increase in strength. The parameter ϵ, which controls the accuracy of the integration, was determined by preliminary test runs and was set to 10^-13.

We study the orbital behaviour of massless bodies representing hypothetical exo-Earths using the maximum eccentricity (ME) and the Lyapunov characteristic indicator (LCI). The ME is defined as $\mathrm{ME}\mathrm{=}\underset{\mathit{t}\mathrm{\in }\mathrm{\left[}\mathrm{0}\mathit{,}{\mathrm{10}}^{\mathrm{5}}\mathrm{\right]}{\mathit{T}}_{\mathrm{1}}}{\max }\mathrm{\left(}\mathit{e}\mathrm{\right(}\mathit{t}\mathrm{\left)}\mathrm{\right)}\mathit{,}$(1)where T₁ is the period of the primary and ME gives the maximum value of eccentricity reached during the whole computation time. For very low values (e < 0.2) and very high values (e > 0.8) of ME, a good agreement can generally be found with the results of chaos indicators (see e.g. Dvorak et al. 2003; Süli et al. 2005; Nagy et al. 2006; Süli et al. 2008). For intermediate values of ME a confident distinction between stable and unstable motion is not possible, since high eccentric orbits (0.2 < e < 0.8) can be stable in the long-term.

As a chaos indicator, we used in this study the LCI, which is a finite time approximation of the largest Lyapunov characteristic exponent (LCE) described in detail in e.g. Froeschlé (1984). The LCI is therefore defined as $\mathrm{LCI}\mathrm{\left(}\mathit{t,}{\mathit{x}}_{\mathrm{0}}\mathit{,}{\mathit{\xi }}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathit{t}}\log \mathrm{\parallel }\mathit{\xi }\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{\parallel }\mathit{,}$(2)where x₀ is the initial condition of the orbit and ξ(t) is the solution of the linearized variational equations. LCI(t,x₀,ξ₀) measures the mean rate of divergence of the orbits.

3. Results

Fig. 1 ME plots in the (a − i) plane. The color code corresponds to the maximum eccentricity – red means low eccentricity and therefore stable motion, green and blue correspond to high eccentricities and therefore unstable motion (see color code). The plots are for different mass ratios and eccentricities of the gas giant: first column: μ = 0.0005, e1 = 0,0.1,0.3 (up to down); second column: μ = 0.001, e1 = 0,0.1,0.3 (up to down); third column: μ = 0.003, e1 = 0,0.1,0.3 (up to down). The black dots indicate the orbits for which we also calculated LCI values (see Fig. 3).

In Fig. 1, we show the results for the three chosen mass ratios and three different eccentricities of the gas giant. All graphs are ME plots, where red means low eccentricity hence stable motion, and blue corresponds to high eccentricities hence to possibly unstable motion. We note that the step-like behaviour in the graphs of Fig. 1 occurs because of the grid of used initial conditions – in particular because of the quite low resolution in the inclination (see Table 1). All plots show the region between the star and the gas giant, where we plotted the inclination of the test-planet versus its semi-major axis. The first row gives the results for e₁ = 0.0 (circular case), and the second and the third for 0.1 and 0.3, respectively (elliptic case). In all cases, the three plots correspond to mass ratios of 0.0005, 0.001, and 0.003 (from left to right). The circular problem is described in Sect. 3.1, while the elliptic problem is detailed in Sect. 3.2.

3.1. Circular gas giant

Fig. 2 Level curves of constant averaged Hamiltonian on the phase space (ecosω,esinω) for the circular restricted three-body problem. The following initial conditions for the test-planets were used for both plots: a = 0.15, i = 15° for e = 0 (upper graph), i = 50° for e = 0 (lower graph).

When the outer giant planet moves on a circular orbit (Fig. 1, panels a1−3), the global dynamics are regular with low maximum eccentricities (e ≤ 0.01) for a test-planet at small inclination values (i ≤ 35°). This stable red-colored region shrinks as the mass ratio increases. For highly inclined orbits, larger values of the test-planets’ eccentricity are observed. Only close-in test-planets show low ME values (during 100 000 periods) for highly inclined orbits. But we note that no additional effect such as tidal effects or general relativity corrections are considered here, which could modify the behaviour of close-in test-planets. Figure 1 also shows the influence of some mean motion resonances, e.g. at 0.48 (3:1), 0.63 (2:1), 0.76 (3:2), and 0.83 (4:3). The dynamical behaviour described above is mainly due to the influence of the Kozai mechanism (Kozai 1962; Lidov 1962). Circular orbits are indeed a well-known secular equilibrium of the restricted three-body problem, which is reduced to two degrees of freedom after short-period averaging (i.e. averaging with respect to the mean anomalies of both bodies) and node reduction (e.g. Poincaré 1892; Malige et al. 2002). Assuming that the outer giant planet is on a circular orbit, this problem is integrable, and its dynamics can be represented in the phase space (ecosω,esinω) (see Thomas & Morbidelli 1996, for more detail). Such a representation is given in Fig. 2 where we plot the curves of a constant Hamiltonian, which correspond to the trajectories of the test-planet.

For small inclinations (Fig. 2 (upper graph), i = 15°), circular orbits of the massless body correspond to a stable equilibrium and no variation in the eccentricity is possible. This explains the red-colored region in Fig. 1. For higher inclinations (Fig. 2 (lower graph), i = 50°), the point e = 0 becomes an unstable equilibrium point, and a separatrix divides the phase space into three parts: two regions are characterized by the libration of ω around 90° or 270° and the third one by the circulation of this angle. The two stable equilibria created by the bifurcation of the equilibrium in the circular orbit are referred to as Kozai equilibria. This change in the stability of the central equilibrium induces large variations in the eccentricity of a test-planet initially on a circular orbit, since its real motion (short periods included) remains close to the separatrix of the reduced problem. Moreover, by increasing the inclinations of the test-planet, the separatrix will reach higher values of eccentricity. The increasing amplitudes of these variations are shown by yellow, green, and light blue colors in Fig. 1.

Fig. 3 LCI-values for some single orbits in the systems with a mass ratio μ = 0.001 (indicated by black dots in Fig. 1). The initial conditions for the integrated orbits are summarized in the upper right corner of each plot. LCI-values going to 0 indicates stable motion, while LCI-values going to a value higher than 0 indicates chaotic motion.

As one can see in Fig. 1, the critical inclination corresponding to this change in stability decreases with increasing semi-major axis ratios. The critical values observed in this numerical study correspond to those calculated analytically by Kozai (1962, see Fig. 1 of his paper).

We checked whether the ME results are in a good agreement with the LCI. Test runs showed that the chaotic character of the orbits close to the unstable equilibrium is reliably confirmed by the LCI indicator but only when considering quite long integration times. This is because the Kozai perturbations are secular (which means that they operate on extremely long timescales), and it is very difficult for a chaos indicator associated with a non-symplectic integrator, like Burlish-Stoer, to exhibit these perturbations in a short computation time, as explained in Libert et al. (2010). To establish the correct orbit character in an acceptable computation time, we point out that both the orbit and the variational equations should be computed simultaneously in a symplectic way, as in the case of the global symplectic integrator of Libert et al. (2010). To check the stability border indicated by the ME plots, we decided to follow the evolution of some well-chosen orbits on an adequatly long integration time (10⁷ years). In Fig. 3, we show, as an example, the evolution of the LCI for four (e₁ = 0.0, upper graph) and two (e₁ = 0.3, lower graph) orbits. While for gas giant eccentricity of 0.3 (e₁ = 0.3) the difference between stable and chaotic motion can be seen in an integration timespan of approximately 5 × 10⁵ years, the difference between stable and chaotic motion for e₁ = 0.0 is not obvious before 10⁶ years. Nevertheless, these examples show that the results of the ME and the LCI are in a good agreement when an adequatly long integration time is chosen for the LCIs (and a sufficiently low precision of the variable-time step integrator).

3.2. Elliptic gas giant

Fig. 4 Similar to Fig. 1 for μ = 0.0005 and e1 = 0.5.

As the eccentricities of the gas giant increase (Fig. 1 panels b1−3 and c1−3), the region of nearly circular motion (red color region of the first row) shrinks and almost disappears at an eccentricity of the giant as high as 0.5 (see Fig. 4). In other words, the unstable area (dark blue) caused by the proximity of the test-planet to the gas giant becomes larger as the giant eccentricities increase. Again one can discern several mean motion resonances as well as the disturbing influence of the change in the stability of the stable equilibrium.

For inclinations lower than the critical values mentioned in the previous section, the dynamics appear different in the sense that the maximal eccentricity values are no longer equal to zero but can reach values as high as 0.2, depending on the proximity of the test-planet to the giant planet: the higher the semi-major axes ratio, the higher the variation in eccentricity of the test-planet. In addition, close to 35° there appears to be a region where the eccentricity hardly varies, which is still visible even for an eccentricity of 0.5 (see Fig. 4). To learn more about this area of small variations in eccentricity, we investigated in more detail a single orbit at inclination 35°. The orbital evolution is shown in Fig. 5. The semi-major axis of the test-planet was set to 0.15, the eccentricity of the gas giant to e₁ = 0.5, and the mass ratio to μ = 0.0005. A coupling between the eccentricity and the inclination of the test-planet occurred, and ω started librating alternately around 90° and 270°. Similar behaviour is observed for test-planets inclined 35°, independently of the non-zero eccentricity of the gas giant.

Fig. 5 The left graph shows the eccentricity (red curve) and the inclination (blue curve – note that we substracted 35° from the inclination values for clearer visibility), while the right graph shows the argument of the pericenter evolutions of a typical test-planet orbit around 35°. The eccentricity of the gas giant is set to e1 = 0.5 and the semi-major axis of the test-planet to 0.15.

The behaviour mentioned above has – to our knowledge – not yet been reported in the literature. From the dynamical point of view, it is of special interest because the elliptic restricted three-body problem has not yet been extensively studied. When the eccentricity of the giant body is non-zero, then the restricted three-body problem cannot be reduced to one degree of freedom and the analytical study of the dynamics becomes more complicated. Our first investigations have found that circular orbits are not necessary the most stable equilibria of this two degree of freedom problem and determined the typical variations in eccentricity observed for test-planets at small inclinations. Our preliminary explanation of the structure around 35°, is that it is associated with a secular resonance whose resonant angle is ω − Ω. A detailed analytical study of the elliptic restricted three-body problem in a Hamiltonian formalism is currently in preparation.

Furthermore, the plots b and d of Fig. 3 for a very high eccentricity of e₁ = 0.3 illustrate that the ME and LCI results are in good agreement. From the results described above, it appears that inclined massless test-planets may exist in stable configurations with gas giants on circular or elliptical orbits. The higher the eccentricity of the gas giant, the smaller the possible range in semi-major axis for the test-planet: for a gas giant situated at 1, the stable range for small inclinations, derived with the LCI, is a₀₀ =  [0.1,0.84]  for e₁ = 0, a₀₁ =  [0.1,0.71]  for e₁ = 0.1, and a₀₃ =  [0.1,0.54]  for e₁ = 0.3. Therefore we can assume that the stability border (a_sb) can be calculated with the formula a_sb = a₀₀·(1 − e₁), which gives the results a₀₁ = 0.76 and a₀₃ = 0.59. Both results lie quite close to the actual values, so the above-mentioned formula can be used at least for a first approximation of the stability border. Furthermore, for gas giants on eccentric orbits, test-planets suffer from limited variations in eccentricity, depending on the value of their semi-major axis (the higher the semi-major axes ratio, the higher the variation in eccentricity of the test-planet) and on their inclination (close to 35°, there is a region of small variation in eccentricity). In the next section, we show that these results may be generalized to Earth-mass planets.

4. Application to extrasolar planetary systems

Fig. 6 ME plots (see color code) of the extrasolar planetary system HD 154345 for different eccentricities of the known planet (from left to right: e1 = 0.044,0.09,0.244). The a-i-grid corresponds to the initial conditions of fictitious planets, covering the region of the HZ. In the upper row, we show the calculations obtained by the restricted three-body problem (test-planet is massless) and draw in the width of the HZ (white lines) for e1 = 0.0, 0.1 and 0.2 (from left to right). In the lower row, we show the calculations obtained by the full three-body problem study for e1 = 0.044, 0.09, 0.244 (from left to right) and a fictitious planet of 2 M ⊕ .

The previous results can be applied for example to the extrasolar planetary system HD 154345. All orbital parameters are summarized in Table 2. In the first two columns, we give the name and the mass of the host star. The following columns show the mass, the semi-major axis, and the eccentricity of the discovered planet. The sixth column gives the mass ratio. The two last columns provide the width of the habitable zone (HZ) according to Kaltenegger et al. (2010) and the width of the HZ when the semi-major axis of the gas giant is normalized to 1. In a previous study (Funk et al. 2010), many extrasolar planetary systems were investigated in detail, by means of a parameter study, where we varied the semi-major axis, the eccentricity and the mass of the known extrasolar planet, as well as the inclination of the fictitious planets. All calculations were performed in the full three-body problem with fictitious planets having 2 M _⊕ . The used initial conditions are summarized in Table 3. Three maximal eccentricity plots of the system HD 154345 corresponding to eccentricities 0.044, 0.09, and 0.244 are displayed in Fig. 6 (bottom panels).

To compare these results with those of the restricted three-body problem discussed in previous section, we choose the plots corresponding to the parameters mass ratio μ = 0.001 and eccentricities of 0.0, 0.1, and 0.2, since these values are close to the observed ones. These three plots are given in Fig. 6 (upper row), where we draw, for each one, the width of the normalized HZ (white lines).

As one can see both results agree quite well, which indicates that our results for the restricted three-body problem can be applied to extrasolar planetary systems as long as the mass of the fictitious planet is not too large. Therefore, we find that the possibility of extrasolar giant planets possessing an inclined earth-mass companion is non-negligible. The long-term evolution in eccentricity of this small body would depend on its inclination and semi-major axis, as well as on the eccentricity of the extrasolar giant planet, as described in Sect. 3. Furthermore, the HZ of some extrasolar systems may belong to the range of semi-major axis identified as stable regions in the previous section. At least for the extrasolar planetary system HD 154345, which is likely to posses a 35° inclined, nearly circular, earth-mass companion in the HZ.

5. Conclusions

We have investigated the long-term evolution of inclined orbits in the circular and elliptic restricted three-body problem. We have integrated the evolution of a grid of initial conditions of the massless test-planets (semi-major axis versus inclination) for different eccentricities and masses of the gas giant. Our results have shown that a giant planet on circular or elliptic orbits may have an inclined massless companion. The higher the eccentricity of the gas giant, the smaller the possible range in semi-major axis for the massless test-planet: for a gas giant situated at a₁ = 1, the range is approximately a =  [0.1,0.8]  for e₁ = 0, a =  [0.1,0.5]  for e₁ = 0.3, and a =  [0.1,0.2]  for e₁ = 0.5 for inclinations up to approximately 35°.

For a circular giant, our simulations illustrate the well-known results of a restricted three-body problem: at small inclinations, circular orbits for the test-planet are in a stable equilibrium, which bifurcates at higher inclinations into the Kozai equilibria. For giants on eccentric orbits, our results show that the dynamics of the problem seems to be quite similar. Nevertheless, we have shown that an Earth-mass companion is affected by limited variations in eccentricity, depending on both the value of its semi-major axis (the higher the semi-major axis ratio, the larger the variation in eccentricity of the test-planet) and its inclination.

We have paid particular attention to a stability region close to a fictitious planet inclination of 35°, which – to our knowledge – has not yet been reported in the literature and our results may indicate that close to 35° a region of small variation in eccentricity appears.

We finally checked that our results for massless test-planets can be applied to extrasolar planetary systems as long as the mass of the fictitious planet is not too large. Hence, we have found that extrasolar giant planets may very well possess an inclined Earth-mass companion. Furthermore, the HZ of some extrasolar systems may belong to the range in semi-major axis identified as a stable region. We showed that the extrasolar planetary system HD 154345 at least is likely to possess a 35° inclined, nearly circular, earth-mass companion in the HZ.

Acknowledgments

B. Funk wants to acknowledge the support by the Austrian FWF Erwin Schroedinger grant no. J2892-N16. The work of A.-S. Libert is supported by an FNRS post-doctoral research fellowship. Á. Süli wants to acknowledge the support by the Hungarian Scientific Research Fund PD-75508. E. Pilat-Lohinger wants to acknowledge the support by the Austrian FWF project No. P19569-N16. The European Union and the European Social Fund have provided financial support to the project under the grant agreement No. TÁMOP-4.2.1/B-09/1/KMR-2010-0003. We thank the referee for his constructive suggestions.

References

All Tables

All Figures

Fig. 1 ME plots in the (a − i) plane. The color code corresponds to the maximum eccentricity – red means low eccentricity and therefore stable motion, green and blue correspond to high eccentricities and therefore unstable motion (see color code). The plots are for different mass ratios and eccentricities of the gas giant: first column: μ = 0.0005, e1 = 0,0.1,0.3 (up to down); second column: μ = 0.001, e1 = 0,0.1,0.3 (up to down); third column: μ = 0.003, e1 = 0,0.1,0.3 (up to down). The black dots indicate the orbits for which we also calculated LCI values (see Fig. 3).

In the text

	Fig. 2 Level curves of constant averaged Hamiltonian on the phase space (ecosω,esinω) for the circular restricted three-body problem. The following initial conditions for the test-planets were used for both plots: a = 0.15, i = 15° for e = 0 (upper graph), i = 50° for e = 0 (lower graph).
In the text

	Fig. 3 LCI-values for some single orbits in the systems with a mass ratio μ = 0.001 (indicated by black dots in Fig. 1). The initial conditions for the integrated orbits are summarized in the upper right corner of each plot. LCI-values going to 0 indicates stable motion, while LCI-values going to a value higher than 0 indicates chaotic motion.
In the text

	Fig. 4 Similar to Fig. 1 for μ = 0.0005 and e1 = 0.5.
In the text

	Fig. 5 The left graph shows the eccentricity (red curve) and the inclination (blue curve – note that we substracted 35° from the inclination values for clearer visibility), while the right graph shows the argument of the pericenter evolutions of a typical test-planet orbit around 35°. The eccentricity of the gas giant is set to e1 = 0.5 and the semi-major axis of the test-planet to 0.15.
In the text

Fig. 6 ME plots (see color code) of the extrasolar planetary system HD 154345 for different eccentricities of the known planet (from left to right: e1 = 0.044,0.09,0.244). The a-i-grid corresponds to the initial conditions of fictitious planets, covering the region of the HZ. In the upper row, we show the calculations obtained by the restricted three-body problem (test-planet is massless) and draw in the width of the HZ (white lines) for e1 = 0.0, 0.1 and 0.2 (from left to right). In the lower row, we show the calculations obtained by the full three-body problem study for e1 = 0.044, 0.09, 0.244 (from left to right) and a fictitious planet of 2 M ⊕ .

In the text