A&A 488, 1133-1147 (2008)
DOI: 10.1051/0004-6361:200809822

Dynamics and stability of telluric planets within the habitable zone of extrasolar planetary systems

Numerical simulations of test particles within the HD 4208 and HD 70642 systems

T. C. Hinse^1,2 - R. Michelsen¹ - U. G. Jørgensen¹ - K. Gozdziewski³ - S. Mikkola⁴

1 - Niels Bohr Institute, University of Copenhagen, Juliane Maries Vej 30, 2100 Ø, Denmark
2 - Armagh Observatory, College Hill, BT61 9DG Armagh, Northern Ireland, UK
3 - Nicolaus Copernicus University, Torun Centre for Astronomy, Gagarin Str. 11, 87-100 Torun, Poland
4 - Turku University Observatory, Väisäläntie 20, Piikkiö, Finland

Received 20 March 2008 / Accepted 3 June 2008

Abstract
Aims. We study gravitational perturbation effects of observed giant extrasolar planets on hypothetical Earth-like planets in the context of the three-body problem. This paper considers a large parameter survey of different orbital configuration of two extrasolar giant planets (HD 70642b and HD 4208b) and compares their dynamical effect on Earth-mass planetary orbits initially located within the respective habitable terrestrial region. We are interested in determining giant-planet orbit (and mass) parameters that favor the condition to render an Earth-mass planet to remain on a stable and bounded orbit within the continuous habitable zone.
Methods. We applied symplectic numerical integration techniques to studying the short and long term time evolution of hypothetical Earth-mass planets that are treated as particles. In addition, we adopt the MEGNO technique to obtain a complete dynamical picture of the terrestrial phase space environment. Both multi-particle and single-particle simulations were performed to follow an Earth-mass planet in the habitable region and its subsequent long term evolution.
Results. Our numerical simulations show that giant planets should be on circular orbits to minimize the perturbative effect on terrestrial orbits. The orbit eccentricity (and hence proximity) is the most important orbital parameter of dynamical significance. The most promising candidate for maintaining an Earth-mass planet on a stable and bounded orbit well-confined to the continuous habitable zone is HD 70642b. Even the large planetary mass of HD 70642b renders an Earth-mass planet habitable during the complete lifetime of the host star. The results allow us to extrapolate similar observed systems and points the necessity further constraining the uncertainty range in giant planet orbital eccentricity by future follow-up observations.

Key words: chaos - methods: N-body simulations - celestial mechanics - astrobiology

1 Introduction

This paper considers dynamical aspects and orbital stability properties of telluric (Earth-like) planets within the circumstellar habitable zone (HZ) of two selected single-planet extrasolar systems HD 4208 and HD 70642. The ability of a planet to support life is based on our knowledge of life on Earth and Earth's location within the solar system. The concept of the HZ, or the radial distance around a star in which a telluric planet's atmosphere might be able to maintain liquid water to eventual initiate and support life, has been presented and discussed in the scientific literature for many years (Huang 1959; Hart 1979; Kasting et al. 1993; Kasting & Catling 2003). For liquid water to exist, the average surface temperature must remain between 273 K and 373 K at 1 atm. atmospheric pressure.

The discovery of the first extrasolar giant (Jupiter-like) planet (Mayor & Queloz 1995) raises the natural question of whether Earth-like planets exist within known extrasolar planetary systems. Most of the known extrasolar planets (in orbit around a main-sequence host star) were detected by ground-based radial velocity observations. This technique favors the detection of a massive planet in a close orbit. Detecting Earth-mass planets around solar like host stars is just outside the current detection capability of the radial velocity technique.

In addition to the formation of observed Jupiter-like planets, terrestrial planets are expected to be formed by accretion processes within the inner region of a circumstellar protoplanetary disk (Agnor et al. 1999; O'Brian et al. 2006; Quintana et al. 2007; Raymond et al. 2004; Quintana & Lissauer 2006; Wetherill 1996; Chambers 2001; Chambers & Wetherill 1998; Kokubo & Ida 2000; Chambers 2003; Thébault et al. 2002). Today, major observational efforts focus on the detection of an Earth-like extrasolar planet by various techniques. Recently, the terrestrial planet formation theory is supported by both microlensing and radial velocity observations: Beaulieu et al. (2006) demonstrate the microlensing detection of a 5.5 Earth-mass planet at 2.6 AU orbital distance from the host star. In addition, Udry et al. (2007) report on the detection of an almost equally small mass planet close to or in the habitable region in orbit around an M-dwarf star. The currently on going satellite mission COROT will possibly detect Earth-mass planets by the transit method. Other future satellite missions sheduled for launch include TPF, KEPLER, and SIM - all capable of detecting extrasolar Earth-mass planets orbiting main sequence host stars.

This motivates us to conduct numerical experiments to study dynamics and stability properties of hypothetical telluric planetary orbits within observed extrasolar planetary systems. This subject is very interesting and has already been investigated by a few groups (Menou & Tabachnik 2003; Jones et al. 2001; Noble et al. 2002; Sándor et al. 2007; Asghari et al. 2004; Jones et al. 2005; Süli et al. 2007). In particular, we are interested in the orbital stability properties of a telluric planet considering giant planet perturbations. It is reasonable to expect that the requirement of long term orbit confinement of a telluric planet within the habitable zone is a necessary condition for the development of biological activity. The question is what orbit configurations of the giant planet are preferred to maintain a stable telluric planet that initially was assumed to be formed in the habitable zone. To answer this question, we chose to investigate and compare two systems containing a giant planet with different orbital properties. The difference in these two systems is mainly in the size of the derived semi-major axis. This work relies heavily on merging two approaches - the fast indicator mappings to resolve the overall phase space structure with longterm secular time scale integrations to understand individual orbits of Earth-like planets. In Sect. 2 we present orbital parameters inferred from radial velocity observations of the two systems HD 70642 (this system is a close analogue to the Sun) and HD 4208. In Sect. 3 we derive the zonal boundaries of the habitable region for each system. Sections 4 and 5 describes the applied numerical methods and adopted initial conditions for the parameter survey. Section 6 presents the results and the presented work is summarized with conclusions in Sect. 7.

2 Observations and orbital parameters of HD 4208b and HD 70642b

   
2.1 HD 70642

A planet, HD 70642b, in orbit around the star HD 70642 has been announced by Carter et al. (2003) based on radial velocity observations. Data aquisition has been ongoing during a 5-year time period within the Anglo-Australian Planet Search Program. As a result of the long observing baseline and constant precision measurements, HD 70642b marks the onset of the emergence of planets in long-period near circular orbits. Orbital parameters of the planet are listed in Table 1. The star HD 70642 is a chromosphericlly inactive nearby solar-like G5 dwarf. Physical properties were obtained by photometrical and spectroscopic observations. The mass of HD 70642 is $1.0 \pm 0.05 ~M_{\odot}$ and its age is estimated to be $4 ~\times ~10^9$ yrs with surface temperature $T_{\rm eff} = 5670 \pm 20$ K. The metallicity of HD 70642 is consistent with the majority of observed host stars harboring planets.

   
2.2 HD 4208

The existence of a single planet, HD4208b, in orbit around HD 4208 has been presented by Vogt et al. (2001) based on radial velocity observations obtained within the Keck planet search program. Orbital parameters, as determined from 1-planet Kepler fits, are listed in Table 1. Physical information on host star properties are based on Strømgren photometry, Hipparcos astrometric, and Keck spectral observations. HD 4208 is a Sun-like star of spectral type G5V with an estimated mass of $0.93 ~M_{\odot}$ and surface temperature $T_{\rm eff} = 5572$ K. This estimate agrees with independent published results obtained within the Geneva-Copenhagen star survey program (Nordström et al. 2004). However, no estimate of the stellar age has been published yet in the literature. Preliminary calculations using stellar evolution models indicate an age of $4.5 \pm 0.5$ Gyr (Southworth 2005, private communication). In the case of HD 4208, the process of isochrone fitting age determination is not reliable due to too large an uncertainty range in observed stellar parameters. In this paper, we assume an age of 4.5 Gyr for HD 4208.

Table 1: Derived orbital parameters for HD 4208b and HD 70642b from synthetic 1-planet Kepler fits to observed radial velocity data.

3 The habitable zone in (a, e)-parameter space

No life we know of could exist without liquid water. It has therefore been customary to define the habitable zone (HZ) as the orbital region where liquid water can exist in principle on the surface of a solid planet. In its simplest approach, the HZ can be calculated as the circular orbital radii where the stellar flux absorbed on an atmosphereless planet equals the Planck flux at temperatures between $0~^{\circ} \textnormal{C}$ and $100~^{\circ}\textnormal{C}$. However, we will take a slightly more advanced model based on Kasting et al. (1993) that includes the thermodynamical effects of an atmosphere.

3.1 Zonal boundaries

In this paper, we adopt numerical estimates for the HZ from one-dimensional radiative-equilibrium atmosphere model calculations presented in Kasting et al. (1993).

Figure 1: Zonal boundaries of the habitable zone for the planetary systems around HD 4208 ( left panel) and HD 70642 ( right panel). ZAMS HZ refers to zero-age main sequence boundaries (initial boundaries) and CHZ refers to the continuous habitable zone (inner triangular region). In geometrical (x,y)-space, the contour lines represent circles of radius $r_{\rm in}$ and $r_{\rm out}$ at different times. The difference in the size of the CHZ triangular region reflects differences in host star properties.

Open with DEXTER

Their work provides a good estimate of the position and extent of the HZ of an Earth-like planet, although crucial parameters obviously have to be ignored, which will have additional climatic effects on the location of actual habitability of a planet from the host star. However, for all of these advanced models, we obviously have to keep in mind that essential input for the models are unknown in real life, such as the actual composition of the exoplanetary atmosphere, the planetary albedo, and the question of whether water actually exists on the planets. We disregard these additional uncertainties in the following and adopt the boundary estimate of the HZ from Kasting et al. (1993). Following their lead the HZ is determined by applying several theoretical criteria to define the inner and outer boundaries. The general strategy to find these criteria is to state environmental conditions for a terrestrial-like planet in order to maintain surface water in its liquid form during a substantial part of thehost star's main-sequence period. In the model calculations, an Earth-like planet is assumed to have a 1 atm. $\textnormal{CO}_{2}/\textnormal{H}_{2}\textnormal{O}/\textnormal{N}_{2}$atmosphere.

Kasting et al. provide three criterions. We adopt the intermediate criterion on the zonal boundaries in this work. This estimate of the inner and outer boundaries of the HZ is referred to as the runaway greenhouse and maximum greenhouse limit. The inner boundary is defined by the onset of water evaporation. Water vapor enhances the greenhouse effect and therefore promotes surface warming. Eventually at a critical point all surface water starts to evaporate into the atmosphere and subsequently (in a runaway manner) is lost from the upper atmosphere by hydrogen escape. The outer boundary is defined as the maximum distance at which a surface temperature of 273 K can be maintained by a cloud-free $\textnormal{CO}_{2}$ atmosphere (i.e. the point where there would be enough $\textnormal{CO}_{2}$ and water in the planet's atmosphere to raise surface temperatures to 273 K).

The continous HZ (CHZ) is introduced to account for the gradual change in stellar luminosity. According to stellar evolution models, the stellar luminosity and surface temperature are a function of time, resulting in an outward shift of the zonal boundaries. Adopting the intermediate criterion Kasting et al. (1993) provides estimates of the time evolution of the HZ boundaries for different stellar masses. This is necessary, since the habitable zone for both HD 4208 (with assumed an age of $4.5 \times 10^{9}$ yrs) and HD 70642 is substantially evolved away from the initial ZAMS location. We define the overlap region between the initial and present HZ annulus as the continuous HZ (CHZ). In this paper, the zonal range of habitability for both stars is determined by linear interpolation between the $0.75 ~M_{\odot}$and $1.0~ M_{\odot}$ boundary evolution lines shown in Kasting et al. (1993, Fig. 14).

This approach to adopting the zonal habitable boundaries is similar to previous work (Jones et al. 2001; Jones & Sleep 2002) that considers exoplanetary host stars of different spectral types. The time evolution of the HZ boundaries has been further developed by Jones et al. (2005) in combining the Kasting et al. model with modern stellar evolution models. A slightly different choice in the criterion for the boundary limits is made in Noble et al. (2002). A more advanced geodynamical model was developed to determine the HZ within the Solar system (Franck et al. 2000). Application of this model to the 47 UMa system is found in Cuntz et al. (2003).

3.2 Orbit constraints within an annulus

For a planet to maintain habitability (and eventually initiate biological activity) during the host star's entire main-sequence life time, the telluric planet's orbit needs to remain confined to within the CHZ. The fulfillment of this requirement is not obvious, if external gravitational perturbations are applied.

The following outline is a simple attempt to impose geometrical constrains on (a,e) orbit parameters implying habitability. This method reflects a new approach to assessing a qualitative test of the orbit confinement in (a,e)-space of the telluric planet and differs from previous work within dynamical analysis work of habitable telluric planets.

The confinement of a telluric planet's orbit to a given annulus (i.e. the CHZ), imposes certain constraints on orbital parameters (a,e). If $r_{\rm in}$and $r_{\rm out}$ denotes the inner and outer boundary distance of a given annulus from the central host star, then the constraint condition on the orbital pericenter and apocenter distances is

$$r_{\rm in} \le a(1-e) \le r_{\rm out}~~~\textnormal{and}~~~ r_{\rm in} \le a(1+e) \le r_{\rm out},$$ (1)

where a is the orbit semi-major axis and e the eccentricity. The loci of points with constant radius will define a continuous line within (a,e)-parameter space. Figure 1 shows the ZAMS, CHZ, and present time HZ boundaries in (a,e)-space for HD 4208 (left panel) and HD 70642 (right panel) using present-day stellar parameters along with the climate model described in Kasting et al. (1993). As a consequence of stellar evolution, the ZAMS HZ boundaries moves outward to the currently observed HZ. For the initial and current states of the host stars, the respective lower triangles mark the region of habitability at these two instants of time. The triangular region in Fig. 1, corresponding to the CHZ, represents the set of (a,e) orbits that are permissible for confining the exoplanet within the CHZ.

4 Applied numerical methods and tests

4.1 Model and direct orbit integrations

The underlying model used to study orbital stability properties of a telluric planet adopts point-mass Newtonian mechanics considering gravitational perturbations within the full and restricted three-body problem on different time scales. In the stellarcentric frame, the planetary problem ((n+1)-body problem; with n=2 being the number of planets) to be solved is given by the equations of motion (Morbidelli 2002)

$$\frac{{\rm d}^{2}\boldsymbol{r}_{i}}{{\rm d}t^{2}} = -\frac{k... ...rac{\boldsymbol{r}_{j}}{\vert\boldsymbol{r}_{j}\vert^3}\Bigg),$$ (2)

where (r₁, m₁) and (r₂,m₂) are the positions (relative to the host star) and masses of the telluric and observed extrasolar planets, respectively. The numeral k² denotes Gauss gravitational constant and m₀ is the mass of the star.

To solve the equations of motion we used the symplectic (hybrid-MVS) and extrapolation (Radau) algorithms as implemented in the MERCURY package (Chambers & Migliorini 1997; Chambers 1999). We split our numerical analysis into three sections: i) computation of MEGNO maps (to be discussed in more detail in Sect. 4.2); ii) multi-particle test simulations over intermediate time scales, and iii) single-massive-planet simulations on long time scales.

   
4.2 Aspects and properties of the MEGNO factor

In dynamical theory, orbital instability is usually associated with chaotic dynamics. A powerful numerical method of quantitatively differentiating between quasi-periodic and chaotic dynamics is the MEGNO (Mean Exponenial Growth factor of Nearby Orbits) technique introduced by Cincotta & Simó (2000); Cincotta et al. (2003). Calculating the MEGNO factor $\langle Y\rangle(T)$(sometimes also referred to as the MEGNO indicator) for a given set of initial conditions provides a measure of chaoticity of the system evolution over the time span T. The MEGNO indicator is closely related to the Lyapunov characteristic number (LCN or maximum Lyapunov exponent), which quantitatively measures exponential divergence of nearby orbits (Morbidelli 2002) and belongs to the class known as fast indicators. This technique is very effective in exploring a multi-parameter, dynamical system. In order to explore the global dynamics of a perturbed telluric planet, we numerically computed the MEGNO indicator to construct phase space stability (instability) maps. This technique is computationally efficient and allows one to explore a large phase space in relatively short times. Detailed theoretical and computational aspects of MEGNO and its applications to planet dynamics work are given in Gozdziewski (2001). The MEGNO technique has been successfully applied in a variety of dynamical analysis studies (Breiter et al. 2005; Gozdziewski 2004,2003; Gozdziewski & Maciejewski 2001), including stability studies of a telluric planet within 47 UMa (Gozdziewski 2002). In the following we briefly review important properties of the MEGNO indicator.

The numerical calculation of the MEGNO indicator is based on the ODEX package using the Gragg-Bulirsch-Stoer (GBS) extrapolation method of Hairer et al. (1993) solving the variational equations (using the technique outlined in Mikkola & Innanen 1999) within the framework of the n-body problem. The variational equations describe the dynamical propagation of a small change in initial conditions. Finding a solution to the system of variational equations allows computation of the time-averaged MEGNO indicator. Because in the chaotic system the variations can grow without bounds, during the course of the integration of the system, it is necessary to renormalize the variational vector after a certain time $\tau$. The renormalization procedure was originally introduced by Benettin et al. (1976) ensuring the prevention of numerical overflow and saturation. In practice, renormalization is a simple numerical operation due to the linear nature of the variational equations. Tancredi et al. (2001) points out that numerical errors are introduced at renormalisation time steps in weak chaotic domains of phase space resulting in a false computation of the LCN. However, this is only true for the 2-particle shadow method (Benettin et al. 1976) in which two adjacent orbits are integrated separately in phase space. Applying periodic renormalization in the framework of the variational equations does not lead to spurious numercial artifacts in the computation of the LCN as derived from the MEGNO technique. However, this statement is only partially true. In Sect. 4.3, we demonstrate the numerical effect of applying a renormalization time step in our effort to find a solution of the variational equations when computing the MEGNO factor and the derived LCN.

Following Cincotta & Simó (2000), Cincotta et al. (2003) and Gozdziewski (2001), if the phase space trajectory is a regular or quasi-periodic solution (characterized by linear divergence of nearby orbits), then asymptotically

$$\lim_{t~\to~\infty} \langle Y\rangle (t) \to 2.$$ (3)

If trajectories are chaotic, then

$$\lim_{t~\to~\infty} \langle Y\rangle (t) \approx \frac{\gamma}{2}t$$ (4)

where $\gamma$ is the LCN of the probed phase space solution with given initial conditions. In general, the MEGNO factor measures the dynamical behavior of the system's phase space point. However, in this work the orbit of the telluric planet is expected to be significantly disturbed by external giant planet perturbations. Hence, for given initial conditions, orbits of a quasi-periodic or chaotic nature will be attributed to the less massive telluric planet. The giant planet's orbit remains unaffected by small-planet perturbations on short time scales. As an example within small-body dynamics Wisdom (1983) demonstrated the presence of chaotic (erratic) behavior of a particle within the 3:1 mean-motion resonance with Jupiter.

   
4.3 The choice of MEGNO renormalization time step

In the literature, rigorous tests and accuracy checks have been performed in the calculation of the MEGNO indicator to ensure numerical stability of results, and possible caveats of this technique have been discussed (Gozdziewski 2001). In the following, we report on an additional test supplementing Gozdziewski's accuracy checks for the reliable computation of $\langle Y\rangle(t)$. This test concerns the choice of the renormalization time step $\tau$: a too long $\tau$ allows the separation vector to grow long enough to introduce numerical errors. Hence, $\tau$ should be kept small enough to maintain numerical accuracy. The Lyapunov characterstic number can be approximately recovered from the time-running evolution law $\gamma = 2\langle Y\rangle(T)/T$ (Cincotta & Simó 2000; see Eq. (4)). The goal of this test is to numerically reproduce the Lyapunov time (i.e. $t_{\rm L} = 1/\gamma$) for a particle perturbed by Jupiter for a range of renormalization time steps. Wisdom (1983) provides the initial conditions for numerically modeling asteroid dynamics within the 3:1 mean-motion resonance. In this particular case, Wisdom consistently derived $\log(\gamma) \approx -3.5$ or $t_{\rm L} \approx 3160$ yrs. Figure 2 shows the time evolution of the LCN, as computed from the MEGNO, in the classical $(\log\lambda(t),\log t)$-graph for various choices of the renormalization time step $\tau/T_{\rm Jup} = 1/4, \ldots, 10~000$, where $T_{\rm Jup}$ is Jupiters orbital period. Each numerical integration invokes identical initial conditions for the test particle. The graphs are derived from the $2\langle Y\rangle/t$evolution law. The range for the variable renormalization time step is from $\tau = 0.25~T_{\rm Jup}$ to $\tau = 10^{4}~T_{\rm Jup}$ as indicated in the figure legend. From Fig. 2 one observes a random scatter of $\log \gamma$for $\tau < 4~T_{\rm Jup}$, showing no systematics in the final $\gamma$ as a function of $\tau$ (see the inset figure in Fig. 2). A systematic trend is first observed for $\tau > 4 T_{\rm Jup}$, clearly showing a significant deviation of $\gamma$ from the value published in Wisdom (1983). Furthermore, Fig. 2 suggests that rescaling the variational equations should be done on the order of $T_{\rm Jup}$ or less. For $\tau = 10^{3}{-}10^{4}~T_{\rm Jup}$, the algorithm is no longer stable, yielding false Lyapunov times.This test concludes that the Lyapunov time of an asteroid can be determined from the time evolution of the MEGNO factor and the renormalization time should be a few tens of Jupiter's orbital period. The experiment shows the importance of choosing a proper renormalization time to obtain reliable results. Too large $\tau$ would allow the dynamical system to have enough time between any two renormalizations of the variational vector to produce numerical overflow due to exponential behavior of the orbit (or variational vector). However, it must be mentioned that no ready-made description exists of how to choose the correct renormalization step size and generally extrapolating the results obtained in this test to other problems is not guaranteed, so numerical tests are always encouraged. In general, the proper choice of the renormalization time depends on the nature of the dynamical problem under study. Intuitively, systems which involve strong perturbations generally require a small renormalisztion step size, since large deviations in the variational vector are expected on short time scales. A possibly universal approach to determining a proper renormalization step size would be to constantly monitor the norm of the variational vector in time and impose some constraint of a maximum allowed length in phase space. Whenever this maximum is reached, the algorithm applies the renormalization of the variational vector at a given time, after which the monitoring process is repeated. Within the present work, we were guided by our test and adopted a renormalization step size equivalent to the orbital period of the outer giant planet. Our choice of the renormalisation time proves adequate when comparing the obtained results with previously published independent data.

	Figure 2: Parameter survey of the renormalization time step $\tau$ to reproduce the Lyapunov time for an asteroid within the 3:1 mean-motion resonance as computed by Wisdom (1983). $T_{\rm Jup}$ is the orbital period of Jupiter. The inset figure shows resolved details on the time evolution for low choices of $\tau$.
Open with DEXTER

4.4 MEGNO stability maps

The use of MEGNO maps is an efficient way to explore the system's phase space structure, thereby providing a global picture of interesting dynamical features associated to resonance phenomena. In this work we apply the MEGNO technique to studying the orbital stability (instability) properties of a telluric planet within initial $(a_{\rm tel},e_{\rm giant})$ parameter space (or any other 2D phase space section with a different choice in initial conditions, Gozdziewski (2003), where $a_{\rm tel}$ denotes the semi-major axis of the hypothetical telluric (test particle) planet and $e_{\rm giant}$ the eccentricity of the perturbing giant planet. Results are shown in Figs. 4a, b and 5. In the following we outline the general construction of MEGNO maps. For a given range in eccentricity and semi-major axis, the grid of initial conditions is given by

$$a_{i} = a_{\rm min} + \frac{\Delta a}{N_{x}}~i = a_{\rm min} + \frac{a_{\rm max} - a_{\rm min}}{N_{x}}~i$$ (5)

and

$$e_{j} = e_{\rm min} + \frac{\Delta e}{N_{y}}~j = e_{\rm min} + \frac{e_{\rm max} - e_{\rm min}}{N_{y}}~j,$$ (6)

where $i = 0, \ldots, N_{x}$ and $j = 0, \ldots, N_{y}$ are integers measuring the scan resolution. Choosing high N_x,N_y will result in a more detailed mapping of phase-space structures within the considered range of (a_i,e_j) space. The maximum and minimum values define the range of the parameter survey. In this work the survey is mainly focused on the uncertainty range in observed orbital elements of the giant planet and the habitable region of the telluric planet. At each grid point (a_i,e_j), the MEGNO indicator $\langle Y \rangle$ is numerically calculated at integration end by simultanously solving the differential form of $\langle Y \rangle$ along with the variational equations using the FORTRAN cs.f code (Gozdziewski 2001). A typical scan grid of (N_x,N_y)=(200,160) is used. This corresponds to a typical resolution of $(\Delta a, \Delta e)=(0.006~\textnormal{AU},0.002)$. Within the maps, $\langle Y \rangle$ is calculated row-by-row. For a given giant planet eccentricity e_j, the habitable region $a_{0}, \ldots, a_{N_{x}}$ is scanned. The total integration time is chosen within the range $10^3{-}10^{4}~ P_{\rm orb}$, where $P_{\rm orb}$ is the orbital period of the giant planet and reflects the characteristic period of the system. This choice is suggested by Gozdziewski (2001).

4.5 Orbit accuracy and computational aspects

Figure 3: Graphical representation of Eq. (7) showing the dependence of the planet mass on the line-of-sight inclination for different orbit eccentricities for the giant planets HD 4208b ( left panel) and HD 70642b ( right panel). Numerical values for K and P within the mass function are taken from Table 1 for each system. The inset in each plot represents the $(i,m_{\rm pl})$ dependencies for e = 0.0, 0.1, 0.2, 0.3, 0.4 with $i \in [50^{\circ }, 90^{\circ }]$. The horizontal line corresponds to $M_{\rm Jup} = 16~M_{\rm Jup}$, which is the approximate transition mass between planets and brown dwarfs.

Open with DEXTER

Accuracy checks were performed in all numerical simulations presented in this work. To capture the correct physics (near the pericenter passage), we monitored the relative energy error dE/E (most sensitive) to probe numerical accuracy. Convergence tests were carried out to determine the minimum accuracy of the conservation of relative energy. Two accuracy parameters are necessary for a GBS-integration, which needs to be specified for a given accuracy requirement. To conserve energy on the order of d $E/E \simeq 10^{-12}$ during a given integration, the absolute $(\delta)$ and relative $(\epsilon)$ error tolerance parameters are both set to $\delta,\epsilon = 10^{-14}$ using double precision arithmetic. In general, the error in relative energy is characterized by a random walk. This is explained by stochastic accumulation of rounding errors during numerical computations. In Radau integrations, the local error accuracy parameter for step size adjustment is set to 10^-14resulting in $({\rm d}E/E)_{\rm max} \simeq 10^{-10}$.

5 Initial conditions and parameter survey

Writing the equations of motion in the stellar-centric frame, the three-body problem has 12 dynamical variables (6 degrees of freedom). In addition, including the masses as free parameters, we then have 15 free parameters to be determined as initial conditions for the numerical integration. Assuming coplanar orbits (eliminating four dynamical variables or two d.o.f.), we are left with $(a,e,\varpi,M)$ for each planet, where a the semi-major axis, e the eccentricity, $\varpi$ is the longitude of pericenter, and Mdenotes the mean anomaly of the orbit. The mass of the host star is held constant during the integration and its nominal value is adopted from the literature (cf. Sects. 2.1 and 2.2). Our parameter survey was done within the observed uncertainty range of the giant planet orbit eccentricity e and giant planet mass $m_{\rm pl}$. These quantities are dynamically interesting. The giant planet mass, orbital element uncertainty range, and nominal values are listed in Table 1.

5.1 Initial conditions for MEGNO and multi-particle simulations

In constructing MEGNO stability maps, the telluric planet was started on a circular orbit with mean longitude $180^{\circ}$. The planar three-body problem was considered, and we set the mass of the telluric planet to be $3 \times 10^{-3}~M_{\rm Jup}$. The survey range in the telluric planets semi-major axis was set to $a \in [0.76~\textnormal{AU}, 1.43~\textnormal{AU}]$ (HD 70642) and $a \in [0.57~\textnormal{AU},1.34~\textnormal{AU}]$ (HD 4208). For these experiments the observed giant planet's eccentricity was surveyed within the observed uncertainty range as shown in Table 1. The remaining orbital elements were set to their nominal published values and the mean anomaly of the giant planet set to $M = 0^{\circ}$ in all numerical experiments.

In multi-particle simulations, we adopt the co-planar restricted three-body model and numerically follow 2000 non-interacting test-particles (i.e without self-gravitation). Initially and for both systems, all particles are started on circular orbits with $M \in [0^{\circ}, 360^{\circ}]$ randomly distributed within the whole range of the HZ.

5.2 Planet mass survey

The derived giant planet mass is a minimum mass (corresponding to an orbit viewed edge-on). An upper limit in this quantity is unknown as long as the angle between the line-of-sight and the normal to the orbit is undetermined. This indeterminacy justifies a i-parameter survey, which directly translates to a survey in the mass parameter. The functional dependence of a planet's real mass on the line-of-sight inclination i and orbit eccentricity e, is obtained by considering the mass-function (Hilditch 2001). The significance of the e-dependency on the planetary mass is unknown and needs some attention in the following.

Figure 4: MEGNO maps of a telluric planet within the dynamical habitable zone for HD 4208 and HD 70642. The MEGNO indicator is color-coded with white ( $\vert\langle Y\rangle - 2.0\vert < 0.001$) representing initial conditions indicating quasi-periodic motion of the telluric planet and with any color different from white representing chaotic motion gradually deviating from quasiperiodic motion (here from $\langle Y \rangle =1.5$ to $\langle Y \rangle = 5$). Contour lines denoting the boundaries of the habitable zone(s) are superimposed, and a total of $32 \times 10^{3}$ initial conditions ( (Nx, Ny) = (200,160)) are considered each running for $5 \times 10^4 \times P$ years. The location of major mean-motion resonances are indicated by an arrow. Left panel: MEGNO maps for HD 4208b for giant planet masses $m_{\rm pl} = 1,4,6,16~M_{\rm Jup}$. Right panel: MEGNO maps for HD 70642 for giant planet masses $m_{\rm pl} = 2,4,12,16~M_{\rm Jup}$.

Open with DEXTER

If $m_{\rm pl}$ measures the giant planet mass, we have the following analytic expression for the minimum mass of the planet companion as a function of orbit eccentricity

$$m_{\rm pl}\sin i = \Bigg ( \frac{m_{0}^{2}}{2\pi~k^{2}}~(1-e^{2})^{3/2}~K^{3}~ P\Bigg)^{1/3},$$ (7)

where m₀ is the host star mass, K the velocity semi-amplitude, P the orbital period, and k is Gauss gravitational constant (in Eq. (7) we assumed $m_{\rm pl} \ll m_{0}$). For different orbit eccentricities the effect on the $(i,m_{\rm pl})$ dependency can be studied. Figure 3 shows this functional relation for orbit eccentricities within the range $e \in [0, 0.9]$. The maximum eccentricity for both HD 4208b and HD 70642b is e = 0.16 as derived from observations (cf. Table 1. From the inset figures we see that the change with the interval $e \in [0.0, 0.2]$ of $m_{\rm pl}$ as a function of i, is nearly indistinguishable. For $i=90^{\circ}$ the difference in mass between circular and high-eccentricity orbits corresponds to a factor two. This consideration shows that the effect of orbit eccentricity on $(m_{\rm pl},i)$dependency is only significant for high-eccentricity orbits.

Based on a statistical analysis on observed minimum masses Zucker & Mazeh (2001) estimates that 5% of solar-type G stars in the solar neighborhood have planets in the range 1-10  $M_{\rm Jup}$. A similar result was obtained by Tabachnik & Tremaine (2002). We decided to consider a giant planet mass survey in the range $m_{\rm pl}\sin i < m_{\rm pl} < 16~M_{\rm Jup}$ for HD 4208b and HD 70642b, respectively.

6 Results and discussion of numerical simulations

6.1 MEGNO maps

In the following, we present results in the form of MEGNO maps considering $(a_{\rm tel},e_{\rm giant})$ sections within the phase space of HD 70642 and HD 4208. In all maps, we superimpose the contour lines of the extent of the habitable zone boundaries within the terrestrial region of each system. The inner triangular region represents the continuous habitable zone. We chose to color-code quasi-periodic motion by white. Initial conditions resulting in quasi-periodic motion have $\vert\langle Y\rangle - 2.0\vert < 0.001$ at the end of the integration. Initial conditions resulting in chaotic orbits of the telluric planet differ from the white color.

6.2 HD 70642

MEGNO maps were computed for masses $m_{\rm pl} = 2,4,6, 8,10,12,14,16~M_{\rm Jup}$of the giant planet. This range in mass corresponds to a line-of-sight inclination within $i \in [7^{\circ},90^{\circ}]$. Figure 4b shows MEGNO maps for four selected masses $(m_{\rm pl} = 2,4,12,16~M_{\rm Jup})$, highlighting the emergence of main dynamical features as a function of increasing giant planet mass. The integration time is the same for each MEGNO scan. For $m_{\rm pl} = 2~M_{\rm Jup}$, we identify the presence of the 4:1, 9:2 and 7:2 mean-motion resonance for a large giant planet orbit eccentricity. If $n_{\rm tel}$ and $n_{\rm giant}$ measure the mean motion of the telluric and giant planet, respectivly, then the approximate nominal location $a_{\rm tel}$, of a mean-motion resonance is obtained from Kepler's 3rd law of motion and given by

$$a_{\rm tel} = \Bigg(\frac{n_{\rm giant}}{n_{\rm tel}}\Bigg)^{... ...rac{M+m_{\rm tel}}{M+m_{\rm giant}}\Bigg )^{1/3}~a_{\rm giant}$$ (8)

where $a_{\rm tel},m_{\rm tel}$ and $a_{\rm giant},m_{\rm giant}$ are the telluric and outer giant planet semi-major axis and mass, and M is the mass of the central star. Increasing the mass of the giant planet mainly results in a strengthening of the effect of mean-motion resonances over a wide range of giant-planet orbit eccentricity. This is particularly the case for the 4:1 mean-motion resonance as observed in the $m_{\rm pl} = 16~M_{\rm Jup}$ case in Fig. 4b. In this case the 4:1 mean-motion resonance starts to act for almost (e < 0.05)circular orbits of the giant perturber. For giant-planet mass $m_{\rm pl} = 16~M_{\rm Jup}$ in Fig. 4b, the 5:1 resonance appears in the center of the present-day habitable zone.

6.3 HD 4208

For HD 4208, Fig. 4a shows the results of MEGNO scans within the HD 4208b parameter space. Four selective scans corresponding to giant planet masses $m_{\rm pl} = 1,4,6,16~M_{\rm Jup}$ are shown. From Fig. 4a the dynamical effect of increasing the outer giant planet's mass is evidently observed. For $m_{\rm pl} = 1~M_{\rm Jup}$ in Fig. 4a, a general instability region occurs on the right side of semi-major axis $a_{\rm tel} \approx 1.30~\textnormal{AU}$ and only dynamically active for high giant-planet eccentricities. This instability region starts to move inward towards the center of the continuous habitable zone for higher masses of the observed giant planet. An interesting feature to note is the decrease in the stability ``island'' at $a_{\rm tel} = 1.26~\textnormal{AU}$extending over a giant planet eccentricity range $e \in [0.0, 0.15]$. Mean-motion resonances are observed and indicated with an arrow in each MEGNO map within the figure panels of Fig. 4a. Of particular interest is the 2:1 mean-motion resonance located at $a_{\rm tel} = 1.05~\textnormal{AU}$ for $m_{\rm pl} = 1~M_{\rm Jup}$. The 2:1 resonance is ``active'' for almost all eccentricities of the outer giant planet. With increasing giant-planet mass, the continuous habitable zone becomes dominated by the 3:1, 5:2, and 2:1 mean-motion resonances. This is clearly evident for $m_{\rm pl} = 16~M_{\rm Jup}$ in Fig. 4a. The dynamical structure of the terrestrial region in HD 4208 has also been studied by Asghari et al. (2004). They present a stability map analysis of hypothetical telluric planets using the Kolmogorov entropy as a quantitative measure to detect and differentiate chaotic/quasi-periodic dynamics. They complement their analysis with (MEM) maximum eccentricity maps. The results in Fig. 4a show an almost one-to-one correlation of the MEGNO technique with the Kolmogorov entropy method. Almost every feature from global (general) chaos to the presence of mean-motion resonances are reproduced by two independent techniques. However, Asghari et al. (2004) did not undertake a mass-parameter survey for the giant planet HD 4208b. The results presented here in the form of MEGNO stability maps considered various masses of the giant planet in HD 4208.

Figure 5: High-resolution MEGNO maps showing the phase-space finestructure of several (high-order) mean-motion resonances within the habitable zone of HD 4208. The mass of HD 4208b is $m_{\rm pl} = 1~M_{\rm Jup}$. The resolution of the scan is (Nx,Ny) = (200, 100), each with integration time $10^{4} \times P_{\rm orb}$. Top panel: finestructure of the 2:1 mean-motion resonances located at 1.06 AU. The solid contour line represents the outer boundary $r_{\rm out} \simeq 1.26~\textnormal{AU}$ of the initial (ZAMS) habitable zone. Lower panel: finestructure of mean-motion resonances within the region $a_{\rm tel} \in [1.1;1.25]$ AU. The two contour lines corresponds to the outer boundary of the ZAMS ( $r_{\rm out} = 1.26$ AU) and present-day ( $r_{\rm out} = 1.34$ AU) habitable zone. Mean motion resonances are indicated by arrows.

Open with DEXTER

6.4 Resonance finestructure

An interesting follow up is to study the phase-space finestructure of mean-motion resonances. Figure 5 (top panel) shows the finestructure of the 2:1 mean-motion resonance within HD 4208. This resonance resembles a wedge-like structure with a well-defined transition region from chaotic to quasi-periodic motions. The pecular V-shape structure is explained by resonance perturbation theory (Murray & Dermott 1999). The time evolution of the telluric semi-major axis captured within the 2:1 mean-motion resonance will exhibit oscillations determined by the resonance libration width.

The lower panel in Fig. 5, shows the finestructure in the region $a_{\rm tel} \in [1.1;1.25]$ AU. The high-resolution MEGNO scan shows a magnification of three prominent mean-motion resonances. In particular they correspond to the 9:5 $(a \simeq 1.145~\textnormal{AU})$, 5:3 $(a \simeq 1.183~\textnormal{AU})$, and 8:5 $(a \simeq 1.22~\textnormal{AU})$ commensurabilities. Each resonance is located at the outer edge of the continuous, habitable zone and dynamically active for near-circular orbits of the giant perturber. The characteristic form of each resonance is V-shaped and reflects the dependence of the libration width as a function of giant-planet eccentricity. In both high-resolution scans several high-order, mean-motion resonances are identified between the discussed resonances. It is interesting to note the quasi-periodic region in the 5:3 mean-motion resonance (lower panel of Fig. 5). This region corresponds to regular motion in the libration zone of the 5:3 mean-motion resonance and is characteristic of all mean-motion resonances. However, the exact dynamical nature (quasi-periodic or chaotic) of a given mean-motion resonance depends on the choice of initial conditions. A more detailed study of resonance finestructure will be conducted in the future.

6.5 Multi-particle simulations

MEGNO stability maps only represent a portion of the total system's phase space. No direct information on the time evolution of the telluric planet's eccentricity and/or semi-major axis is provided by MEGNO maps. In the following, we study the effect of giant planet (repeated) perturbations on telluric planet orbital parameters. This is done by numerically following $2 \times 10^{3}$ particles and studying their time evolution within (a,e)space. Since orbital energy and angular momentum are related to semi-major axis and orbit eccentricity, this approach should give some detailed understanding of the effect of resonance perturbations. Each particle is assumed to model a telluric planet initially formed within the HZ. In Figs. 8a, b and 7a, b, we show selected results in the form of (a,e)-snapshots of multi-particle simulations for HD 4208 and HD 70642. The particles were randomly distributed within $a_{\rm tel} \in [0.76~\textnormal{AU}; 1.35~\textnormal{AU}]$ (HD 70642) and $a_{\rm tel} \in [0.60~\textnormal{AU}, 1.34~\textnormal{AU}]$ (HD 4208). The HZ boundaries are superimposed in each snapshot and the CHZ is the inner triangular region. In particular, we focus our discussion on particle dynamics with the CHZ. Each figure shows simulation snapshots for two choices in initial parameters. The total integration time is 10⁶ years. Snapshot times were chosen so as to demonstrate the most prominent dynamical transitions and features. Learning from results obtained by MEGNO stability analysis, a total of four giant planet masses were considered $(m_{\rm pl} = [2,4,12,16]~M_{\rm Jup})$. Each mass parameter was then paired with either two values in the orbit eccentricity $(e_{\rm giant} = [0.0,0.16])$ corresponding to the minimum and maximum values of the observed eccentricity uncertainty range.

Figure 6: Simulation snapshots for particles (telluric planets) within HD 70642 considering different giant planet parameters. For each panel the time is indicated within each snapshot. A total of $2 \times 10^{3}$ particles are randomly distributed within the HZ. The continuous HZ is marked by the inner (0.92 AU) and outer (1.17 AU) boundaries. The innermost and outermost boundaries are located at 0.80 AU and 1.35 AU respectively. The total integration time is 106 years. Left 4-figure panel: $m_{\rm pl} = 2~M_{\rm Jup}$ and e= 0.0. Right 4-figure panel: $m_{\rm pl} = 4~M_{\rm Jup}$ and e=0.16.

Open with DEXTER

	Figure 7: Same as Fig. 6, but for different choices of $m_{\rm pl}$. Left 4-figure panel: $m_{\rm pl} = 12~M_{\rm Jup}$ and e = 0.0. Right 4-figure panel: $m_{\rm pl} = 16~M_{\rm Jup}$ and e=0.16.
Open with DEXTER

6.6 HD 70642

For HD 70642 we present simulation snapshots of a swarm of particles under the perturbation of different giant planet parameters. In general, within Figs. 6 and 7, we observe the development of oscillations in particle eccentricities induced by giant-planet perturbations. The main dynamics are observed within a temporal change in orbit eccentricity at a constant semi-major axis. This is best seen by examining single-particle dynamics, to be discussed later. This dynamical behavior is best demonstrated for the giant planet parameters shown in Figs. 6b and 7b. Initially the particle swarm is excited in eccentricity giving rise to a eccentricity gradient throughout the HZ after a fast relaxation time (10⁴ to 10⁵ yrs). This e-excitement is greater for particles in the proximity of the giant planet. The e-oscillations of neighboring telluric orbits are phased at the beginning of each simulation except at mean-motion resonance locations. These resonance perturbations quickly introduce dephasing of orbit eccentricity oscillations. Comparing Fig. 6b with MEGNO stability maps for HD 70642 (the $m_{\rm pl} = 4~M_{\rm Jup}$ case), we identified the prominent 3:1 (at $\simeq 1.3~\textnormal{AU}$) mean-motion resonances. Furthermore, the simulations in Figs. 6a, b and Figs. 7a, b shows that the oscillation amplitude at a given semi-major axis is constant (except at mean-motion resonances). However, the oscillation frequency increases with time, and high-frequency oscillations are observed for particles close to the giant perturber at the beginning. Animated MPEG simulation movies can be obtained upon request. Progressing in time, these frequencies increase and the e-oscillations begin to narrow. In addition, we observe that telluric planets in general are confined to within the CHZ for a low- $e_{\rm giant}$ giant planet orbit. This is also true if HD 70642b is a massive $m_{\rm pl} = 12{-}16~M_{\rm Jup}$ giant planet on a circular orbit. Increasing the giant planet eccentricity results in a migration of the periastron distance closer to the HZ region, introducing strong concurrent perturbations on a telluric planet and resulting in large oscillatory eccentricity excitations. Considering the case $e_{\rm giant} = 0.16$, we observe that nearly all particles exceed the CH boundaries. Only a telluric planet initially formed and located at $a_{\rm tel} \simeq 1.05~\textnormal{AU}$will marginally stay within the CHZ. This result is nearly independent of giant-planet mass (compare b-panel in Figs. 6 and 7).

6.7 HD 4208

Figure 8: Simulation snapshots for particles (telluric planets) orbiting HD 4208 considering different giant planet parameters. For each panel the time is indicated within each snapshot. A total of $2 \times 10^{3}$ particles are randomly distributed within the HZ with initial e=0. The continuous HZ is marked by the inner (0.74 AU) and outer (1.26 AU) boundaries. The innermost and outermost boundaries are located at 0.63 AU and 1.34 AU respectively. The total integration time is 106 years. Left 4-figure panel: $m_{\rm pl} = 1~M_{\rm Jup}$ and e=0.0. Right 4-figure panel: $m_{\rm pl} = 1~M_{\rm Jup}$ and e=0.16.

Open with DEXTER

	Figure 9: Same as Fig. 8, but for different choices of $m_{\rm pl}$. Left 4-figure panel: $m_{\rm pl} = 12~M_{\rm Jup}$ and e=0.0. Right 4-figure panel: $m_{\rm pl} = 16~M_{\rm Jup}$ and e=0.16.
Open with DEXTER

Similar dynamical simulations were carried out for HD 4208. In Figs. 8a, 8b, we observe telluric e-excitations at $a_{\rm tel} \simeq 1.06~\textnormal{AU}$ and $a_{\rm tel} \simeq 1.2~\textnormal{AU}$, corresponding to the 2:1 and 4:3 mean-motion commensurabilities. Particles within the 2:1 resonance are observed to be excited in eccentricity up to $e \simeq 0.18$ within the CHZ. Again the dynamics are mainly oscillatory in eccentricity taken at the same semi-major axis. Particles initially located with $a_{\rm tel} < 1.0~\textnormal{AU}$ are all on low-e orbits during the integration time span and confined to within the CHZ. In Fig. 8b the dynamical picture is more dramatic when considering HD 4208b in a high-e orbit. First, we observe a quick excitation of particle eccentricity at simulation start. The e-oscillations are characterized by large amplitudes with an initial transient phasing of nearby orbits. With time the frequency increases and e-oscillations tend to become narrower leading to a crossing of nearby orbits. This observation is useful when considering relative particle velocities and related accretion possibilities within the terrestrial region. This is discussed in more detail later. At mean-motion resonances, particle orbits are out of phase nearly instantly. Particles initially located within a narrow band $(\Delta a \simeq 0.06~\textnormal{AU})$ centered on $a \simeq 1.0~\textnormal{AU}$ (corresponds to 12% of the total range of the CHZ) have orbital eccentricities well confined within the CHZ exhibiting e-oscillations between $e_{\rm min} = 0.0$ and $e_{\rm max} \simeq 0.23$. For $a > 1.1~\textnormal{AU}$, particles are removed from the system by ejection (or accreted onto the giant planet). At the end of the 10⁶ year integration, the formation of a gap is observed with a cleared region within $a_{\rm min} \simeq 1.1~\textnormal{AU}$ and $a_{\rm max} \simeq 1.3~\textnormal{AU}$. From the MEGNO stability map (Fig. 4a), we observe that this region is dominated by chaotic particle motion resulting in orbital radial mixing with subsequent ejection or accretion. The survival of a few particles at $a < 1.3~\textnormal{AU}$ (3:2 mean-motion resonance) exhibits regular motion (oscillatory in eccentricity), and this location is indicated as being quasi-periodic within the corresponding MEGNO map (cf. Fig. 4a). In Figs. 9a, b we show additional simulation snapshots for various combinations in giant-planet parameters. The dynamics as shown in Fig. 9a are characterized by the 3:1 mean-motion resonance around $a \simeq 0.8~\textnormal{AU}$ and the removal of particles with initial $a_{\rm tel} > 1.1~\textnormal{AU}$. At $a_{\rm tel} \simeq 1.0~\textnormal{AU}$ we observe a large e-excitation of particles. This excitation is best explained by the close proximity of the giant planet. Furthermore, it is interesting to note the ``missing effect'' of the 5:2 mean-motion resonance at $a \simeq 0.9~\textnormal{AU}$. The corresponding MEGNO stability map (Fig. 4a) indicates the presence of this resonance. However, its effect is at most only at the 10%-level in eccentricity excitation and less effective when compared with the dynamical effect of the 3:1 mean-motion resonance. Finally, in Fig. 9b, we show the worst-case scenario in possible orbital parameters for HD 4208b considering $m_{\rm pl} = 16~M_{\rm Jup}, e=0.16$. In this case, almost all particles (with a > 0.9 AU) have been removed from the CHZ after 10⁶ years and roughly 10-15% of the total initial population of test particles survived. In summary the dynamical picture within the HZ of HD 4208 is characterized by mean-motion resonances and large-scale chaos leading to particle removal by either ejection or collision. This behavior is mainly explained by the proximity of the giant planet to the terrestrial habitable region for either eccentric giant planet orbits and/or high giant-planet mass.

   
6.8 Long-term and resonance dynamics of single particles

Short-term multi-particle integrations allowed us to perform a large parameter survey of various initial conditions, providing a dynamical picture of the short-term topology of the habitable phase space region. In general, nothing can be concluded on the subsequent evolution of single-particle orbital elements. In the three-body problem stable and confined orbits are only proven up to the integration time. It is possible that particles apparently exhibiting quasi-periodic regular motion (confined within the CHZ in orbit eccentricity, Figs. 6a, 7a) within the first 10⁶ years may suddenly undergo dynamical changes with large excursions in orbital parameters (sticky orbits). In order to learn more about the long term stability of the system, we are left with the question of the subsequent dynamical evolution of the particles. Since it is computationally time-consuming to probe each and every particle over billions of years, we integrated a few spot-test-objects over a time period of 1 billion years. For such a long time scale, it is no longer valid to model the telluric planet as a test particle. In the following we assign the telluric planet one Earth mass.

Figure 10: Long-term integrations of an Earth mass planet under the gravitational perturbation of HD 70642b for three locations in the HZ. The mass of the giant planet is $m_{\rm pl} = 16~M_{\rm Jup}$ with initial orbit eccentricity e=0.16. Each figure shows the semi-major axis ( left panel), eccentricity ( middle panel), and the polar representation ( right panel). All three initial conditions show bounded motion over 109 years. For practical reasons, we plot data points every 700 000 years in all panels. It is important to note that the apparent chaotic variation is a sampling effect and that aliased frequncies are introduced. In fact, all three orbits have a quasi-periodic MEGNO signature with $\vert\langle Y\rangle - 2.0\vert < 0.001$, as can be seen from Fig. 4b.

Open with DEXTER

Figure 11: Demonstrating $(\omega ,e)$-correlation for a massive telluric planet (in the HD 70642 system) with initial semi-major axis a) a0=1.045 AU and b) a0 = 1.123 AU. The latter corresponds to the nominal location of the 5:1 mean-motion resonance. We applied an on line low-pass filter with decimation (Quinn et al. 1991) to remove frequencies comparable to the telluric planets orbital period. The smoothed curves are overplotted on top of the unsmooted data points in the a) panel. It is seen that in the secular system the apsidal line is librating around the giant planets apsidal line. In panel b), we only plot the smoothed elements. The secular dynamics in the 5:1 mean-motion resonance is characterized by the alternation between libration and circulation of the apsidal line indicating strong chaotic behavior identified by the MEGNO technique as shown in panel c).

Open with DEXTER

We selected three initial conditions for the semi-major axis of a telluric planet. Comparing Fig. 6 with Fig. 7, it is seen that the maximum eccentricity excitation of a telluric planet (at a given semi-major axis) is only weakly dependent on the mass of HD 70642b, but strongly depends on the orbit eccentricity of the perturbing giant planet. In the following we adopt $m_{\rm pl} = 16~M_{\rm Jup}$ and e=0.16 for HD 70642b and refer to this setup as the ``worst-case-scenario'' to probe quasi-periodic motion (cf. Fig. 7b) during a 1 billion year integration. Considering this parameter setup will enable us to judge the long term dynamical evolution of a telluric planet perturbed by a giant planet on a low-eccentric orbit. For giant planet circular orbits, the accumulation of giant planet perturbations are expected to be less over time. If quasi-periodic motion is demonstrated for the worst-case-scenario setup, then it would provide a high degree of confidence in the conclusion that the bulk of particles in the terrestrial habitable region are rendered quasi-periodic (at least) when extrapolated to consider low-eccentric giant planet perturbations. In contrast, if the worst-case-scenario exhibits erratic and/or unbounded motion of the telluric planet, then we are left with the tedious (and almost impossible) task of checking the long-term dynamical evolution of individual planets initially started in the (C)HZ. In Fig. 10 we show the time evolution of orbital Kepler elements for three telluric planets orbiting HD 70642. The planets are followed over a total integration time of 10⁹ years using the 14th-order Radau algorithm as implemented in the MERCURY package. Comparing with the Burlisch-Stoer algorithm, we performed numerical tests indicating that Radau reduces the relative energy error by 1-2 orders of magnitudes with only little expense of extra CPU time.

All three telluric planets were initially started on a circular orbit at three different initial locations in the HZ. Figure 10 shows the time evolution of a,e along with the polar representation $(e_{\rm tel}~\cos~(\omega_{\rm tel}-\omega_{\rm pl}), e_{\rm tel}~\sin~(\omega_{\rm tel}-\omega_{\rm pl}))$ for each integration. Here $\omega_{\rm tel},\omega_{\rm pl}$ denotes the argument of pericenter (or apsidal longitude) of the tellus and the giant planet, respectively. For all three cases, it is demonstrated that the telluric orbit is bounded and exhibits quasi-periodic motion during the 10⁹ years, and we can conclude that a hypothetical telluric planet remains on a stable orbit out to a semi-major axis of at least a = 1.25 AU. This result allows us to conclude with high confidence the following. If a telluric planet was formed in the HZ of HD 70642, it would remain on a stable orbit for at least one billion years, as long as the giant-planet orbit eccentricity and mass remain lower than the worst-case scenario parameter setup. However, how likely terrestrial planets form by accretion processes under strong giant planet perturbations is not addressed in this work. Also we demonstrate instability. In an additional integration with $a_{0} \sim 1.30$ AU the telluric planet collided with the host star after $4.5 \times 10^{6}$ years.

In the following, we discuss the time evolution of the telluric planets eccentricity on short time scales outside and inside (next subsection) a mean-motion resonance. Figure 11a shows the time evolution of the eccentricity for a massive telluric planet initially started with a₀ = 1.045 AU over a 8000-year period. Data points are shown every 36 days. In the same figure, we overplot the time variation of the argument of periastron of the telluric planet at the same sampling rate. The figure not only shows a 2000-year period in the eccentricity variation, but also demonstrates a correlation between the apsidal line and eccentricity variation. In Fig. 11a the horizontal line indicates the initial argument of pericenter of HD 70642b and is equal to 277 degrees. From the figure we observe that, whenever the apsidal difference of the two planets is $(\omega_{\rm tel}-\omega_{\rm pl}) \sim 0$ (moment of apsidal alignment), the telluric orbit takes on its most elongated form $(e_{\rm tel} \sim 0.12)$. In Fig. 11a this is indicated by the crossing of the saw-tooth (blue) curve with the horizontal line. In addition, whenever the orbit is close to circular, the apsidal line of the telluric planet is unlocked and free to rotate (precess) for a short period of time and the apsidal difference is at maximum with $\omega_{\rm tel}-\omega_{\rm pl} \sim \pm 75^{\circ}$. To confirm this result, we repeated the simulations with the same initial conditions using the SWIFT (BS + RMVS) algorithms. It is important to stress that the apsidal alignment geometry between the two planets is not causing chaotic evolution or resulting in global orbital instability of the telluric planets orbit. The orbit in Fig. 11a is characterized by a quasi-periodic time evolution as indicated by the MEGNO factor (dashed line) in Fig. 11c. After having noticed this $(\omega ,e)$-correlation, we studied the nature of this correlation for different initial argument of pericenters (all other elements are similar to the simulation presented in Fig. 11a) of the giant planet HD 70642b. In particular, we looked at 3 additional, different values of $\omega_{\rm pl}=0,60,180$ degrees. In all three cases (results can be obtained upon request), the maximum eccentricity of the telluric planet occurred whenever the apses of the two planets were aligned for a short period of time. In addition, we repeated the simulation of Fig. 11a, but changed the initial eccentricity of the telluric planet to e = 0.15 and also changed the initial argument of pericenter of HD 70642b to $\omega = 75$degrees. The time evolution of the telluric planet's eccentricity and argument of periastron is shown in Fig. 12. Again, we observe that, when the apses of the planets are aligned (at 75 degrees), the eccentricity of the telluric (massive) planet is at maximum. In addition, we also observe that the minimum eccentricity of tellus occurs whenever the apse difference is 180 degrees. This event occurs in Fig. 12 when the argument of pariastron of tellus intersects the upper horizontal line (75 deg + 180 deg = 255 degrees).

6.9 Dynamics inside the 5:1 mean-motion resonance

The dynamical picture is changed when the two planets are in mean-motion resonance. In the following, we consider the5:1 mean-motion resonance of a massive telluric planet with HD 70642b. For now, we postpone any detailed dynamical analysis of the mean-motion dynamics and focus on the long term evolution inside this resonance. Initial conditions are similar to Fig. 11a, except that the telluric planet is started with semi-major axis a₀ = 1.123 AU, which corresponds to the 5:1 mean-motion resonance. The result of a 8000-year integration (15th order Radau) is shown in Fig. 11b. Compared to the dynamical behavior of the eccentricity in Fig. 11a, we now observe a more erratic time evolution of this element when the two planets participate in the 5:1 commensurability. Also, we observe an increase in eccentricity to $e \sim 0.18$ of the telluric orbit. We explain this increase in eccentricity because of the more frequent mutual planetary encounters at periapse alignment. However, the erratic zig-zag behavior of the eccentricity is a clear hint of chaotic time evolution, and it is mostly related to motion near the separatrix separating phase space in quasi-periodic (regular) and chaotic dynamics. Interestingly, and despite being chaotic, the initial conditions have been integrated for 10⁹ years (not shown in this work, but similar to Fig. 10), and during this period, the system remains bounded and hence dynamically stable. However, it should be noted that the $(\omega ,e)$-correlation is a well known dynamical phenomenon in solar-system, secular asteroid dynamics at mean-motion resonance (Lecar et al. 2001, p. 588). Here, we demonstrate that the same mechanism also applies to gravitationally interacting planets. In a later study, we plan to study the details of the various mean-motion resonances of this system by considering the time evolution of the corresponding critical angle around mean-motion resonances.

7 Summary, discussion, and conclusion

We have performed orbit stability simulations of telluric planets considering orbital perturbations of giant planets in two extrasolar systems. The planet HD 4208b is on a close orbit with $a \sim 1.7$ AU, and HD 70642b is on a wide orbit with $a \sim 3.3$ AU. Additional orbit parameters for both systems are shown in Table 1. For the two systems we studied the effect of gravitational perturbations on hypothetical telluric planets over a wide parameter range in giant-planet orbit eccentricities and mass. The simulations performed are divided into i) short term (10⁶ years) and ii) long term simulations (10⁹ years). Dynamical effects on short time scales were studied by direct multi-particle simulations in combination with the MEGNO technique. Short time scale simulations make it computationally feasible to study the time evolution of a large quantity of particles in the framework of a parameter survey study. The random distribution of 2000 telluric planets within the habitable region enabled us to obtain an instant picture of the dynamics in this region as a consequence of giant-planet perturbations. Each particle is considered to be isolated in the context of the three-body problem. To get an idea of the long-term dynamical evolution of a massive Earth-mass planet in HD 70642, we studied three different cases of the initial semi-major axis considering the dynamically ``worst-case'' combination of orbital parameters for HD 70642b. In the following we denote an unstable orbit to be a short-lived orbit (i.e. a particle surviving over 10⁶ years). Particles on moderate-to-high eccentric orbits are generally unstable. The reason for this instability is radial mixing: objects on increasingly eccentric orbits become crossing orbits, thereby visiting a larger region of phase space. This implies an increased probability of close encounters or even direct collisions with the giant planet/host star. In effect, mean-motion resonances provide an effective mechanism for particle removal/ejection.

Figure 12: Time evolution of eccentricity (sine curve) and argument of periastron of tellus with identical initial condition as given in Fig. 11a, except that the telluric planet was initially started on an orbit with e0=0.15 and $\omega _{\rm pl} = 75$ degrees. The lower horizontal line corresponds to $\omega = 75$ degrees and the upper line corresponds to $\omega = 75~\textnormal{deg} + 180~\textnormal{deg} = 255$ degrees.

Open with DEXTER

7.1 MEGNO technique

Results from the calculation of the MEGNO factor are in close agreement with previous results (Asghari et al. 2004) published independently and obtained by the Kolmogorov entropy method of quantifying chaotic dynamics. This is the first time a direct comparison with a different fast-indicator technique has been done establishing some confidence in their use. Combining the calculation of MEGNO maps with direct particle simulations enabled us to correlate dynamical features in the MEGNO maps with eccentricity excitations of particle orbits within the habitable zone. Most of these features are identified as corresponding to the location of orbital mean-motion resonances with $\langle Y \rangle \gg 2.0$. The main dynamical effect of a mean-motion resonance is additional pumping of orbital eccentricities and the perturbation effects at resonances generally starts to operate at high giant-planet orbit eccentricity and mass. The identification of chaotic regions within a MEGNO map implies not necessarily dynamically unstable orbits for a given set of initial conditions. By dynamical instability, we mean the event of large excursions of the orbit eventually resulting in orbit crossing and/or close encounters (one object entering the Hill sphere of another object) with other planets. Based on the results of this work, we would not interpret every chaotic feature observed in a MEGNO map as indicative of or synonymous with being dynamically unstable. The problem is that no direct link exists between the degree of instability and $\langle Y \rangle$, which is related to the Lyapunov time. As an example, the solar system's Lyapunov time is estimated as 5 Myr, indicating that it is chaotic, but still all planetary orbits are dynamically stable over a period of at least 4 Gyr (Lecar et al. 2001) - a property of the system termed stable chaos (Milani et al. 1997; Milani & Nobili 1992) to indicate a chaotic system (short Lyapunov time) exhibiting stable motion on long time scales. From theory, a dynamical system has a spectrum of Lyapunov exponents (possibly complex eigenvalues) each associated to a given eigenvector of the system. Each Lyapunov exponent describes the rate of change of its corresponding eigenvector in time, and the number of positive (vanishing) Lyapunov exponents indicates the number of independent directions in phase space along which the orbit exhibits chaotic (quasi-periodic) behavior (Morbidelli 2002). In the general case, the time evolution of a Kepler element is given by a linear combination of all Lyapunov exponents. To probe whether an initial condition is in a chaotic region of phase space, it has proven sufficient to only determine the largest (maximum) Lyapunov exponent (MLE) in the spectrum of eigenvalues. The MEGNO factor provides an estimate of the MLE characterizing the rate of exponential divergence in the system. However, from the MEGNO factor alone, it is not obvious which Kepler element of the telluric planets orbit will exhibit chaotic behavior and the present work (co-planar orbits) suggests that the typical chaotic variation is mainly found in the time variation of the eccentricity and apsidal line. Hence a stability analysis should include a discussion of these parameters. It is important to note that the apparent chaotic variation in the semi-major axis in Fig. 10 is an artificial effect due to the sampling frequency of data output. After pointing out these remarks, the construction of MEGNO maps is a powerful way to obtain a quick overview of dynamically interesting features indicating chaotic motion of a given phase space region. However, we recommend that a MEGNO map should always be supplemented with additional dynamical information showing the change in other orbital quantities (i.e. semi-major axis, eccentricity and apsidal line variation) for a given problem to enhance the physical information contained in $\langle Y \rangle$ on the orbital time evolution.

7.2 HD 70642

We conclude the following about the stability of telluric planets in the habitable region of HD 70642. We demonstrated the long term (10⁹ years) stability of 3 telluric orbits with different initial semi-major axis. The stability was probed by adopting the ``worst-case-scenario'' combination of observed giant planet orbit parameters. However, it was found that the telluric planet's orbit eccentricity exhibits large variations on a 2000-year period, and it is speculative at this point whether a telluric planet would be able to develop some form of biology under such dynamical variations. From the long term simulations we then concluded that global stability in the habitable region is obtained for HD 70642b on a circular orbit. In the case of a telluric planet in the CHZ, the mutual distance between the orbits of the planets is only 2.26 AU, and we observed only small variations in the telluric orbit eccentricity. Furthermore, it was found that the giant-planet mass is the last important parameter when studying its dynamical influence on the terrestrial region (when circular orbits were considered). Even a mass of $m_{\rm gp} = 16~M_{\rm Jup}$ orbiting the host star on a circular orbit (at distance 3.3 AU) would render the habitable region stable.

7.3 HD 4208

The dynamical picture of the terrestrial region in HD 4208 is completely different. Because of the giant planets' close proximity to the habitable zone, stronger gravitational perturbations are introduced when compared to HD 70642. Simulations suggest that the only permissable parameter combination for HD 4208b is a minimum mass giant planet in a low-e orbit. This dynamical setup still introduces a strong mean-motion resonance operating well inside the CHZ (Fig. 8a) at $a \simeq 1.05$ AU. However, the mass of HD 4208b seems to be dynamically significant now. Retaining the giant planet on a circular orbit, but increasing the mass by an order of magnitude has the effect of particle removal with $a \simeq 1.10$ AU and onwards. Considering the dramatic cases in which HD 4208b is on an eccentric orbit with mass parameters either $m_{\rm pl} = 1~M_{\rm Jup}$ or $16~M_{\rm Jup}$, shows that telluric planets crosses the CHZ inner and outer boundaries. For the high-mass range, no telluric planet survives over 10⁶ years (Fig. 9b). Numerical simulations performed in this work suggest that extrasolar giant planets should be on circular orbits with a semi-major axis that not much smaller than a = 3.3 AU. If no additional perturbations are present, such a dynamical configuration would cause a terrestrial Earth-mass planet to remain on a stable and bounded orbit within the CHZ. At the time of writing, half a dozen of the 230 exoplanetary systems have a > 3.3 AU with moderate-to-high orbital eccentricities spanning the range e = 0.16 to 0.70. In view of this work, these systems seem to be hostile to the dynamical environment covering the habitable terrestrial region. We suggest that future observations of already known exoplanets should aim to further constrain and minimize the uncertainty range of the orbital eccentricity. This parameter seems to have the most dynamical significance when considering the dynamical stability properties of telluric planets in the habitable terrestrial region.

Acknowledgements

This work benefited greatly from suggestions by the referee who pointed out and clarified important misleading aspects in the original submitted manuscript! This research project was supported by the Danish Natural Science Research Council, FNU. Numerical simulations were performed at the supercomputing facility at the University of Copenhagen (DCSC-KU) run on behalf of the Danish Center for Scientific Computing. T.C.H. would like to acknowledge Åke Nordlund for providing access to DCSC. K.G. is supported by the Polish Ministry of Science, Grant 1P03D 021 29.

References

• Agnor, C., Canup, R., & Levison, H. 1999, Icarus, 142, 219 [NASA ADS] [CrossRef]
• Asghari, N., Broeg, C., Carone, L., et al. 2004, A&A, 426, 353 [NASA ADS] [CrossRef] [EDP Sciences]
• Beaulieu, J.-P., Bennett, D. P., Fouque, P., et al. 2006, Nature, 439, 437 [NASA ADS] [CrossRef] (In the text)
• Benettin, G., Galgani, L., & Strelcyn, J.-M. 1976, Phys. Rev. A, 14, 2338 [NASA ADS] [CrossRef] (In the text)
• Breiter, S., Melendo, B., Bartczak, P., & Wytrzyszczak, I. 2005, A&A, 437, 753 [NASA ADS] [CrossRef] [EDP Sciences]
• Carter, B., Butler, R. P., Tinney, C., et al. 2003, ApJ, 593, L43 [NASA ADS] [CrossRef] (In the text)
• Chambers, J. 1999, MNRAS, 304, 793 [NASA ADS] [CrossRef]
• Chambers, J. 2001, Icarus, 152, 205 [NASA ADS] [CrossRef]
• Chambers, J. 2003, Lunar and Planetary Science, 34, 2000C [NASA ADS]
• Chambers, J., & Migliorini, F. 1997, BAAS, 27-06-P
• Chambers, J. E., & Wetherill, G. W. 1998, Icarus, 136, 304 [NASA ADS] [CrossRef]
• Cincotta, P., & Simó, C. 2000, ApJS, 147, 205 [NASA ADS]
• Cincotta, P., Giordano, C., & Simó, C. 2003, Physica D, 182, 151 [NASA ADS] [CrossRef]
• Cuntz, M., von Bloh, W., Bounama, C., & Franck, S. 2003, Icarus, 162, 214 [NASA ADS] [CrossRef] (In the text)
• Franck, S., Block, A., von Bloh, W., et al. 2000, Planet. Space Sci., 48, 1099 [NASA ADS] [CrossRef] (In the text)
• Gozdziewski, K. 2001, A&A, 378, 569 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Gozdziewski, K. 2002, A&A, 393, 997 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Gozdziewski, K. 2003, A&A, 398, 1151 [NASA ADS] [CrossRef] [EDP Sciences]
• Gozdziewski, K. 2004, ApJ, 610, 1093 [NASA ADS] [CrossRef]
• Gozdziewski, K., & Maciejewski, A. 2001, ApJ, 563, L81 [NASA ADS] [CrossRef]
• Hairer, E., Norsett, S., & Wanner, G. 1993, Solving ordinary differential equations I, nonstiff problems, 2nd Ed., Springer Series in Computational Mathematics (Springer-Verlag), see also http://www.unige.ch/math/folks/hairer/ (In the text)
• Hart, M. 1979, Icarus, 37, 351 [NASA ADS] [CrossRef]
• Hilditch, R. 2001, An Introduction to Close Binary Stars (Cambridge University Press) (In the text)
• Huang, S.-S. 1959, Pub. A. S. P., 71, 421H [NASA ADS] [CrossRef]
• Jones, B., & Sleep, P. 2002, A&A, 393, 1015 [NASA ADS] [CrossRef] [EDP Sciences]
• Jones, B., Sleep, P., & Chambers, J. 2001, A&A, 366, 254 [NASA ADS] [CrossRef] [EDP Sciences]
• Jones, B., Underwood, D., & Sleep, P. 2005, ApJ, 622, 1091 [NASA ADS] [CrossRef]
• Kasting, J., Whitmire, D., & Reynolds, R. 1993, Icarus, 101, 108 [NASA ADS] [CrossRef]
• Kasting, J. F., & Catling, D. 2003, ARA&A, 41, 429 [NASA ADS] [CrossRef]
• Kokubo, E., & Ida, S. 2000, Icarus, 143, 15 [NASA ADS] [CrossRef]
• Lecar, M., Franklin, F., Holman, M., & Murray, N. 2001, ARA&A, 39, 581 [NASA ADS] [CrossRef] (In the text)
• Mayor, M., & Queloz, D. 1995, Nature, 378, 355 [NASA ADS] [CrossRef] (In the text)
• Menou, K., & Tabachnik, S. 2003, ApJ, 583, 473 [NASA ADS] [CrossRef]
• Milani, A., & Nobili, A. M. 1992, Nature, 357, 569 [NASA ADS] [CrossRef]
• Milani, A., Nobili, A. M., & Knezevic, Z. 1997, Icarus, 125, 13 [NASA ADS] [CrossRef]
• Morbidelli, A. 2002, Modern Celestial Mechanics - Aspects of Solar System Dynamics (Taylor & Francis) (In the text)
• Murray, C., & Dermott, S. 1999, Solar System Dynamics (Cambridge University Press) (In the text)
• Noble, M., Musielack, Z., & Cuntz, M. 2002, ApJ, 572, 1024 [NASA ADS] [CrossRef]
• Nordström, B., Mayor, M., Andersen, J., et al. 2004, A&A, 418, 989 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• O'Brian, D., Morbidelli, A., & Levison, H. 2006, Icarus, 184, 39 [NASA ADS] [CrossRef]
• Quinn, T. R., Tremaine, S., & Duncan, M. 1991, AJ, 101 (6), 2287 (In the text)
• Quintana, E., & Lissauer, J. 2006, Icarus, 185, 1 [NASA ADS] [CrossRef]
• Quintana, E., Adams, F., Lissauer, J., & Chambers, J. 2007, ApJ, 660, 807 [NASA ADS] [CrossRef]
• Raymond, S., Quinn, T., & Lunine, J. 2004, Icarus, 168, 1 [NASA ADS] [CrossRef]
• Sándor, Z., Süli, Á., Érdi, B., Pilat-Lohinger, E., & Dvorak, R. 2007, MNRAS, 375, 1495 [NASA ADS] [CrossRef]
• Süli, Á., Dvorak, R., & Érdi, B. 2007, AN, 328, 781 [NASA ADS]
• Tabachnik, S., & Tremaine, S. 2002, MNRAS, 335, 151 [NASA ADS] [CrossRef] (In the text)
• Tancredi, G., Sánchez, A., & Roig, F. 2001, AJ, 121, 1171 [NASA ADS] [CrossRef] (In the text)
• Thébault, P., Marzari, F., & Scholl, H. 2002, A&A, 384, 594 [NASA ADS] [CrossRef] [EDP Sciences]
• Udry, S., Bonfils, X., Delfosse, X., et al. 2007, A&A, 469, 43 [NASA ADS] [CrossRef] (In the text)
• Vogt, S., Butler, P., Marcy, G., et al. 2001, ApJ, 568, 352 [NASA ADS] [CrossRef] (In the text)
• Wetherill, G. 1996, Icarus, 119, 219 [NASA ADS] [CrossRef]
• Wisdom, J. 1983, Icarus, 56, 51 [NASA ADS] [CrossRef] (In the text)
• Zucker, S., & Mazeh, T. 2001, ApJ, 562, 1038 [NASA ADS] [CrossRef] (In the text)