Effects of Extreme Obliquity Variations on the Habitability of Exoplanets

1. Introduction

The habitability of a world depends on a host of properties, from observable quantities like its mass and distance from the parent star to those that are difficult if not impossible to measure: atmospheric composition, surface reflectivity, ice, water distribution, and so on. In the case of stars as massive as our Sun, detecting Earth-mass planets in any orbit is difficult with modern technology. In the last decade, attention has turned primarily to the discovery of rocky planets orbiting in the habitable zone (HZ), a shell around a luminous object in which an Earth-like planet could support liquid water on its surface (Dole, 1964; Kasting et al., 1993; Kopparapu et al., 2013), as these worlds are best suited for the development and sustainment of life as we know it. In its latest revision (Kopparapu et al., 2013), the HZ is calculated for a highly idealized case, in which many properties of the star, planet, and planetary system are ignored. Following the identification of possible processes that can impact habitability (Heller and Armstrong, 2014), we explore how gravitational perturbations from additional planets can affect the climate. We find that, in many cases, these perturbations can extend the outer edge of the HZ, thereby increasing the number of planets in the Galaxy that are potentially habitable.

While the vast majority of work on habitability has used a replica of Earth to determine orbits that are potentially habitable, there are notable exceptions. Some studies have explored the habitability of synchronously rotating planets (Joshi et al., 1997; Joshi, 2003; Edson et al., 2011; Pierrehumbert, 2011; Wordsworth et al., 2011; Yang et al., 2013). Others (Abe et al., 2011; Zsom et al., 2013) considered planets that are much drier than Earth. Several studies (Williams and Pollard, 2002, 2003; Dressing et al., 2010; Spiegel et al., 2010) varied the eccentricity and obliquity of an Earth-like planet and found that larger values tend to increase the globally averaged temperature on a planet, while holding the semimajor axis constant.

While these studies made great strides in understanding Earth's climate sensitivity to rotation rate, obliquity, and eccentricity, aside from Spiegel et al. (2010), they largely ignored that the latter two properties evolve with time due to gravitational perturbations from other bodies. Earth maintains a relatively constant axial tilt either due to the presence of the Moon, as suggested by Laskar et al. (1993), or due to the inherent stability of Earth's axis, as indicated by Lissauer et al. (2012). However, as illustrated by Williams (1998a, 1998b), changes in the architecture of our solar system—such as moving Jupiter inward—can result in dramatic variations in the obliquity of Earth even with the presence of a Moon. Still, it has been suggested that the relatively small variations in Earth's obliquity result in a stable climate conducive to the development of life. Adding to this stability is the fact that the orbital eccentricity remains smaller than about 0.05 due to the approximately circular orbits of the large planets of the Solar System.

It is possible that small obliquities and circular orbits are not a requirement for habitability. Williams and Pollard (2003) used general circulation models (GCMs) to determine how different obliquity variations affect Earth's climate. They found that Earth-like planets with high obliquities were no more likely to experience extreme runaway greenhouse or snowball Earth events, making them just as habitable as Earth. Later, Spiegel et al. (2010) examined how large eccentricity oscillations affect the climates of rocky exoplanets. They found that, in some cases, planets could break out of a snowball event during periods of high eccentricity. It was the goal of the current study to build on these previous results and explore self-consistent models of the climates of planets that experience rapid, large-amplitude, and possibly chaotic oscillations of eccentricity and obliquity.

Orbit-induced seasonal effects like the ice-albedo feedback determine the limit of the outer edge of the HZ. As the surface temperature drops, volatile ices such as CO₂ and water can condense on the surface. The high albedos of their solid phases inhibit a planet's ability to absorb solar radiation, which reduces the temperature further. Mars, if it possessed sufficient surface water, would have been in danger of falling prey to these snowball episodes. Geological evidence exists for these episodes in Earth's past (Hoffman et al., 1998). For Earth, a dynamic CO₂ recycling system works to offset the negative effects of these events on timescales of millions of years (Walker et al., 1981). As the planet cools, weathering rates slow down and lock CO₂ in the atmosphere. As the CO₂ builds up, the greenhouse effect increases, eventually melting the ice. On early Mars, as now, such robust CO₂ cycling could not have been enough to resurrect the planet from these snowball events. From seasonally resolved modeling, it seems that other factors, in particular orbital and obliquity variations, would have to play a role (Armstrong et al., 2004). Mars' obliquity has probably undergone significant evolution in the past (Laskar and Robutel, 1993) due to gravitational torques by the other planets in the Solar System, a property that may have permitted at least episodic liquid water at the surface.

With this Solar System context in mind, we turn our attention to predicting the habitability of exoplanets. We do not anticipate the photometric and spectroscopic data needed to characterize exoplanets until the launch of the James Webb Space Telescope in ∼2017. Spacecraft like Kepler and ground-based radial velocity surveys are already returning data that can pin down, through computational analysis, important orbital parameters that impact the climate, such as eccentricity, timing of perihelion passage, and the evolution of the spin axis of the planet. Through the coupling of these data to N-body simulations and simple, fast climate codes, a more comprehensive picture of the HZ of a system can be obtained. This type of analysis can be used to prioritize future characterization observations, which are likely to be challenging, expensive, and based on precious little information. In the present study, we found that planetary system architecture, that is, the distribution of masses and orbits of the other planets in the system, can play a significant role in defining the extent of the HZ.

The rotational evolution of the bodies in our solar system can be accurately modeled because the masses and orbits of the planets are extremely well measured. For exoplanets, the situation is more difficult. Radial velocity surveys (e.g., Butler et al., 2006) are only able to place a lower bound on mass and cannot measure the relative inclination between orbital planes. Kepler can constrain inclinations (Fabrycky et al., 2012) but is heavily biased toward the discovery of planets in coplanar configurations. While these systems are analogous to our solar system, we note that for the one system for which astrometry (which is not biased toward any particular inclination) has measured a mutual inclination, υ Andromedae (McArthur et al., 2010), the relative inclination is 30 degrees. Moreover, studies that predict the large eccentricities of exoplanetary orbits simultaneously predict large inclinations (Marzari and Weidenschilling, 2002; Chatterjee et al., 2008; Raymond et al., 2010; Barnes et al., 2011). It is possible, perhaps likely, that there exists a population of planetary systems with large inclinations and a potentially habitable planet. These architectures will induce much larger changes in orbital inclination, which in turn induces large obliquity oscillations. The Gaia mission may be able to determine the range of architectures for giant planets (Casertano et al., 2008; Sozzetti et al., 2013).

As no rocky planet is currently known to orbit with sibling planets with high mutual inclinations, we explored the phenomena with 17 hypothetical, dynamically stable systems. We found that there is a direct link between the orbital architecture of a planetary system and the possible range of climate conditions on a potentially habitable planet. We used these models to constrain the orbital conditions of a hypothetical planet and found that orbital and rotational evolution tends to push the outer edge of the HZ out, relative to planets where no evolution occurs.

Below, we outline a model that links physically realistic orbital architectures to the spin evolution of a hypothetical Earth-like planet and finally to its climate. In Section 2, we discuss the motivation behind the systems we have modeled in an effort to obtain a set that spans a range of orbital elements. In Section 3, we outline the model used to evolve the spin axis of the planet. In Section 4, we present a simplified energy balance model (EBM) designed to be robust across wide variations of these orbital changes and fast enough that million-year integrations require only a few hundred seconds of computational time. We then present the results of these models in Section 5 followed by a discussion in Section 6.

2. Orbital Simulations

To make a first assessment of the potential of similar architectures to support habitable worlds, we created 17 hypothetical systems with moderate inclinations, always initially including an Earth-like planet on a circular orbit 1 AU from a solar-mass star. The orbital architectures are arbitrary but consistent with the distribution of orbital elements of known planets. While the potentially habitable world is always the same, its siblings have a wide array of properties. Systems consist of 2–3 planets, with eccentricities ranging from 0 to 0.3 and mutual inclinations from 10 to 30 degrees. The orbital properties of these cases, in astrocentric coordinates, are presented in Table 1.

These systems were selected after careful consideration of orbital stability. We numerically integrated each case for 100 Myr in “hybrid” mode with Mercury (Chambers, 1999) and confirmed that the evolution of every orbital element appeared periodic. Additionally, energy was conserved to better than 1 part in 10⁶, which is sufficient for numerical accuracy (Barnes and Quinn, 2004). We then rotated each system such that the reference plane corresponds to the fundamental plane. For those systems with super-Jupiter planets, this conversion can change some orbital elements significantly. We then reran each case at very high resolution for ∼1 Myr, conserving orbital angular momentum to 1 part in 10¹². We are therefore confident that no numerical inaccuracies are propagated into the rotational calculations.

3. Obliquity Modeling

Using the results from the orbital runs described above, we employed the obliquity model of Laskar (1986a,b) as used in previous orbit coupled modeling by Armstrong et al. (2004). There are two primary factors that influence the evolution of the obliquity. First, variations in the geometry of the orbit that are governed by the overall system architecture are present in the orbital elements derived from the N-body simulations. The gravitational influence of the other massive bodies in the system affects the eccentricity, e, the inclination, i, the argument of perihelion, ω, and the longitude of the acceding node, Ω, of the Earth-sized planet, here cast in terms of the eccentricity-inclination variables

where is longitude of perihelion. Changes in parameters like the inclination result in changes in the spin axis orientation relative to the fundamental plane, that is, the obliquity, as illustrated in Fig. 1. In other words, the rotational angular momentum is decoupled from the orbital angular momentum.

In addition to these geometric factors, the direct torques from the central star are included as a term, R(ψ), in the precession, p_A, and obliquity, ψ, evolution equations,

with

Here, the semimajor axis, a, is measured in AU; the planet's angular velocity, ν, is measured in rad day⁻¹; M=1.0 in solar units; k²=GM/4π², with G, the gravitational constant, measured in units of AU³ M⁻¹ day²; and E_D is the dynamical ellipticity, a measure of the non-sphericity of the planet. Finally, the relativistic precession is accounted for with

where the value of κ_r is

where c is the speed of light, M_p is the planet's mass in units of solar masses, and n and a are related via Kepler's law

as outlined in Laskar (1986a, 1986b). The values used for these model parameters are listed in Table 2.

Table 2.

Parameters for the Obliquity Calculation

Parameter	Symbol	Value	Units
Gravity constant	k	0.01720209895
Angular spin rate	ν	2π	rad day−1
Dynamical ellipticity	ED	0.00328005	unitless
Initial obliquity	ψ0	23.44	degrees
Initial precession	pA	0.0	degrees
Planet mass	Mp	1.0	Earth masses

Open in a new tab

For each of the 17 orbital runs, the Earth-mass planet in the simulation has a 1-day rotation with an Earth-like dynamical ellipticity and orbits a Sun-like star. The obliquity, precession angle, and rates of change in these parameters are computed each time step from the input orbital parameters, and those results are used as the initial conditions for the following time step in which a fourth-order Runge-Kutta integrator is used. These results are then available for further coupling to the EBM described below. The primary limitation to the method is that only direct torques from the central body are included. Also, there is an implicit assumption that the planet is rapidly rotating and that it can be accurately described by the dynamical ellipticity relative to the central body. Since none of our models involve tidally locked worlds or close-passing planets or moons, these assumptions are justified.

Of the 17 systems studied, we selected seven to analyze that represented the full spectrum of outcomes from the models:

• (1) Earth-Jupiter-Saturn system for comparison (System 1, Baseline—Fig. 2).

• (2) Two systems with high mean obliquity but modest variations (Systems 2 and 3—Figs. 3 and 4).

• (3) Two systems with wide and rapid variations in obliquity (Systems 4 and 5—Figs. 5 and 6).

• (4) Two systems with wide and slow variations in obliquity (Systems 6 and 7—Figs. 7 and 8).

Figure 2 illustrates the baseline run with (Moon-less) Earth, Jupiter, and Saturn. The panels on the left represent the relevant orbital parameters used in the obliquity calculations, and the right-hand panels are the parameters of interest for the climate calculations. These solutions are qualitatively similar to recent work (e.g., Lissauer et al., 2012) but differ in the detailed periodicities and magnitudes of the obliquity and precession due to the specific orbital solution used in the model. However, the climatically important parameters have relatively low amplitude and are slowly varying compared to other models.

Figure 3 shows a model with the Earth-mass planet in a system with two other planets 10 times its size, with large eccentricity variations and high mutual inclination. This architecture results in wide swings in obliquity from the starting value to about 80 degrees but has relatively slow variations.

Figure 4 shows a model with two 10 Jupiter-mass planets orbiting in the system with the Earth-mass planet, again with high eccentricity, high mutual inclinations, and an obliquity that oscillates around 85 degrees with a period of approximately 200,000 years.

Figures 5 and 6 illustrate some of the cases where the obliquity variations become rapid (periods of <16,000 years) at relatively high amplitudes between 10 and 60 degrees. In each of these cases, there is high mutual inclination and large eccentricities due to two 10 Jupiter-mass planets in Fig. 4 and one Jupiter-mass planet in Fig. 5.

Lastly, Figs. 7 and 8 show an Earth-like planet that has a Jupiter-mass planet in the system with a mutual inclination of 10 degrees (Fig. 7) and 30 degrees (Fig. 8), resulting in a wide range in obliquity that varies more slowly than the previous cases.

The precession rates illustrated in the figures are positive, as expected, except near ψ=0 or ψ>90 degrees. These are numerical artifacts of the geometry of the semianalytical model. As the obliquity goes beyond 90 degrees, the precession rate reverses due to the fact that the “south” pole is now the “north” pole. Additionally, the cotangent term in Eq. 5 makes the precession poorly defined near 0 degrees. From a climate perspective, the precession is unimportant at low obliquities, and this is only relevant for “bookkeeping.”

4. A Simplified Energy Balance Model

With the orbital variations and obliquity calculations in hand, we assessed the surface conditions using a simplified EBM. This calculation requires a model that can simulate planets on the timescales for glacier growth/retreat: 10⁴ to 10⁵ years. GCMs are generally too complex to be practical on these timescales. And both GCMs and detailed EBMs, to a greater or lesser extent, require detailed information about the planet or assumptions based on relevant observations. In the case of exoplanets, especially those with widely varying orbital parameters, making these assumptions is difficult if not impossible. To include enough detailed physics to allow for direct calculation of climate effects in absence of observation would render our code too cumbersome. To act as a first-order comparison between the other systems and the baseline model, with as few input parameters as possible, we have developed a simplified EBM to examine the conditions on the surface. While the lack of detail makes specific climate predictions difficult for individual planets, general comparisons should be valid.

We modeled the surface as a one-dimensional latitude grid of 90 bands from −90 S to +90 N and the atmosphere as a single slab gray absorbing layer. We modified a simplified model outlined by Armstrong et al. (2004) to include basic deposition and evaporation of water ice over the seasonal cycle.

For each time step in the million-year orbital-obliquity model, we modeled three years of the seasonally resolved climate, one year to act as a spin-up from the initial conditions, and two years to compute global surface properties. Each model year starts with t=0 at the time of perihelion passage, T_peri=0. From this, we can compute the eccentric anomaly, E, from the implicit equation

using Newton's method (where all the values above are measured in MKS units). The eccentric anomaly is related to the true anomaly, f, by

With the true anomaly, we compute the instantaneous distance to the star,

The instantaneous stellar distance allows us to compute the incident stellar radiation,

where L_* is the stellar luminosity. From this, we can compute the daily mean top-of-atmosphere instellation at any point on the globe,

where η is the half-angle of daylight—a measure of the length of the day—at a given latitude δ; and the substellar latitude, δ_*=ψ sin(L_s), is determined by the obliquity, ψ, and the stellar longitude, L_s=L_sp+f, where L_s=0 is the northern hemisphere vernal equinox. L_sp=L_sp0 – p_A – ω is the solar longitude of perihelion, determined by the arbitrary constant L_sp0 (in our case, 90 degrees represents northern hemisphere summer solstice), the spin precession, p_A, and argument of perihelion, ω. The half-angle of daylight, η, is given by

To estimate the surface temperature, we first determine the energy balance at the surface at each latitude bin,

where A is the surface albedo of either ground or ice depending on local conditions, ɛ_s is the surface emissivity, T_e is the planet's atmosphere-free equilibrium temperature, and F_surf is the heat flux from the surface. The atmosphere is modeled as a single slab with opacity τ that is equivalent to the number of absorbing layers required to achieve a surface temperature of T_s, related to T_e by

The τ parameter is essentially a one-dimensional model of the greenhouse effect. That is, adjusting τ as a free parameter allows us to include a contribution to greenhouse warming, with τ=0.095 reproducing the Earth-like global surface temperature in our baseline model. To estimate the surface temperature, we set ΔE=0, solve first for the equilibrium temperature, and then solve for the surface temperature.

Ice deposition is handled parametrically in the model by choosing a global ice deposition rate that comes into effect when the surface temperatures dip below 273 K. At this point, the surface albedo is set to the value for ice, and ice is allowed to accumulate. If the surface temperature exceeds 273 K in the presence of ice, the surface temperature is held at 273 while the deposits are evaporated according to the difference between what the temperature of the surface would be in the absence of ice and the freezing point of water, T_ice=273. In this case, the mass loss rate in units of kg s⁻¹ m⁻² is given by

where L_h=3.34×10⁵ is the latent heat of fusion for ice in J kg⁻¹. Once the ice is gone, the albedo is reset to the surface value, and the surface temperature evolves normally. The ice deposition rate is set such that the thickness and extent of the ice-snow regions are roughly comparable with values for present-day Earth in the baseline model.

In summary, the only “parameters” are the albedo of land/ice, the atmospheric opacity (which is varied as a free parameter in the study), the ice deposition rate, and a parameterized version of thermal exchange with the surface. To calibrate the model, we chose reasonable values for the albedo, tuned the deposition rate to get reasonable ice caps, and modified the atmospheric opacity and surface flux to hit a global mean temperature of ∼288 K. The parameters of the model are outlined in Table 3.

This approach has several important features:

• (1) The model is extremely fast, allowing us to run thousands of cases over millions of years to explore the full parameter space of the model.

• (2) The physics is extremely simplified and coupled tightly to the orbital parameters. Therefore, we can see the first-order effects of the orbital variations, the main goal of this study.

• (3) The model is determined by only three free parameters: the opacity (or number of absorbing layers), the surface heat flux, and the ice deposition rate.

• (4) The final results are normalized to the baseline run, which, while limiting what we can say about the specific climate, allows for a quantitative comparison of the outer edge of the HZ for the hypothetical planets and modern Earth.

This model reproduces the first-order climate effects and the ice-albedo feedback (see Fig. 9). Each panel plots the two model years after the spin-up year. The top plot is the baseline run, showing the surface temperatures for our nominal “Earth.” The discontinuities result from the abrupt change in surface albedo. The second model is identical to the first, except the eccentricity has been increased to 0.15 to show the asymmetric effect on the climate, in this case creating a more persistent polar cap in the southern hemisphere. The third model is the same as the first but shows the effect of increasing the obliquity to 90 degrees.

Using this model, we computed approximately 4400 test cases for the 17 models over a range of semimajor axes from 0.80 to 3.0 and opacities ranging from 0.095, obtained from tuning the Earth baseline model to 288 K, to 0.26, nearly 3 times the amount of absorbing material in the atmosphere. For each time step in the orbital model, we ran an instance of the EBM, computing the three-year climate run, and averaged the last two years to compute the global mean temperature. In addition, we computed two additional measures of the habitability of the surface, the Temperature Habitability Index (THI), and the Ice Habitability Index (IHI). THI is defined here as the time-averaged fraction of the surface that is between 273 and 373 K. IHI is defined as the time-averaged ice-free regions of the planet regardless of surface temperature. We then normalized the quantities to the THI and IHI for the t=0 baseline Earth, which has a THI=0.766 and IHI=0.485. If the THI (or IHI) is greater than 1.0, this planet is “more habitable” than the baseline Earth. If the THI or IHI is less than 1.0, it is less habitable. A THI of 0 indicates a world that is either frozen (in which case the IHI will also be 0) or has no temperatures below 373 K. The maximum values for THI and IHI by this metric are 1.3 and 2.1, respectively.

This method has some distinct advantages for exploring the effects of the orbit on climate. Since the model is seasonally resolved, even some planets that might have an annually averaged global mean surface temperature below freezing may still experience significant periods of time in the “habitable” range and, due to the orbital effects, avoid a planetwide snowball.

5. Results

Figure 10 shows the aggregate results for the seven systems studied in detail. The first column is the THI, the second column is the IHI, and the third column is the global mean temperature, all as a function of both the distance from the host star and the opacity of the atmosphere for Systems 1–7. Since our simplified EBM ignores complications due to the runaway greenhouse, we will restrict our discussion to the outer edge of the HZ.

The baseline model shows that the outer edge of the HZ is approximately 1.4 AU for Earth at τ=0.095. Each of the other high-obliquity systems and/or high-eccentricity systems shows a systematic—and sometimes dramatic—increase in the outer edge of the HZ. Those increases are tabulated in Table 4. The maximum increase is System 3, which has both high eccentricity and large obliquity, with large variations in both parameters.

Table 4.

Summary of the Systems Used in the Complete Climate Comparison

N	ē	Pe, Myr	, deg	Pψ, Myr	lout,AU	lstatic,AU	ES, %	EV, %
1	0.05±0.023	0.25	22.6±1.1	0.17	1.4	1.4	n/a	n/a
2	0.14±0.042	0.20	60.6±20.2	0.33	2.1	2.0	50	5
3	0.33±0.025	0.04	78.6±18.7	0.17	2.7	2.5	93	8
4	0.08±0.036	0.02	40.7±12.1	0.02	1.9	1.7	36	12
5	0.08±0.038	0.05	40.6±19.1	0.02	1.9	1.7	36	12
6	0.10±0.048	0.14	58.7±21.4	0.33	2.1	1.9	50	11
7	0.08±0.041	0.07	28.8±11.4	0.17	1.7	1.4	21	21
8	0.0001±0.00006	0.007	37.6±8.6	0.23	1.7	1.6	21	6
9	0.09±0.038	0.01	18.4±4.1	0.02	1.4	1.3	−7	8
10	0.02±0.001	0.08	23.9±0.7	0.04	1.4	1.4	0	0
11	0.08±0.035	0.01	34.3±7.1	0.02	1.7	1.6	21	6
12	0.19±0.064	0.33	34.3±7.4	0.20	1.9	1.6	36	19
13	0.03±0.010	0.33	23.7±0.6	0.09	1.4	1.4	0	0
14	0.02±0.008	0.25	24.2±1.0	0.09	1.4	1.4	0	0
15	0.17±0.026	0.25	22.6±1.2	0.08	1.5	1.4	7	7
16	0.02±0.008	0.33	47.9±19.8	0.20	1.9	1.8	36	6
17	0.08±0.043	0.17	24.6±4.7	0.10	1.5	1.3	−7	15

Open in a new tab

Table includes the system number, N; the mean eccentricity, ē, along with its dominate period of oscillation, P_e, in millions of years; mean obliquity, , with its period, P_ψ; the calculated outer edge of the habitable zone, l_out, for temperature based on the habitability index; the outer edge of the habitable zone for the static runs, l_static; and the orbit enhancement factor (E_S) and the variable enhancement factor (E_V); see text.

To separate the effects due to the large amplitudes of the eccentricity and obliquity and their variations, we performed static calculations using the mean values of e and ψ to compute the static outer edge. By comparing this value to the outer edge of the HZ in the variable runs, l_out, we can determine how much of the expansion of the outer edge is due to the variability and how much is due to the large values of e and ψ.

Table 4 lists this comparison as the HZ enhancement factor, as a percentage, due to the static orbital properties of the simulations,

and the HZ enhancement factor due to the variability of those parameters,

where l_baseline is the outer edge of the baseline system, and l_static is the outer edge in the static case. From this analysis, we see that the increase in the outer edge is dominated by the large values of e and ψ. However, a non-negligible component—in one case the sole component, as in System 7—is due to the variability of those values.

The values for E_S and E_V are listed in Table 4 and shown schematically in Fig. 11. In the figure, the height of the bar is equal to E_S for each system. The green portion of the bar indicates the percentage of the enhancement that is due to the variability alone. For example, our most enhanced system, System 3, had a 93% increase in the HZ compared to the baseline model. Of that, 8% of the enhancement is due to the variability of the system. In one case, System 7, the enhancement is entirely due to the variability of the parameters.

In some cases, the nonvariable systems move the HZ inward from the baseline value of 1.4. System 9 shows no enhancement in the full simulation, but the HZ moves inward when variability is removed (hence the negative “enhancement”). In System 17, the static run produces an outer limit that is 7% smaller than the baseline system (again causing negative enhancement), but dynamics allow a 15% increase compared to its static value.

Figure 12 compares the outer edge of the HZ for the variable cases (top panels) and static cases (bottom panels) as a function of the eccentricity and obliquity. The error bars on Systems 1–17 are derived from the standard deviation of the eccentricity and obliquity from the complete simulations. The points in the top panels lie to right of those in the bottom, showing that the variability increases the HZ. However, the effect is most strongly correlated with the obliquity.

In an effort to quantify the relationship, Table 5 lists the linear regression slope and intercept, along with an error-weighted goodness of fit indicator,

where y_i is the “observed” value of either ψ or e, f (x_i) is the linear fit, and σ_i is the uncertainty computed for the given parameter. In the nonvariable cases, since we have no estimate of the uncertainty, the errors in the goodness of fit calculation are taken as 1% of the value of the data point. The results show little correlation with eccentricity but a very strong relationship between obliquity and the outer edge of the HZ. Removing the variability reduces the intercept by 5 degrees and steepens the slope by 7 degrees per astronomical unit, which means the outer edge systematically moves inward when the variability is removed.

6. Discussion

Our simulations show that the evolution of planets' orbit and rotation can increase the maximum separation between a star and a habitable planet by up to 93%. By controlling for the natural extension due to larger eccentricity and obliquity, we found that their oscillations can extend the outer edge by up to 20% and never decrease it. Thus, the number of potentially habitable planets in the Galaxy may be larger than previously thought.

We interpret our results to mean that planets with large and rapid obliquity oscillations are more likely to be habitable than those with negligible oscillations, such as Earth. This perspective is at odds with the notion that the stability of Earth's obliquity is important to the development of life. While it still may be true that rapid oscillations can be detrimental, and certainly at some point obliquity cycles could be too large and rapid, our results clearly show that rapid obliquity evolution can be a boon for habitability. At the least, one should not rule out life on planets with rapid obliquity cycles.

Our results are important for future telescopic searches for life, such as the Terrestrial Planet Finder (TPF). Although a final design has yet to be selected, TPF's mission is to directly image potentially habitable planets. These observations depend critically on large star-planet separations in order to disentangle stellar light from reflected planetary light. Our results show that potentially habitable planets can exist at larger star-planet separations than previously appreciated, improving the odds that TPF can discover an inhabited planet.

By necessity, our study was based on hypothetical planets. While our model systems are extreme from a Solar System point of view, they are relatively tame by exoplanet standards. Eccentricities larger than 0.9 have been discovered (Naef et al., 2001; Jones et al., 2006; Tamuz et al., 2008), very large mutual inclinations are implied by the misalignment between some stars' rotation axis and the orbital planet of a companion planet (Triaud et al., 2010; Naoz et al., 2011), and eccentricity-inclination coupling can drive large and rapid oscillations, if both properties are large (Kozai, 1962; Takeda et al., 2007; Barnes et al., 2011). Thus, our results should not be viewed as an extreme possibility but rather as in the middle of a spectrum of spin-orbit coupling.

Our approach has been simplified in several important ways. While our N-body simulations are the best representation of possible orbits, our rotational model is simplified. A better approach would be to calculate the direct torques from all the bodies in the system and adjust the angular momentum distribution accordingly. Such a model is much more computationally expensive but could be incorporated into an N-body model without too much extra computational cost. Our EBM is also highly idealized, and future improvements could include ocean/land dichotomies, the physics of glacier advancement and retreat, cloud physics, and ultimately even a three-dimensional global circulation model. Each of these additions, however, adds free parameters to a model with very few constraints. While these features will improve realism, we will continue to suffer from a dearth of observational constraints. Nonetheless, as we move toward identifying planets worthy of detailed spectroscopic follow-up, such modeling could provide additional insight for prioritization.

While our results demonstrate that planetary system architecture can influence the position of the HZ, it remains unclear how robust this influence is. For example, our planets all began with a spin rate of 24 h and an obliquity of 23.5 degrees. How do different choices change the picture presented here? Future work should explore a range of initial conditions and determine whether certain architecture always drives the planet into a particular obliquity cycle. If true, then we may be able to characterize a planet's obliquity without direct measurements. While such a study was beyond the scope of this study, the possibility of tightly constraining obliquity is tantalizing and certainly worthy of a follow-up investigation.

Our study suggests that rapid changes in obliquity and eccentricity increase the outer edge of the HZ. We quantify that relationship with linear trends in the enhancement factor with obliquity, but we did not find a threshold to achieve a specific quality that permits significant expansion. We blame the small number of systems we studied for this ambiguity and leave its identification for future work. We note that prior to running a simulation, it is very difficult to know how the orbital and rotational angular momenta will evolve; thus it could take considerable effort to produce a suite of architectures that suitably cover parameter space.

Our study has shown how orbital architecture is a crucial factor when assessing planetary habitability. While previous work has mostly focused on static planetary properties, planets are expected to lie in multiplanet systems; hence the sequence of states must be considered. For the foreseeable future, we will have very few constraints on the properties of potentially habitable planets, and we therefore must leverage any information we have.

Acknowledgments

This work was supported by the NASA Astrobiology Institute's Virtual Planetary Lab lead team and an award from the NAI Director's Discretionary Fund.

Abbreviations

EBM, energy balance model; GCMs, general circulation models; HZ, habitable zone; IHI, Ice Habitability Index; THI, Temperature Habitability Index; TPF, Terrestrial Planet Finder.