A&A 426, 353-365 (2004)
DOI: 10.1051/0004-6361:20040390

Stability of terrestrial planets in the habitable zone of Gl 777 A, HD 72659, Gl 614, 47 Uma and HD 4208

N. Asghari¹ - C. Broeg¹ - L. Carone¹ - R. Casas-Miranda¹ - J. C. Castro Palacio¹ - I. Csillik¹ - R. Dvorak² - F. Freistetter² - G. Hadjivantsides¹ - H. Hussmann¹ - A. Khramova¹ - M. Khristoforova¹ - I. Khromova¹ - I. Kitiashivilli¹ - S. Kozlowski¹ - T. Laakso¹ - T. Laczkowski¹ - D. Lytvinenko¹ - O. Miloni¹ - R. Morishima¹ - A. Moro-Martin¹ - V. Paksyutov¹ - A. Pal¹ - V. Patidar¹ - B. Pecnik¹ - O. Peles¹ - J. Pyo¹ - T. Quinn⁵ - A. Rodriguez¹ - C. Romano¹³ - E. Saikia¹ - J. Stadel⁴ - M. Thiel³ - N. Todorovic¹ - D. Veras¹ - E. Vieira Neto¹ - J. Vilagi¹ - W. von Bloh¹ - R. Zechner¹ - E. Zhuchkova¹

1 - Helmholtz-Institute for Supercomputational Physics, Am Neuen Palais 10, 14469 Potsdam, Germany
2 - Universitätssternwarte Wien, Türkenschanzstr. 17, 1180 Wien, Austria
3 - Institut für Physik, Universität Potsdam, Am Neuen Palais 10, 14469 Potsdam, Germany
4 - Institut für Theoretische Physik, Universität Zürich, Winterthurerstr. 190, 8057 Zürich, Switzerland
5 - Department of Astronomy, University of Washington, Seattle 98195-1580, USA

Received 4 March 2004 / Accepted 16 June 2004

Abstract
We have undertaken a thorough dynamical investigation of five extrasolar planetary systems using extensive numerical experiments. The systems Gl 777 A, HD 72659, Gl 614, 47 Uma and HD 4208 were examined concerning the question of whether they could host terrestrial-like planets in their habitable zones (HZ).

First we investigated the mean motion resonances between fictitious terrestrial planets and the existing gas giants in these five extrasolar systems. Then a fine grid of initial conditions for a potential terrestrial planet within the HZ was chosen for each system, from which the stability of orbits was then assessed by direct integrations over a time interval of 1 million years. For each of the five systems the 2-dimensional grid of initial conditions contained 80 eccentricity points for the Jovian planet and up to 160 semimajor axis points for the fictitious planet. The computations were carried out using a Lie-series integration method with an adaptive step size control. This integration method achieves machine precision accuracy in a highly efficient and robust way, requiring no special adjustments when the orbits have large eccentricities.

The stability of orbits was examined with a determination of the Rényi entropy, estimated from recurrence plots, and with a more straightforward method based on the maximum eccentricity achieved by the planet over the 1 million year integration. Additionally, the eccentricity is an indication of the habitability of a terrestrial planet in the HZ; any value of e>0.2 produces a significant temperature difference on a planet's surface between apoapse and periapse. The results for possible stable orbits for terrestrial planets in habitable zones for the five systems are: for Gl 777 A nearly the entire HZ is stable, for 47 Uma, HD 72659 and HD 4208 terrestrial planets can survive for a sufficiently long time, while for Gl 614 our results exclude terrestrial planets moving in stable orbits within the HZ.

Studies such as this one are of primary interest to future space missions dedicated to finding habitable terrestrial planets in other stellar systems. Assessing the likelihood of other habitable planets, and more generally the possibility of other life, is the central question of astrobiology today. Our investigation indicates that, from the dynamical point of view, habitable terrestrial planets seem to be compatible with many of the currently discovered extrasolar systems.

Key words: stars: individual: Gl 777A - stars: individual: HD 72659 - stars: individual: Gl 614 - stars: individual: 47 Uma - stars: individual: HD 4208 - stars: planetary systems

1 Introduction

In this study the extrasolar systems Gl 777 A, HD 72659, Gl 614, 47 Uma and HD 4208 were examined to assess whether they could host additional terrestrial-like planets in their habitable zones (=HZ).

Since the discovery of the first extrasolar planetary system about 10 years ago (Mayor & Queloz 1995), a major point of dynamical investigations has been the determination of stable regions in extrasolar planetary systems, where additional planets on stable orbits could exist. Today we know about 105 planetary systems with 120 planets, where 13 systems have more than one planet (both confirmed and unconfirmed cases).

Because of the observational methods in use, our knowledge of extrasolar planets is highly biased: almost all planets were detected via ground-based measurements of the central star's radial velocity. This method favors the detection of very large (Jupiter-sized) planets that move very close to the main star. Several space missions are planned (i.e. COROT, KEPLER, TPF, DARWIN) that will be able to find planets with a much smaller mass (comparable to that of the Earth) by using other detection methods. The first results of these missions will only be available after more than a decade, however, theoretical studies can currently make predictions on the existence of stable terrestial planets within the known systems. This serves to establish testable theories of planet formation and evolution, as well as aiding the future searches for such habitable terrestial planets.

When one looks at the distribution of semimajor-axes of all known exosolar planets, one can see that they are located from $\sim$0.02 AU to $\sim$6.5 AU. These distances cover the so-called "habitable zone'' (=HZ) of a typical main sequence star, ranging from $\sim$0.7 AU to $\sim$1.5 AU (Kasting et al. 1993). Some of these large planets move in the region around 1 AU that overlaps with the HZ, and leave no room for possible terrestrial planets orbiting the host star within the HZ - only a large satellite (like Titan) or a Trojan-like planet could have all the properties which are necessary for orbital stability and also the development of a stable atmosphere in such cases. When the large planets move well outside (or inside) the HZ they will perturb any additional undetected planet that may move within this zone.

For extrasolar systems with only one planet the possible existence of additional planets has been investigated (e.g. Dvorak et al. 2003b; Pál & Sándor 2003). In an extensive paper (Menou & Tabachnik 2003) the stability of possible terrestrial planets was studied for all extrasolar systems and a classification according to their stability was established, but there the role of resonances was underestimated. Resonances - mean-motion resonances and secular resonances - can either stabilize the motion of a planet and protect it from close encounters and collisions or, in other cases, intensify the gravitational perturbations and therefore destroy the stability of the orbits. Resonances between the orbits of the observed large planets have been studied (Callegari et al. 2002; Hadjidemetriou 2002; Érdi & Pál 2003), but a detailed study of the perturbations on possible additional planets due to resonances is still lacking. Thus, any study of the dynamics inside the HZ of an extrasolar system has also to include an investigation of the resonances.

Determinations of the dynamical stability of multiple planetary systems (Kiseleva et al. 2002a,b; Beaugé & Michtchenko 2003) and of planets in double stars (Pilat-Lohinger & Dvorak 2002a; Pilat-Lohinger et al. 2002b, 2003; Dvorak et al. 2003a; Lammer et al. 2003) have been done recently.

In this article we investigate the dynamics inside the habitable zone of five different extrasolar systems: Gl 777 A, 47 Uma, HD 72659, Gl 614 and HD 4208. Except 47 Uma, all these systems consist of the central star and a Jovian planet. According to the spectral type of the main star, the habitable zone was determined; a detailed analysis of the mean motion resonances inside the HZ was carried out and the full width of the HZ was investigated according to the stability of possible additional terrestrial planets. This article is organized as follows: Sect. 2 describes the five extrasolar systems for which we performed a dynamical study for additional planets, in Sect. 3 we introduce the dynamical models and the methods used to analyze the time series data resulting from our numerical experiments. Section 4 shows how we dealt with the terrestrial planet orbits in the major mean-motion resonances of the Jovian planet(s). Subsequent sections focus on each of the 5 systems individually, detailing our investigations and results. We conclude with a summary of these results and their implications, including comparison to earlier, less detailed studies by Menou & Tabachnik (2003). The appendicies include concise descriptions of the analytical and numerical methods used, and may serve as a useful reference to the reader.

   
2 Description of the extrasolar systems

Table 1 shows the main characteristics of the five extrasolar planetary systems which we investigated in this study. Except the system 47 Uma (which hosts two planets), all consist of a giant planet moving outside the HZ of the star; this is a situation similar to our Solar System (SS) when we take into account only Jupiter. Thus, an appropriate dynamical model is the restricted three body problem (star + planet + massless terrestrial planet). A different dynamical model was taken for the system 47 Uma, which looks quite similar to our SS when only Jupiter and Saturn are considered. Both 47 Uma and the SS have the more massive planet orbiting closer to the star with the second planet at approximately twice the distance. Furthermore, the 47 Uma planets have roughly a 3:1 mass ratio and exhibit low eccentricity orbits very much like Jupiter and Saturn in our own SS. The restricted four body problem, known to be a good model for an asteroid in the SS (taken as a massless body), additionally considers the secular resonances acting on the motion of a massless test body and will also serve as a good model for terrestrial planets within 47 Uma.

Table 1: Parameters of the investigated systems - obtained from http://www.obspm.fr/planets (1. Sep. 2003).

In a wider sense we can say that all the central stars are solar like stars with masses slightly smaller or larger than the Sun. The spectral types for these main sequence stars (G0 to K0) allow HZs in the range between 0.5 AU and 1.5 AU. With the exception of HD 4208 and 47 Uma all the Jovian planets' orbits have significantly larger eccentricities than that of Jupiter with sometimes very large uncertainties. On the other hand, we expect the determined semimajor axes to be rather precise. Our parameter study will take into account this discrepancy in the precision of the measured orbital elements of the exoplanets under study and will be detailed in the following section. The minimum masses of the Jovian planets lie exactly in the range we expect for gas giants. With regard to the unknown inclinations of these extrasolar systems, we always took the masses of Table 1 as correct, assuming that the orbital plane of the system is seen edge-on. It should be kept in mind that the true masses of these planets could be significantly larger. Statistically we can say that in 5 of 6 cases the inclination will change the mass by a factor smaller than four compared to the one given in the tables. Recent computations have shown that with these larger masses the strenghts of the perturbations will not significantly change (Sandor et al. 2004).

   
3 Simulation method and stability analysis

The availability of a supercomputer with 128 processors for this investigation enabled the direct computation of orbits to assess stability. Futhermore, the extent of the computational resources favoured the use of a very precise numerical integration scheme, the Lie-integration method, which is free from numerical difficulties experienced by other (lower order) techniques, particularly in the case of highly eccentric orbits. The Lie-integration method uses an adaptive stepsize and is quite precise and fast, as has been shown in many comparative test computations with other integrators such as Runge-Kutta, Bulirsch-Stoer or symplectic integrators. Although symplectic integrators are very effective when eccentricities remain small, the Lie integrator is a better choice in studies such as this one, where very large eccentricity orbits are explored. We detail this integration method in Appendix B, and further details can be found in Hanslmeier & Dvorak (1984) as well as Lichtenegger (1984).

The uncertainties in the orbital elements derived from observations make it sufficient to consider only Newtonian forces in a dynamical model of two massive bodies (the central star and the discovered planet) and a massless test body, representing a potential terrestrial planet. As mentioned previously, this system corresponds to the three dimensional restricted 3-body problem. For the 47 Uma system where two Jovian planets are known, we followed the dynamics of a restricted 4-body problem with a central star, both discovered Jovian planets, and a massless test body representing a potential terrestrial planet. This increased the computing time considerably for 47 Uma compared to the other restricted 3-body systems.

For the integrations we use a fine 2-dimensional grid of initial conditions; for the terrestrial planets varying the distance to the central star, exploring the HZ, and for the (inner - 47 Uma) Jovian planet varying the eccentricity of its orbit. We vary the eccentricity of the Jovian planet because the values derived from observations have large uncertainties, and yet this is an important parameter for stability of other planets within the system. For the nodes and the pericenter of the planets we took the ones published by the observers; in the cases where these values are unknown we set them to zero. Table 2 outlines the ranges for our initial conditions. We always started the initial orbits of the terrestrial planets as being circular; they turned out to deform quite quickly, exploring a range of more eccentric orbits so that varying their initial eccentricities would seem to be unnecessary. The eccentricity of the fictitious planet's orbit, guided by the formation scenario for rocky terrestrial planets (Richardson et al. 2000) was set to zero; because of technical reasons the inclination was set to almost zero ( $i=0\hbox{$.\!\!^\circ$ }001$). Nevertheless the computations were carried out in three-dimensions.

The integration of the terrestrial planet's orbit was always started on the connecting line between the star and the Jovian planet. An integration time of 1 million years was chosen, which is long enough to reveal the stability character of the planets for the 5 systems of this study. The default parameter grid of 80 by 80 which we used for the 5 systems leads to a total integration time of several $10^{{\rm 10}}$ years, when we also include test computations and the fact that some of the orbits were integrated for an extended $10^{{\rm 7}}$ years.

Table 2: Initial semimajor axes of the fictitious terrestrial planets $a_{{\rm T}}$ and the eccentricities of the Jovian planets $e_{{\rm G}}$.

For the analysis of stability we used a straightforward check based on the eccentricities. For this we examined the behavior of the eccentricity of the terrestrial planets along their orbit and used the largest value as a stability criterion; in the following we call it the maximum eccentricity method (=MEM). This simple check has already been used in other studies of this kind and was found to be quite a good indicator of the stability character of an orbit (Dvorak et al. 2003a). An orbit was deemed unstable when the eccentricity exceeded a value of e=0.5, after which we stopped further computation. In all our former studies this stability limit turned out to be an appropriate tool because all the terrestrial planet orbits with e=0.5 turned out to suffer, sooner or later, from a close encounter with the large planet, causing the terrestrial planet to escape (Dvorak et al. 2003a). Although some orbits in multiple exosolar systems have eccentricities larger than 0.5 their special configuration allows non-crossing orbits. By placing additional fictitious planets between their orbits with eccentricities larger than 0.5 they always suffer in the long run from close encounters leading to unstable orbits.

For the habitability of a planet we also used an additional criterion based directly on the eccentricity of the orbit within the HZ. This was done in order to take into account the variations in the "solar'' insolation on the surface of the terrestrial planet. A good approximation (Lammer 2004; private communication), requiring e<0.2, is sufficient to keep this variation in insolation small enough during an orbit.

On the other hand we computed the Rényi entropy, a measure which is often used in nonlinear dynamics to determine how predictable an orbit is. These values are comparable to the Fast Lyapunov Indicators (=FLI) introduced by Froeschlé et al. (1997). It is slower than determining the FLIs but the advantage of the entropy method is that it can be used a posteriori. We estimated the Rényi entropy by means of Recurrence Plots (RPs), a tool of data analysis that has found numerous applications in many different fields in recent years (Webber et al. 1994; Marwan et al. 2002; Thiel et al. 2003). RPs were initially introduced to simply visualize the behaviour of trajectories in phase space. The distances of every pair of points of the system's trajectory are represented in a 2-dimensional binary matrix. In this way RPs yield different patterns depending on the system's character. One can introduce different measures that quantify the obtained structures. It was shown in Thiel et al. (2003) that it is also possible to estimate the Rényi entropy based on the distribution of diagonal lines obtained in the RP. The details of this second approach are covered in Appendix A. The two methods were used complementary because on one side the MEM is the right tool to assure the stability of an orbit in the sense of being habitable (the eccentricities stay small) whereas on the other side the entropy method gives a direct measure of chaos and unveals the resonance structure of phase space in more detail (compare e.g. Figs. 10 and 11 for HD 72659 and Figs. 15 and 16 for HD 4208).

   
4 The stability within resonances

We know from our planetary system about the importance of resonances between the mean motions of two bodies. Because all our planets have only small orbital eccentricities there are different simplified models available which can provide interesting results concerning the structure of phase space inside these resonance (e.g. Engels & Henrard 1994; Malhotra 1998). Most of these models can be used for the natural satellites (e.g. Malhotra & Dermott 1990) but in the case of extrasolar planets with large eccentricities theoretical considerations are not yet obtained. Therefore we choose to use a numerical approach for each system separately which we describe in the following.

For the investigation of the resonances, we choose initial conditions placed in the most relevant mean-motion resonances (=MMRs) of the fictitious planet with the Jovian planet inside but also outside the HZ. These resonances were checked for stability in 8 different positions of the terrestrial planet (corresponding to $M=0^{{\rm\circ}},45^{{\rm\circ}},90^{{\rm\circ}}, 135^{{\rm\circ}},180^{{\rm\circ}},225^{{\rm\circ}},270^{{\rm\circ}},315^{{\rm\circ}})$. Additionally the computations were carried out with the Jovian planet initially placed at the apoastron and periastron. For a detailed list of the resonant positions that were investigated for each system see Table 3.

Table 3: Stability of orbits in mean motion resonances. The numbers give the stable orbits according to the 8 different initial conditions.

	Figure 1: Schematic view of the stability of orbits in the resonances of 1st and 2nd order in HD 4208. Full (empty) circles stand for stable (unstable) orbits in apoastron and periastron position. When the stability is different we mark the apoastron by a triangle, the periastron with a square.
Open with DEXTER

	Figure 2: Schematic view of the stability of orbits in the resonances of 3rd order in HD 4208. Description like in Fig. 1.
Open with DEXTER

	Figure 3: Schematic view of the stability of orbits in the resonances of 4th order in HD 4208. Description like in Fig. 1.
Open with DEXTER

As an example we discuss the results of the investigation of the MMRs for the system HD 4208. We studied the following mean motion resonances up to the fourth order: 2:1, 3:2 and 4:3 (first order); 3:1, 5:3, 7:5 (second order); 4:1, 5:2, 7:4, 8:5, 10:7 and 11:8 (third order); 5:1, 7:3, 11:7, 13:9, 15:11 (fourth order). As shown in the figure for the first set of resonances (1st and 2nd order) the orbits close to the central star, which move well inside the HZ, are all stable (Fig. 1). For the MMRs close to the Jovian planet we can see a preference for stable orbits for the initial conditions $M=0^{\rm\circ}$, and $180^{\rm\circ}$. For the 3rd order resonance (Fig. 2) the picture is very inhomogeneous; for the Jovian planet in the apoastron position the orbits are stable even for the resonances close to the giant planet. The 4th order resonances are not destabilizing an orbit as we can see from Fig. 3; most of them are stable! The percentage of stable orbits in resonances is very large for HD 4208.

Details of the results from the investigation of resonances for all systems can be found in Table 3. For the resonances acting in three of the systems far from the perturbing planet almost all of them are stable; closer to the perturbing planet, they are more and more unstable (47 Uma, HD 72659 and HD 4208). Two systems are very much dominated by unstable motion in resonances: Gl 777 A and Gl 614.

The results of the orbital computations which we started exactly in the resonances (with the properly chosen semimajor axis which corresponds to the MMR we were investigating) are in a certain sense redundant. For some resonances (e.g. the 3:1 MMR) the starting positions of the giant planet in the pericenter or the apocenter and the starting position of the terrestrial planet on the connecting line in between - on a circular orbit - should give the same results. Also there is a symmetry for the initial mean anomalies M and $M+180^{\rm\circ}$ for the fictitious planet for some of the resonances (e.g. the 2:1 MMR). The reason that we have undertaken the computations for all positions is the following: when the results for the stability analysis were not the same - even when they should be equal - we took it as a sign of unstable motion. What we observed in our computations, which is in fact a known property, is that the unavoidable numerical errors can reveal the stability character of an orbit. Thus the results presented in Table 3 can be regarded as a good estimate for the stability of motions in a MMR. We emphasize that a detailed theoretical study for larger eccentricities is highly desirable.

5 Gl 777 A

The first discovery of a planet in Gl 777 A (=HD 190360) was reported by Naef et al. (2003) from the Geneva group of observers. This extrasolar planetary system is a wide binary with a very large separation (3000 AU); for our dynamical investigations of motions close to one star there was no need to take into account the perturbations of the very far companion. The central star is of spectral type G6 IV with $0.9 M_{\odot}$ and has a planet of minimum mass 1.33  ${{M_{\rm J}}}$ with a semimajor-axis of 4.8 AU. Because of the large eccentricity (e=0.48) the possible region of motion for additional planets is confined to a<2.4 AU (=periastron). Nevertheless, to have a global stability picture of possible additional planets, we investigated the stability in the region of the MMR from the 4:3 to the 5:1 resonance located at a=1.64 AU. From Table 3 one can see that only a few percent of the orbits started in the MMRs are stable.

	Figure 4: Main characteristics of the extrasolar system Gl 777 A. The light grey region shows approximately the position of the HZ; the dark grey bar indicates how closely the planet approaches the central star in its orbit.
Open with DEXTER

The interesting region of habitability (see Fig. 4), where planets could have temperature conditions to allow liquid water on the surface, corresponds roughly to 0.7< a < 1.3 AU, where we ignore the eccentricity of the terrestrial planet. We have started our computations in a larger region ( 0.5 < a < 1.3) with a grid spacing of $\Delta a = 0.01$ AU and changed also the eccentricity of the known planet between 0.4 < e < 0.5 with a gridsize of $\Delta e = 0.01$. The results of the two methods of analysis of the orbital behavior are shown in Figs. 5 and 6. In the first plot we show the results of the MEM, where two features are immediately visible: 1) strong vertical lines due to high order resonances; and 2) unstable orbits due to high eccentricity and high semimajor axes values (red or yellow colors). The latter feature is easy to understand because closer to the existing planet the perturbations are larger. The two methods complement each other in the information they convey; the MEM tells us about the variable distance to the central star and consequently it is a direct measure of the differential energy flux (insolation) on the planet. We can therefore determine where the variation of this distance does not exceed 50 percent, corresponding to an eccentricity of e=0.2. The Rényi entropy is a more sensitive probe of the dynamical character of the orbit, giving us a measure of the degree of chaos. In particular high order resonance features are made very clear using this second method, and we can even see the resonances acting when the eccentricity of the planet is as low as e=0.4 (the bottom of Fig. 6).

As a result for habitability of a terrestrial planet inside the orbit of the Jovian planet, we find that for the system Gl 777 A there is quite a good chance that planets will last long enough in the HZ to acquire the necessary conditions for life in the region with a about 1 AU.

6 47 Uma

The first planet in 47 Uma was discovered by Butler & Marcy (1996), and the second, the outer one, was reported in a paper by Fischer et al. (2002). Former studies investigated the zones inside the orbit of the more massive inner planet (Jones & Sleep 2002; Cuntz et al. 2003) which may host terrestrial planets. The Jovian planet has an almost circular orbit with very small errors in the eccentricity. Because of the errors for the eccentricity of the inner Jovian planet's orbit, we varied its eccentricity from 0< e < 0.12 for our computations. It turned out that the system became unstable when $e_{{\rm inner}} >0.12$ and consequently also no inner terrestrial planet would survive. Almost all the MMRs of the fictitious planet with the inner planet are in the HZ (Fig. 7).

The respective computations, summarized in Table 3, show that almost 1/3 of planets located in these resonances would survive with moderate eccentricities, which can be explained by the small eccentricity of the Jovian planet. The situation is similar to Jupiter and Saturn but with a scaling factor of approximately 0.4 in a. Also in our Solar System, asteroids are unstable in the outer main belt at and beyond the 2:1 resonance, although some islands of stability exist, like the Hilda asteroids in the 3:2 MMR. The vertical line in Fig. 8, located near 1 AU, indicates unstable orbits that mark the 3:1 resonance. The thin vertical line of large eccentricities close to a=0.82 AU corresponds to the 4:1 mean motion resonance of the terrestrial planet with the inner large planet. We can observe a broad curved "unstable'' line between 0.8 AU and 0.92 AU, which is caused by secular resonances (=SR) of the perihelia; this was confirmed by test computations, where we omitted the outer planet. It is remarkable that these SR act very strong when the initial eccentricity of the inner planet is close to zero, then it becomes weaker and disappears for 0.08 < e < 0.1, but it is again visible for e>0.1. It is interesting to note how the period of the periastron of the giant planet changes ( $4000~ {\rm yr} < \dot \omega < 20~000 ~{\rm yr}$) with increasing eccentricity. Between the SR and the 3:1 resonances a small region of orbits with small enough eccentricities allowing habitability (e<0.2) may survive. On the outer part of the HZ the very strong 2:1 resonance (broad line close to 1.3 AU) limits the region of stable orbits, but still leaves a large region of stable orbits between the 3:1 and the 2:1 resonance. This is in good agreement with the results of the previously mentioned studies (e.g. Cuntz et al. 2003), which did not cover such a large volume of phase space of possible motions in the HZ. The entropy plot does not show any features other than the ones we can see from the MEM.

	Figure 5: Initial condition diagram for fictitious planets in the system Gl 777 A: initial semimajor axes of the planet versus the eccentricity of the Jovian planet. The maximum eccentricity of an orbit during its dynamical evolution is marked with different colors.
Open with DEXTER

	Figure 6: Initial condition diagram for fictitious planets in the system Gl 777 A: initial semimajor axes of the planet versus the eccentricity of the Jovian planet. The value of the entropy (=entropy plot) of an orbit is marked in different colors.
Open with DEXTER

	Figure 7: Main characteristics of the extrasolar system 47 Uma. Specifications like in Fig. 4.
Open with DEXTER

7 HD 72659

The G0V star HD 72659 was found to have a companion from the Keck Precision Doppler survey (Butler et al. 2002). The Jovian planet (2.55 ${{M_{\rm J}}}$) has an orbit with a semimajor axis a=3.24 AU and an eccentricity of e=0.18. The MMRs are located from 2.47 AU (3:2) to 1.1081 AU (5:1); the 5:1, 4:1 and 3:1 are well inside the periastron position of 2.657 AU and lie in the HZ (around 1 AU, see Fig. 9).

The resonances turned out to be stable in more than $40 \%$ of the orbits; especially the high order resonances close to the HZ are stable in both initial conditions (periastron and apoastron position). As a consequence we expected that these planetary systems may host additional terrestrial planets in stable orbits.

	Figure 8: Results of the MEM for 47 Uma.
Open with DEXTER

	Figure 9: Main characteristics of the extrasolar planetary system HD 72659. Specifications like in Fig. 4.
Open with DEXTER

	Figure 10: Results of the MEM for HD 72659.
Open with DEXTER

	Figure 11: Entropy plot for HD 72659.
Open with DEXTER

	Figure 12: Main characteristics of the extrasolar planetary system Gl 614. Specifications like in Fig. 4.
Open with DEXTER

	Figure 13: Results of the MEM for Gl 614. This plot shows that for this system it is very unlikely for a terrestrial planet to survive in a stable orbit which also permits liquid water to exist on its surface.
Open with DEXTER

	Figure 14: Main characteristics of the extrasolar planetary system HD 4208. Specifications like in Fig. 4.
Open with DEXTER

	Figure 15: Results from the MEM for HD 4208.
Open with DEXTER

	Figure 16: Entropy plot for HD 4208.
Open with DEXTER

Because of the uncertainties of the observed Jovian planet's eccentricity we varied it from 0.08 to 0.30 with a stepsize of 0.22/80 = 0.00275 and chose the initial semimajor axis of the fictitious terrestrial planet to satisfy 0.4 AU < a < 1.2 AU. The results are shown in Fig. 10 (MEM) and Fig. 11 (entropy plot). We can identify quite well in these plots the resonances up to the 7:1 resonance (only in the entropy plot). Again one can see that the dynamics of a single orbit can be determined quite well with this method; it not only confirms what is depicted in Fig. 10, it also shows many more details especially for the motions in resonances. On the contrary the MEM is the appropriate tool for determining the eccentricity, which is - together with the semimajor axes - the crucial parameter for our research of determining planets in habitable zones.

Globally we can see a quite stable HZ in this extrasolar system which allows planets on orbits with small eccentricities. The strong unstable line close to 1.6 AU corresponds to the 3:1 resonance, while the other resonances, although giving rise to large perturbations in the eccentricities, are confined to the center of the resonance. In Fig. 10 we can also see that for the most probable eccentricity value of the Jovian planet (e=0.18) all orbits up to the 3:1 resonance are stable with low eccentricities (e < 0.2) and as a consequence the HZ could be populated by a terrestrial planet (or even more planets depending on their masses).

8 Gl 614

A planet around the K0V star Gl 614 (=14 Herculis) was discovered by Naef et al (2004) which is quite massive (4.89 ${{M_{\rm J}}}$) and orbits the central star with a semimajor axis a=2.85 AU and an eccentricity of e=0.38. The computation of orbits in MMRs showed that none of them are stable up to the 6:1 resonance in apoastron and periastron of the existing planet (for details see Table 3); this is not surprising when one looks at the large mass and the semimajor axis combined with the large eccentricity: the apoastron of the planet is 1.76 AU. The characteristics of this systems are shown in Fig. 12.

Looking at the results from the MEM (Fig. 13) one can see that for the most probable value of the eccentricity of the Jovian planet only the orbits with a<0.7stay in the stability and habitability (e<0.2) limits within the HZ.

9 HD 4208

Using the data of the Keck Precision Velocity Survey, Vogt et al. (2002) discovered that the system HD 4208 contains a planet with a slightly smaller mass than Jupiter on an almost perfectly circular orbit at 1.67 AU. The G5V star has a mass almost like our Sun. Our simulations of the dynamical evolution of terrestrial planets were carried out in the same way as for the other four systems but with special emphasis on the resonances in this system (see Sect. 3). Unfortunately, we made an error in the mass of the central star (we took it to be 14% smaller than the one determined by the observers) for the computations in the HZ (not for the stability in resonances). The habitable zone will be shifted outwards only a few percent because of this mass error; motions in the resonances (in the HZ) will be, for all practical purposes, unaffected.

The MMRs of first, second, third and forth order were investigated for this system. For a total of 18 resonances in two different positions of the existing planet (apoastron and periastron) for 8 different values of the mean anomaly of the fictitious planets we computed the stability (for a discussion see Sect. 3).

We expected unstable motion for this system from a =1.40 AU outwards, a limit which we computed using estimations derived by Wisdom (1980) and Duncan et al. (1989) for the onset of global chaos. We took the mean value between Wisdom's and Duncan's estimate (namely $\Delta a = 1.27 \mu ^{2/7}$ - where $\mu$ is the mass ratio of the primaries - which leads to a $\Delta a=0.28$for this system). In the two plots (Figs. 15 and 16) one can see that in fact the orbits are unstable from a=1.3 AU.

From the 2:1 resonance inward, for all initial conditions of the Jovian planet, we can expect stable orbits with low enough eccentricities (e<0.2) to allow habitable terrestrial planets with liquid water, the only exception being the resonant orbits (the vertical lines in both figures).

10 Conclusions

We have undertaken a dynamical study of five extrasolar planetary systems using extensive numerical experiments to answer the question of whether they could host terrestrial-like planets in habitable zones. For the single-planet systems we used the elliptic restricted three body problem, and for the system 47 Uma the mutual perturbations of the two Jovian planets were taken into account and therefore the restricted four-body problem served as a dynamical model. Because of the dependence of the stability of an orbit in a mean motion resonance on angular position we have selected 8 different positions for the mean anomaly of the planet with a step of $45^{\rm\circ}$ in the apoastron and the periastron position of the Jovian planet.

The characteristics of each system dictated the initial conditions, chosen in a fine 80 by 80 grid within the habitable zone, from which the orbits were computed using a robust numerical method (Lie-series integration) for 1 million years. The grid of the initial conditions of the fictitious terrestrial planets was chosen to cover the whole habitable zone of the system, and also to model the uncertainties in the elements of the observed planet(s). The stability of orbits was assessed with two methods, namely the computation of the Rényi entropy as measure of the chaoticity of an orbit and the determination of the maximum eccentricity of the orbit of a fictitious planet during its orbital evolution of 1 million years.

We can say that our computations for such a fine grid, taking into account also the essential role of the MMRs, lead to a deeper insight concerning the dynamics of the five systems which we studied. We also give the percentage of orbits which survived in the paper (=MT) of Menou & Tabachnik (2003) where they investigated all known extrasolar planetary systems with respect to possible additional terrestrial planets. We note that a direct comparison of MT with the percentages of "our'' survivors is not useful here because of the different approaches used; we have emphasized the role of the MMRs and neglected possible inclinations. However, we know that terrestrial planets will form within a protoplanetary disk thus staying with small orbital inclinations (Richardson et al. 2000; Lissauer 1993); additionally in a recent publication (Dvorak et al. 2003a) it was shown that the inclinations of the fictitious planets up to 15 degrees do not change the stability of orbits in the HZ. The results for possible stable orbits may be summarized as follows:

• In the system Gl 777 A the stability zone for the motion of terrestrial planets is well inside of the HZ and suggests any planets residing there will survive for a sufficiently long time (in MT, 86.8% of the orbits were found to be stable).

• In the system 47 Uma we can also say that there is a good chance for planets to move inside the HZ with small eccentricities between the main resonances; this is a result which is consistent with others but not with MT where only 28% remained.

• HD 72659 turned out to be a very good candidate for hosting planets in the HZ; again this does not confirm the results of MT - they found that only 40.2% of the orbits were stable.

• The results of the computations for Gl 614 show that it is very unlikely that there is an additional planet moving in the HZ - these results are consistent with MT (9.2% stable orbits)

• For HD 4208 in the HZ there is enough room left for terrestrial planets and that they could survive for a sufficiently long time; these results are more or less consistent with those of MT where 50.2% of the orbits were stable.

New observational possibilities provided by missions like COROT, DARWIN or the TPF make the first search for terrestrial exoplanets possible in the next decades. In an ESA study the goal of the missions is summarized as follows: "To detect and study Earth-type planets and characterize them as possible abodes of life''. In this sense dynamical studies like the one we present here should help to define promising targets for observations.

Acknowledgements

We all deeply thank the Helmholtz Institute for hosting the summerschool 2003 and for providing the necessary financial support. We also acknowledge the generosity of the Max Planck Albert Einstein Institut where we were able to perform these massive computations on their new high performance cluster, PEYOTE. Special thanks go to Profs. Kurths and Ruediger, the directors of the HISP in Potsdam for their hospitality.

   
Appendix A: Recurrence Plots

A.1 Introduction

In his seminal paper Henri Poincaré (1890) introduced the concept of recurrences in phase space, when he studied the stability of the solar system. Recurrence Plots (RPs) were introduced to visualize the recurrences of trajectories of dynamical systems in phase space (Eckmann et al. 1987). These plots have proved to be rather useful in the analysis of time series, as they give a first impression about the behavior of the system under study. However, in order to go beyond mere visual inspection, which depends on the person doing the analysis and hence is always to some extent subjective, measures that quantify the structures found in RPs were introduced (Webber & Zbilut 1994). This quantification of RPs has found numerous applications in many fields, such as Geology (Marwan et al. 2002), Physiology (Marwan et al. 2002; Zbilut et al. 2002), Climatology (Marwan & Kurths et al. 2002), etc. Especially, it has been shown that it is possible to estimate dynamical invariants based on RPs, like the Rényi entropy and the correlation dimension (Thiel et al. 2003). The estimation of these invariants by means of RPs has some advantages with respect to other algorithms usually applied to the analysis of time series. The method of RPs avoids, for example, the problem of the choice of embedding parameters, and it can also be applied to non-stationary data (Thiel et al. 2004). We concentrate on the estimation of the Rényi entropy, as it can be interpreted as the inverse of the prediction time (or Lyapunov timescale), i.e., it can be used to measure the prediction horizon of the system.

A.2 Definition of Recurrence Plots and examples

RPs were introduced to simply visualize the behavior of trajectories in phase space (Eckmann et al. 1987). Suppose we have a dynamical system represented by the trajectory $\{\vec x_i\}$ for $i=1,\ldots,N$ in a d-dimensional phase space. Then we compute the recurrence matrix

$$R_{i,~j} = \Theta(\varepsilon-\left\vert\vert\vec x_{i} - \vec x_{j}\right\vert\vert), \quad ~ i, j=1\dots N,$$ (A.1)

where $\varepsilon$ is a predefined threshold and $\Theta(\cdot)$ is the Heaviside function. The graphical representation of R_i, j, called a Recurrence Plot, is obtained by encoding the value one for a "black'' point and zero for a "white'' point. The recurrence rate RR is defined as the percentage of black points in the RP, i.e.

$$RR= \frac{1}{N^2}\sum\limits_{i,j=1}^{N}R_{i,j}.$$

Figure A.1: Prototypical examples of RPs: a) RP of a sine function; b) RP of the Rössler system in chaotic regime; c) RP of white noise.

Figure A.1a shows the RP of a sine function, i.e a circle in phase space. The plot is characterized by non-interrupted diagonal lines. Figure A.1b is the RP of the Rössler system in a chaotic regime. In this case the predominant structure are diagonal lines, which are interrupted. Figure A.1c represents the RP of white noise. It is homogeneous with mainly single points, indicating a stochastic system, as the state at time t+1 is independent of the one at t. From these plots, we see that there is a certain connection between the length of diagonal lines and the ratio of determinism or predictability inherent to the system. This connection can be explained as follows: suppose that the states at t=i and t=j are neighboring, i.e. R_i,j=1. If the system behaves predictably, then similar situations lead to a similar future, i.e. the probability for R_i+1,j+1=1 is high. For perfectly predictable systems, this leads to infinitely long diagonal lines (like in the RP of the sine function). In contrast, if the system is stochastic, the probability for R_i+1,j+1=1 is small, i.e., we find single points or short lines. If the system is chaotic, initially neighboring states diverge exponentially. The faster the divergence, i.e. the higher the Lyapunov exponent, the shorter the diagonals. In the next section we develop this observation and explain formally the relationship between the length of the diagonal lines in the RP and the predictability of the system.

Figure A.2: Left panel: number of diagonal lines of at least length l versus l in the RP of the Bernoulli map for different values of the threshold $\varepsilon$. The mean slope of the curves is equal to 0.6917. Right panel: Rényi entropy vs. $\varepsilon$ for the Bernoulli map.

  
A.3 Estimation of the Rényi entropy based on the RP

In this section we recall first the definition of the Rényi entropy of second order and then present the mathematical relation between this information measure and the RPs. Let us consider a trajectory $\vec x(t)$ in a bounded d-dimensional phase space and the state of the system is measured at time intervals $\tau$. Let $\{1,2,...,M(\varepsilon)\}$ be a partition of the attractor in boxes of size $\varepsilon$. Then p(i₁,...,i_l) denotes the joint probability that $\vec x(t=\tau)$ is in the box i₁, $\vec x(t=2\tau)$ is in the box i₂, ..., and $\vec x(t=l\tau)$ is in the box i_l. The order-2 Rényi entropy (Rényi 1970; Grassberger 1983) is then defined as

$$K_2=-\lim_{\tau\to 0}\lim_{\varepsilon \to 0}\lim_{l \to \infty}\frac{1}{l\tau}\ln \sum_{i_1,...,i_l} p^2(i_1,\dots,i_l).$$ (A.2)

Roughly speaking, this measure is directly related to the number of possible trajectories that the system can take for l time steps in the future. If the system is perfectly deterministic, in the classical sense, there is only one possibility for the trajectory to evolve and hence K₂=0. In contrast, one can easily show that for purely stochastic systems, the number of possible future trajectories increases to infinity so fast, that $K_2 \to \infty$. Chaotic systems are characterized by a finite value of K₂, as they belong to a category between pure deterministic and pure stochastic systems. Also in this case the number of possible trajectories diverges but not as fast as in the stochastic case. The inverse of K₂ has units of time and can be interpreted as the mean prediction horizon/time of the system.

In the paper of Thiel et al. (2003) it was shown that there exist a direct relationship between K₂ and RPs, and it is as follows

$$\ln P^c_{\varepsilon}(l) \sim \varepsilon^{D_2}\exp(-\hat{K}_2(\varepsilon)\tau),$$ (A.3)

where $P^c_{\varepsilon}(l)$ is the cumulative distribution of diagonal lines in the RP, i.e. the probability of finding a diagonal in the RP of at least length l . D₂ is the correlation dimension of the system under consideration (Grassberger & Procaccia 1983). Therefore, if we represent $P^c_{\varepsilon}(l)$ in a logarithmic scale versus l we should obtain a straight line with slope $-\hat{K}_2(\varepsilon)\tau$ for large l's, which is independent of $\varepsilon$ for a rather large range in $\varepsilon$. This is shown in the left panel of Fig. A.2 for the chaotic Bernoulli map x_n+1=2 x_n mod(1). One obtains $\hat K_2=0.6733$for 10 000 data points, in good accordance with the values found in the literature (Ott 1993).

Finally one represents the slope of the curves for large l in dependence on $\varepsilon$, where a plateau is found for chaotic systems. The value of this plateau determines $\hat{K}_2$ (Fig. A.2, right panel). If the system is not chaotic, one has to consider the value of the slope for a sufficiently small value of $\varepsilon$.

A.4 Automatization of the algorithm

To compute the stability diagrams presented in this paper, which consist of about 13000 grid points, the estimation of K₂ had to be automated. The crucial point in the automatization is the estimation of the scaling region of $\ln P^c_{\varepsilon}(l)$ vs. l and the plateau in $K_2(\varepsilon)$vs. $\varepsilon$. In both cases we applied a cluster dissection algorithm (Spaeth 1992). The algorithm divides the set of points into distinct clusters and performs a linear regression in each cluster. Then the sum of all square residuals is minimized in order to determine the scaling region and the plateau. To find both regions automatically, we used the following parameters:

• We considered only diagonal lines up to length $l_{\max}=400$. Longer lines were excluded because of finite size effects.
• We considered only values of $P^c_{\varepsilon}(l)$ with $P^c_{\varepsilon}(l)> 500$ for the same reason as in the last item.
• We used 40 different values for $\varepsilon$, corresponding to 40 equally spaced recurrence rates RR between 1% and 95%, to have a well defined plateau in $K_2(\varepsilon)$ vs. $\varepsilon$.
• We used 10 000 data points of each simulated trajectory. The more data points one uses, the more pronounced the scaling regions. Note that the computation time increases approximately with N².
• For the applied cluster dissection algorithm we had to specify the number of clusters in each run. For the detection of the scaling region in $\ln P^c_{\varepsilon}(l)$ vs. l, we chose 2 different clusters and used the slope of the largest cluster. For the detection of the plateau in $K_2(\varepsilon)$ vs. $\varepsilon$, we chose 3 clusters and used the value of the cluster with the minimum absolute slope.

These choices have proven to be the most appropriate ones for the estimation of the scaling regions. All these parameters are defaults of the computer program. The only needed input is the file with the trajectory data.

   
Appendix B: The Lie-series method

B.1 General properties of the Lie-series

Gröbner (1967) defined the Lie-operator D as follows:

$$D=\theta_1(z)\frac{\partial}{\partial z_1} + \theta_2(z)\frac... ...l z_2} + \dots + \theta_n(z)\frac{\partial}{\partial z_n}\cdot$$ (B.1)

D is a linear differential operator; the point $z=(z_1,z_2,\dots,z_n)$ lies in the n-dimensional z-space, the functions $\theta_i(z)$ are holomorphic within a certain domain G, e.g. they can be expanded in a converging power series. Let the function f(z) be holomorphic in the same region as $\theta_i(z)$. Then D can be applied to f(z):

$$Df=\theta_1(z)\frac{\partial f}{\partial z_1} + \theta_2(z)\f... ...z_2} + \dots + \theta_n(z)\frac{\partial f}{\partial z_n}\cdot$$ (B.2)

If we proceed applying D to f we get

         D²f = D(Df)

$$\vdots$$

Dⁿf = D(D^n-1f).

The Lie-series will be defined in the following way:

$$L(z,t)=\sum_{\nu=0}^{\infty} \frac{t^\nu}{\nu!} D^\nu f(z) = f(z)+ tDf(z) + \frac{t^2}{2!}D^2f(z)+\dots$$

Because we can write the Taylor-expansion of the exponential function

$${\rm e}^{tD} = 1 + tD^1 + \frac{t^2}{2!}D^2+\frac{t^3}{3!}D^3+\dots$$ (B.3)

L(z,t) can be written in the symbolic form

$$L(z,t)={\rm e}^{tD}f(z).$$ (B.4)

The convergence proof of L(z,t) is given in detail in Gröbner (1967).

One of the most useful properties of Lie-series is the Exchange Theorem:
Let F(z) be a holomorphic function in the neighborhood of $(z_1,z_2,\dots,z_n)$ where the corresponding power series expansion converges at the point $(Z_1,Z_2,\dots,Z_n)$; then we have:

$$F(Z)= \sum_{\nu=0}^{\infty} \frac{t^\nu}{\nu!} D^\nu F(Z)$$ (B.5)

or

$$F({\rm e}^{tD})z={\rm e}^{tD}F(z).$$ (B.6)

Making use of it we can demonstrate how Lie-series solve DEs. Let us give the system of DEs:

$$\frac{{\rm d}z_i}{{\rm d}t}=\theta_i(z)$$ (B.7)

with $(z_1,z_2,\dots,z_n)$. The solution of (B.7) can be written as

$$z_i={\rm e}^{tD}\xi_i$$ (B.8)

where the $\xi_i$ are the initial conditions z_i(t=0) and D is the Lie-operator as defined in (B.1). In order to prove (B.8) we differentiate it with respect to the time t:

$$\frac{{\rm d}z_i}{{\rm d}t}=D{\rm e}^{tD}\xi_i={\rm e}^{tD}D\xi_i.$$ (B.9)

Because of

$$D\xi_i=\theta_i(\xi_i)$$ (B.10)

we obtain the following result which turns out to be the original DE (B.7):

$$\frac{{\rm d}z_i}{{\rm d}t}={\rm e}^{tD}\theta_i(\xi_i)=\theta_i({\rm e}^{tD}\xi_i)=\theta_i(z_i).$$ (B.11)

B.2 A simple example

To demonstrate the principle of the Lie-integration, we will show how one proceeds in the simple case of the harmonic oscillator, a 2nd order DE where the solution is known:

$$\frac{{\rm d}^2x}{{\rm d}t^2}+\alpha^2 x = 0.$$ (B.12)

The first step consists in separating into two 1st order DEs such that

$$\frac{{\rm d}x}{{\rm d}t}=y=\theta_1(x,y)$$

$$\frac{{\rm d}y}{{\rm d}t}=-\alpha^2x=\theta_2(x,y)$$

with the initial conditions $z(t=0)=\xi$ and $y(t=0)=\eta$. With this notation we find the Lie-operator of the form

$$D=\theta_1\frac{\partial}{\partial \xi} + \theta_2\frac{\part... ...artial \xi} - \alpha^2\xi\frac{\partial}{\partial \eta}\cdot$$ (B.13)

The solution can now be written as a Lie-series

$$x={\rm e}^{\tau D}\xi \hskip0.7cm {\rm and} \hskip0.7cm y={\rm e}^{\tau D}\eta$$ (B.14)

where $t-t_0=\tau$ and for $\tau=0$ this are the initial conditions. Being aware of the symbolic development of ${\rm e}^{\tau D}$ we can compute the first terms:

$$D^1\xi=\eta=\theta_1$$

$$D^2\xi=D\eta=-\alpha^2\xi=\theta_2$$

$$D^3\xi=-\alpha^2 D\xi=-\alpha^2 \eta$$

$$D^4\xi=-\alpha^2 D\eta=\alpha^4 \xi$$

$$D^5\xi=\alpha^4 D\xi=\alpha^4 \eta$$

$$D^6\xi=\alpha^4 D\eta=-\alpha^6 \xi$$

$$\vdots$$

For the Lie-terms even respectively odd in the order consequently one can find

$$D^{2n}\xi=(-1)^n\alpha^{2n}\xi$$

$$D^{2n+1}\xi=(-1)^n\alpha^{2n}\eta ~$$

which leads to the solution for z:

$$z=\xi+\tau\eta-\frac{\tau^2}{2!}\alpha^2\xi-\frac{\tau^3}{3!}\alpha^2\eta + \frac{\tau^4}{4!}\alpha^4\xi \dots$$

Finally we get after the factorization of $\xi$ and of $\eta$

$$z=\xi\left(1-\frac{\tau^2}{2!}\alpha^2+\frac{\tau^4}{4!}\alpha^4-\frac{\tau^6}{6!}\alpha^6+\dots\right)$$

$$\hspace*{5.5mm}+\frac{\eta}{\alpha}\left(\tau\alpha-\frac{\ta... ...c{\tau^5}{5!}\alpha^5-\frac{\tau^7}{7!}\alpha^7+\dots\right) ~$$

which is exactly the known solution of the harmonic oscillator:

$$z(t)=\xi\cos\alpha\tau + \frac{\eta}{\alpha}\sin{\alpha\tau}.$$

References

• Beaugé, C., & Michtchenko, T. A. 2003, MNRAS, 341, 760 [NASA ADS] [CrossRef] (In the text)
• Butler, R., & Marcy, G. 1996, AJ, 464, L153 [NASA ADS] [CrossRef] (In the text)
• Butler, R., Marcy, G., Vogt, S., et al. 2002, AJ, 582, 455 [NASA ADS] (In the text)
• Callegari, N., Michtchenko, T. A., & Ferraz-Mello, S. 2002, A&AS, 34, 30 [NASA ADS].01 (In the text)
• Cuntz, M., von Bloh, W., Bounama, C., & Franck, S. 2003, Icarus, 162, 214 [NASA ADS] [CrossRef] (In the text)
• Duncan, M., Quinn, T., & Tremaine, S. 1989, Icarus, 82, 402 [NASA ADS] [CrossRef] (In the text)
• Dvorak, R., Pilat-Lohinger, E., Funk, B., & Freistetter, F. 2003a, A&A, 398, L1 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Dvorak, R., Pilat-Lohinger, E., Funk, B., & Freistetter, F. 2003b, A&A, 410, L13 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Eckmann, J. P., Kamphorst, S. O., & Ruelle, D. 1987, Europhys. Lett., 4, 973 [NASA ADS] (In the text)
• Engels, J., & Henrard, J. 1994, CMDA, 58, 215 [NASA ADS] (In the text)
• Érdi, B., & Pál, A. 2003, Dynamics of resonant exoplanetary systems, Proc. 3rd Austrian-Hungarion Workshop on Trojans and related topics, ed. F. Freistetter, R. Dvorak, & B. Érdi, 3 (In the text)
• Fischer, D., Marcy, G., Butler, R., Laughlin, G., & Vogt, S. 2002, AJ, 564, 1028 [NASA ADS] [CrossRef] (In the text)
• Froeschlé, C., Lega, E., & Gonczi, R. 1997, CMDA, 67, 41 [NASA ADS] (In the text)
• Grassberger, P. 1983, Phys. Lett., 97A, 227 [NASA ADS] (In the text)
• Grassberger, P., & Procaccia, I. 1983, Physica D, 9, 189 [NASA ADS]
• Groebner, W. 1967, Die Lie-Reihen und ihre Anwendungen (Berlin: VEB Deutscher Verlag der Wissenschaften) (In the text)
• Hadjidemetriou, J. D. 2002, CMDA, 83, 141 [NASA ADS] (In the text)
• Hanslmeier, A., & Dvorak, R. 1984, A&A, 132, 203 [NASA ADS] (In the text)
• Jones, B. W., & Sleep, P. N. 2002, A&A, 393, 1015 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, Icarus, 101, 108 [NASA ADS] [CrossRef] (In the text)
• Kiseleva-Eggleton, L., Bois, E., Rambaux, N., & Dvorak, R. 2002a, AJ, 578, L145 [NASA ADS] [CrossRef] (In the text)
• Kiseleva-Eggleton, L., Bois, E., Rambaux, N., Dvorak, R., & Rivera, E. J. 2002b, AAS, 201, #24.04
• Lammer, H., Dvorak, R., Pilat-Lohinger, E., et al. 2003, Atmosphere and orbital stability of exosolar planets orbiting gamma Cephei. In Proceedings EGS-AGU-EUG Joint Assembly, Nice, in press (In the text)
• Lichtenegger, H. 1984, CMDA, 34, 357 [NASA ADS] (In the text)
• Lissauer, J. J. 1993, ARA&A, 31, 129 [NASA ADS] (In the text)
• Malhotra, R. 1998, Orbital Resonances and Chaos in the Solar System. In Solar System Formation and Evolution, ed. D. Lazzaro, R. Vieira Martins, S. Ferraz-Mello, J. Fernandez, & C. Beauge, ASP Conf. Ser. 149, 37 (In the text)
• Malhotra, R., & Dermott, S. F. 1990, Icarus, 85, 444 [NASA ADS] [CrossRef] (In the text)
• Marwan, N., & Kurths, J. 2002, Phys. Lett. A, 302, 299 [NASA ADS] [CrossRef] (In the text)
• Marwan, N., Thiel, M., & Nowaczyk, N. R. 2002, Nonlinear Processes in Geophysics, 9, 325 [NASA ADS] (In the text)
• Marwan, N., Wessel, N., Meyerfeldt, U., Schirdewan, A., & Kurths, J. 2002, Phys. Rev. E, 66, 026702 [NASA ADS] [CrossRef]. (In the text)
• Mayor, M., & Queloz, D. 1995, Nature, 378, 355 [CrossRef] (In the text)
• Menou, K., & Tabachnik, S. 2003, AJ, 583, 473 [NASA ADS] [CrossRef] (In the text)
• Naef, D., Mayor, M., Korzennik, S. G., et al. 2003, A&A, 410, 1051 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Naef, D., Mayor, M., Beuzit, J. L., et al. 2004, A&A, 414, 351 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Ott, E. 1993, Chaos in Dynamical Systems (Cambridge University Press) (In the text)
• Pál, A., & Sándor, Zs. 2003, Dynamical stability of the habitable zones of exoplanetary systems. In Proceedings of the 3rd Austrian-Hungarion Workshop on Trojans and related topics, ed. F. Freistetter, R. Dvorak, & B. Érdi, 25 (In the text)
• Pilat-Lohinger, E., & Dvorak, R. 2002a, CMDA, 82, 143 [NASA ADS] (In the text)
• Pilat-Lohinger, E., Funk, B., Freistetter, F., & Dvorak, R. 2002b, Stability of planetary orbits in double stars. In Proceedings of the First European Workshop on Exo-Astrobiology, ed. H. Lacoste, 547 (In the text)
• Pilat-Lohinger, E., Funk, B., & Dvorak, R. 2003, A&A, 400, 1085 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Poincarè, H. 1890, Acta Math., 13, 1 [NASA ADS] (In the text)
• Rényi, A. 1970, Probability theory (North-Holland) (appendix) (In the text)
• Richardson, D. C., Quinn, T., Stadel, J., & Lake, G. 2000, Icarus, 143, 45 [NASA ADS] [CrossRef] (In the text)
• Rössler, O. E. 1976, Phys. Lett., 57A, 397 [NASA ADS]
• Sandor, Zs., Pilat-Lohinger, E., Dvorak, R., Érdi, B., & Süli, A. 2004, CMDA, in preparation (In the text)
• Späth, H. 1992, Cluster Dissection and Analysis (Chichester: Horwood) (In the text)
• Thiel, M., Romano, M. C., & Kurths, J. 2003, Applied Nonlinear Dynamics, 11, N3 [NASA ADS] (In the text)
• Thiel, M., Romano, M. C., Kurths, J., & Read, P. 2004, Chaos, in press (In the text)
• Vogt, S., Butler, R., Marcy, G., et al. 2002, AJ, 568, 352 [NASA ADS] [CrossRef] (In the text)
• Webber, C. L., & Zbilut, J. P. 1994, J. Appl. Physiol., 76, 965 [NASA ADS] (In the text)
• Wisdom, J. 1980, AJ, 85, 1122 [NASA ADS] [CrossRef] (In the text)
• Zbilut, J. P., Thomasson, N., & Webber, C. L. 2002, Med. Eng. Phys., 24, 53 [NASA ADS] [CrossRef] (In the text)