Obliquity variations of terrestrial planets in habitable zones

Icarus 168 (2004) 223–236 www.elsevier.com/locate/icarus

Obliquity variations of terrestrial planets in habitable zones Keiko Atobe,a,∗ Shigeru Ida,a and Takashi Ito b a Department of Earth and Planetary Sciences, Faculty of Science, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo 152-8551, Japan b Astronomical Data Analysis Computer Center, National Astronomical Observatory, Mitaka, Tokyo 181-8588, Japan

Received 24 March 2003; revised 15 October 2003

Abstract We have investigated obliquity variations of possible terrestrial planets in habitable zones (HZs) perturbed by a giant planet(s) in extrasolar planetary systems. All the extrasolar planets so far discovered are inferred to be jovian-type gas giants. However, terrestrial planets could also exist in extrasolar planetary systems. In order for life, in particular for land-based life, to evolve and survive on a possible terrestrial planet in an HZ, small obliquity variations of the planet may be required in addition to its orbital stability, because large obliquity variations would cause significant climate change. It is known that large obliquity variations are caused by spin–orbit resonances where the precession frequency of the planet’s spin nearly coincides with one of the precession frequencies of the ascending node of the planet’s orbit. Using analytical expressions, we evaluated the obliquity variations of terrestrial planets with prograde spins in HZs. We found that the obliquity of terrestrial planets suffers large variations when the giant planet’s orbit is separated by several Hill radii from an edge of the HZ, in which the orbits of the terrestrial planets in the HZ are marginally stable. Applying these results to the known extrasolar planetary systems, we found that about half of these systems can have terrestrial planets with small obliquity variations (smaller than 10◦ ) over their entire HZs. However, the systems with both small obliquity variations and stable orbits in their HZs are only 1/5 of known systems. Most such systems are comprised of short-period giant planets. If additional planets are found in the known planetary systems, they generally tend to enhance the obliquity variations. On the other hand, if a large/close satellite exists, it significantly enhances the precession rate of the spin axis of a terrestrial planet and is likely to reduce the obliquity variations of the planet. Moreover, if a terrestrial planet is in a retrograde spin state, the spin–orbit resonance does not occur. Retrograde spin, or a large/close satellite might be essential for land-based life to survive on a terrestrial planet in an HZ.  2003 Elsevier Inc. All rights reserved. Keywords: Celestial mechanics; Earth; Extrasolar planets; Resonances; Terrestrial planets

1. Introduction More than 110 extrasolar planets around nearly 100 Sunlike stars are now known to exist (see, e.g., http://cfa-www. harvard.edu/planets/). The main technique used to detect these planets is to observe the reflex motion of host stars due to planetary motion through the Doppler shift of stellar spectral lines (e.g., Marcy et al., 2000). The planets so far found are as massive as Jupiter, so all the known extrasolar planets are considered to be Jupiter-like giant planets. However, it is reasonable to expect that terrestrial planets (relatively small rocky planets) could be present in some of the extrasolar planetary systems where a giant planet(s) has

* Corresponding author.

E-mail address: [email protected] (K. Atobe). 0019-1035/$ – see front matter  2003 Elsevier Inc. All rights reserved. doi:10.1016/j.icarus.2003.11.017

been detected. Are there terrestrial planets like the Earth that could habor life among extrasolar planetary systems? In general, organisms which we are familiar with on Earth require liquid water during at least part of their life cycle. In order for life to emerge on a terrestrial planet, it is also necessary that its orbit remains confined to the habitable zone (HZ) of its host star over the length of time required for biological evolution. A stellar HZ is the range of distances from a star allowing the existence of liquid water at the surface of a terrestrial planet (Abe, 1993; Kasting et al., 1993). The orbital stability of hypothetical Earth-like planets inside the HZs of extrasolar planetary systems has already been investigated through numerical simulations by many authors (Gehman et al., 1996; Jones et al., 2001; Jones and Sleep, 2002; Noble et al., 2002; Menou and Tabachnik, 2003). Here, orbital stability means that the orbits of the hypothetical terrestrial planets remain confined to stellar HZs over the time required

224

K. Atobe et al. / Icarus 168 (2004) 223–236

for the evolution of life. In general, if the orbits of giant planets are located in the vicinity of the HZs, terrestrial planets cannot exist long enough on stable orbits inside the HZs because of the strong perturbing influence of the giant planets. Menou and Tabachnik (2003) studied the orbital stability of 85 known extrasolar planetary systems. Their results indicate that more than half of the 85 systems are unlikely to harbor habitable terrestrial planets. Not only long-term orbital stability of terrestrial planets in stellar HZs but also moderate climate on these planets may be necessary for the evolution of life, in particular for a land-based life. Planetary climate is greatly influenced by insolation distribution and intensity, which depend on the orbital parameters and spin motion of the planet, especially eccentricity, obliquity, and precession motion (Milankovitch, 1941; Berger, 1984; Berger, 1989). Obliquity is the angle of a spin axis of a planet relative to its orbital plane. Because of their equatorial bulge, planets are subject to the torque arising from the gravitational force of a host star (and possibly of their satellites and of other planets). This torque causes precessional motion of the spin axis of the planet about its orbit normal. On the other hand, gravitational interactions with other planets cause the orbit normals to precess and nutate about the normal to the invariable plane of the system. So the obliquity of a planet generally changes periodically. We study the amplitude of obliquity variations of terrestrial planets in HZs perturbed by gas giants. Obliquity changes the latitudinal distribution of yearly insolation on the planet, which could be an important driver of climate variations. Ward (1974) showed that the total fluctuation in the Earth’s obliquity is about 1.5◦ –2.0◦ in a period of about 4 × 104 years. According to the more detailed calculation, the Earth’s obliquity varies by ± 1.3◦ around its mean value of 23.3◦ (Laskar and Robutel, 1993; Laskar, 1996). The glacial cycle of the Earth has been triggered by the periodic change of obliquity to some extent, through the variation of insolation at the edge of a glacier in the northern hemisphere (e.g., Crowley and North, 1991). Even such small variation of obliquity has produced significant variations in the Earth’s climate, at least during the Quaternary. Mars’ obliquity is believed to have suffered from a largescale oscillation with an amplitude of more than ∼ 10◦ around its average of ∼ 25◦ on a time scale of ∼ 105 –106 years (Ward, 1973, 1974, 1979). Such a large variation is the result of a coupling between the precession rate of Mars’ spin axis and one of the frequencies of its orbital precession. This is called the spin–orbit resonance. Laskar and Robutel (1993) showed that the obliquity of terrestrial planets other than the Earth could have experienced large and chaotic variations throughout their histories. As a result of chaotic variations on a longer timescale, the range of Mars’ obliquity change can be as large as from 0◦ to 60◦ . Furthermore, obliquity of an Earth without the Moon would vary radically in the range from 0◦ to more than 85◦ (Laskar and Robutel, 1993;

Laskar et al., 1993; Laskar, 1996). If the Earth’s obliquity had such a large variation, its climate would drastically change, providing a significant obstacle to the development of complex life on Earth. The previous analyses of the behavior of planetary obliquity, however, have been concerned mostly with the present orbital configuration of the planets in our Solar System. We are interested in the evaluation of obliquity of possible terrestrial planets in extrasolar planetary systems with various orbital configurations of planets. For this purpose, it is useful to introduce generalized analytical formulae which approximately express the amplitude of the obliquity variations as functions of planetary dynamical parameters. The oscillation amplitude of obliquity in nonresonant regions has been derived through a linear analysis of the equation of spin motion (e.g., Sharaf and Budnikova, 1969a, 1969b; Vernekar, 1972; Ward, 1974; Berger, 1978). With the analytical formulae for non-resonant regions, Ward et al. (2002) showed the amplitude for a mass less body in the place of Mars as a function of the rotational period and the semimajor axis of the particle. However, the nonresonant formulae cannot describe the oscillation amplitude of obliquity in the regions near spin–orbit resonances where the oscillation amplitude is not small. The variation of planetary obliquity is also written in Hamiltonian form (e.g., Kinoshita, 1977; Laskar and Robutel, 1993; Laskar, 1996), and analyses of the spin–orbit resonance in the Solar System are available (e.g., Colombo, 1966; Henrard and Murigande, 1987; Laskar and Robutel, 1993; Touma and Wisdom, 1993; Laskar, 1996). Analytical expression of obliquity in resonance is given by Ward (1982) and Correia and Laskar (2003). In this paper, we reanalyze the spin–orbit resonance, assuming a simple planetary system to derive the generalized formulae for extrasolar planetary systems. Calculating the obliquity variations using the analytical formulae, we investigate how strongly possible terrestrial planets in the known extrasolar planetary systems are affected by the spin–orbit resonances. We find that possible terrestrial planets in only 1/5 of the known extrasolar planetary systems exhibit both small obliquity variations (
10 ◦ , where  is the variation amplitude of obliquity) and orbital stability over their entire HZs and most of them are systems with close-in giant planets. Our results show that terrestrial planets on stable orbits usually suffer large obliquity variations unless they have retrograde spin or a large/close satellite. In Section 2, we briefly summarize the spin–orbit resonance and introduce the analytical expressions of obliquity variations. Using the analytical formulae, we discuss dynamical properties of giant planets causing large obliquity variations of terrestrial planets in HZs, in Section 3. In Section 4, we apply our results to the known extrasolar planetary systems. Sections 5 and 6 are devoted to discussion and conclusion.

Obliquity of extrasolar planets

225

2. Obliquity variation In this section, we briefly summarize basic equations and analytical expressions of the amplitude of obliquity variations which we use in this paper. We investigate obliquity variations of a hypothetical massless terrestrial planet in a planetary system with a given giant planet(s). Spin axis sˆ of a terrestrial planet precesses around its orbit normal (the vector normal to the orbital plane) nˆ because of the gravitational torque raised by other objects, following (e.g., Ward, 1974) d sˆ = α(ˆs · n)(ˆ ˆ s × n), ˆ (1) dt where sˆ and nˆ are scaled to be unit vectors. We assume that the spin axis coincides with the principle axis and neglect the wobble, because the planet’s rotation is usually sufficiently faster than the precession and the nutation. The precession constant α is 3Gmc ED α= 4πa 3 −3

mc a −5 8.6 × 10 1M 1 AU −1
D × (rad yr−1 ), (2) 24 h where G is the gravitational constant, mc is the mass of the central star, D and a are the spin rotation period and the semimajor axis of the terrestrial planet, and E is its dynamical ellipticity, assuming E ∝ D −2 (e.g., Laskar and Robutel, 1993). In Eq. (2), we only consider the torque produced by a central star. The torque caused by the giant planets can be neglected. The torque produced by satellites can be taken into account by appropriately increasing α (e.g., Ito et al., 1995). Actually, the Earth’s α is increased by a factor of 3 by the Moon. Obliquity  of the terrestrial planet (the angle between sˆ and n; ˆ sˆ · nˆ = cos ) changes through change in nˆ caused by the gravitational perturbations exerted by other planets, which are the giant planet(s) here. Let I and Ω be the orbital inclination and the longitude of the ascending node of the terrestrial planet referred to the invariant plane (see Fig. 1). Because we are interested in terrestrial planets on nearly circular orbits that have not been largely affected by close scattering or orbital resonances, their orbits are nearly on the equatorial plane of the original protoplanetary disks, which may be similar to the invariant plane. Hence we assume that is I 1 and use the secular perturbation theories (e.g., Brouwer and Clemence, 1961; Murray and Dermott, 1999) to describe variations of I and Ω: I sin Ω = I0 sin(Bt + γ ) +

Ne 

Ii sin(fi t + γi ),

(3)

i=1

I cos Ω = I0 cos(Bt + γ ) +

Ne  i=1

Ii cos(fi t + γi ),

(4)

Fig. 1. Relationship between the reference plane, orbital and equatorial planes of the terrestrial planet. sˆ is the unit vector in the direction of the spin axis, nˆ is that normal to the orbital plane,  is the obliquity, θ is the angle of the equatorial plane relative to the reference plane, ψ is the longitude of the equatorial plane, I is the orbital inclination, and Ω is the longitude of the ascending node relative to the reference plane.

where fi are the eigenfrequencies of this system (i = 1, 2, . . . , Ne ). Free oscillation frequency B and forced inclinations Ii are given by Eqs. (A.1) and (A.2) in Appendix A. Free inclination I0 , and phase parameters γ and γi are determined by initial conditions. Substitution of Eqs. (3) and (4) into (1) through sˆ = (sin θ sin ψ, − sin θ cos ψ, cos θ ) where ψ and θ are the longitude and inclination of the equatorial plane of the terrestrial planet to the invariant plane (Fig. 1), and nˆ = (sin I sin Ω, − sin I cos Ω, cos I ) yields (see Appendix A)  dψ −α cos  + α cos  cot θ I0 cos(ψ − Ω0 ) dt  Ne  + (5) Ii cos(ψ − Ωi ) + O(I 2 ), i=1

d(cos ) I0 B sin θ sin(ψ − Ω0 ) dt Ne  + Ii fi sin θ sin(ψ − Ωi ) + O(I 2 ),

(6)

i=1

where Ω0 ≡ Bt + γ and Ωi ≡ fi t + γi . Since I 1, the ˙ is −α cos . When it approaches precession frequency (ψ) one of the orbital precession frequencies, i.e., ψ˙ − Ω˙ 0 0 or ψ˙ − Ω˙ i 0, d(cos )/dt constant and  changes secularly to have a large amplitude. This is the spin–orbit resonance. Since Ω˙ 0 and Ω˙ i are always negative (e.g., Eq. (A.1)), spin– orbit resonance occurs only for the planets with prograde spin ( < 90◦ ). Far from the spin–orbit resonance, oscillation amplitude of obliquity, , is obtained by expanding Eq. (1) around the

226

K. Atobe et al. / Icarus 168 (2004) 223–236

initial obliquity ini with I 1 and  1. The maximum  is written as       Ne    fi Ii BI0     + ||max =  (7)   α cos ini + B α cos ini + fi  i=1

(e.g., Ward, 1974; Ward et al., 2002). The above formula is no longer applicable at α cos ini + B 0 (spin–orbit resonance of ψ˙ − Ω˙ 0 0), and at α cos ini + fi 0 (that of ψ˙ − Ω˙ i 0). Note that B − fi 0 is also a singularity because Ii diverges (secular resonance; see Eq. (A.2)). Next, we consider  near the spin–orbit resonance. We consider a system with a giant planet and a terrestrial planet. In this case, the invariant plane of the system (reference plane) coincides with the orbital plane of the giant planet, so that there is no forced motion (fi = 0, Ii = 0) in this frame. Then, I = I0 , which does not change with t, and Ω = Bt +γ (= Ω0 ) (Eqs. (3) and (4)). √ Let x = ψ − Ω0 and t˜ = αI0 |B| sin r t where r is the resonant (ψ˙ Ω˙ 0 = B) obliquity given by cos r = |B|/α. Taking only the first order term in I , Eq. (5) is reduced to dx B + α cos  . = y = −√ dt˜ αI0 |B| sin r

(8)

Since I 1, θ . Then Eq. (6) is d2 x sin  (9) = sin x sin x. sin r dt˜ 2 Integrating Eq. (9) with Eq. (8), we obtain an integration constant H as 1 H = cos x + y 2 , (10) 2 H should be essentially equivalent to a pendulum-type Hamiltonian that has been discussed in the previous literature, such as in Laskar and Robutel (1993). Trajectories with H < 1 show libration around the stable point (x = π and y = 0). The maximum libration in y is −2 to 2 when H → 1. If −2 < y < 2, there is x with which a trajectory has the maximum libration. We call the region of −2 < y < 2 the “resonant region.” Given α, B and I0 ,  uniquely corresponds to some particular value of y (Eq. (8)). The maximum amplitude of libration in obliquity is   | cos |max = 2 I0 |B| sin r /α = 2 I0 cos r sin r . (11) The range of  in the resonant region is cos r − | cos |max
cos 
cos r + | cos |max . When the variation is small enough, | cos |max sin r ||max , so that ||max = √ 2 I0 / tan r . These expressions are consistent with the formulae given by Ward (1982) or Correia and Laskar (2003). Note that the resonant half width | cos |max (||max ) has direct dependence only on I0 and r . It depends indirectly on the other parameters (the semimajor axis a, the rotation period D of the terrestrial planet, the semimajor axis aj , and mass mj of the giant planet) through r . Also note that Eq. (7) is immediately reproduced by a similar argument for nonresonant region (H > 1) where x circulates

Fig. 2. Examples of obliquity variation of a terrestrial planet in the system with one giant planet. We calculated the case with mc = 1M , m1 = 1MJ (MJ is the mass of Jupiter), a1 = 5.2 AU, a = 1.1 AU, D = 24 h and I0 = 1.3◦ . (a) The maximum and minimum values of  during the integration on 3 × 106 years as a function of initial obliquity ini . (b) The maximum oscillation amplitude of  of the terrestrial planet as a function of ini . Crosses show the numerical values. The solid line shows analytical expressions given by Eq. (11), and dashed lines are those given by Eq. (7).

with small oscillation amplitude in y. If r is replaced by ini , |y| |y| (y ≡ y − y(x = π)), so that y − cos x/y (see Eq. (10)) and ||max 1/|y|| sin ini |, resulting in Eq. (7). We numerically integrate Eq. (1) with Eqs. (3) and (4). Figures 2 show maximum and minimum values of  of a hypothetical terrestrial planet with various ini from 0◦ to 90◦ integrated over 3 × 106 years, which is longer than a few times the typical timescale of spin–orbit resonance. For each ini , the numerical results with initial ψ − Ω0 from 0◦ to 360◦ at every 10◦ are plotted together. The other parameters are explained in the caption. The resonant region is between − = 36◦ and + = 66◦ , centered on r = 51.7◦ . The square region centered on r in Fig. 2a is explained by the Hamiltonian argument (also see Laskar and Robutel (1993)). In the resonant region −
ini
+ , the maximum variation

Obliquity of extrasolar planets

227

range is the same for any ini . In Fig. 2b, we plot the maximum amplitude of  as a function of ini . The crosses are the numerical values, the solid lines are analytical calculations using Eq. (11), and the dashed lines are those of Eq. (7). Combining Eqs. (11) and (7), we can analytically evaluate obliquity amplitude within 10% accuracy. If other giant planets exist, in addition to the resonance associated with free oscillation (ψ˙ Ω˙ 0 ), that with forced oscillation (ψ˙ Ω˙ i ) also exist. In the same way as Eq. (11), |fi | , α  | cos |max,i = 2 Ii cos r,i sin r,i ,

cos r,i =

(12) (13)

where i = 1, 2, . . . , Ne . If they are isolated, total  is summation of  of individual resonances (e.g., Correia and Laskar (2003), see also Fig. 3a). However, if the resonant regions [−0

+0 for free oscillation and −i

+i for forced one (i = 1, 2, . . .)] overlap,  is not the simple summation but is described by the extended merged resonance of mini (−i )

maxi (+i ) (i = 0, 1, 2, . . .) (Fig. 3b).

3. Obliquity variation in a habitable zone in extrasolar planetary systems In Section 2, we have introduced analytical expressions for the amplitude of obliquity variation (||max ). Hereafter, we denote ||max by . Using analytical formulae, we study dynamical properties of giant planets causing large obliquity variations of terrestrial planets in a habitable zone (HZ). As mentioned in Section 1, for land-based life to evolve and survive on a terrestrial planet in an HZ, small or limited obliquity variation of the planet may be required, in addition to the planet’s orbital stability. Because we do not know the specific locations of terrestrial planets, their spin periods, and critical  which has a severe effect on landbased life in extrasolar planetary systems, we calculate the probability that a terrestrial planet has  larger than some given values (crit ), averaged over an entire HZ and likely spin periods. As mentioned in Section 2, spin–orbit resonance occurs when ψ˙ B or ψ˙ fi . Corresponding obliquity variation amplitude  in resonance is estimated by Eqs. (11) and (13), and the amplitude in nonresonance is estimated by Eq. (7). Given the initial orbital elements and masses of the host star and the giant planet(s) in the system, B depends only on the semimajor axis a of the terrestrial planet (Eq. (A.1) in Appendix A) and α depends on a and its spin period D (Eq. (2)). In general, there are multiple fi which depends on the masses mj and the semimajor axis aj of the planets but are independent of a in multiplanet systems. Figures 4 show the examples of obliquity variation amplitudes for the terrestrial planet, as a function of a and D: • system 1: m1 = 1MJ , a1 = 0.1 AU, I0 = 0.05 (Fig. 4a);

Fig. 3. Maximum oscillation amplitude of obliquity of a terrestrial planet in the system with two giant planet, m1 = 1MJ , a1 = 5.2 AU, m2 = 1MS (MS is the mass of Saturn), a2 = 9.5 AU, (a) a = 0.8 AU, and D = 24 h and (b) a = 1.4 AU, and D = 2.4 h. Crosses show the numerical values. Dashed lines show analytical values by the combination of Eqs. (11), (13) and (7) if each resonance (ψ˙ B and ψ˙ fi ) is isolated. Solid lines show the merged values.

• system 2: m1 = 1MJ , a1 = 1.0 AU, I0 = 0.05 (Fig. 4b); • system 3: m1 = 1MJ , a1 = 0.1 AU, I12 = 0.05; m2 = 1MJ , a2 = 1.0 AU, I22 = 0.05 (Fig. 4c). Here we chose the parameters such that the resonances (see below) occur near a = 1 AU. In Figs. 4, we adopt the mass of the host star mc = 1M and initial obliquity ini = 25◦ . Figures 4 show that  depends on a and D. The regions where obliquity has a large amplitude can be categorized into 3 types: (I) The regions of spin–orbit resonances resulting from ψ˙ Ω˙ 0 ( B). (II) The regions of spin–orbit resonances resulting from ψ˙ Ω˙ i ( fi ) (i = 1, 2, . . . , Ne ). (III) The regions of secular resonances, B = fi (Ii has large values according to Eq. (A.2) in Appendix A).

228

K. Atobe et al. / Icarus 168 (2004) 223–236

(a)

proximate expression of the type I resonant region as follows: the dependence of B on a comes from the a(1) (α1 ), where dependence of the Laplace coefficients b3/2 α1 = min(a, a1 )/max(a, a1) (see Eqs. (A.1) and (A.3) in (1) Appendix A). In general, b3/2 (α1 ) cannot be expressed with a simple function of a and a1 . However, when the orbits of terrestrial and giant planets are sufficiently separated, i.e., α1 1, the numerical value of the Laplace coefficient approaches 3(a1 /a) at a > a1 and 3(a/a1) at a < a1 . Substituting these expressions into Eq. (A.1), B can be written as


 mc −1/2 m1 a1 p −3 B −4.5 × 10 1M 1MJ 1 AU
q a × (rad yr−1 ), (14) 1 AU where p = 2 and q = −7/2 for a > a1 and p = −3 and q = 3/2 for a < a1 . Since ψ˙ B at the spin–orbit resonance, we obtain the spin period at the resonance, Dr , as



  cos ini mc 3/2 m1 −1 a1 p Dr 0.42 cos 25◦ 1M 1MJ 1 AU q 
a × (h), (15) 1 AU

(b)

(c) Fig. 4. (a) Obliquity variation amplitude () for a terrestrial planet as a function of a and D, in the case of m1 = 1MJ , a1 = 0.1 AU and I0 = 0.05 (system 1). The color bar indicates  in degree. (b) Obliquity variation amplitude for a terrestrial planet as a function of a and D, in the case of m1 = 1MJ , a1 = 1.0 AU and I0 = 0.05 (system 2). At a = 1.0 AU (the location of the giant planet), we set  = 0◦ . (c) Obliquity variation amplitude for a terrestrial planet as a function of a and D, in the case of m1 = 1MJ , a1 = 0.1 AU, I12 = 0.05, m2 = 1MJ , a1 = 1.0 AU, and I22 = 0.05 (system 3). At a = 1.0 AU (the location of the giant planet), we set  = 0◦ .

It is noted that regions II and III do not appear in the systems with one giant planet. In the case of a system with one giant planet, B can be expressed in a simple form. We derive a simple ap-

where p = −2 and q  = 1/2 for a > a1 and p = 3 and q  = −9/2 for a < a1 . Although this formula is valid only when the orbits of the giant and the terrestrial planets are well separated, it provides a good account of the characteristics of resonant regions on the a–D plane. We can see that Dr approximately ∝ a 1/2 in Fig. 4a and that Dr approximately ∝ a −9/2 in Fig. 4b. When a system contains multiple giant planets, B and Dr cannot be expressed in a simple form. However, we can approximate the values of Dr , because B in a multigiant planet system is given by the sum of B in single-giant planet systems. Superposing Fig. 4a on Fig. 4b, type I resonant regions crosses at a ∼ 0.4 AU and D ∼ 15 h. Removing the portion of type I resonant regions with D > 15 h in Figs. 4a and 4b and connecting them, we obtain the resonant region seen at  0.5 AU in Fig. 4c. A type II resonant region is found in the region at 0.5–1 AU in Fig. 4c. Since the eigenfrequencies fi (i = 1, 2, . . . , Ne ) are independent on a, Dr ∝ a −3 (see Eq. (2)). When B = fi , obliquity has large amplitude irrespective of the spin–orbit relation because Ii given in Eq. (A.2) is large (type III resonant region). A type III region appears at ∼ 2.7 AU in Fig. 4c. Since B and fi (i = 1, 2, . . . , Ne ) do not depend on D, this region is parallel to the D-axis. In order to quantify the effect of spin–orbit resonances on the terrestrial planets in an HZ, we did a set of calculations as follows: Considering the limited area on the a–D plane with appropriate D and a for the terrestrial planets in an HZ, we computed the fraction of the area where  is larger than some given values (cirt ) to evaluate the probability of the terrestrial planets having  > crit . Here, we

Obliquity of extrasolar planets

229

Table 1 Examples of area fractions System System 1 System 2 System 3

  10◦*1 (%)

  20◦*2 (%)

  40◦*3 (%)

45.86 0.00 16.38

31.39 0.00 7.50

6.15 0.00 0.12

Note. * The percentages indicate the fraction of the area where  is larger than (1) 10◦ , (2) 20◦ , and (3) 40◦ in the area with D = 8–72 h and a = 0.7–1.3 AU (HZ).

adopt crit = 10◦ , 20◦ , and 40◦ . Since we consider terrestrial planets like the Earth, we here take 8–72 h as the appropriate range of D, and a = 0.7–1.3 AU as an HZ (the HZ around a zero-age star with 1M ). We also adopt ini = 25◦ . The different choices of the range of D, HZs, and ini do not significantly affect the global features of the probability distributions we show below. The choice of I does not affect the global features either, but as shown in Eqs. (11), (13) and (7), larger I results in a larger variation of  and enhances area fractions in what follows. We discuss the choice of I for the known extrasolar systems in the next section. Table 1 shows area fractions with   10◦ ,  20 ◦ and  40◦ in systems 1, 2, and 3 in Figs. 4. The area fractions in this table indicate the degree of overlapping of resonant regions with HZs. The type I resonant region overlaps with the HZ in system 1, but not in system 2. In system 3, the type II resonant region overlaps with the HZ. Figure 5 shows the results of similar procedures for systems with one giant planet with a different semimajor axis a1 , in the case of m1 = 1MJ , I = 0.05 and initial obliquity ini = 25◦ . Spin–orbit resonance in the HZ is most effective when a1 ∼ 0.1 and 5 AU, and it is least effective when a1 ∼ 1 AU. We explain this dependence on a1 , using Eq. (15). The resonant spin period Dr increases with decreasing a1 for a > a1 or with increasing a1 for a < a1 . Equation (15) is equivalent to  cos ini 1/2 mc 3/4 m1 −1/2  0.13 cos  25◦ 1M 1MJ  

   × a 1/4 D −1/2 (AU) for a > a1 , 1 AU 24 h a1r cos ini −1/3 mc −1/2 m1 1/3 a 3/2 (16)    25◦ 1M 1MJ 1 AU  3.9 cos  D 1/3  × 24 h (AU) for a < a1 , where a1r is the semimajor axis of the giant planet that causes a spin–orbit resonance at a. This equation implies that the obliquity variation of terrestrial planets with a 1 AU and D 24 h is in a resonance when a1 ∼ 0.13 and ∼ 3.9 AU in the cases with ini = 25◦ , mc = 1M , and m1 = 1MJ , which is consistent with Fig. 5. When the giant planet is too close to the terrestrial planets, a1 ∼ a (∼ 1 AU), B is much larger than α cos r with D = 8–72 h. In this case, very small D (rapid rotation) is necessary for the terrestrial planets to be in resonance. However, such a rapid rotation might break up the planet. On the other hand, when the giant planet is too far from terrestrial planets (a1  10 AU or a1  0.05 AU), B is too small for the spin–orbit resonance.

Fig. 5. The area fraction with   10◦ (solid line),  20◦ (dashed line), and  30◦ (dot-dashed line) as a function of the semimajor axis of a giant planet a1 , in the case of mc = 1M , m1 = 1MJ and I0 = 0.05.

For the spin–orbit resonance to happen in this case, the spin period must be much longer than 24 h, which is unrealistic unless tidal force exerted by a host star or thick atmosphere significantly damps the spin. Thus, Eq. (16) explains the dependence on a1 in Fig. 5. Figures 6 show the area fraction as a function of the mass of a giant planet m1 , in the cases of I = 0.05 and a1 = 0.1 AU (Fig. 6a), and a1 = 5.0 AU (Fig. 6b). In these cases, m1 ∼ 1MJ is the most effective. Equation (15) is rewritten as

m1r

 cos ini mc 3/2 a1 −2 a 1/2   1.7 cos 25◦ 1M 0.1 AU 1 AU  

  × D −1 (MJ ) for a > a1 , 24 h cos ini mc 3/2 a1 3 a −9/2   2.2 cos 25◦ 1M  5 AU 1 AU   D −1  × 24 h (MJ ) for a < a1 ,

(17)

where m1r is the mass of the giant planet at a1 that causes spin–orbit resonances at a. When a ∼ 1 AU, the most effective m1 is estimated as 1.7MJ for at a1 = 0.1 AU, and as 2.2MJ for a1 = 5 AU. These are consistent with Figs. 6. Figure 7 shows the area fraction with   20◦ as a function of a1 and m1 . The giant planets with m1  1MJ affect the obliquity of terrestrial planets in HZs when they are located at 0.05–0.1 AU or at 3–5 AU. The obliquity of terrestrial planets at ∼ 1 AU is hardly affected when the orbit of a giant planet is close enough to the HZ when m1  1MJ . On the contrary, the orbits of terrestrial planets become unstable when the giant planet is located in the vicinity of the HZ. The giant planets which cause the spin–orbit resonance are different from those which destabilize the orbits of the terrestrial planets. The above estimation is valid for systems with a single giant planet. When other giant planets are in the system and other resonant regions appear, the probability of terrestrial planets in an HZ having a large  would generally increase, although modulation of resonant regions sometimes decreases the probability as seen in Fig. 4c and Table 1.

230

K. Atobe et al. / Icarus 168 (2004) 223–236

Fig. 6. The area fraction with   10◦ (solid line),  20◦ (dashed line), and  30◦ (dot-dashed line) as a function of mass of a giant planet m1 , in the case of mc = 1M and I0 = 0.05. The semimajor axis of the giant planet a1 is (a) 0.1 AU and (b) 5.0 AU.

Fig. 7. The area fraction with   20◦ as a function of a1 and m1 . The color bar indicates the area fraction in percent.

4. Applications to known extrasolar planetary systems We apply the results in the previous sections to known extrasolar planetary systems to evaluate obliquity variations of possible terrestrial planets in their HZs. The samples of extrasolar planetary systems used here are listed in Tables 2, 3, and 4 with relevant data on the stellar masses, and the planets’ masses, semimajor axes and eccentricities. These data were taken from the Extrasolar Planets Encyclopedia (http://cfa-www.harvard.edu/planets/). We categorized extrasolar planetary systems into 4 groups (A, B1, B2, and C) according to the degree of the obliquity variation of hypothetical terrestrial planets in the HZs (see below). In the current radial velocity surveys, only mp sin i (mp is the planet’s mass, and i is the inclination between the lineof-sight of observers and a normal to the orbital plane of planets) is measured. Since we have few observational constraints on i, we adopt planetary masses corresponding to i = 45◦ , which means 1.4 times as large as the minimum

values assuming i = 90◦ . As mentioned in Section 2, we assume that possible terrestrial planets have nearly circular orbits on their protoplanetary disks’ planes. The inclinations between the orbital plane of the giant planets and the original disks are not known. If orbital eccentricities ep of the giant planets are produced by close scattering with other planets, the inclinations of the giant planets I would be ∼ ep /2 (Ida et al., 1993; Marzari and Weidenschilling, 2002) on average. The eccentricities of giant planets can also be pumped up by disk-planet interactions or distant resonant perturbations by other planets. Then Ip ∼ ep /2 does not necessarily hold. However, here we adopt the simple assumption Ip ∼ ep /2. We will discuss the variety of choices of Ip later. In Section 4, we adopted the region between 0.7 and 1.3 AU as an HZ for simplicity. A variety of criteria have been used to define the boundaries of an HZ. According to Kasting et al. (1993), for the inner boundary of the HZ, we use the minimum distance from the host star where a runaway greenhouse effect occurs leading to the evaporation of all surface water. As for the outer boundary of the HZ, the maximum distance at which a cloud-free massive CO2 atmosphere of a terrestrial planet can maintain a surface temperature of 273 K. The location of the HZ depends on the stellar mass and moves with stellar age because stellar luminosity changes with age. With the estimated stellar mass, we use zero-age main sequence HZs for simplicity. The inner and outer radii of HZs are listed in the 6th column of Tables 2, 3, and 4 for each system. As stated in Section 1, orbital stability of possible terrestrial planets has been investigated in previous research. Menou and Tabachnik (2003) investigated the orbital dynamics of 85 extrasolar planetary systems, integrating 100 terrestrial planets (treated as test particles) in HZs for 106 years. They found that orbital stability is correlated with the degree of overlap between an HZ and the gravitational zone of influence of a discovered giant planet, defined by the region extending from Rin = (1 − ep )ap − 3RHill to

Obliquity of extrasolar planets

231

Table 2 Extrasolar planetary systems in group A System

M∗ (M )

Mp sin i (MJup )

a (AU)

e

HZ (AU)

  10◦*1 (%)

  20◦*2 (%)

  40◦*3 (%)

Class

HD46375 HD187123 HD209458 HD179949 HD75289 BD-10 3166 HD76700 51Peg HD49674 HD168746 GI86 HD195019 rhoCrB HD114762 HD121504 HD178911 HD52265 HD73526 HD134987 HD40979 HD169830 HR810 HD28185 HD128311 HD108874 HD27442 HD114783 HD20367 HD19994 HD23079 HD4028 HD10697 HD196050 HD216435 HD30177 HD23596 HD72659 GI777A HD83443

1.00 1.06 1.05 1.24 1.05 1.10 1.00 0.95 1.00 0.92 0.79 1.02 0.95 0.82 1.00 0.87 1.13 1.02 1.05 1.08 1.40 1.03 0.99 0.80 1.00 1.20 0.92 1.05 1.35 1.10 0.93 1.10 1.10 1.25 0.95 1.10 0.95 0.90 0.79

0.000 0.030 0.000 0.050 0.054 0.000 0.000 0.000 0.000 0.081 0.046 0.050 0.028 0.334 0.130 0.124 0.290 0.340 0.240 0.250 0.340 0.161 0.070 0.210 0.200 0.020 0.100 0.230 0.200 0.060 0.040 0.120 0.280 0.140 0.220 0.314 0.180 0.000 0.080 0.420 0.020 0.340 0.160

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.02 0.00 0.00 0.00 6.61 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.25 0.00 1.34 0.00 0.00 1.72 7.48 0.00 9.77

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.16 0.00 0.00 0.00 0.00 0.00 3.37 0.00 2.99

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.06 0.00 0.00

1.03

0.041 0.042 0.045 0.045 0.046 0.046 0.049 0.050 0.057 0.065 0.110 0.140 0.230 0.300 0.320 0.320 0.490 0.660 0.780 0.818 0.823 0.925 1.030 1.060 1.070 1.180 1.200 1.250 1.300 1.480 1.690 2.000 2.500 2.600 2.600 2.720 3.240 3.650 0.038 0.174 0.115 0.241 5.900

0.70/1.30 0.75/1.40 0.75/1.40 1.10/2.25 0.75/1.40 0.85/1.60 0.70/1.30 0.70/1.30 0.70/1.30 0.65/1.25 0.50/1.00 0.70/1.30 0.70/1.30 0.50/1.00 0.70/1.30 0.60/1.20 0.85/1.60 0.70/1.30 0.75/1.40 0.85/1.60 1.40/3.00 0.70/1.30 0.70/1.30 0.50/1.00 0.70/1.30 0.93/1.80 0.65/1.25 0.75/1.40 1.33/2.85 0.85/1.60 0.65/1.25 0.85/1.60 0.85/1.60 1.10/2.25 0.70/1.30 0.85/1.60 0.70/1.30 0.60/1.20 0.50/1.00

55Cnc

0.249 0.520 0.690 0.840 0.420 0.480 0.197 0.470 0.120 0.230 4.000 3.430 1.100 11.000 0.890 6.292 1.130 3.000 1.580 3.160 2.960 2.250 5.700 2.630 1.650 1.430 0.990 1.070 2.000 2.540 0.810 6.590 3.000 1.230 7.700 7.190 2.550 1.150 0.350 0.170 0.840 0.210 4.050

0.70/1.30

7.13

3.56

0.69

I I I I I I I I I I I I I III II II–III II III–IV III III II–III IV IV IV IV III III–IV IV III III–IV II–III III–IV III–IV III–IV III–IV IV II–IV II I I I I III

Note. * The percentages indicate the fraction of the area where  is larger than (1) 10◦ , (2) 20◦ , and (3) 40◦ in the area with D = 8–72 h and a = 0.7–1.3 AU (HZ).

Rout = (1 + ep )ap + 3RHill , where ap and ep are the semimajor axis, eccentricity of the giant planet and RHill is the planet’s Hill radius defined by RHill ≡ ap (mp /3mc )1/3 . Menou and Tabachnik (2003) categorized orbital stability of the terrestrial planets in an HZ, based on the degree of overlap between the gas giant’s zone of influence (ZI) and the system’s habitable zone (HZ): class I (ap
0.25 AU), class II (ap > 0.25 AU and no overlap between HZ and ZI), class III (ap > 0.25 AU and partial overlap between HZ and ZI), and class IV (ap > 0.25 AU and HZ is fully inside ZI). Menou and Tabachnik (2003) also showed the probabilities of terrestrial planets remaining stable to be 73.9±8.3% for class I,

29.9 ± 17.3% for class II, 0.92 ± 2.31% for class III, and 0% for class IV. More orbits of terrestrial planets are unstable in systems with larger overlap. The orbital classification (I–IV) are listed in the last column of Tables 2, 3, and 4 in order to show the relation between orbital stability and obliquity variation. Tables 2, 3, and 4 show the probability with   10 ◦ ,  20◦ and  40 ◦ for terrestrial planets in HZs in 87 extrasolar planetary systems. Here, we consider spin periods of the terrestrial planets as 8–72 h and the HZs listed in the 6th column in the tables for each system. The extrasolar planetary systems can be classified into four distinct groups

232

K. Atobe et al. / Icarus 168 (2004) 223–236

Table 3 Extrasolar planetary systems in group B M∗ (M )

Mp sin i (MJup )

a (AU)

e

HZ (AU)

  10◦*1 (%)

  20◦*2 (%)

  40◦*3 (%)

Class

B1 Jupiter system Tau Boo HD130322 HD114386

1.00 1.30 0.79 0.75

1.00 3.870 1.080 0.990

5.203 0.046 0.088 1.620

0.048 0.018 0.048 0.280

0.70/1.30 1.25/2.70 0.50/1.00 0.50/1.00

13.14 11.74 12.74 20.14

11.41 3.99 10.81 3.89

0.00 0.00 0.00 1.35

II I I III

B2 70Vir HD223084 GJ3021 HD8574 HD150706 HD92788 HD142 HD177830 HD4203 HD210277 HD147513 HD141937 HD213240 HD190228 HD136118 HD50554 HD106252 HD33636 HD39091

1.10 1.05 0.90 1.10 0.98 1.06 1.10 1.17 1.06 0.99 0.92 0.16 1.22 1.30 1.24 1.10 1.05 0.99 1.10

6.600 1.210 3.320 2.230 1.000 3.800 1.360 1.280 1.640 1.240 1.000 9.700 4.500 4.990 11.900 4.900 6.810 7.800 0.370

0.430 0.440 0.490 0.760 0.820 0.940 0.980 1.000 1.090 1.097 1.260 1.520 2.030 2.310 2.335 2.380 2.610 2.700 3.340

0.400 0.480 0.505 0.400 0.380 0.360 0.370 0.430 0.530 0.450 0.520 0.410 0.450 0.430 0.366 0.420 0.540 0.410 0.620

0.85/1.60 0.75/1.40 0.60/1.20 0.85/1.60 0.70/1.30 0.75/1.40 0.85/1.60 0.93/1.80 0.75/1.40 0.70/1.30 0.65/1.25 0.70/1.30 0.93/1.80 1.25/2.70 1.10/2.25 0.85/1.60 0.75/1.40 0.70/1.30 0.85/1.60

100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.11 0.00 0.00 0.00 0.00 0.00 0.68 0.12 2.81

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

II–III II–III III III IV IV III–IV III–IV IV IV IV IV IV IV IV IV IV IV IV

System

based on , although the detailed values of the probability changes according to different choices of the ranges of D and a, I , and sin i. About half of the systems show small obliquity variations (
10 ◦ ) over their entire HZs (group A). This group is mostly comprised of systems containing giant planets with small eccentricities. A corresponding small I makes resonant regions narrow, resulting in small probability. About 1/4 of the systems show large changes of obliquity with more than 20◦ in their HZs (group C). This group is mostly comprised of systems containing giant planets with ∼ 1MJ and semimajor axis ∼ 0.1 AU, or planets with large eccentricities. In the systems containing giant planets with mp ∼ 1MJ and ap ∼ 0.1 AU, resonant regions considerably overlap HZs, as shown in Fig. 6a. In the systems with giant planets on eccentric orbits, resonant regions are extended by large inclinations which are inferred from the large eccentricities. All the systems with two giant planets except for the HD83443 system, and a system with three giant planets, the υ Andromedae system, belong to group C. The rest of them, which we call groups B1 and B2, show intermediate features between group A and group C. In the systems of group B1, resonant regions and HZs partially overlap. We also list the result for the system which only contains Jupiter around a host star with 1M except for massless terrestrial planets (we call this the “Jupiter system” hereafter). The Jupiter system belongs to this group. In group B2, resonant regions do not overlap with HZs, but the obliquity variation amplitude is not small enough: it is larger

than 10◦ but smaller than 20◦ over the entire HZs. The systems in group B2 contain giant planets with e ∼ 0.3–0.6, which causes large obliquity variations in nonresonant regions. Both group A and group C include all the orbital stability classes. As expected from the discussion in Section 3, a large amplitude of obliquity variation is not necessarily associated with an unstable orbit. In general, a giant planet that is closer to the HZ or massive destabilizes the orbits of terrestrial planets. However, obliquity is less affected by spin–orbit resonances if the giant planet is closer to the HZ or more massive. More than 1/3 of the systems in which terrestrial planets are unstable in the HZs (class III or IV) show small obliquity variations. On the other hand, nearly 1/3 of the systems in which the terrestrial planets have stable orbits in HZs (class I or II) show large obliquity variations. The systems that have both stable orbits and small obliquity variations are only 1/5 of all the systems. Note that 2/3 of such systems are systems with the giant planet(s) with ap  0.1 AU. In most of the multiplanet systems so far discovered,  of terrestrial planets in HZs would be more than 20◦ . More resonant regions exist in multiplanet systems than in single-planet systems, because spin–orbit resonances associated with ψ˙ fi occur in addition to the resonance with ψ˙ B. Moreover, secular resonance also occurs (see Section 2). In the HD74156 system, the spin–orbit resonance with an eigenfrequency fi is responsible for the obliquity variation  40◦ . In the 47UMa system, which contains gi-

Obliquity of extrasolar planets

233

Table 4 Extrasolar planetary systems in group C System

M∗ (M )

Mp sin i (MJup )

a (AU)

e

HZ (AU)

  10◦*1 (%)

  20◦*2 (%)

  40◦*3 (%)

Class

HD68988 HD217107 HD108147 HD6434 HD16141 HD114729 14Her HD216437 Epsilon Eridani HD80606 HD89744 HIP75458 HD222582 16CygB HD2039 Gliese 876

1.20 0.98 1.27 1.00 1.00 0.93 0.79 1.07 0.80 0.90 1.40 1.05 1.00 1.01 0.98 0.32

46.33 44.54 53.52 51.10 12.98 5.09 7.94 1.24 82.37 100.00 100.00 100.00 100.00 23.09 50.84 1.34

13.08 12.56 53.52 21.68 4.27 2.25 3.68 0.31 47.07 0.00 0.00 0.00 0.00 0.02 0.00 0.00

0.60/1.20

98.69

39.45

1.15

HD38529

1.39

1.40/3.00

100.00

99.98

2.60

HD74156

1.05

0.75/1.40

100.00

99.98

36.68

HD168443

1.01

0.70/1.30

100.00

71.83

1.33

HD82943

1.05

0.75/1.40

100.00

100.00

0.00

HD12661

1.07

0.80/1.45

100.00

100.00

0.00

HD160691

1.08

0.85/1.60

100.00

62.67

26.67

47UMa

1.03

0.70/1.30

79.10

38.39

18.50

υ Andromedae

1.30

0.140 0.140 0.498 0.300 0.280 0.330 0.326 0.340 0.608 0.927 0.700 0.710 0.710 0.670 0.690 0.100 0.270 0.100 0.690 0.290 0.360 0.649 0.395 0.529 0.228 0.540 0.410 0.350 0.200 0.310 0.800 0.096 0.100 0.012 0.280 0.270

64.45 65.18 53.52 90.31 57.51 55.49 60.09 63.38 83.62 100.00 100.00 100.00 100.00 100.00 100.00 80.00

0.91

0.071 0.070 0.104 0.150 0.350 2.080 2.500 2.700 3.30 0.439 0.880 1.340 1.350 1.700 2.200 0.210 0.130 0.540 2.500 0.129 3.680 0.276 3.470 0.290 2.850 0.730 1.160 0.830 2.560 1.500 2.300 2.100 3.730 0.059 0.829 2.530

0.93/1.80 0.70/1.30 0.75/1.40 0.70/1.30 0.70/1.30 0.65/1.25 0.50/1.00 0.80/1.45 0.50/1.00 0.60/1.20 1.40/3.00 0.75/1.40 0.70/1.30 0.70/1.30 0.70/1.30 0.10/0.20

HD37124

1.900 1.280 0.410 0.480 0.215 0.900 3.300 2.100 0.860 3.410 7.200 8.640 5.400 1.500 5.100 1.890 0.560 0.750 1.200 0.780 12.70 1.560 7.500 7.700 16.90 0.880 1.630 2.300 1.570 1.700 1.000 2.410 0.760 0.690 1.890 3.750

1.25/2.70

100.00

47.28

0.00

I I I I II III–IV III–IV III–IV III–IV II–III III IV IV IV IV IV III–IV III III–IV I IV II IV II IV II–IV IV III–IV III–IV IV IV III–IV II I II IV

ant planets relatively distant from the host star with small eccentricities, a secular resonance is located within its HZ, resulting in the large variation of obliquity. In the HD82943 and HD160691 systems, HZs are bounded by secular resonances. Only the HD83443 and 55Cnc systems show small obliquity variations (
10◦ ). In these systems, the orbits of terrestrial planets in HZs are also stable. It should be noted that multiplanet systems in group C contain giant planets with large eccentricities, which cause the large obliquity variations of terrestrial planets in HZs even in non-resonant regions. Current surveys of the extrasolar planets are limited to the detection of planets located within 3–4 AU from their host stars. Planets on more distant orbits with small eccentricities will be detected. As shown in Fig. 7, distant giant planets have a weaker effect on the obliquity of terrestrial planets in HZs than the giant planets at 2–4 AU. Therefore in future surveys we may find multiple giant planet systems having terrestrial planets with small obliquity variations in HZs.

5. Discussion In Section 4, we evaluated the probability of terrestrial planets in HZs having large obliquity variations in 87 known extrasolar planetary systems. The results suggest that the terrestrial planets with small obliquity variations in HZs may be uncommon. In this section, we discuss whether the resonances are avoidable. In future surveys, giant planets at  5–10 AU could be detected as additional planets in known planetary systems. As suggested by the results in Fig. 7, such distant giant planets have less influence on the HZs. The results shown in Section 4 do not suffer from significant modification if additional giant planets are detected in future. However, the obliquity variations of a terrestrial planet in an HZ may be affected by other terrestrial planets that could exist near or within the HZ. As Fig. 7 suggests, small planets with masses 1MJ can cause spin–orbit resonances. These effects are sometimes more severe than those caused by giant

234

K. Atobe et al. / Icarus 168 (2004) 223–236

planets. For example, in our solar system, terrestrial planets contribute to the change of Mars’ obliquity, as do Jupiter and Saturn (Ward, 1982). The effect of a satellite(s) cannot be neglected for the obliquity variation, if the satellite is sufficiently massive 3 (m is the and/or close to the terrestrial planet so that ms /aps s satellite mass, and aps is the distance between the planet and the satellite) is  mc /a 3 (e.g., Eq. (2)). The present Earth owes the stability of its spin axis to the Moon (e.g., Ward, 1982; Laskar et al., 1993; Laskar, 1996). Because of the lunar torque on the Earth’s equatorial bulge which is twice as large as that of the Sun, the Earth’s precession constant is increased by a factor of 3. This is equivalent to reducing the rotation period by a factor of 3 in the results in previous sections. As a result of this effect, Earth is out of resonant regions. In the past, the Moon was closer to the Earth, so that the lunar torque was stronger than the current one. Thus, the Moon has stabilized the Earth’s obliquity throughout the history of orbital evolution of the Moon. Orbital separation between a terrestrial planet and other planets must be larger than a critical value to maintain the orbital stability of the terrestrial planet, which imposes an upper limit on orbital precession frequencies of the terrestrial planet. A giant planet would marginally maintain the orbital stability of a terrestrial planet in an HZ at least during the lifetime of the host star, if the terrestrial planet is separated from the giant planet by several RHill (RHill is the Hill radius of the giant planet). When the giant planet is separated from the HZ by 5RHill , the orbital precession frequency of the terrestrial planet at 1 AU is 50–100 yr−1 . Hence, 50–100 yr−1 may be the upper limit of the orbital precession frequency. Note that this value is almost of the same order as that of the spin precession frequency of the terrestrial planet at 1 AU that is subject to the torque caused by the central star. If the torque produced by a satellite increases the spin precession frequency over the upper limit of the orbital precession frequency, spin–orbit resonance does not occur. 3  m /a 3 ) satellite would A sufficiently large/close (ms /aps c generally tend to stabilize the obliquity of its host planet. Not only the obliquity of terrestrial planets, but also that of giant planets is affected by spin–orbit resonance. Saturn’s primordial obliquity may have been smaller and because modified to the present value (∼ 26◦ ) by spin– orbit resonance with Neptune (Hamilton and Ward, 2002a, 2002b, 2002c). We also estimated the obliquity variations of giant planets in the known multiple-planet systems assuming rotation periods similar to Jupiter or Saturn (∼ 10 h), and found that no giant planets in multiple-planet systems found so far would be in spin–orbit resonances. As mentioned in Section 2, the obliquity of planets in a retrograde spin state is not affected by the spin–orbit resonance. In the Solar System, however, three terrestrial planets have prograde spins. The primordial spin of Venus, which currently has a retrograde spin, could have been prograde until a large obliquity variation and the dissipative effects of gravitational and atmospheric tides and core-mantle friction

overturned the spin (Correia and Laskar, 2001, 2003). However, it is possible that planets are formed in the retrograde spin state and maintain their primordial obliquity. Planetary spin caused by planetesimal accretion has two components: the systematic component produced by accumulation of a large number of small planetesimals, and the random component produced by a small number of large impacts (e.g., Lissauer et al., 2000, and references therein). The systematic component tends to result in a relatively slow prograde spin. On the other hand, the random component can produce a relatively rapid spin and does not have any preferred spin direction. For the terrestrial planets, it is expected that the random component is dominant (e.g., Lissauer et al., 2000). If terrestrial planets could have formed in a retrograde spin state, they would maintain their primordial obliquity without large variations. As shown in Section 5, among the known extrasolar planetary systems, those with short-period giant planets with semimajor axes  0.1 AU are favorable sites for terrestrial planets to exist with orbital stability and small obliquity variations in HZs. However, the short-period giant planets may have been formed in the outer regions beyond the HZs and migrated inward to their present locations by disk-planet interactions (Lin et al., 1996). It would not be easy for terrestrial planets in the HZs to survive during the passage of the giant planets. If terrestrial planets could have form in an HZ from a residual protoplanetary disk after the passage, such terrestrial planets would stay in stable orbits with small obliquity variations.

6. Conclusion We have investigated obliquity variations of possible terrestrial planets in HZs in extrasolar planetary systems. The amplitude of obliquity variations may be one of the most important factors for survival of land-based life on terrestrial planets. Using the analytical expressions for the obliquity variations of terrestrial planets, we statistically studied general properties of obliquity variations of terrestrial planets in an HZ perturbed by a giant planet(s). We derived the existing probability of terrestrial planets with obliquity variation amplitude larger than given critical values. The results are as follows: (1) Terrestrial planets in the HZs are most affected by spin– orbit resonance when the giant planets are separated from the HZ so that the orbits of terrestrial planets in the HZ are marginally stable. (2) About half of the known extrasolar planetary systems can have hypothetical terrestrial planets with small obliquity variations (
10 ◦ ) over the entire HZs. (3) The systems that can have hypothetical terrestrial planets with both small obliquity variations and stable orbits in the HZs are only 1/5 of known extrasolar planetary

Obliquity of extrasolar planets

systems. Since most of them are comprised of systems with close-in giant planets, actual existence of such terrestrial planets may be not be common. (4) In general, more planets (giant planets or other terrestrial planets) result in larger obliquity variations of terrestrial planets. Almost all of the known multiple-planet systems (except for HD83443 and 55Cnc systems) show large obliquity variations (  20◦ ). (5) The precession by the torque of a large/close satellite tends to reduce the obliquity variation of terrestrial planets in an HZ. Planets’ retrograde spin removes spin– orbit resonance. Both effects can increase the possibility of the existence of terrestrial planets with stable obliquity and orbits. Masses and orbital radii of the most effective giant planets for enhancing obliquity variations are plotted in Fig. 7. For example, in the system including a giant planet with 1MJ and with a HZ of 0.7–1.3 AU, the obliquity of terrestrial planets in the HZ is most radically changed (  20◦ ) if the giant planet orbits at ∼ 0.05–0.1 or 3–5 AU. The orbits of the terrestrial planets are marginally stable (they are destabilized if the giant planet is closer to the HZ or is more massive). In systems where planetary orbital migration and damping of eccentricity have not occurred significantly, orbital configuration of planets may be in a marginally stable state. In these situations, our results suggest that obliquity of terrestrial planets tends to be near the spin–orbit resonances, and often have large variations on timescales of 105 –106 years, as long as the spin of the terrestrial planets is prograde. Therefore, retrograde spin of the terrestrial planets or the existence of a large/close satellite might be essential for land-based life to survive on a terrestrial planet in an HZ.

235

where n is the mean motion of the terrestrial planet, Fj =

1 (1) αj b3/2 (αj ), 8a ∗

(A.3)

a ∗ ≡ max(a, aj ), αj ≡ min(a, aj )/max(a, aj ), a and aj are the semimajor axis of the terrestrial planet and the giant planet j , bs(m)(αj ) are Laplace coefficients (e.g., Brouwer and Clemence, 1961), Ng is number of giant planets, and µj is the mass of giant planet j normalized by the host star mass, mj /mc . Ij i are eigenvectors for fi . We expand Eq. (1) in terms of I 1, in the Cartesian coordinates where the invariant plane is on the x–y plane (see Fig. 1). In these coordinates, the spin vector sˆ and the orbit normal nˆ are given by     sin θ sin ψ sin I sin Ω sˆ = − sin θ cos ψ ; nˆ = − sin I cos Ω . (A.4) cos θ cos I For θ , ψ, I and Ω, see Fig. 1. Substituting Eq. (A.4) into the z component of Eq. (1), we obtain θ˙ α cos I (sin ψ cos Ω − cos ψ sin Ω) + O(I 2 ),

(A.5)

with Eqs. (3) and (4), the above equation read as   Ne  Ii sin(ψ − Ωi ) θ˙ α cos  I0 sin(ψ − Ω0 ) + i=1

+ O(I ), 2

(A.6)

where Ω0 = Bt + γ and Ωi = fi t + γi . From the x and the y components of Eq. (1), we obtain   ψ˙ α cos  −1 + I cot θ (cos ψ cos Ω + sin ψ sin Ω) + O(I 2 )

Acknowledgments We thank William Ward for useful comments concerning the discussion about the known extrasolar planetary systems. We also thank Hidekazu Tanaka for valuable advice on the analytical formulation. The referees also suggested directions which bettered the quality of this paper a great deal.

Appendix A. Derivation of Eqs. (5) and (6) Spin axis sˆ of a terrestrial planet precesses according to Eq. (1). Free oscillation frequency B and forced inclination Ii of the orbital motion of the terrestrial planet in Eqs. (3) and (4) are given by B = −2n

Ng 

µj aFj ,

(A.1)

j =1

Ii =

2n

Ng

j =1 µj aFj Ij i

B − fi

,

(A.2)

−α cos  + α cos  cot θ   Ne  Ii cos(ψ − Ωi ) + O(I 2 ), × I0 cos(ψ − Ω0 ) + i=1

(A.7) where Eqs. (3) and (4) are used in the last equation. Since cos  = sˆ · nˆ = cos θ cos I + sin θ sin I cos(ψ − Ω), d(cos ) = −θ˙ sin θ cos I + θ˙ cos θ sin I cos(ψ − Ω) dt d − (ψ − Ω) sin θ sin I sin(ψ − Ω). (A.8) dt Substituting Eqs. (A.5) and (A.7) into the above equation, and neglecting terms higher than the second order of I , we finally obtain d(cos ) = I0 B sin θ sin(ψ − Ω0 ) dt Ne  Ii fi sin θ sin(ψ − Ωi ) + O(I 2 ). + i=1

(A.9)

236

K. Atobe et al. / Icarus 168 (2004) 223–236

References Abe, Y., 1993. Physical state of very early Earth. Lithos 30, 223–235. Berger, A.L., 1978. Long-term variations of daily insolation and Quaternary climatic changes. J. Atmos. Sci. 35, 2362–2367. Berger, A. (Ed.), 1989. Climate and Geo-Science—A Challenge for Science and Society in the 21st Century. Kluwer Academic, Dordrecht. Berger, A.L., Imbrie, J., Hays, J., Kukla, G., Saltzman, B. (Eds.), 1984. Milankovitch and Climate—Understanding the Response to Astronomical Forcing. Reidel, Norwell, MA. Brouwer, D., Clemence, G.M., 1961. Methods of Celestial Mechanics. Academic Press, New York. Colombo, G., 1966. Cassini’s second and third laws. Astron. J. 71, 891–896. Correia, A.C.M., Laskar, J., 2001. The four final rotation states of Venus. Nature 411, 767–770. Correia, A.C.M., Laskar, J., 2003. Long-term evolution of the spin of Venus. II. Numerical simulations. Icarus 163, 24–45. Crowley, J.T., North, G.R., 1991. Paleoclimatology. Oxford Univ. Press, New York. Gehman, C.S., Adams, F.C., Laughlin, G., 1996. The prospects for Earthlike planets within known extrasolar planetary systems. Publ. Astron. Soc. Pac. 108, 1018–1023. Hamilton, D.P., Ward, W.R., 2002a. The obliquities of the giant planets. Bull. Am. Astron. Soc. 34, 944. Hamilton, D.P., Ward, W.R., 2002b. Obliquity of Saturn: analytical model. Bull. Am. Astron. Soc. 34, 890. Hamilton, D.P., Ward, W.R., 2002c. Obliquity of Saturn: numerical model. Bull. Am. Astron. Soc. 34, 2809. Abstract. Henrard, J., Murigande, C., 1987. Colombo’s top. Celest. Mech. Dyn. Astron. 40, 345–366. Ida, S., Kokubo, E., Makino, J., 1993. The origin of anisotropic velocity dispersion of particles in a disc potential. Mon. Not. R. Astron. Soc. 108, 875–889. Ito, T., Masuda, K., Hamano, Y., Matsui, T., 1995. Climate friction: a possible cause for secular drift of the Earth’s obliquity. J. Geophys. Res. 100, 15147–15161. Jones, B.W., Sleep, P.N., 2002. The stability of the orbits of Earth-mass planets in the habitable zone of 47 Ursae Majoris. Astron. Astrophys. 393, 1015–1026. Jones, B.W., Sleep, P.N., Chambers, J.E., 2001. The stability of the orbits of terrestrial planets in the habitable zones of known exoplanetary systems. Astron. Astrophys. 366, 254–262. Kasting, J.F., Whitmire, D.P., Reynolds, R.T., 1993. Habitable zones around main sequence stars. Icarus 101, 108–128. Kinoshita, H., 1977. Theory of the rotation of the rigid Earth. Celest. Mech. 15, 277–326. Laskar, J., 1996. Large scale chaos and marginal stability in the Solar System. Celest. Mech. Dyn. Astron. 64, 115–162.

Laskar, J., Robutel, P., 1993. The chaotic obliquity of the planets. Nature 361, 608–612. Laskar, J., Joutel, F., Robutel, P., 1993. Stabilization of the Earth’s obliquity by the Moon. Nature 361, 615–617. Lin, D.N.C., Bodenheimer, P., Richardson, D.C., 1996. Orbital migration of the planetary companion of 51 Pegasi to its present location. Nature 380, 606–607. Lissauer, J.J., Dones, L., Ohtsuki, K., 2000. Origin and evolution of terrestrial planet rotation. In: Canup, R.M., Righter, K. (Eds.), Origin of Earth and Moon. Univ. of Arizona Press, Tucson, pp. 101–112. Marcy, G.W., Cochran, W.D., Mayor, M., 2000. Extrasolar planets around main-sequence stars. In: Mannings, V., Boss, A., Russell, S.S. (Eds.), Protostars and Planets, fourth ed. Univ. of Arizona Press, Tucson, pp. 1285–1311. Marzari, F., Weidenschilling, S.J., 2002. Eccentric extrasolar planets: the jumping Jupiter model. Icarus 156, 570–579. Menou, K., Tabachnik, S., 2003. Dynamical habitability of known extrasolar planetary systems. Astrophys. J. 583, 473–488. Milankovitch, M., 1941. Kanon der Erdbestrahlung und seine Anwendung auf das Eiszeitproblem. In: Königlich Serbische Academie Publication, vol. 133. Königlich Serbische Academie, Belgrad, 633 p. Murray, C.D., Dermott, S.F., 1999. Solar System Dynamics. Cambridge Univ. Press, Cambridge. Noble, M., Musielak, Z.E., Cuntz, M., 2002. Orbital stability of terrestrial planets inside the habitable zones of extrasolar planetary systems. Astrophys. J. 572, 1024–1030. Sharaf, S.G., Budnikova, N.A., 1969a. On secular perturbations in the elements of the Earth’s orbit and their influence on the climates in the geological past. Bull. Inst. Theor. Astron. 11, 231–261. In Russian. Sharaf, S.G., Budnikova, N.A., 1969b. Secular perturbations in the elements of the Earth’s orbit and the astronomical theory of climate variations. Trudy Inst. Theor. Astron. 14, 48–84. In Russian. Touma, J., Wisdom, J., 1993. The chaotic obliquity of Mars. Science 259, 1294–1297. Vernekar, A.D., 1972. Long-period global variations of incoming solar radiation. Meteor. Monogr. 12, 1–21. Ward, W.R., 1973. Large-scale variations in the obliquity of Mars. Science 181, 260–262. Ward, W.R., 1974. Climate variations on Mars 1. Astronomical theory of insolation. J. Geophys. Res. 79, 3375–3386. Ward, W.R., 1979. Present obliquity oscillations of Mars—fourth-order accuracy in orbital E and I. J. Geophys. Res. 84, 237–241. Ward, W.R., 1982. Comments on the long-term stability of the Earth’s obliquity. Icarus 50, 444–448. Ward, W.R., Agnor, C.B., Canup, R.M., 2002. Obliquity variations in planetary systems. Proc. Lunar Planet. Sci. Conf. 33rd. Abstract 2017.