1 - Astronomy Group, The Open University, Milton Keynes, MK7 6AA, UK
2 - Armagh Observatory, College Hill, Armagh, BT61 9DG, UK

Abstract
We have investigated whether terrestrial planets can exist in orbits in known exoplanetary systems such that life could have emerged on those planets. Four contrasting systems have been examined in which giant planets have been detected. Mixed-variable symplectic numerical integration has been used to investigate the orbits of putative terrestrial planets within the habitable zone of each system (the range of distances from the star within which water at the surface of a terrestrial planet would be in the liquid phase). We have shown that Rho CrB and 47 UMa could have terrestrial planets in orbits that remain confined to their habitable zones for biologically significant lengths of time. We have also shown that the Gliese 876 and Ups And systems are very unlikely to have such orbits.

The question we have addressed is: "can terrestrial planets exist in orbits in known exoplanetary systems such that life could have emerged on those planets?''

Over 50 exoplanetary systems (exosystems) are now known. Main sequence stars of about a solar mass predominate because these have been given preference in searches. The semimajor axes of the planets' orbits do not exceed 3.3 AU, and in most cases are less than 1 AU. The planets have been discovered by observing cyclic Doppler shifts in the stellar spectral lines. This gives $m\sin(i_{\rm o})$ , where m is the mass of the planet and $i_{\rm o}$ is the inclination of the orbital plane of the planet with respect to the plane of the sky (Mayor & Queloz 1995). Values of $m\sin(i_{\rm o})$ range from 0.16 $m_{\rm J}$ (HD 83443b) to 11 $m_{\rm J}$ (HD 114762b), where $m_{\rm J}$ is the mass of Jupiter. There are low mass companions that exceed 13 $m_{\rm J}$ , but the likelihood that some of these objects are brown dwarfs rather than planets is rather high (Burrows et al. 1997). Recently HD 209458b has been observed in transit (Charbonneau et al. 2000), leading to a mass estimate of 0.69 $m_{\rm J}$ . The absence of masses less than 0.16 $m_{\rm J}$ and the preponderance of small semimajor axes are presumably observational selection effects - large planets close to a star give a large Doppler amplitude and a short Doppler cycle.

It is believed (Lissauer 1987; Pollack et al. 1996) that all the known exoplanets are rich in hydrogen and helium, and therefore resemble the giant planet Jupiter rather than a class unknown in our Solar System, namely, the super-massive terrestrial planet (the Earth has a mass of only $3.15 \; 10^{-3}~ m_{\rm J})$ . This belief is supported by the radius of HD 209458b, about 1.5 $R_{\rm J}$ , where $R_{\rm J}$ is the radius of Jupiter. Therefore, the known exoplanets probably formed not where they are now, but further from the star, and then migrated inwards (Boss 1995). Various migration schemes have been proposed. It is expected that migration would have left the zone traversed by the giant, and much of the space interior to this zone, devoid of terrestrial bodies. However, given that migration is caused either by planetesimal scattering by the giant (Murray et al. 1998) or by the interaction of the giant with the nebular disc (Goldreich & Tremaine 1980; Ward 1996; Lin et al. 1998), it is possible that there was enough material left over after giant migration to create terrestrial planets.

The detection of terrestrial planets is beyond present capabilities, and so one can only theorise about their existence. We assume that terrestrial planets could be present in at least some of the exosystems, and we have concentrated on four contrasting systems - Gliese 876, Ups And, and particularly Rho CrB and 47 UMa. Whether life could have emerged on terrestrial planets in any of these systems depends on whether the terrestrial orbits remain confined long enough to the habitable zone of each system.

2 The habitable zone

All life on Earth requires liquid water during at least part of its life cycle. Consequently, it is usual to define the habitable zone (HZ) as the range of distances from a star within which any water at the surface of a terrestrial planet would be in the liquid phase (Kasting et al. 1993). A variety of criteria have been used to define the boundaries of the HZ. For the inner boundary we use the maximum distance from the star where a runaway greenhouse effect occurs leading to the evaporation of all surface water. For the outer boundary we use the maximum distance at which a cloud-free ${\rm CO}_{2}$ atmosphere can maintain a surface temperature of 273 K. Kasting et al. (1993) have used these criteria in conjunction with climate models to obtain values for the boundary distances for various stars, and we have used these values. The values are conservative because of simplifying features in the models. Notably, at the inner boundary H₂O cloud formation is neglected, which would cause this boundary to move inwards. At the outer boundary ${\rm CO}_{2}$ cloud formation is neglected. Forget & Pierrehumbert (1997) have shown that the net effect of the formation of ${\rm CO}_{2}$ clouds is to warm the surface, through a scattering greenhouse effect from the small particles that make up the cloud, and this moves the outer boundary of the HZ outwards.

For zero-age main-sequence stars (ZAMS stars) the boundaries of the HZ lie closer to the star the later its spectral type. This is because of the combined effects of the star's lower luminosity and the shift in its spectrum to longer wavelengths. As a star ages its luminosity increases and the inner boundary moves outwards. The outer boundary also moves outwards unless "cold starts'' are prohibited, in which case it remains fixed at the ZAMS value. The estimated ages of Rho CrB, Ups And, and 47 UMa are 6 Gyr, 3.3 Gyr, and 7 Gyr respectively (Gonzalez 1998; Gonzalez & Laws 2000), and so the outward movement of the boundaries will have been significant. For Gliese 876 the main-sequence lifetime is so long (of the order of 10³ Gyr) that the boundary movement has been negligible.

3 The method of orbital investigation

To establish whether terrestrial orbits could remain confined to the HZ of a system we must investigate the stability of orbits launched in the HZ. Analytical integration is restricted to three bodies in (near) coplanar, (near) circular orbits. By contrast, the giants are in eccentric orbits, the terrestrial bodies would acquire eccentric orbits, and we cannot be restricted to three bodies or to coplanar orbits. We therefore use numerical integration. In particular we use the second-order mixed-variable symplectic (MVS) integrator contained within the Mercury integrator package (Chambers 1999). This integrator originated with Levison & Duncan (1994). It has since been extensively tested by one of us (JEC) and by others. The version in Mercury is accurate in integrating the orbits of planetary bodies. MVS integrators are about ten times faster than other integrators when, as in exosystems, one body (the star) is the dominant gravitational influence. The symplectic property is that there is no build-up of errors in the total energy and total angular momentum of the system.

MVS integrators cannot handle close encounters between planets accurately, because the part of the Hamiltonian that describes the planetary interaction is then comparable with the star-planet Hamiltonian. The smallest distance at which it is safe to use the integrator is about three times the Hill radius of the planet in the encounter with the larger Hill radius. The Hill radius $R_{\rm H}$ is defined as

For there to be a stable orbit we require that the > $3 R_{\rm H}$ criterion be met, but this alone is not sufficient. We also require that the orbit remains confined to the HZ, otherwise it is unlikely that life could evolve. We take confinement to mean that the semimajor axis remains in the HZ at all times in an integration that is not halted by a close encounter. An even tighter criterion would additionally restrict the orbits in the HZ to some upper limit of eccentricity, but we have not adopted such a criterion.

**Table 1:** The exosystems selected for study
star	planet(s)	comment
name	type	mass/ $M_{\odot}$	m sin $(i_{\rm o})/m_{\rm J}$	$a/{\rm AU}$	e
Gliese 876	M4V	0.336	$2.1\pm0.2$	$0.21\pm0.1$	$0.27\pm0.03$	low mass star
Rho CrB	G0Va	1.00	1.1	0.23	$0.028\pm0.04$	small a and e
Ups And	F8V	1.3	0.71	0.059	$0.034\pm0.15$	a family of giants, one
			2.11	0.83	$0.18\pm0.11$	with high e
			4.61	2.50	$0.41\pm0.11$
47 UMa	G1V	1.03	2.41	2.10	$0.096\pm0.03$	most like Solar System

To perform an integration we set up one or two terrestrial planets in orbits with zero eccentricity and in most cases with zero inclination with respect to the plane of the giant's orbit. At least one of the terrestrial planets is placed in the HZ. If there is just one terrestrial planet then it is called EM and in most cases its mass is equal to that of the Earth plus the Moon; in a few cases its mass is greater. If there is a second terrestrial planet it is called V, is given in most cases the mass of Venus (in a few cases it is given greater mass), and it is placed interior to EM. The mass of the giant planet is set at multiples of its minimum mass, from $1\times$ , which corresponds to the orbit viewed edge-on i.e. $i_{\rm o} = 90^\circ$ , up to $8\times$ , which corresponds to the almost face-on $i_{\rm o} = 7.2^\circ$ , or to a smaller multiple if this is necessary to keep the giant from exceeding the 13 $m_{\rm J}$ threshold for becoming a brown dwarf. In the initial-value input file the argument of the pericentre of the giant is set to $180^\circ$ , and that of the terrestrial planets to zero (except for Gliese 876 where these values were reversed). For all planets the longitude of the ascending node and the mean anomaly are each set to zero. (Note that with zero eccentricity or zero inclination some of these angles are undefined, but the input file requires some value, and we chose zero.)

We then integrate for simulated times between 10⁸ and 10⁹ years, unless the integration is halted automatically by a close encounter $(3 R_{\rm H})$ . Ideally we would have liked to integrate all systems for 10⁹ years, because for the Earth the biosphere was well established at this age (Chyba 1993). However, in order to avoid integration instabilities, the integration time-step needs to be less than one-twelfth of the orbital period of the planet with the shortest period (see below), and so for some systems an integration for 10⁹ years would then consume several thousand hours of CPU time on the Compaq Alpha-based workstations used. Therefore, in order always to explore reasonable ranges of orbital parameters we had to set the integration for less than 10⁹ years in some cases.

4 Results

Details of the four contrasting exosystems selected for study are given in Table 1, where m is the mass of the giant (in Jupiter masses $m_{\rm J}$ ), a is the semimajor axis (sma) of the orbit, and e is its orbital eccentricity. Figure 1 displays each of the systems in Table 1, showing the ZAMS habitable zone HZ(0) (shaded), the boundaries of the habitable zone today HZ(t) (vertical dashed lines), and the giant(s) (black discs). The solid line shows the total excursion $2 \Delta{r}$ of the giant planet due to its eccentricity, and the horizontal dashed line extends to $(3 {R}_{\rm H} + \Delta{r})$ each side of the giant for the case when the giant has its minimum mass. In all the integrations the changes in a and e of the giant were very small, typically of order 0.01% in a and 0.1% in e.

4.1 Rho CrB

Inspection of Fig. 1 shows that Rho CrB could well have stable orbits confined to the HZs. We have performed integrations for simulated times 1-3.5 $\;10^{8}$ years. Table 2 summarises the results. The first column is the line number in the table. The second column shows the mass of the giant in terms of $m_{\rm J}$ , its minimum mass being 1.1 $m_{\rm J}$ . The third column shows the factor by which the mass of Venus and the Earth-Moon was multiplied to get the masses of V and EM respectively in the simulation. The next four columns show the starting values of the sma a and the inclination i of V and EM - the starting eccentricities e are always zero. Columns 1-7 define what we call a configuration. The next four columns show the extreme ranges of a and e seen in the simulation. The final column shows the outcome, where, for example, "EM-G 8.6'' in line 1 means that EM came within $3 R_{\rm H}$ of the giant after 8.6 simulated years, and where, for example in line 3 "> $1 \ 10^{8}$ '' means that there were no close encounters up to the end of a simulation of 10⁸ years.

Rho CrB is a time-consuming system to explore, because the orbital period of the giant is only 39.645 days. The time step in the integration therefore needs to be less than 3.3 days, which leads to more than about 30 hours CPU time per 10⁸ years simulated time on the work stations used. We chose 3.0 days as the time step, but repeated two of the configurations with a time step of 1.5 days, to test whether 3.0 days was short enough. These two configurations are indicated in Table 2 (lines 4 and 20). In each case the ranges of the sma and eccentricity of the terrestrial planet(s) are the same to better than 0.2%. The inclinations started at zero, and therefore remained at zero. Figure 2 shows the details of the two simulations of the configuration in line 4 of Table 2. (Note that the eccentricity variation is probably just the sort of secular variation commonly seen in the Solar System.) As a further test of the 3.0 day time step in this configuration we evaluated the Jacobi constant every couple of hundred years over the first 10⁶ years, and found that changes in it were confined to $(+0.8/-1.6) \; 10^{-4}\%$ , with no discernible secular trend at the $10^{-5}\%$ level.

Table 2 and Fig. 1 show that stable orbits are readily found in the HZs of Rho CrB. Consider the case when only EM and the giant are present, with the giant at eight times its minimum mass (lines 1-5). In this case $(3 {R}_{\rm H} + \Delta{r})$ extends to 0.33 AU. The inner edge of HZ(0) is just inside the mean motion resonance with a giant: EM orbital period ratio of 1:6, corresponding to $a_{\rm EM} = 0.7594$ AU, yet when EM is launched here the orbit is stable for at least 10⁸ years. The orbit also remains confined to HZ(0). When the start sma is 1.32 AU the orbit is stable for at least $3 \; 10^{8}$ years, and only just strays outside HZ(0). Only in lines 1 and 2 do we see instability. These orbits start well inside the inner boundary of HZ(0), and at these rather small distances from the star instability is expected. Note that in line 1 the final evaluations of $a_{\rm EM}$ and $e_{\rm EM}$ are not included in the ranges shown, because at the encounter the orbit was then hyperbolic. There is probably little physical significance to such an orbit because at such an encounter the integrator is then working outside its range of applicability.

Table 2: Orbital integrations for Rho CrB
giant V, start start start start range range range range outcome/

mass EM $a_{\rm V}$ $i_{\rm V}/^\circ$ a $_{\rm EM}$ $i_{\rm EM}/^\circ$ $a_{\rm V}$ /AU $e_{\rm V}$ $a_{\rm EM}$ /AU $e_{\rm EM}$ years

$m_{\rm J}$ mass /AU /AU

factor

1 8.8 1 no V no V 0.35   0 no V no V 0.350-0.386 0-0.285 EM-G 8.6

2 0.3651 0.365-3.34 0-0.931 EM-G 7.0

3 0.7594 0.759-0.820 0-0.089 $>1\; 10^{8}$

4(a) 1.32 1.32-1.45 0-0.098 $>3 \; 10^{8}$

5 8 1.00 1.00-1.09 0-0.092 $>1\; 10^{8}$

6 1.1 1 0.723   0 1.00   0 0.723-0.730 0-0.030 1.00-1.01 0-0.026 $>1\; 10^{8}$

7 10 10 0.723-0.730 0-0.033 1.00-1.01 0-0.025 $>1\; 10^{8}$

8 20 20 0.723-0.729 0-0.029 1.00-1.01 0-0.024 $>1\; 10^{8}$

9 0.35   0   0 0.350-0.354 0-0.066 1.00-1.01 0-0.030 $>1\; 10^{8}$

10 2.2 0.723 0.723-0.736 0-0.038 1.00-1.02 0-0.035 $>1\; 10^{8}$

11 4.4 0.723-0.750 0-0.054 1.00-1.04 0-0.054 $>1\; 10^{8}$

12(b) 8.8 0.77 0.723-0.780 0-0.086 0.770-0.833 0-0.093 $>3.5\;10^{8}$

13(b) 0.80 0.723-0.780 0-0.086 0.800-0.865 0-0.092 $>1\; 10^{8}$

14 0.83 0.723-0.780 0-0.103 0.830-0.90 0-0.093 $>2\;10^{8}$

15 0.85 0.723-0.780 0-0.087 0.850-0.92 0-0.093 $>1\; 10^{8}$

16 1.00 0.723-0.780 0-0.088 1.00-1.088 0-0.093 $>3 \; 10^{8}$

17 1.00 1.40 1.00-1.09 0-0.092 1.40-1.54 0-0.101 $>1\; 10^{8}$

18 1.20 1.20-1.31 0-0.105 1.40-1.54 0-0.108 $>1\; 10^{8}$

19 1.25 1.25-1.38 0-0.132 1.40-1.59 0-0.121 EM-V $8.2 \; 10^{6}$

20 (a) 1.35 1.34-1.48 0-0.097 1.40-1.54 0-0.086 EM-V $1.3 \; 10^{2}$

21 2 0.723 1.00 0.723-0.780 0-0.087 1.00-1.09 0-0.093 $>1\; 10^{8}$

22 4 0.723-0.780 0-0.088 1.00-1.09 0-0.093 $>1\; 10^{8}$

23 8 0.723-0.780 0-0.087 1.00-1.09 0-0.095 $>2\;10^{8}$

24 10 10 0.723-0.780 0-0.096 1.00-1.09 0-0.093 $>3 \; 10^{8}$

25 1.00   0 1.524   0 1.00-1.09 0-0.093 1.52-1.69 0-0.105 $>1\; 10^{8}$

26 1.30 1.40 1.30-1.43 0-0.098 1.40-1.54 0-0.102 $>1\; 10^{8}$

(a) Same outcome with a time step of 1.5 days and 3.0 days.

**Table 2:** Orbital integrations for Rho CrB
	giant	V,	start	start	start	start	range	range	range	range	outcome/
	mass	EM	$a_{\rm V}$	$i_{\rm V}/^\circ$	a $_{\rm EM}$	$i_{\rm EM}/^\circ$	$a_{\rm V}$ /AU	$e_{\rm V}$	$a_{\rm EM}$ /AU	$e_{\rm EM}$	years
	$m_{\rm J}$	mass	/AU		/AU
		factor
1	8.8	1	no V	no V	0.35	0	no V	no V	0.350-0.386	0-0.285	EM-G 8.6
2					0.3651				0.365-3.34	0-0.931	EM-G 7.0
3					0.7594				0.759-0.820	0-0.089	$>1\; 10^{8}$
4(a)					1.32				1.32-1.45	0-0.098	$>3 \; 10^{8}$
5		8			1.00				1.00-1.09	0-0.092	$>1\; 10^{8}$

6	1.1	1	0.723	0	1.00	0	0.723-0.730	0-0.030	1.00-1.01	0-0.026	$>1\; 10^{8}$
7				10		10	0.723-0.730	0-0.033	1.00-1.01	0-0.025	$>1\; 10^{8}$
8				20		20	0.723-0.729	0-0.029	1.00-1.01	0-0.024	$>1\; 10^{8}$
9			0.35	0		0	0.350-0.354	0-0.066	1.00-1.01	0-0.030	$>1\; 10^{8}$
10	2.2		0.723				0.723-0.736	0-0.038	1.00-1.02	0-0.035	$>1\; 10^{8}$
11	4.4						0.723-0.750	0-0.054	1.00-1.04	0-0.054	$>1\; 10^{8}$

12(b)	8.8				0.77		0.723-0.780	0-0.086	0.770-0.833	0-0.093	$>3.5\;10^{8}$
13(b)					0.80		0.723-0.780	0-0.086	0.800-0.865	0-0.092	$>1\; 10^{8}$
14					0.83		0.723-0.780	0-0.103	0.830-0.90	0-0.093	$>2\;10^{8}$
15					0.85		0.723-0.780	0-0.087	0.850-0.92	0-0.093	$>1\; 10^{8}$
16					1.00		0.723-0.780	0-0.088	1.00-1.088	0-0.093	$>3 \; 10^{8}$
17			1.00		1.40		1.00-1.09	0-0.092	1.40-1.54	0-0.101	$>1\; 10^{8}$
18			1.20				1.20-1.31	0-0.105	1.40-1.54	0-0.108	$>1\; 10^{8}$
19			1.25				1.25-1.38	0-0.132	1.40-1.59	0-0.121	EM-V $8.2 \; 10^{6}$
20 (a)			1.35				1.34-1.48	0-0.097	1.40-1.54	0-0.086	EM-V $1.3 \; 10^{2}$
21		2	0.723		1.00		0.723-0.780	0-0.087	1.00-1.09	0-0.093	$>1\; 10^{8}$
22		4					0.723-0.780	0-0.088	1.00-1.09	0-0.093	$>1\; 10^{8}$
23		8					0.723-0.780	0-0.087	1.00-1.09	0-0.095	$>2\;10^{8}$
24				10		10	0.723-0.780	0-0.096	1.00-1.09	0-0.093	$>3 \; 10^{8}$
25			1.00	0	1.524	0	1.00-1.09	0-0.093	1.52-1.69	0-0.105	$>1\; 10^{8}$
26			1.30		1.40		1.30-1.43	0-0.098	1.40-1.54	0-0.102	$>1\; 10^{8}$

(b) Time step 1.5 days instead of 3.0 days.

In the rest of Table 2 the terrestrial planet V is also present. This provides a more stringent condition because close encounters between EM and V are now possible. Lines 6-8 have the giant with minimum mass and V and EM launched in their Solar System positions. Stability and confinement are found in all three cases for at least 10⁸ years. The configurations differ in the start inclinations, but this makes very little difference to the outcome; this is also apparent in lines 23 and 24. Lines 6, 10, 11, and 16 show the effect of increasing the giant's mass by factors of two to 8 times the minimum. As expected, there is an increase in the ranges of a and e, but the ranges remain small.

In lines 12-26 the giant has 8 times its minimum mass. Near the inner edge of HZ(0), V and EM can be put surprisingly close together and yet avoid close encounters for at least $3.5 \; 10^{8}$ years (line 12); to be sure of this a time step of 1.5 days was used. Study of the details of the orbits reveals that when the eccentricities are large the longitudes of the pericentres of V and EM differ by less than about $30^{\rm o}$ . This is an example of e-omega coupling, which in the Solar System protects some asteroids on planet-crossing orbits. This deserves further investigation, but with respect to the present work we need only note that it aids stability and confinement in the HZs. With EM near the outer edge of HZ(0), V has close encounters with EM if it is launched beyond 1.2 AU (lines 19, 20).

Lines 21-23 show that when the masses of V and EM are increased this has very little effect on the outcome. The dominant effect on a and e of V and EM is the giant.

HZ(t) lies further from the giant than does HZ(0) (Fig. 1). It is therefore to be expected that stability and confinement will be even more prevalent in HZ(t) than it is in HZ(0).

4.2 47 UMa

We have investigated 47 UMa more extensively than the other three systems. This is a solar-type star with the giant outside the HZ, and so of the four systems studied it most resembles the Solar System. The orbital periods of interest are sufficiently long for us to have been able to adopt 10⁹ years as the standard integration time, with a time step 1/18 of the orbital period of the body with the shortest period. Figure 3 shows the outcome with just EM present, and with the giant at 1, 1.5, 2, and 4 times its minimum mass, corresponding to $14.5^\circ \leq i_{\rm o} \leq 90^\circ$ .

At the minimum giant mass $(2.41~ m_{\rm J})$ stable orbits are found inward of 1.3 AU. Just beyond this distance the mean-motion resonance is encountered with a giant: EM orbital period ratio of 2:1, corresponding to ${a}_{\rm EM} = 1.32292$ AU. Table 3 shows that this results in EM-G close encounters on a time scale of 10⁸ years. At 1.4 AU the close encounter is after only a few years, and is presumably because $(3 {R}_{\rm H} + \Delta{r})$ extends inwards of 1.4 AU (Fig. 3). Another mean motion resonance is encountered at 1.00957 AU - this is the giant: EM orbital period ratio of 3:1, and it results in a narrow zone of instability. Figure 4 shows the end of the orbital evolution of EM when it was launched at this resonance; it is the rather sudden increase in $e_{\rm EM}$ that leads to the encounter.

Table 3: Orbital integrations for 47UMa with EM only
giant start start range $a_{\rm EM}/{\rm AU}$ range $e_{\rm EM}$ outcome/

mass/ $m_{\rm J}$ $a_{\rm EM}$ $i_{\rm EM}/^\circ$ years

/AU

2.41 1.00   0 $1.000\pm0.002$ 0-0.032 $>1\;10^{9}$

1.00957 0.855-1.155 0-0.659 EM-G $9.6\; 10^{4}$

1.02 1.015-1.023 0-0.181 $>1\;10^{9}$

1.20 1.196-1.205 0-0.134 $>1\;10^{9}$

1.30 1.289-1.305 0-0.087 $>1\;10^{9}$

10 1.288-1.305 0-0.081 $>1\;10^{9}$

1.32292   0 1.274-1.367 0-0.315 EM-G $3.0\;10^{8}$

1.35 1.273-1.373 0-0.370 EM-G $4.5\;10^{7}$

10 1.272-1.370 0-0.382 EM-G $3.3\;10^{5}$

1.40   0 1.400-1.411 0-0.053 EM-G 8.9

3.615 0.70   0 $0.7000\pm0.0004$ 0-0.095 $>1\;10^{9}$

0.83339 0.829-0.836 0-0.085 $>1\;10^{9}$

1.00 0.996-1.004 0-0.049 $>1\;10^{9}$

1.00957 0.986-1.021 0-0.603 EM-G $1.30\;10^{4}$

1.20 1.194-1.208 0-0.132 $>1\;10^{9}$

1.30 1.281-1.310 0-0.110 $>1\;10^{9}$

1.32292 1.281-1.347 0-0.183 EM-G $8.0\;10^{2}$

1.35 1.297-1.357 0-0.236 EM-G $3.2\;10^{4}$

1.36 1.342-1.360 0-0.154 EM-G 39

1.40 1.400-1.409 0-0.091 EM-G 8.8

1.6026 ... ... EM-G 3.0

4.82 0.83339   0 0.828-0.835 0-0.113 $>1\;10^{9}$

1.00 0.994-1.006 0-0.071 $>1\;10^{9}$

1.00957 0.994-1.040 0-0.623 EM-G $7.3\;10^{3}$

1.02 0.955-1.156 0-0.722 EM-G $1.6\;10^{3}$

1.03 1.030-1.034 0-0.189 $>1\;10^{9}$

1.10 1.095-1.106 0-0.144 $>1\;10^{9}$

1.15 1.146-1.159 0-0.140 $>1\;10^{9}$

10 1.146-1.159 0-0.140 $>1\;10^{9}$

1.20   0 1.192-1.209 0-0.132 EM-G $6.1\;10^{2}$

10 1.191-1.208 0-0.135 EM-G $2.8\;10^{3}$

1.40   0 1.400-1.406 0-0.064 EM-G 5.3

9.64 0.70   0 0.699-0.701 0-0.140 $>1\;10^{9}$

0.80 0.798-0.802 0-0.091 $>1\;10^{9}$

0.83339 0.83339-0.852 0-0.540 EM-G $3.9\;10^{8}$

0.90 0.897-0.903 0-0.086 $>1\;10^{9}$

0.95 0.946-0.954 0-0.058 $>1\;10^{9}$

1.00 1.00-1.01 0-0.055 EM-G 95

1.03 1.030-1.034 0-0.044 EM-G 12

1.04 1.040-1.046 0-0.025 EM-G 8.8

1.05 1.050-1.057 0-0.009 EM-G 5.7

1.10 1.10-1.109 0-0.014 EM-G 2.7

**Table 3:** Orbital integrations for 47UMa with EM only
giant	start	start	range $a_{\rm EM}/{\rm AU}$	range $e_{\rm EM}$	outcome/
mass/ $m_{\rm J}$	$a_{\rm EM}$	$i_{\rm EM}/^\circ$			years
	/AU
2.41	1.00	0	$1.000\pm0.002$	0-0.032	$>1\;10^{9}$
	1.00957		0.855-1.155	0-0.659	EM-G $9.6\; 10^{4}$
	1.02		1.015-1.023	0-0.181	$>1\;10^{9}$
	1.20		1.196-1.205	0-0.134	$>1\;10^{9}$
	1.30		1.289-1.305	0-0.087	$>1\;10^{9}$
		10	1.288-1.305	0-0.081	$>1\;10^{9}$
	1.32292	0	1.274-1.367	0-0.315	EM-G $3.0\;10^{8}$
	1.35		1.273-1.373	0-0.370	EM-G $4.5\;10^{7}$
		10	1.272-1.370	0-0.382	EM-G $3.3\;10^{5}$
	1.40	0	1.400-1.411	0-0.053	EM-G 8.9
3.615	0.70	0	$0.7000\pm0.0004$	0-0.095	$>1\;10^{9}$
	0.83339		0.829-0.836	0-0.085	$>1\;10^{9}$
	1.00		0.996-1.004	0-0.049	$>1\;10^{9}$
	1.00957		0.986-1.021	0-0.603	EM-G $1.30\;10^{4}$
	1.20		1.194-1.208	0-0.132	$>1\;10^{9}$
	1.30		1.281-1.310	0-0.110	$>1\;10^{9}$
	1.32292		1.281-1.347	0-0.183	EM-G $8.0\;10^{2}$
	1.35		1.297-1.357	0-0.236	EM-G $3.2\;10^{4}$
	1.36		1.342-1.360	0-0.154	EM-G 39
	1.40		1.400-1.409	0-0.091	EM-G 8.8
	1.6026		...	...	EM-G 3.0
4.82	0.83339	0	0.828-0.835	0-0.113	$>1\;10^{9}$
	1.00		0.994-1.006	0-0.071	$>1\;10^{9}$
	1.00957		0.994-1.040	0-0.623	EM-G $7.3\;10^{3}$
	1.02		0.955-1.156	0-0.722	EM-G $1.6\;10^{3}$
	1.03		1.030-1.034	0-0.189	$>1\;10^{9}$
	1.10		1.095-1.106	0-0.144	$>1\;10^{9}$
	1.15		1.146-1.159	0-0.140	$>1\;10^{9}$
		10	1.146-1.159	0-0.140	$>1\;10^{9}$
	1.20	0	1.192-1.209	0-0.132	EM-G $6.1\;10^{2}$
		10	1.191-1.208	0-0.135	EM-G $2.8\;10^{3}$
	1.40	0	1.400-1.406	0-0.064	EM-G 5.3
9.64	0.70	0	0.699-0.701	0-0.140	$>1\;10^{9}$
	0.80		0.798-0.802	0-0.091	$>1\;10^{9}$
	0.83339		0.83339-0.852	0-0.540	EM-G $3.9\;10^{8}$
	0.90		0.897-0.903	0-0.086	$>1\;10^{9}$
	0.95		0.946-0.954	0-0.058	$>1\;10^{9}$
	1.00		1.00-1.01	0-0.055	EM-G 95
	1.03		1.030-1.034	0-0.044	EM-G 12
	1.04		1.040-1.046	0-0.025	EM-G 8.8
	1.05		1.050-1.057	0-0.009	EM-G 5.7
	1.10		1.10-1.109	0-0.014	EM-G 2.7

At 1.5 times the minimum giant mass $(3.62~ m_{\rm J})$ the outcome is much the same as at the minimum mass, with somewhat larger ranges of $a_{\rm EM}$ and $e_{\rm EM}$ as expected. The 4:1 mean motion resonance at 0.83339 AU was now included and EM is stable there (Table 3). This orbit has not yet been investigated for libration. At 1.30 AU, the orbit is stable for 10⁹ years, and Fig. 3 shows that it is within $(3 \it {R}_{\rm H} + \Delta {r})$ of the giant. However, with $\Delta r = 0.20$ AU, it is not within 3 $R_{\rm H}$ of the giant, and it is possible for a small proportion of terrestrial orbits, by chance, to survive for long periods, perhaps aided by e-omega coupling. In summary, taking the results at 2.41 $m_{\rm J}$ and 3.62 $m_{\rm J}$ together, Fig. 3 shows that we have found stable orbits across most of HZ(0), and in the inner region of HZ $(\it {t})$ .

At twice the minimum giant mass $(4.82~{m}_{\rm J})$ the increasing reach of $(3 {R}_{\rm H} + \Delta{r})$ has pushed the stable zone to within 1.2 AU of the star. At four times the minimum giant mass $(9.64~ {m}_{\rm J}, {i}_{\rm o} = 14.5^\circ)$ stable orbits were found only inward of 1.0 AU, and the 4:1 mean motion resonance at 0.83339 AU is now unstable. At this high mass we have found stable orbits only in the inner region of HZ(0), and nowhere in HZ(t).

From Table 3 it can be seen that the range $\Delta{a}$ in sma is generally smaller than the excursion 2ae associated with the maximum eccentricity. This is also apparent in lines 3-5 of Table 2. This means that the larger effect of the giant is to pump up the orbital eccentricity of EM. This is in accord with analytical studies of systems of two planets of low mass in near-circular, low-inclination orbits (Hénon & Petit 1986).

Figure 3 shows that, except near mean motion resonances, the orbit of EM is stable provided its start sma is (nearly) outside the distance $(3 {R}_{\rm H} + \Delta{r})$ from the giant. This is a reasonable outcome, but it is not obvious because the orbit of EM evolves. Numerical integration is needed to see if this expectation is met, and also to explore the changes in the sma of EM and the increases in its eccentricity. Exploration of stability at resonances also requires numerical integration.

When V is added to EM in the 47 UMa system there is further instability due to collisions between EM and V. The results are shown in Fig. 5 for the case of a giant planet with twice the minimum mass $({i}_{\rm o} = 30^\circ)$ . The shaded range of distances along the x-axis corresponds to HZ(0) (Fig. 1). The start sma of EM is plotted along the x-axis and the initial distance between EM and V along the y-axis. Remember that V is always interior to EM. HZ(0) has been divided into two parts by the diagonal dot-dashed line. Below the line only EM is launched in HZ(0), and above the line both EM and V are so launched. The inner boundary of HZ(t) is at 1.0 AU. Open circles show the start sma of stable orbits, crosses show those of orbits terminated by an EM-G close encounter, and asterisks those terminated by an E-V close encounter. Table 4 gives the details for all final outcomes, but there is only a representative sample of ranges in a and e (the ranges for the other runs were not kept because each run overwrote the files of the earlier run).

Table 4: Orbital integrations for 47 UMa with EM and V, and with giant mass $= 4.82\, m_{\rm J}$
start start start start range $a_{\rm v}$ /AU range range range outcome/

$a_{\rm v}$ /AU $i_{\rm v}/^{\circ}$ $a_{\rm EM}$ /AU $i_{\rm EM}/^\circ$ $e_{\rm v}$ $a_{\rm EM}$ /AU $e_{\rm EM}$ years

0.40 0 1.0 0 0.40-0.4004 0-0.005 0.994-1.006 0-0.070 $>1\;10^{9}$

0.41 1.01 ... ... ... ... EM-V $3.5\;10^{3}$

0.42 1.02 ... ... ... ... EM-V $1.1\;10^{3}$

0.43 1.03 0.43-0.4304 0-0.016 1.023-1.034 0-0.186 $>1\;10^{9}$

0.50 1.00 ... ... ... ... $>1\; 10^{8}$

0.51 1.01 ... ... ... ... EM-V $5.8\;10^{3}$

0.52 1.02 ... ... ... ... EM-V $8.1\;10^{2}$

0.53 1.03 ... ... ... ... $>1\; 10^{8}$

0.60 1.00 $0.60\pm0.0006$ 0-0.150 $1.00\pm0.006$ 0-0.070 $>1\;10^{9}$

1.00957 $0.60\pm0.0006$ 0-0.151 0.995-1.018 0-0.481 EM-V $4.7\;10^{4}$

1.20 $0.60\pm0.0004$ 0-0.112 1.192-1.210 0-0.132 EM-G $3.0\;10^{3}$

0.62 1.02 ... ... ... ... EM-V $7.5\; 10^{3}$

0.63 1.03 0.629-0.631 0-0.140 1.023-1.035 0-0.196 $>1\;10^{9}$

0.65 0.95 ... ... ... ... $>1\; 10^{8}$

0.675 0.675-0.676 0-0.105 0.948-0.952 0-0.082 $>1\;10^{9}$

1.00 $0.675\pm0.0006$ 0-0.110 0.994-1.006 0-0.072 $>1\;10^{9}$

0.70 0.95 ... ... ... ... $>1\; 10^{8}$

1.00 0.699-0.7015 0-0.100 0.993-1.006 0-0.074 $>1\;10^{9}$

1.00957 0.699-0.7005 0-0.101 0.994-1.018 0-0.315 EM-V $9.4\;10^{3}$

1.02 ... ... ... ... EM-V $7.1\;10^{2}$

1.03 ... ... ... ... EM-V $5.9\;10^{7}$

1.04 0.697-0.701 0-0.20 1.033-1.047 0-0.181 $>1\;10^{9}$

1.15 $0.70\pm0.0006$ 0-0.101 1.146-1.159 0-0.138 $>1\;10^{9}$

0.723 1.00 0.720-0.723 0-0.105 0.994-1.008 0-0.307 EM-V $1.3\;10^{6}$

10 0.721-0.724 0-0.101 0.994-1.007 0-0.481 EM-V $6.0\;10^{6}$

1.05 0 0.713-0.724 0-0.247 1.035-1.063 0-0.204 $>1\;10^{9}$

1.15 0.723-0.724 0-0.098 1.145-1.158 0-0.139 $>1\;10^{9}$

0.725 0.95 ... ... ... ... $>1\; 10^{8}$

0.74 1.04 ... ... ... ... EM-V $6.3 \; 10^{7}$

0.75 0.95 ... ... ... ... $>1\; 10^{8}$

1.05 0.742-0.769 0-0.220 1.015-1.066 0-0.178 EM-V $5.1\;10^{8}$

1.15 $0.75\pm0.0007$ 0-0.097 1.145-1.159 0-0.138 $>1\;10^{9}$

0.775 0.95 0.774-0.777 0-0.099 0.947-0.953 0-0.097 $>1\;10^{9}$

1.05 ... ... ... ... EM-V $3.4\;10^{7}$

1.10 0.774-0.776 0-0.139 1.095-1.106 0-0.143 $>1\;10^{9}$

1.15 ... ... ... ... $>1\; 10^{8}$

0.80 0.95 ... ... ... ... EM-V $6.8\;10^{7}$

1.05 ... ... ... ... EM-V $2.6\;10^{7}$

1.10 0.798-0.802 0-0.107 1.094-1.107 0-0.148 $>1\;10^{9}$

1.15 0.799-0.801 0-0.096 1.145-1.159 0-0.138 $>1\;10^{9}$

0.825 1.03 0.823-0.827 0-0.092 1.022-1.034 0-0.186 EM-V $2.3\;10^{3}$

1.10 ... ... ... ... EM-V $4.5\;10^{7}$

1.15 0.824-0.829 0-0.087 1.140-1.159 0-0.261 EM-V $1.2\;10^{8}$

0.85 1.00 0.0847-0.851 0-0.077 0.997-1.005 0-0.074 EM-V $6.6\;10^{2}$

1.15 0.846-0.856 0-0.150 1.129-1.163 0-0.244 EM-V $4.8\;10^{8}$

0.875 0.874-0.882 0-0.114 1.123-1.158 0-0.195 EM-V $1.2\;10^{6}$

0.90 1.20 ... ... ... ... EM-G $2.4\;10^{3}$

0.925 ... ... ... ... EM-G $3.0\;10^{3}$

**Table 4:** Orbital integrations for 47 UMa with EM and V, and with giant mass $= 4.82\, m_{\rm J}$
start	start	start	start	range $a_{\rm v}$ /AU	range	range	range	outcome/
$a_{\rm v}$ /AU	$i_{\rm v}/^{\circ}$	$a_{\rm EM}$ /AU	$i_{\rm EM}/^\circ$		$e_{\rm v}$	$a_{\rm EM}$ /AU	$e_{\rm EM}$	years
0.40	0	1.0	0	0.40-0.4004	0-0.005	0.994-1.006	0-0.070	$>1\;10^{9}$
0.41		1.01		...	...	...	...	EM-V $3.5\;10^{3}$
0.42		1.02		...	...	...	...	EM-V $1.1\;10^{3}$
0.43		1.03		0.43-0.4304	0-0.016	1.023-1.034	0-0.186	$>1\;10^{9}$
0.50		1.00		...	...	...	...	$>1\; 10^{8}$
0.51		1.01		...	...	...	...	EM-V $5.8\;10^{3}$
0.52		1.02		...	...	...	...	EM-V $8.1\;10^{2}$
0.53		1.03		...	...	...	...	$>1\; 10^{8}$
0.60		1.00		$0.60\pm0.0006$	0-0.150	$1.00\pm0.006$	0-0.070	$>1\;10^{9}$
		1.00957		$0.60\pm0.0006$	0-0.151	0.995-1.018	0-0.481	EM-V $4.7\;10^{4}$
		1.20		$0.60\pm0.0004$	0-0.112	1.192-1.210	0-0.132	EM-G $3.0\;10^{3}$
0.62		1.02		...	...	...	...	EM-V $7.5\; 10^{3}$
0.63		1.03		0.629-0.631	0-0.140	1.023-1.035	0-0.196	$>1\;10^{9}$
0.65		0.95		...	...	...	...	$>1\; 10^{8}$
0.675				0.675-0.676	0-0.105	0.948-0.952	0-0.082	$>1\;10^{9}$
		1.00		$0.675\pm0.0006$	0-0.110	0.994-1.006	0-0.072	$>1\;10^{9}$
0.70		0.95		...	...	...	...	$>1\; 10^{8}$
		1.00		0.699-0.7015	0-0.100	0.993-1.006	0-0.074	$>1\;10^{9}$
		1.00957		0.699-0.7005	0-0.101	0.994-1.018	0-0.315	EM-V $9.4\;10^{3}$
		1.02		...	...	...	...	EM-V $7.1\;10^{2}$
		1.03		...	...	...	...	EM-V $5.9\;10^{7}$
		1.04		0.697-0.701	0-0.20	1.033-1.047	0-0.181	$>1\;10^{9}$
		1.15		$0.70\pm0.0006$	0-0.101	1.146-1.159	0-0.138	$>1\;10^{9}$

0.723		1.00		0.720-0.723	0-0.105	0.994-1.008	0-0.307	EM-V $1.3\;10^{6}$
			10	0.721-0.724	0-0.101	0.994-1.007	0-0.481	EM-V $6.0\;10^{6}$
		1.05	0	0.713-0.724	0-0.247	1.035-1.063	0-0.204	$>1\;10^{9}$
		1.15		0.723-0.724	0-0.098	1.145-1.158	0-0.139	$>1\;10^{9}$
0.725		0.95		...	...	...	...	$>1\; 10^{8}$
0.74		1.04		...	...	...	...	EM-V $6.3 \; 10^{7}$
0.75		0.95		...	...	...	...	$>1\; 10^{8}$
		1.05		0.742-0.769	0-0.220	1.015-1.066	0-0.178	EM-V $5.1\;10^{8}$
		1.15		$0.75\pm0.0007$	0-0.097	1.145-1.159	0-0.138	$>1\;10^{9}$
0.775		0.95		0.774-0.777	0-0.099	0.947-0.953	0-0.097	$>1\;10^{9}$
		1.05		...	...	...	...	EM-V $3.4\;10^{7}$
		1.10		0.774-0.776	0-0.139	1.095-1.106	0-0.143	$>1\;10^{9}$
		1.15		...	...	...	...	$>1\; 10^{8}$
0.80		0.95		...	...	...	...	EM-V $6.8\;10^{7}$
		1.05		...	...	...	...	EM-V $2.6\;10^{7}$
		1.10		0.798-0.802	0-0.107	1.094-1.107	0-0.148	$>1\;10^{9}$
		1.15		0.799-0.801	0-0.096	1.145-1.159	0-0.138	$>1\;10^{9}$
0.825		1.03		0.823-0.827	0-0.092	1.022-1.034	0-0.186	EM-V $2.3\;10^{3}$
		1.10		...	...	...	...	EM-V $4.5\;10^{7}$
		1.15		0.824-0.829	0-0.087	1.140-1.159	0-0.261	EM-V $1.2\;10^{8}$
0.85		1.00		0.0847-0.851	0-0.077	0.997-1.005	0-0.074	EM-V $6.6\;10^{2}$
		1.15		0.846-0.856	0-0.150	1.129-1.163	0-0.244	EM-V $4.8\;10^{8}$
0.875				0.874-0.882	0-0.114	1.123-1.158	0-0.195	EM-V $1.2\;10^{6}$
0.90		1.20		...	...	...	...	EM-G $2.4\;10^{3}$
0.925				...	...	...	...	EM-G $3.0\;10^{3}$

In Fig. 5 the instability at the giant: EM 3:1 mean motion resonance is clear at ${a}_{\rm EM} = 1.00957$ AU, though when V is present it is due to an EM-V close encounter in all the cases investigated, i.e. $({a}_{\rm EM} - a_{\rm V}) \leq 0.6$ AU. Away from this resonance, the presence of V has made an inroad into the stability at $({a}_{\rm EM} - a_{\rm V}) \leq 0.325$ AU, through close encounters between EM and V.

Tables 3 and 4 show that in the integrations that we have done on 47 UMa, only a few close encounters occurred between 10⁸ and 10⁹ years, and these were either at resonances or near the boundary between stable and unstable zones. Therefore, except near such boundaries, we have seen a tendency for instability in terms of the $3 R_{\rm H}$ criterion to show itself within 10⁸ years.

4.3 Gliese 876 and Ups And

We have also had a brief look at the Gliese 876 and Ups And systems assuming that in the latter the three giants have the same orbital inclination. From Fig. 1 it can be seen that in both of these systems, even if the giant has its minimum mass, the distance $(3 {R}_{\rm H} + \Delta{r})$ from a giant always spans the whole of HZ(0) and HZ(t). It might therefore be expected that there are no stable orbits in the HZs. Numerical integration confirms this expectation, as shown in Table 5, where the time step was 1/12 of the orbital period of the body with the shortest orbital period. For Gliese 876 the $3 R_{\rm H}$ criterion was violated after no more than a year even with the minimum giant mass and with EM launched inside the inner boundary of the HZ. Ranges for the eccentricities and semimajor axes of EM after launch are not given because of the very short time before the close encounter.

Table 5: Orbital integrations for Gliese 876 and Ups And
giant
V, start start start start range range range range outcome/

mass/ $m_{\rm J}$
EM $a_{\rm V}$ $i_{\rm V}/^\circ$ a $_{\rm EM}$ $i_{\rm EM}/^\circ$ $a_{\rm V}$ $i_{\rm V}/^\circ$ $a_{\rm EM}/{\rm AU}$ $e_{\rm EM}$ years

mass /AU /AU

factor

Gliese 876

2.1 1 no V no V 0.10 0 no V no V ... ... EM-G 0.8

10 ... ... EM-G 1.0

20 ... ... EM-G 1.0

30 ... ... EM-G 1.0

40 ... ... EM-G 1.0

Ups And

G1 0.71 1 no V no V 1.30 0 no V no V 1.27-1.34 0-0.340 EM-G3 $1.4\;10^{2}$

G2 2.11 1.35 1.32-1.39 0-0.234 EM-G3 $1.2\;10^{2}$

G3 4.61 1.40 1.38-1.47 0-0.132 EM-G3 $1.2\;10^{7}$

1.45 1.45-1.50 0-0.118 EM-G3 $7.7\;10^{4}$

1.50 1.49-1.56 0-0.276 EM-G2 $2.1\;10^{5}$

3.40 3.40-3.63 0-0.105 EM-G3 $1.4\;10^{6}$

The giants in Ups And are labelled G1, G2, G3 out from the star. The minimum giant masses are used.

Table 5 shows a similar confirmation of expectations for Ups And. Even with the minimum giant masses none of the orbits in HZ(0) into which we launched EM are stable. Rivera & Lissauer (2000) have studied small objects in Ups And and though they used somewhat different orbital parameters for the giants, they found instabilities in the HZ on time scales comparable to and rather shorter than ours. Studies of the Ups And system are made difficult by the very short orbital period of the inner giant, only 4.617 days. Therefore extremely short time steps are required, less than about 0.38 days, and consequently several hundred hours of CPU time per 10⁸ years are needed.

5 Comparison with other work

To our knowledge, the work of Gehman et al. (1996) is the closest to ours. They made an analytical study of systems consisting of a giant planet in a circular orbit around a star, plus a single terrestrial planet in a circular coplanar orbit. They obtained conditions for Hill stability. One of the systems studied was 47 UMa. At the time of their work the giant's mass was estimated as $2.4~ {m}_{\rm J} \sin(i_{\rm o})$ and its orbital eccentricity as 0.03, not very different from the zero they need to assume in their approach. The current values are $2.41~ {m}_{\rm J} \sin(i_{\rm o})$ , and 0.096. For the minimum mass $(2.4~ {m}_{\rm J})$ they found that the Hill stable zone interior to the giant's orbit extends inwards from about 1.6 AU.

We have integrated their system and found that the stable zone extends inwards from about 1.4 AU. However, our $3 R_{\rm H}$ criterion is rather stiffer than their Hill stability criterion and so in this case only we have used a hybrid integrator (Chambers 1999) that is stable and accurate well inside the $3 R_{\rm H}$ limit of the "pure'' MVS integrator. When we use it with a stability criterion similar to Gehman et al. then we too find that the stable zone extends inwards from about 1.6 AU.

6 Conclusions

We conclude that Rho CrB and 47 UMa could have terrestrial planets in orbits that remain confined to their habitable zones for biologically significant lengths of time - for at least a few hundred million years. We also conclude that the Gliese 876 and Ups And systems are very unlikely to have planets in such orbits. To a first order the masses and initial orbital inclinations of the terrestrial planets have little effect on the outcome; it is the mass of the giant, and the initial semimajor axes of all the planetary orbits that matter. When a second terrestrial planet is present, close encounters between the two terrestrial bodies reduce the range of initial semimajor axes for which the orbits are stable.

In 47 UMa mean motion resonances cut strips of instability in the habitable zone. Away from resonances and boundaries between stable and unstable zones there is a tendency for instability in terms of the $3 R_{\rm H}$ criterion to show itself within 10⁸ years.

The stability of the orbits of terrestrial planets in the habitable zones of known exoplanetary systems

1 Introduction