Introduction
About one and a half decades have passed since the detection of the first planet hosted by a solar-type star other than the Sun. The field of research on extra-solar planets and associated aspects of exobiology has considerably matured in the meantime. One important aspect with significant ramifications for the search of life outside the Solar System concerns the identification of Earth-mass planets in stellar habitable zones (HZs), including detailed investigations regarding their orbital stability.
In principle, there are two categories of studies available with respect to this type of research. The first category encompasses studies of theoretical star–planet systems, especially binary and multiple stellar systems, in order to decipher the orbital stability of the planets (e.g., Pilat-Lohinger & Dvorak Reference Pilat-Lohinger and Dvorak2002; Musielak et al. Reference Musielak, Cuntz, Marshall and Stuit2005; Cuntz et al. Reference Cuntz, Eberle and Musielak2007; Szenkovits & Makó Reference Szenkovits and Makó2008, among a large body of other literature). The second category deals with observed star planet systems, typically consisting of a star in the centre and a small number of extra-solar giant planets (EGPs), with (a) hypothetical Earth-mass planet(s) taken into consideration. Detailed studies in the framework of this latter category were given by Gehman et al. (Reference Gehman, Adams and Laughlin1996), Jones et al. (Reference Jones, Sleep and Chambers2001), Noble et al. (Reference Noble, Musielak and Cuntz2002), Cuntz et al. (Reference Cuntz, von Bloh, Bounama and Franck2003), Menou & Tabachnik (Reference Menou and Tabachnik2003) and many others. These papers explore which systems are in principle able to host terrestrial planets, especially planets in the stellar HZs. Clearly, this latter aspect is of extreme importance to the flourishing field of exobiology, although some complications need to be taken into account, especially if stellar evolutionary timescales become part of the overall picture (e.g., Jones et al. Reference Jones, Underwood and Sleep2005, Reference Jones, Sleep and Underwood2006; Lopez et al. Reference Lopez, Schneider and Danchi2005; von Bloh et al. Reference von Bloh, Cuntz, Schröder, Bounama and Franck2009).
There is one aspect, however, that has received little attention from previous studies, which is the dependency of the ejection time on the starting position of the planet, especially in systems where the orbital stability is not guaranteed. An existing example includes the work by Noble et al. (Reference Noble, Musielak and Cuntz2002). They investigated the orbital stability of terrestrial planets inside the HZs of three stellar systems, which are 51 Peg, 47 UMa and HD 210277. The centre stars are relatively similar to the Sun concerning mass, spectral type and effective temperature. Orbital stability was attained for terrestrial planets in the HZ of 51 Peg and a significant portion of the HZ of 47 UMa; however, no orbital stability was found for hypothetical Earth-mass planets in the HZ of HD 210277. In this case, a Jupiter-type planet crosses into the stellar HZ, thus effectively thwarting habitability for this system.
In regard to HD 210277, Noble et al. (Reference Noble, Musielak and Cuntz2002) varied the starting position of the terrestrial planet, originally placed at the centre of the HZ, in 30° increments. They found that the ejection time of the terrestrial planet changed between 0.81 yr and 61.97 yr with little systematics in terms of the attained variation for the set of orbital stability simulations. It is the main motivation of the present work to extend this type of study to other systems, namely HD 20782 and HD 188015. We intend to identify statistical properties of the ejection times of the terrestrial planet, depending (if possible) on the assumed masses of the giant planets and on the starting positions of the giant and terrestrial planets, including the relative angles (phases) between these objects.
Note that our current study exhibits inherent similarities to the previous work by Fatuzzo et al. (Reference Fatuzzo, Adams, Gauvin and Proszkow2006). They provided a detailed statistical analysis of ejection times, encompassing about 500 orbital configurations, for Earth-mass planets in binary systems. This study included variations of the companion mass, the companion eccentricity, the companion periastron, and the planet's inclination angle i relative to the stellar binary plane. The authors identified distributions of ejection times τs (in years), referred to as T Ej in the present paper. With respect to log τs they typically turned out to be Gaussian in width, although the longest-lived systems displayed substantial (non-Gaussian) tails. This latter result indicates that seemingly unstable systems may still hold some merit in the study and realization of planetary habitability.
In the following, we describe our theoretical approach. We will discuss the adopted methods and the system parameters concerning HD 20782 and HD 188015, and we will outline the adopted definition of habitability. Thereafter, we will describe our results with a focus on the statistical analysis of the planetary ejections. In addition, we will discuss selected examples of simulations for HD 20782 and HD 188015 with different starting conditions for the giant and Earth-mass planets. Finally, we will present our conclusions.
Theoretical approach
Methods and system parameters
Our study is based on the two star–planet systems HD 20782 and HD 188015. Each system contains a Jupiter-type planet with a minimum mass of 1.78 and 1.50 M J, respectively; see Table 1 for information on the adopted stellar and planetary parameters with all symbols having their usual meaning. Moreover, the central stars of the two systems are very similar to our Sun. The masses of the two stars are 0.98±0.05 and 1.09±0.05 M ⊙ and the stellar effective temperatures are 5758±44 and 5746±44 K, respectively. Note that some orbital properties of this system were studied previously by Cuntz & Yeager (Reference Cuntz and Yeager2009) in their effort to analyse the validity of the Hill radius criterion for planetary ejections from the system.
Table 1. Stellar and planetary parameters

Note. Data as given by Butler et al. (Reference Butler2006) and references therein, notably Valenti & Fischer (Reference Valenti and Fischer2005).
Note that the two star–planet systems HD 20782 and HD 188015 are fundamentally different compared to the solar case. In these systems the giant planet either remains in or crosses into the stellar HZ, thus effectively thwarting the possibility of habitable Earth-type planets; see Noble et al. (Reference Noble, Musielak and Cuntz2002) for previous simulations of the HD 210277 system also showing this type of behaviour. In case of HD 188015, the orbit of the giant planet is essentially circular, with a perigee and apogee of 1.04 and 1.37 AU, respectively, both located in the stellar HZ (see Fig. 1). On the other hand, the orbit of the giant planet in the system of HD 20782 is extremely elliptical, noting that its perigee is at 0.10 AU and its apogee is at 2.63 AU.

Fig. 1. Orbits of the giant planets for the two represented systems, which are HD 188015 (dashed line) and HD 20782 (narrow dotted line). We also give the HZs (grey area) for the two systems. The two HZs are assumed to have identical extents (see text).
For our simulations concerning each system, we consider both the observed giant planet and a hypothetical terrestrial planet of one Earth-mass, i.e. 3.005×10−6M ⊙, which allows us to execute a grid of model simulations. The method of integration uses a fourth-order Runge–Kutta integration scheme (Press et al. Reference Press, Flannery, Teukolsky and Vetterling1989). The code has been extensively tested against known analytical solutions, including the two-body and restricted three-body problem (see Noble et al. Reference Noble, Musielak and Cuntz2002; Cuntz et al. Reference Cuntz, Eberle and Musielak2007; Eberle et al. Reference Eberle, Cuntz and Musielak2008, for detailed results). In the framework of our simulations that are limited to about 103 yr due to the ejection of the Earth-mass planet from the HZ, we apply a time-step of 10−4 yr for the integration scheme. This value is found to be fully appropriate. A case study comparing the planetary orbits based on two different integration time-steps, i.e. 10−4 and 10−5 yr, is given in Fig. 2 with no discernable differences found. The initial conditions for the orbit of the Earth-mass planet is chosen such that it is assumed to be in a circular orbit about the star, although it is evident that it will be significantly affected immediately by the giant planet preventing its circular motion. The initial coordinates and velocities of all objects, i.e. x 0, y 0, v x0 and v y0, as well as their masses are also set with an absolute precision (decimal places) of 10−4.

Fig. 2. Orbital stability simulations for an Earth-mass planet in the system HD 188015 using two different precisions, which are 10−4 (top) and 10−5 (bottom), with a simulation time of 103 yr. The giant planet has initially been placed at perigee. The starting position of the Earth-mass planet is given as δ=0.5 at a starting angle of θ0=180°.
Stellar habitable zones
Since the effective temperatures and luminosities of HD 20782 and HD 188015 agree with the Sun within the uncertainty bars, we will assume that the HZs of the two stars are identical to the solar HZ. Therefore, following Kasting et al. (Reference Kasting, Whitmire and Reynolds1993), the inner and outer edges of the HZs are assumed to be located at 0.84 and 1.67 AU, respectively, from the star. The underlying definition of habitability is footed on the assumption that liquid surface water is a prerequisite for life, a key concept that is also the basis of ongoing and future searches for extrasolar habitable planets (e.g., Catanzarite et al. Reference Catanzarite, Shao, Tanner, Unwin and Yu2006; Cockell et al. Reference Cockell2009). Note that the numerical evaluation of these limits is based on an Earth-type planet with a CO2/H2O/N2 atmosphere.
The inner limit of habitability is set by the loss of water from the upper planetary atmosphere through photodissociation and subsequent escape of hydrogen to space associated with a run-away greenhouse effect. The outer limit of habitability is given by the maximum greenhouse effect (Kasting et al. Reference Kasting, Whitmire and Reynolds1993; Underwood et al. Reference Underwood, Jones and Sleep2003); here a surface temperature of 273 K can be maintained by a cloud-free CO2 atmosphere. We would also like to point out that concerning the outer edge of habitability, even less conservative limits have been proposed in the meantime (e.g., Forget & Pierrehumbert Reference Forget and Pierrehumbert1997; Mischna et al. Reference Mischna, Kasting, Pavlov and Freedman2000). They are based on the assumption of relatively thick planetary CO2 atmospheres and invoke strong backwarming that is further enhanced by the presence of CO2 crystals and clouds. However, as these limits, which can be as large as 2.4 AU for solar-like stars, depend on distinct properties of the planetary atmosphere, they are not relevant for our study. Moreover, the significance of these limits has recently been challenged based on detailed radiative transfer simulations (Halevy et al. Reference Halevy, Pierrehumbert and Schrag2009).
Results and discussion
Statistics of terrestrial planet ejections
As part of our study we performed various sets of simulations both for HD 20782 and HD 188015 aimed at exploring the onset of orbital instability for the terrestrial planet and its ejection from the stellar HZ. At the start of each of our simulations, we place the giant planet either at perigee (3 o'clock position) or apogee (9 o'clock position) relative to the host star. We pursue a set of simulations for each system with different initial positions for the terrestrial planet inside the stellar HZ also noting that its velocity components are set up to orbit the star on a circular orbit. The distances of the terrestrial planets from the host star are chosen according to r j=r HZ,i+δj⋅(r HZ,o−r HZ,i) with r HZ,i as inner radius of the HZ, r HZ,o as outer radius of the HZ and δj=0.1, 0.5 and 0.9 (j=1, 2, 3) as distance parameter.
For a given starting position of the giant planet (perigee or apogee) and starting distance (defined by δ) of the terrestrial planet, we perform a total of 36 simulations that are distinguished by the starting angle θ0 of the terrestrial planet. Its position is varied counterclockwise in 10° increments with 0° corresponding to the perigee of the giant planet. As expected, the orbit of the Earth-mass planet in many of our simulations is almost immediately considerably disturbed by the gravitational pull of the giant planet, leading to its ejection from the HZ on a relatively short timescale compared to timescales needed for establishing astrobiological settings.
The main focus of our study is to provide a detailed investigation of the ejection times of an Earth-mass planet from the HZs of HD 20782 and HD 188015 to elucidate the associated statistical properties. Tables 2–5 show the arithmetic mean and the median of the ejection times, i.e. T Ejmean and T Ejmedian (see the Appendix), assuming N*=N for the different sets of models. It is found that T Ejmean is usually noticeably larger than T Ejmedian; see Table 6 for a logarithmic evaluation of the ejection times. This behaviour is due to the skewness in the distribution. Note that the skewness of the T Ej distributions is considerably positive in 9 out of 16 sets of models, and considerably negative in only one set of models (see Table 7).
Table 2. Ejection times TEj mean (in years) for HD 20782

Note. f(m) denotes the factor by which the minimum mass of the EGP is multiplied.
Table 3. Ejection times TEj mean (in years) for HD 188015

Note. f(m) denotes the factor by which the minimum mass of the EGP is multiplied.
Table 4. Statistical properties of ejection times (in years) for HD 20782

Note. f(m) denotes the factor by which the minimum mass of the EGP is multiplied.
Table 5. Statistical properties of ejection times (in years) for HD 188015

Note. f(m) denotes the factor by which the minimum mass of the EGP is multiplied.
Table 6. Logarithmic evaluation of ejection times

Note. f(m) denotes the factor by which the minimum mass of the EGP is multiplied. The unit of the ejection times is years.
Table 7. Summary of skewness

Note. If the skewness is evaluated based on log T Ej, Skew is negative for all simulations but one, which is the model series for HD 20782 with the Earth-mass planet placed at δ=0.5 and the giant planet with f(m)=1.3 originally at the apogee position.
Moreover, it is found that the ejection times are consistently larger if the starting position of the terrestrial planet is placed at δ=0.5, i.e. in the middle of the HZ, instead of δ=0.1 or 0.9, as expected. The difference is typically a factor of 2 for HD 20782 in regard to the plurality of models and a factor of 10 for HD 188015. This astonishingly large difference in the factor is caused by the variation in strength of gravitational interaction between the EGP and the Earth-mass planet in the two systems, noting that the orbit of the EGP in the HD 20782 system is highly elliptical, whereas in the HD 188015 system it is not. Apparently, for large eccentricities the variation in T Ej is less dependent on the planet's starting position.
It is also found that the standard deviation SD(T Ej) is relatively large for the different sets of models. In fact for most sets of models, the immediate ejection of the Earth-mass planet is consistent with 1.5 SD or less. On the other hand, for many of the simulations there appears to be a significant ‘positive tail’ concerning the ejection times, indicating longer durations of orbital stability. Concerning δ=0.5, the case of optimal orbital stability, the ejection times (as gauged by T Ejmean) are generally longer for HD 188015 than for HD 20782. The likely reason is that in the latter star–planet system, the extremely elliptical orbit of the EGP leads to a larger number of relatively close encounters with the Earth-mass planet, resulting in its ejection from the HZ. We also investigated the behaviour of the excess kurtosis XKur (see the Appendix) for our models, but we were unable to identify a well-pronounced pattern.
We also explored how the increments concerning the starting angle θ0 of the Earth-mass planets affect our results, particularly the attained values of T Ejmean (see Tables 2 and 3). Here the two extremes are: increments of 90° (four models) and increments of 10° (36 models) as used both in regard to the HD 20782 and HD 188015 simulations. Surprisingly, we found that the choice of increment is not very important for the manifestation of T Ejmean, as well as T Ejmedian. Nevertheless, model evaluations based on increments of 90° will miss out on almost all details of the statistical analysis for the planetary ejection times.
One of these properties is that a small change in the starting angle θ0 can result in a drastically different outcome. This behaviour is depicted in Figs 3 and 4 for HD 20782 and HD 188015, respectively, which consider model simulations where the EGP is placed at perigee or apogee, as well as model simulations where the mass of the EGP is increased by 30%. A typical example for the star–planet system HD 20782 occurs if the giant planet is placed at perigee and the Earth-mass planet is placed at δ=0.5. If the Earth-mass planet has a starting angle of θ0=110°, the ejection time is 96.2 yr, but for a starting angle of θ0=120°, it is only 0.48 yr (i.e. a difference of a factor of 200). This phenomenon is akin to the mathematical behaviour of phase space that contains ‘sensitive dependence on initial conditions’ (e.g., Smith & Szebehely Reference Smith and Szebehely1993); see the Case studies section for a detailed depiction of examples.

Fig. 3. Display of log T Ej for an Earth-mass planet dependent on its starting position, angle θ0, for HD 20782. The top figure shows the comparison between the models with starting positions δ=0.1 (solid line; squares), δ=0.5 (dashed line; circles) and δ=0.9 (dotted lines; triangles) with the giant planet originally placed at perigee. The middle figure refers to the same setting, except that the giant planet is now originally placed at apogee. The bottom figure depicts models with the Earth-mass planet at a starting position of δ=0.5. The mass of the giant planet has now been increased by 30% and is originally placed at perigee (solid line; squares) or at apogee (dashed line; circles).

Fig. 4. Same as Fig. 3, but for HD 188015.
Another feature visible in Fig. 3, depicting the dynamics of the HD 20782 star–planet system, is that, with the EGP originally placed at perigee, the ejection time for the hypothetical Earth-mass planet with the starting position of the EGP at perigee is especially short for various particular starting angles. This effect is caused by the periodicity of the EGP orbit that is notably elliptical in this system. Conversely, for the HD 188015 star–planet system, the essentially circular case, there appears to be a maximum in the ejection time distributions for the Earth-mass planet near ±60° apart from the position of the Jupiter-type planet. This appears to indicate enhanced stability for Trojan analogues in this system, associated with the Lagrangian points L4 and L5 (e.g., Roy Reference Roy2005); note that a recent study on the stability region of hypothetical habitable Trojan planets for a broad variety of systems has been given by Schwarz et al. (Reference Schwarz, Dvorak, Sűli and Érdi2007).
In Table 6 we provide a logarithmic evaluation of the ejection times, which is also motivated to better compare our findings with the previous work by Fatuzzo et al. (Reference Fatuzzo, Adams, Gauvin and Proszkow2006), who obtained a detailed statistical analysis of ejection times for terrestrial planets from binary systems while using a logarithmic scale. In our studies, there is no qualitative difference between evaluations based on log T Ejmean or on (log T Ej)mean for either star–planet system. However, there is a significant difference concerning the skewness. It is found that if the skewness analysis is applied to log T Ej rather than T Ej, the value of Skew (see the Appendix) is negative for all models of HD 20782, but one, and for all models of HD 188015 (see Table 7). This type of result is actually to be expected for an analysis provided on a logarithmic scale which, for example, for a uniform distribution would provide a negative skew (i.e. a tail to the left with the bulk of the distribution to the right). It is noteworthy that the study of ejection times in binary systems by Fatuzzo et al. (Reference Fatuzzo, Adams, Gauvin and Proszkow2006) arrives at a positive skew even though this type of study occurs within a logarithmic frame. Thus, in our case the positive tail of ejection times is considerably less pronounced than in the Fatuzzo et al. investigation.
Finally, our study shows that the ejection times for the Earth-mass planets from the HZs of HD 20782 and HD 188015, originally placed at the centre of the respective HZ, vary by a factor of ~200 and ~1500, respectively, depending on the starting positions of the giant and terrestrial planets. This result is featured by Table 8, which depicts minimum and maximum ejection times log T Ej for HD 20782 and HD 188015 concerning the different sets of simulations. It is also found that if the mass of the EGP is increased by 30%, the variation in ejection times for HD 188015 increases to a factor of ~6000. The reason for this latter behaviour is two-fold. First, in some cases, especially if the EGP is outside the stellar HZ, the higher mass of the EGP results in a faster ejection of the Earth-mass planet from the HZ owing to the increased gravitational attraction between the two objects. Secondly, however, in other cases the Earth-mass planet appears to be bound more closely to the star–planet system if the increased gravitational attraction between the EGP and the Earth-mass planet is able to initiate such a behaviour.
Table 8. Minimum and maximum ejection times log TEj for HD 20782 and HD 188015

Note. f(m) denotes the factor by which the minimum mass of the EGP is multiplied. The unit of the ejection times is years.
Case studies
As previously discussed, we performed different sets of numerical computations to identify the influence of the starting positions of the giant planet and Earth-mass planet on the outcome of the simulations. In the following, we will showcase various different simulation runs to provide comparisons of the behaviour of the Earth-mass planet for various initial conditions. By default, most of these runs are depicted for a simulation time of 103 yr. Figures 5–7 are for the star–planet system HD 20782. In this system, the Jupiter-type planet is on an extremely elliptical orbit, which typically results in large variations of the gravitational attraction toward the Earth-mass planet. First, we will examine two examples for the case δ=0.1 and then additional examples for δ=0.5 and 0.9. In order to demonstrate how the outcomes vary based on altered initial conditions, we will only change one initial condition between the two comparison figures.

Fig. 5. Case studies of orbits for a terrestrial planet (dotted line) in the HZ (grey area) of HD 20782. The starting positions of the Earth-mass planets are given as δ=0.1 at an angle of θ0=20°. The giant planet (solid line) was initially placed at perigee (top figure) and apogee (bottom figure). The Earth-mass planets were ejected after 23.9 and 7.6 yr, respectively.

Fig. 6. Case studies of orbits for a terrestrial planet (dotted line) in the HZ (grey area) of HD 20782. The starting positions of the Earth-mass planets are given as δ=0.5 at an angle of θ0=10°. The simulations also refer to different cases for the giant planet (solid line), originally placed at perigee (left column) and apogee (right column) and with minimum mass (top row) and its mass increased by 30% (bottom row), respectively. The Earth-mass planets were ejected after 8.0, 20.7, 34.9 and 18.7 yr, respectively (top left, top right, bottom left, bottom right).

Fig. 7. Case studies of orbits for a terrestrial planet (dotted line) in the HZ (grey area) of HD 20782. The starting positions of the Earth-mass planets are given as δ=0.9 at an angle of θ0=160°. The giant planet (solid line) was initially placed at perigee (top figure) and apogee (bottom figure). The Earth-mass planets were ejected after 51.5 and 2.5 yr, respectively.
The first set that we have chosen to compare is for the apogee and perigee starting positions of the giant planet (Fig. 5). As can be seen from the set of plots, this change in starting position by 180° for the giant planet produces a radically different outcome. All of the initial conditions for the terrestrial planet placed at a starting angle θ0 of 20° remain the same for both figures. When the giant planet is placed at perigee the terrestrial planet is ejected from the stellar HZ after only 23.9 yr. If we compare this to the case of the giant planet starting at apogee we see the terrestrial planet being ejected from the HZ in even less time, 7.6 yr.
This is interesting due to the fact that in the second case the terrestrial planet is starting much further away from the giant planet than in the first case, thus being exposed to a significantly smaller amount of gravitational attraction, and yet remains within the HZ for a much shorter time than when it starts closer to the giant planet. On the other hand, if the starting angle of the terrestrial planet is selected as 30°, 40° and 50°, the encountered ejection times are 12.1 yr, 14.2 yr and 2.3 yr, respectively, for the giant planet placed at perigee position. However, if the giant planet is originally placed at apogee position, the ejection times of the Earth-mass planets are 10.9 yr, 141.1 yr and 14.1 yr, respectively. This is a highly intriguing example that small changes in the initial conditions of a non-dissipative system can lead to significantly different outcomes (see also Fig. 3).
The next set of comparisons shown deals with the comparison between initial mass and increased mass as well as variation in the starting position of the giant planet, as depicted in Fig. 6. Each chart shows the orbits of a single simulation in each case. They compare the initial mass and the increased mass for both the perigee and apogee starting positions of the giant planet. In all four cases, the starting position of the Earth-mass planet has been chosen as θ0=10°, a choice which is admittedly highly arbitrary. It can be seen again from these plots that changing the starting position from perigee to apogee of the minimum-mass giant planet increases the amount of time the terrestrial planet is within the HZ by almost a factor of three.
However, this is not the case for the increased mass. While the increased mass does increase the amount of time that the terrestrial planet resides within the HZ for the perigee starting position of the giant planet, it reduces the amount of time spent within the HZ when the giant planet starts at apogee. It can also be seen from the figure that the changes in starting position and mass of the giant planet produce extremely different outcomes for the orbit of the terrestrial planet. In three of the four cases it can be seen that the terrestrial planet is ejected from the HZ as well as from the system in a very short amount of time. However, in the case of the giant planet starting at perigee and with an increased mass, it is seen that the terrestrial planet is ejected from the HZ in 34.9 yr but still remains to be part of the system. The orbit of the terrestrial planet appears to take on a more elliptical path with each close encounter with the giant planet, pushing it farther and farther away from the host star and ultimately the stellar HZ.
Figure 7 is relatively similar to Fig. 5 in respect to what is being compared with the exception of the initial placement of the terrestrial planet within the HZ. For this set of simulations we are now looking at a starting placement of δ=0.9 and a starting angle of 160°. It can be seen that by changing the starting position of the giant planet the orbit of the terrestrial planet is again significantly altered. In the case of the giant planet at apogee the terrestrial planet is very close to the starting position of the giant planet as can be seen from the early interaction. Owing to this interaction, the terrestrial planet is almost immediately ejected from the system, whereas for the case of perigee it takes a longer period of time before the interaction between the giant planet and the terrestrial planet causes the terrestrial planet to be ejected from the system. This was not the case for the case studies depicted in Figs 5 and 6. Although the terrestrial planet started much closer to the perigee starting position of the giant planet, it was seen that when the giant planet started at apogee that the terrestrial planet was ejected from the HZ much more quickly, which is again very surprising. In the cases depicted in Fig. 7, the difference between the length of stay within the HZ is about a factor of 20.
The next three figures (Figs 8–10) are all associated with the star system HD 188015. In this system, the giant planet has a semimajor axis of a p=1.203 and a relatively small eccentricity of e p=0.137 (see Table 1), and thus spends all of its time in the stellar HZ. Thus, there is typically a stronger gravitational interaction between the Jupiter-type and the Earth-mass planet compared to the system HD 20782, although strong variations in the gravitational force exist owing to the changes in orbital positions of the objects.

Fig. 8. Case studies of orbits for a terrestrial planet (dotted line) in the HZ (grey area) of HD 188015. The starting positions of the Earth-mass planets are given as δ=0.1 at an angle of θ0=250°. The giant planet (solid line) was initially placed at perigee (top figure) and apogee (bottom figure). The Earth-mass planets were ejected after 13.8 and 8.3 yr, respectively.

Fig. 9. Case studies of orbits for a terrestrial planet (dotted line) in the HZ (grey area) of HD 188015. The starting positions of the Earth-mass planets are given as δ=0.5 at an angle of θ0=160°. The simulations also refer to different cases for the giant planet (solid line), originally placed at perigee (left column) and apogee (right column) and with minimum mass (top row) and its mass increased by 30% (bottom row), respectively. The Earth-mass planets were ejected after 53.8, 11.9, 5.7 and 5.9 yr, respectively (top left, top right, bottom left, bottom right).

Fig. 10. Case studies of orbits for a terrestrial planet (dotted line) in the HZ (grey area) of HD 188015. The starting positions of the Earth-mass planets are given as δ=0.9 at an angle of θ0=40°. The giant planet (solid line) was initially placed at perigee (top figure) and apogee (bottom figure). The Earth-mass planets were ejected after 1.3 and 3.3 yr, respectively.
We again examine the various starting positions of the terrestrial planet ranging from δ=0.1, 0.5 and 0.9 to compare the different outcomes of the simulation runs. The first figure that we examine in greater detail (Fig. 8) is somewhat similar to Fig. 5, albeit with noticeable differences in the star system and the initial starting position of the terrestrial planet. Here we see that for the case of perigee the terrestrial planet's starting position is far enough away from the position of the giant planet that it is able to remain within the system. The interactions with the giant planet do, however, alter the orbit of the terrestrial planet allowing it to fluctuate between the central part of the HZ and being just near the inner border of the HZ, which starts to occur after only 13.8 yr. For the case of apogee, we see that the terrestrial planet and the giant planet start very near to each other causing an interaction to happen almost immediately and ejecting the terrestrial planet from the system, after only 8.3 yr.
Figure 9 is again a comparison between the initial mass and the increased mass of the giant planet with the terrestrial planet originally placed at 160°, and the giant planet placed at either the perigee or the apogee position. For this particular star–planet system we are able to identify a completely different outcome from what we saw with the previous system. The terrestrial planet spends a noticeable amount of time, i.e. 53.8 yr, within the HZ for the perigee starting position of the giant planet with the original mass. However, when the mass is increased, the terrestrial planet is ejected from the HZ after only 5.7 yr. This means that it is now ejected about nine times more quickly than in the simulation with the original mass. For the case of apogee, we see that the terrestrial planet resides within the HZ for 11.9 yr, but is then ejected more quickly when the mass of the giant planet has been increased. In this case the terrestrial planet is ejected twice as fast once the mass of the giant planet is increased. It can also be identified from this set of plots that the Earth-mass planet is spending the most time within the HZ when it is starting near the perigee starting position, which is also what we found for the first two case studies of the previous star–planet system.
The final set of orbits that we show are for that of a starting case of 40° with an initial starting position of δ=0.9 within the HZ. It can be seen in Fig. 10 that the terrestrial planet is ejected from the HZ when the initial position of the giant planet is at perigee, whereas the planet is not necessarily ejected from the system when the giant planet's initial position is at apogee. Nevertheless, in both cases the Earth-mass planets are still almost immediately ejected from the stellar HZ. For the starting position of the giant planet at perigee the terrestrial planet remains within the HZ for a length of 1.3 yr, while for the case of the giant planet starting at apogee the terrestrial planet only remains within the HZ for 3.3 yr. The durations of the planets in the HZ are extremely short but nonetheless correspond to drastically different orbits for the Earth-mass planet.
Conclusions
The goal of our study was to provide a detailed statistical analysis of the ejection of Earth-mass planets from the HZs of the solar twins HD 20782 and HD 188015 and to assess various statistical properties associated with the ejection. The systems of HD 20782 and HD 188015 each possess a giant planet that crosses into the stellar HZ, thus disallowing the existence of habitable terrestrial planets due to the initiation of orbital instabilities; see Noble et al. (Reference Noble, Musielak and Cuntz2002) for previous results for the system HD 210277, which shows a similar type of behaviour. Note that in the case of HD 188015, the orbit of the giant planet is essentially circular, i.e. e p=0.137, whereas in the case of HD 20782, the planetary orbit is extremely elliptical with e p=0.925. As starting positions for the giant planets, we considered both the apogee and perigee positions, whereas the starting positions of the Earth-mass planets were widely varied.
We identified the following statistical properties pertinent to the ejection of the Earth-mass planets from the stellar HZs.
(1) The ejection times are largest if the starting position of the terrestrial planet is placed at δ=0.5, i.e. in the middle of the HZ, instead of δ=0.1 or 0.9, as was expected. The difference is typically a factor of 2 for HD 20782 and a factor of 10 for HD 188015. This notable difference in the factor is caused by the different orbit-dependent strength of the gravitational interaction between the EGP and the Earth-mass planet in the two systems due to the different eccentricities.
(2) Concerning δ=0.5, the most interesting case, the ejection times (as gauged by T Ejmean) are generally longer for HD 188015 than for HD 20782. The reason is that in the latter star–planet system the extremely elliptical orbit of the EGP leads to a larger number of relatively close encounters with the Earth-mass planet, resulting in its ejection from the HZ.
(3) The ejection times for the Earth-mass planets from the HZs of HD 20782 and HD 188015, originally placed at the centre of the respective HZ, vary by a factor of ~200 and ~1500, respectively, depending on the starting positions of the giant and terrestrial planets. It was also found that if the mass of the EGP is increased by 30%, the variation in ejection times for HD 188015 increases by a factor of ~6000.
(4) Note that the reason for this latter behaviour is two-fold. First, in some cases, especially if the EGP is outside the stellar HZ, the higher mass of the EGP results in a faster ejection of the Earth-mass planet from the HZ owing to the increased gravitational attraction between the two objects. Secondly, however, in other cases the Earth-mass planet appears to be bound more closely to the star–planet system if the increased gravitational attraction between the EGP and the Earth-mass is able to enforce such a behaviour due to their orbital positions.
(5) We also studied the manifestation of T Ejmean based on different increments for the starting positions of the Earth-mass planet, which are 90° (four models) up to 10° (36 models). It was found that for both stars, for most sets of models, there were no discernable differences with regard to the selected increment. Surely, assessments based on 90° model increments will miss out on almost all details of the statistical analysis for the planetary ejection times including the skewness analysis (see below).
(6) The analysis of skewness of the ejection times T Ej (but not log T Ej) reveals a significantly positive skew for the majority of models. The skewness is considerably positive in nine out of 16 sets of models, and considerably negative in only one set of models. This result is also evident if T Ejmean is compared to T Ejmedian for the individual models. Thus there is a considerable number of models with relatively large ejection times for the Earth-mass planet compared to the respective means, although all models are orbitally unstable as discussed.
It should be noted, however, that there is no way an ‘Earth’ could have formed in the classical HZ of either of the two systems studied. The short survival times of an ‘Earth’ forced into the HZ is thus no surprise considering the giant planets' proximity. On the other hand, it is noteworthy that considerable differences in the survival times of the Earth-mass planets exist, which may be relevant for establishing guidelines of stability for systems with less intrusive giant planets. In this regard, our results might entail interesting astrobiological properties for systems, which at first sight do not appear to allow for stable terrestrial planets. The point is that due to the skewness of the distribution of ejection times (T Ej), also attained in the orbital stability study of Fatuzzo et al. (Reference Fatuzzo, Adams, Gauvin and Proszkow2006) for binary systems, there is the possibility of ‘pockets of relative stability’, often related to resonances, in systems that are apparently unstable.
In these systems the realization of stability may then depend on the starting condition of the terrestrial planet. In physical systems rather than synthetic systems, this condition may possibly be associated with its formation history. This prospect should serve as a stark motivation to continue exploring the orbital stability of planets in multiple stellar systems including those which appear to be at the brink of instability. An interesting recent example in this respect may be the system ν Octantis (Ramm et al. Reference Ramm, Pourbaix, Hearnshaw and Komonjinda2009), which consists of two stellar components with M 1≃1.4M ⊙ and M 2≃0.5M ⊙ on a noticeably elliptical orbit (e=0.2358). The radial velocity technique seems to indicate the presence of a planet about midway between the stars supposedly in resonance with the binary system that, however, appears highly unlikely based on current planet formation and orbit stability expectations.
Acknowledgements
This work has been supported by the Department of Physics, University of Texas at Arlington. We also would like to thank an anonymous referee whose comments added to the quality of the paper.
Appendix. Summary of statistical measures
In the following we summarize the statistical measures relevant to our study for the sake of completeness and accessibility; see e.g. Upton & Cook (Reference Upton and Cook2008) for additional information. The statistical measures are discussed in detail in most standard textbooks of statistics. The computation of the mean of a distribution or data set is based on the well-established equation

We also make use of the median. The median is described as the number separating the higher half of a distribution or data set from the lower half. If there is an odd number of data, the median can be found by sorting and picking the middle one. If there is an even number of data (as in our studies), the median is usually not unique, and it is customary to take the mean of the two middle values.
In addition, we utilize the definition of variance and standard deviation. The variance of a distribution or data set is a measure of its statistical dispersion. It is defined as

with N*=N or N*=N−1 depending on the adopted statistical model. The standard deviation (SD) is the square root of the variance.
There are two other statistical measures of importance, which are skewness and excess kurtosis. The skewness is defined as

and is a measure of symmetry. A distribution or data set that is symmetric about the centre has a skewness of zero as, for example, the normal distribution. A negative value of Skew indicates that there is a noticeable (i.e. assymetric) tail to the left, whereas the bulk of the distribution is positioned to the right. Conversely, a positive value of Skew indicates an assymetric tail to the right with the bulk of the distribution positioned to the left.
The excess kurtosis is defined as

and indicates whether a distribution or data set has a highly acute peak and has thick tails to the left and right of the peak (positive excess kurtosis) or whether it has a broad, i.e. lower and wider peak with thin tails to the left and right (negative excess kurtosis). Note that the excess kurtosis of the normal distribution is zero.