<figure><img src="image_1.png" /></figure>

HAL
open science

Exo-Milankovitch Cycles. I. Orbits and Rotation States

Russell Deitrick, Rory Barnes, Thomas Quinn, John Armstrong, Benjamin

Charnay, Caitlyn Wilhelm

► To cite this version:

Russell Deitrick, Rory Barnes, Thomas Quinn, John Armstrong, Benjamin Charnay, et al.. Exo-Milankovitch Cycles. I. Orbits and Rotation States. The Astronomical Journal, 2018, 155 (2), pp.60. 10.3847/1538-3881/aaa301 . hal-02304197

HAL Id: hal-02304197

https://hal.science/hal-02304197v1

Submitted on 20 Aug 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Exo-Milankovitch Cycles. I. Orbits and Rotation States

Russell Deitrick<fnref n="1" />,<fnref n="2" />, Rory Barnes<fnref n="1" />,<fnref n="2" />, Thomas R. Quinn<fnref n="1" />,<fnref n="2" />, John Armstrong<fnref n="2" />,<fnref n="3" />, Benjamin Charnay<fnref n="2" />,<fnref n="4" />, and Caitlyn Wilhelm<fnref n="1" />,<fnref n="2" />

Received 2017 September 8; revised 2017 December 12; accepted 2017 December 17; published 2018 January 16

## Abstract

The obliquity of the Earth, which controls our seasons, varies by only ~2°.5 over ~40,000 years, and its eccentricity varies by only ~0.05 over 100,000 years. Nonetheless, these small variations influence Earth's ice ages. For exoplanets, however, variations can be significantly larger. Previous studies of the habitability of moonless Earth-like exoplanets have found that high obliquities, high eccentricities, and dynamical variations can extend the outer edge of the habitable zone by preventing runaway glaciation (snowball states). We expand upon these studies by exploring the orbital dynamics with a semianalytic model that allows us to map broad regions of parameter space. We find that, in general, the largest drivers of obliquity variations are secular spin-orbit resonances. We show how the obliquity varies in several test cases, including Kepler-62 f, across a wide range of orbital and spin parameters. These obliquity variations, alongside orbital variations, will have a dramatic impact on the climates of such planets.

**Key words:** planetary systems – planets and satellites: atmospheres – planets and satellites: dynamical evolution and stability

## 1. Introduction

The habitable zone (HZ) is the region around a star in which a planet with a roughly Earth-like atmosphere (i.e., dominated by N₂, CO₂, and water) can maintain liquid water on its surface for billions of years (Huang 1959). Though initially this concept was discussed only for the Sun and Sun-like stars, it was expanded to include stellar types F–M by Kasting et al. (1993). Their study used a one-dimensional radiation-convective model to simulate the climates of Earth-like planets orbiting stars of types F, G, K, and M at different orbital distances to determine where the surface temperature was in the correct range for stable liquid water. This definition of the HZ does not include the possibility of subsurface life, such as life that might exist on Mars or Europa; however, currently this definition is appropriate for exoplanets because it is likely to be much more difficult to detect subsurface life across interstellar distances. Many studies have since used one-dimensional radiative-convective models to expand on and refine this work (Selsis et al. 2007; Goldblatt et al. 2013; Kopparapu et al. 2013, 2014).

These HZ definitions are based on the average stellar flux received by a planet and do not fully account for the effects of eccentricity and obliquity. To understand the importance of these parameters, it is worth reviewing the research done on Earth's ice ages and similar theories applied to exoplanets.

The connection between Earth's ice ages and its orbital/ obliquity evolution was first suggested by Joseph Alphonse Adhémar in 1842. This theory was continued by James Croll in the second half of the nineteenth century, but it was not until the twentieth century, when Milutin Milanković performed the first highly accurate calculations of Earth's insolation history and subsequently collaborated with climate scientist Wladimir Köppen, that the theory achieved much success. Hence the theory became commonly known as “Milankovitch theory,”

and the insolation variations were dubbed “Milankovitch cycles.”

Milanković and Köppen identified what is still believed to be the primary connection between orbital variations and the ice ages: summertime insolation at high latitude. During summer, if high-latitude insolation is strong enough to raise temperatures high enough, for a long enough period of time, to melt all snow and ice, then from year to year, ice sheets on land will shrink. Conversely, if the insolation is weaker and temperatures are cooler in summer, then ice will tend to grow from year to year, leading to larger ice sheets.

After Milanković, calculations of Earth's orbital and insolation history continued to improve (some examples are Brouwer & van Woerkom 1950; Sharaf & Boudnikova 1969; Bretagnon 1974; Berger 1978; Berger & Loutre 1991). The development of the dynamical theory of the Earth–Moon–Sun system by Kinoshita (1975, 1977) facilitated perhaps the most accurate determination of the Earth's orbital and obliquity evolution by Laskar (1986) and Laskar et al. (1993a).

There is, however, a great deal of controversy surrounding Milankovitch theory (see Wunsch 2004; Maslin 2016). While it is clear that eccentricity, obliquity, and precession can affect Earth's climate, it is unclear how the rather small variations in these parameters (for Earth) trigger or influence such dramatic changes. Several theories have proven tentatively successful at explaining this (Clark & Pollard 1998; Roe 2006; Tziperman et al. 2006; Huybers & Tziperman 2008), but a detailed review is beyond the scope of this work. For now, we will move on to exoplanets, which can have much more dramatic orbital variations than Earth.

Habitable exoplanets may not, for example, have a large satellite like Earth's Moon. The importance of the Moon to Earth's obliquity was pointed out by Ward (1974). Laskar et al. (1993a, 1993b; and later, Lissauer et al. 2012) showed that, without the Moon, Earth's obliquity could oscillate by a much

<footnote n="1">¹Department of Astronomy, University of Washington, Seattle, WA 98195-1580, USA; deitrr@astro.washington.edu</footnote>

<footnote n="2">²Virtual Planetary Laboratory, University of Washington, Seattle, WA 98195-1580, USA</footnote>

<footnote n="3">³Department of Physics, Weber State University, Ogden, UT 84408-2508, USA</footnote>

<footnote n="4">⁴LESIA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités, UPMC Univ. Paris 06, Univ. Paris Diderot, Sorbonne, Paris Cité, 5 Place Jules Janssen, F-92195 Meudon, France</footnote>
larger amount (~20° over ~1 Myr, rather than the ~2.5° we find with the Moon). As the small oscillation Earth experiences over ~40 kyr seems to influence the ice ages, it was expected, by extrapolation, that this much larger oscillation would be catastrophic.

Habitable exoplanets should not be expected to have large moons like Earth's, as we do not yet understand the probability of the Moon-forming impact. Thus it is important to try to quantify the effects of obliquities very different from Earth's.

One of the first studies to explore the effect of obliquity on habitability was Williams & Kasting (1997). These authors used a 1D energy balance model (EBM) with a parameterized carbon cycle to account for effects obliquity can have on weathering rates and thus changes in atmospheric CO2 pressures. This paper addressed the problem that was proposed by Laskar et al. (1993b), that Earth might be uninhabitable without the Moon to stabilize its obliquity. Williams & Kasting (1997) concluded that Earth could remain habitable with a wide range of obliquities (23°5–90°), but suggested that large regions at the poles would be frozen for extended periods of time at obliquity = 90°.

Williams & Pollard (2003) furthered this study of the obliquity effects on climate using a 3D global climate model (GCM). Though they did not model a large range of solar fluxes, they did find a number of interesting results. Like Williams & Kasting (1997), they found that an Earth-like planet at the current solar flux can remain habitable with any obliquity—none of their simulations experienced a snowball or runaway greenhouse event. Additionally, they noted that at high obliquity, the planet can develop ice belts, that is, ice sheets and sea ice at the equator, rather than ice caps at the poles. In fact, there is some evidence indicating the formation of tropical and midlatitude glaciers on Mars (Forget et al. 2006; Head et al. 2006), which might have occurred during high-obliquity times.

Further studies have investigated the climate response to different obliquities. Spiegel et al. (2009) used a 1D EBM to find the boundaries of the HZ across the range of obliquity 0°–90°. They found that the edges of the HZ can change by as much as 10%, depending on the obliquity and the efficiency of heat diffusion, and showed that high-obliquity planets have high latitudes that are warmer on average than the equatorial regions, but experience extreme seasonal variations. Rose et al. (2017) examined glaciation under different obliquities using an analytical EBM, showing that equatorial ice is much less stable and less likely than polar ice caps.

To lowest order, the effect of eccentricity on climate can be understood by examining the annually averaged stellar flux. The mean stellar flux goes as $S \propto (1 - e^2)^{-1/2}$, where e is the eccentricity (Laskar et al. 1993a). Despite the fact that a planet spends more time at apocenter than at pericenter (Kepler's second law), increasing the orbital eccentricity actually *increases* the amount of flux the planet receives over an orbit. In reality, the effect of eccentricity on climate is more complicated than this simple picture, because of seasonal variations. For example, the duration and intensity of seasons will depend on the shape of the planet's orbit, potentially leading to nonlinear effects on the climate. Further, while the average stellar flux is informative, the primary quantity of interest is temperature: we would like to know whether liquid water can be maintained on the surface of the planet. One informative value is the equilibrium temperature, that is, the

temperature of the pressure level at which the atmosphere becomes optically thick. As noted by Selsis et al. (2007), above $T_{eq} \approx 270$ K, a planet with an Earth-like atmosphere would enter a runaway greenhouse. However, the equilibrium temperature can sometimes be deceptive, for example, Venus's $T_{eq} \approx 227$ K. Rather, the desired quantity is the surface temperature of the planet, and to estimate that, a climate model is necessary, as the relationship with the incoming flux and surface temperature is extremely complicated because of the greenhouse effect, atmospheric circulation, precipitation, and so on. We turn now to several studies that have quantified the relationship between eccentricity and surface temperature.

The HZ is usually conceived as a circular annulus surrounding the star. A planet can have a semimajor axis firmly within the limits, yet can make excursions beyond those limits because of nonzero eccentricity. It is not immediately obvious whether or not such excursions render the planet uninhabitable, because the atmosphere and ocean buffer the heat loss or gain. Williams & Pollard (2002) used a 3D GCM and a 1D EBM to look at the effects of large eccentricity and concluded that seasonal excursions beyond the limits of the HZ are survivable even up to an eccentricity of ~0.7. More recent GCM work by Bolmont et al. (2016) similarly found that excursions beyond the HZ are not typically fatal; however, they found that at $e \sim 0.6$ in the HZs of hotter stars, planets spend so much time near apoastron that snowball events occur.

Using the same EBM as Spiegel et al. (2009), Dressing et al. (2010) showed the stronger dependence of the HZ boundaries on eccentricity, finding that the outer edge of the HZ can be pushed out by up to 0.8 au at $e = 0.5$, for example, for a planet orbiting a solar-type star. These authors did not see snowball events occurring at apoastron, even up to eccentricity of 0.9.

The first study to look at how orbital *variations*, or exo-Milankovitch cycles, affect climate for exoplanets was Spiegel et al. (2010). These authors focused on the eccentricity cycles induced by having a giant planet perturbing a rocky HZ planet. Notably, they found that a frozen planet could be deglaciated if the eccentricity is excited to large enough values.

To further explore the dynamical influence of climate, Armstrong et al. (2014) constructed various planetary systems that included an Earth-mass HZ planet and one or two companions of various masses. To our knowledge, the first study to bring orbital dynamics, obliquity evolution, and climate modeling together for exoplanets, they modeled the evolution of these systems over 1 Myr and then used the resulting orbital and obliquity output as initial conditions for a simple 1D EBM. By also including the growth and decay of ice sheets (i.e., quasi-permanent ice on land), they found that the outer edge of the HZ depends on the orbital evolution. One shortcoming of the study was that the EBM used did not include latitudinal heat diffusion, which was shown by Spiegel et al. (2008, 2009) to be very important. Thus their results are best viewed as a proof of concept, rather than an exact determination of the outer limit of the HZ. That study also did not include a detailed examination of how the large-amplitude obliquity oscillations arise.

Atobe et al. (2004) looked at obliquity evolution of hypothetical Earth-sized planets in systems with known giant planets. They used a second-order orbital theory to map the HZs in terms of the severity of obliquity oscillations, assuming coplanarity. In particular, they noted the importance of secular spin-orbit resonances in triggering large-amplitude variations.
However, Atobe et al. (2004) did not consider the effects of higher order orbital theory and of mutual inclinations.

One of the first studies to apply the concept of exo-Milankovitch cycles to known exoplanets, Brasser et al. (2014) modeled the dynamics of the planetary system orbiting HD 40307 and calculated the resulting insolation that would be received by planet g. This outer planet, if it exists (see Díaz et al. 2016), is probably a gaseous mini-Neptune (Rogers 2015), rather than a rocky super-Earth, but nevertheless Brasser et al. (2014) remains an important study of the potential effects that dynamics can have on habitability. They noted a secular resonance that, depending on planet g's initial obliquity and rotation rate, could strongly influence the obliquity evolution and thus the insolation received at each latitude. These authors did not apply a climate model, however.

More recently, Forgan (2016) looked at Milankovitch cycles for several planets in binary systems with an EBM similar to that used by Spiegel et al. (2009). The presence of a second star can induce large changes in the surface temperatures, particularly for S-type planets (those that orbit a single star within a binary). This problem warrants further study, particularly because surface albedos depend very strongly on the spectral type of the host stars (Shields et al. 2013), so that the stellar fluxes of each star may need to be treated differently in the climate model.

Finally, Shields et al. (2016) modeled the dynamics, including tidal forces, of the Kepler-62 planets and the possible climate of planet f. Kepler-62 f lies in the outer part of the star's HZ, according to the calculations of Kopparapu et al. (2013), so an increased abundance of CO₂ or some other greenhouse gas is necessary to keep the planet from freezing. In addition to demonstrating the stability of the system, these authors also modeled the climate of planet f using two different 3D GCMs and found that with its mean eccentricity, the planet can be habitable but that this depends strongly on its obliquity and atmospheric CO₂ pressure.

In this work, we study the orbital and obliquity evolution of Kepler-62 f, HD 40307 g (briefly), and a mutually inclined Earth-Jovian system. Previous studies have looked at the orbital evolution and obliquity damping (due to tides) of Kepler-62 f (Bolmont et al. 2015; Shields et al. 2016), but have not explored the obliquity evolution due to gravitational perturbations. And while Armstrong et al. (2014) explored the effect of mutual inclinations in a handful of configurations, an in-depth mapping of parameter space for a mutually inclined system, including different initial rotation states, has not yet been done.

The goals of this study are (1) to introduce our secular (orbit-averaged) orbital framework, which allows us to map wide regions of parameter space; (2) to demonstrate the complexity of planetary systems and the potential ways orbital evolution can influence climate, including large eccentricities and mutual inclinations, in preparation for the discovery of dynamically “hot” systems; and (3) to generate initial conditions for climate modeling in a follow-up study (R. Deitrick et al. 2017, in preparation). Our formulations allow us to understand the fundamental causes of large obliquity variations, and to identify higher-order secular resonances that depend on eccentricity and inclination.

Section 2 contains a description of our model, which we validate against an N-body model in Section 3. In Section 4, we lay out results from the model for several example systems

(Kepler-62, HD 40307, and a test system) and discuss secular resonances and Cassini states. In Section 5, we offer some interpretation and broader context, and we conclude in Section 6. Finally, the Appendix contains the disturbing function used in the orbital model, written in terms of our integration variables.

## 2. Methods

### 2.1. Orbital Model

Our model for the orbital evolution, called DISTORB (for “Disturbing function Orbits”), uses the literal, fourth-order disturbing function developed in Murray & Dermott (1999) and Ellis & Murray (2000). We use only the secular terms, meaning that the rapidly varying terms that depend on the mean longitudes of the planets are ignored on the assumption that these terms will average to zero over long timescales. This is a valid assumption as long as no planets are in the proximity of mean-motion resonances. There are, in fact, two orbital models in DISTORB: the first is simply a direct Runge-Kutta integration of the fourth-order equations of motion with variable time-stepping; the second is the Laplace-Lagrange eigenvalue solution, which reduces the accuracy in the disturbing function to second order, but returns a solution that is explicit in time and thus provides a solution in much less computation time. Most of the work we present here takes advantage of the fourth-order solution for its accuracy, but we discuss some results using the Laplace-Lagrange method in Sections 4.1 and 4.2.

The equations of motion are Lagrange's equations (see Murray & Dermott 1999). In the secular approximation, the equations for semimajor axis and mean longitude, and any disturbing function derivative with respect to these variables, are ignored. Additionally, to avoid singularities in the equations for the longitudes of pericenter and ascending node that occur at zero eccentricity and inclination, respectively, we rewrite Lagrange's equations and the disturbing function in terms of the variables (a form of Poincaré coordinates)

$$h = e \sin \varpi \quad (1)$$

$$k = e \cos \varpi \quad (2)$$

$$p = \sin \frac{i}{2} \sin \Omega \quad (3)$$

$$q = \sin \frac{i}{2} \cos \Omega, \quad (4)$$

where e is the orbital eccentricity, i is the inclination, Ω is the longitude of the ascending node, and ϖ = ω + Ω is the longitude of periastron (see Figure 1).

Lagrange's equations for secular theory are then

$$\frac{dh}{dt} = \frac{\sqrt{1 - e^2}}{na^2} \frac{\partial \mathcal{R}}{\partial k} + \frac{kp}{2na^2\sqrt{1 - e^2}} \frac{\partial \mathcal{R}}{\partial p} + \frac{kq}{2na^2\sqrt{1 - e^2}} \frac{\partial \mathcal{R}}{\partial q} \quad (5)$$

$$\frac{dk}{dt} = -\frac{\sqrt{1 - e^2}}{na^2} \frac{\partial \mathcal{R}}{\partial h} - \frac{hp}{2na^2\sqrt{1 - e^2}} \frac{\partial \mathcal{R}}{\partial p} - \frac{hq}{2na^2\sqrt{1 - e^2}} \frac{\partial \mathcal{R}}{\partial q} \quad (6)$$
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Geometry used in the obliquity model DISTROT. The light gray represents the planet's orbital plane, while the darker gray represents a plane of reference. The important orbital angles are the inclination, i, the longitude of ascending node, Ω, and the argument of pericenter, ω. The longitude of pericenter is a “dog-leg” angle, ϖ = Ω + ω. The angle Λ is measured from the vernal point Υ at time t to the ascending node, Ω. The precession angle is defined as ψ = Λ − Ω (also a dog-leg angle). The reference point for Ω is usually chosen as the vernal point at some known date for the solar system, but there is probably a more sensible choice for exoplanetary systems.</figcaption></figure>

$$
\begin{equation}
\begin{aligned}
\frac{dp}{dt} = & -\frac{kp}{2na^2\sqrt{1-e^2}} \frac{\partial \mathcal{R}}{\partial h} + \frac{hp}{2na^2\sqrt{1-e^2}} \frac{\partial \mathcal{R}}{\partial k} \\
& + \frac{1}{4na^2\sqrt{1-e^2}} \frac{\partial \mathcal{R}}{\partial q}
\end{aligned}
\tag{7}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
\frac{dq}{dt} = & -\frac{kp}{2na^2\sqrt{1-e^2}} \frac{\partial \mathcal{R}}{\partial h} + \frac{hq}{2na^2\sqrt{1-e^2}} \frac{\partial \mathcal{R}}{\partial k} \\
& - \frac{1}{4na^2\sqrt{1-e^2}} \frac{\partial \mathcal{R}}{\partial p},
\end{aligned}
\tag{8}
\end{equation}
$$

where $\mathcal{R}$ is the disturbing function (see Appendix), and $a, n, h, k, p, q$ and $e$ are the semimajor axis, mean motion, and eccentricity, respectively. See Berger & Loutre (1991) for the complete set of Lagrange's equations in $h, k, p, q$, and $e$, including mean-motion (e.g., resonant) effects.

General relativity is known to affect the apsidal precession (associated with eccentricity) of planetary orbits, so we include a correction to Equations (5) and (6). Following Laskar (1986), the apsidal corrections are

$$
\left. \frac{dh}{dt} \right|_{\text{GR}} = \delta_R k \qquad (9)
$$

$$
\left. \frac{dk}{dt} \right|_{\text{GR}} = -\delta_R h, \qquad (10)
$$

where

$$
\delta_R = \frac{3n^3 a^2}{c^2(1 - e^2)}, \qquad (11)
$$

and $c$ is the speed of light.

The secular approximation allows us to take large time steps
(years to hundreds of years) and thus to run thousands of
simulations quickly and explore parameter space with relatively
minimal computer usage compared with N-body models.

2.2. Obliquity Model

The obliquity model DISTROT (for “Disturbing function Rotation”) is derived from the Hamiltonian for rigid body motion introduced by Kinoshita (1975) and Kinoshita (1977) and later used by Laskar (1986), Laskar et al. (1993a), Armstrong et al. (2004, 2014), and several others. In the absence of large satellites (such as the Moon), the equations of motion for a rigid planet are

$$
\begin{equation}
\begin{split}
\frac{d\psi}{dt} = R(\varepsilon) - \cot(\varepsilon) [A(p, q)\sin\psi \\
\qquad\qquad\qquad\qquad\qquad\qquad + B(p, q)\cos\psi] - 2\Gamma(p, q) - p_g
\end{split}
\tag{12}
\end{equation}
$$

$$
\frac{d\epsilon}{dt} = -B(p, q)\sin\psi + A(p, q)\cos\psi, \quad (13)
$$

where $\psi$ is the “precession angle” (see Figure 1 and the
following paragraph), $\epsilon$ is the obliquity, and

$$
R(\epsilon) = \frac{3\kappa^2 M_*}{a^3 \nu} \frac{J_2 M r^2}{C} S_0 \cos \epsilon \quad (14)
$$

$$
S_0 = \frac{1}{2} (1 - e^2)^{-3/2} \qquad (15)
$$

$$
A(p, q) = \frac{2}{\sqrt{1 - p^2 - q^2}} [\dot{q} + p\Gamma(p, q)] \quad (16)
$$

$$
B(p, q) = \frac{2}{\sqrt{1 - p^2 - q^2}} [\dot{p} - q\Gamma(p, q)] \quad (17)
$$

$$
\Gamma(p, q) = q\dot{p} - p\dot{q}. \tag{18}
$$

Note the sign error in Equation (8) of Armstrong et al. (2014), corrected in our Equation (16). Our Equation (14) does not contain the lunar constants present in Equation (24) of Laskar (1986). Here, $\epsilon$ represents the obliquity, $p$ and $q$ are as in Equations (3) and (4), $\dot{p}$ and $\dot{q}$ are their time derivatives (Equations (7) and (8)), $\kappa$ is the Gaussian gravitational constant, $M_*$ is the mass of the host star in solar units, $\nu$ is the rotation frequency of the planet in rad day$^{-1}$, $CM^{-1}r^{-2}$ is the specific polar moment of inertia of the planet, and $J_2$ is the gravitational quadrupole of the (oblate) planet.

The final term in Equation (12), $p_g$, accounts for precession due to general relativity and is equal to $\delta_R/2$ (Barker & O'Connell 1970). The symbol $\psi$ refers to the precession angle, defined as $\psi = \Lambda - \Omega$, where $\Lambda$ is the angle between the vernal point $\Upsilon$, the position of the Sun/host star at the planet's northern spring equinox, and the location of the ascending node, $\Omega$, measured from some reference direction $\Upsilon_0$ (often taken to be the direction of the vernal point at some reference date, hence the use of the symbol $\Upsilon$);<fnref n="5" /> see Figure 1. The convention of defining $\Upsilon$ as the location of the Sun at the northern spring equinox is sensible for the solar system since we observe from Earth's surface, but it is a confusing definition to use for exoplanets, for which the direction of the rotation axis is unknown. We adhere to the convention for the sake of consistency with prior literature. In the coming decades, it may become possible to determine the obliquity and orientation of an exoplanet's spin axis; in that event, care should be taken in determining the initial $\psi$ for obliquity modeling. One can

<footnote n="5">⁵ The vernal point occurred in the constellation Aries during Ptolemy’s time, so it is also called the “first point of Aries,” and the “ram’s horn” symbol is used.</footnote>
equivalently refer to $\Upsilon$ as the position of the planet at its northern spring equinox $\pm 180^\circ$. The relevant quantity for determining the insolation, however, is the angle between periastron and the spring equinox, $\omega + \Lambda = \bar{\omega} + \bar{\psi}$.

An additional complication for obliquity evolution is, of course, the presence of a large moon. We do not include the component of the Kinoshita (1975) model that accounts for the lunar torque because the coefficients used are specific to the Earth–Moon–Sun three-body problem and were calculated from the Moon’s orbital evolution (and are therefore not easily generalized). However, we can approximate the effect of the Moon by *forcing* the precession rate (Equation (14)) to be equal to the observed terrestrial value. We do not use this feature of the model in our exploration of exoplanets, but we use it to validate our models by reproducing the Earth’s obliquity evolution.

Equation (12) for the precession angle contains a singularity at $\varepsilon = 0$. To avoid numerical instability, we instead recast Equations (12) and (13) in terms of the rectangular coordinates:

$$ \xi = \sin \varepsilon \sin \psi \quad (19) $$

$$ \zeta = \sin \varepsilon \cos \psi \quad (20) $$

$$ \chi = \cos \varepsilon. \quad (21) $$

The third coordinate, $\chi$, is necessary to preserve sign information when the obliquity crosses $90^\circ$ (see Laskar et al. 1993a). The equations of motion for these variables are then

$$ \frac{d\xi}{dt} = -B(p, q)\sqrt{1 - \xi^2 - \zeta^2} + \zeta[R(\varepsilon) - 2\Gamma(p, q) - p_g] \quad (22) $$

$$ \frac{d\zeta}{dt} = A(p, q)\sqrt{1 - \xi^2 - \zeta^2} - \xi[R(\varepsilon) - 2\Gamma(p, q) - p_g] \quad (23) $$

$$ \frac{d\chi}{dt} = \xi B(p, q) - \zeta A(p, q). \quad (24) $$

The value of $J_2$ is a function of the planet’s rotation rate and density structure. In hydrostatic equilibrium, the Darwin–Radau relation (Cook 1980; Murray & Dermott 1999) shows that the planet’s gravitational quadrupole moment, $J_2$, scales with the rotation rate squared. Earth is very close to hydrostatic equilibrium, while Mars and the Moon are not. For the purposes of this study, we assume that planets are in hydrostatic equilibrium, and we scale $J_2$ according to

$$ J_2 = J_{2\oplus} \left(\frac{\nu}{\nu_{\oplus}}\right)^2 \left(\frac{r}{r_{\oplus}}\right)^3 \left(\frac{M}{M_{\oplus}}\right)^{-1}, \quad (25) $$

where $\nu$ is the rotation rate, $r$ is the mean planetary radius, and $M$ is its mass, and Earth’s measured values of $J_{2\oplus} = 1.08265 \times 10^{-3}$, $\nu_{\oplus} = 7.292115 \times 10^{-5}$ radians s$^{-1}$, $r_{\oplus} = 6.3781 \times 10^6$ m, and $M_{\oplus} = 5.972186 \times 10^{24}$ kg are used for reference ($J_2$ from Cook 1980; $r_{\oplus}$ and $M_{\oplus}$ from Prša et al. 2016). This is identical to the method used by Brasser et al. (2014). Like that study, we assume that $J_2$ cannot go below the measured value of Venus from Yoder (1995), which would ordinarily occur at rotation periods $\gtrsim 13$ days. The assumption then is that there is a limit to hydrostatic equilibrium for a partially rigid body, which is difficult to validate, but it does not affect our results

<table><caption>**Table 1**
Initial Conditions for Kepler-62</caption><thead><tr><th>Planet</th><th>b</th><th>c</th><th>d</th><th>e</th><th>f</th></tr></thead><tbody><tr><td>m (M<sub>⊕</sub>)</td><td>A<br>B</td><td>2.72<br>2.3</td><td>0.136<br>0.1</td><td>14<br>8.2</td><td>6.324<br>4.4</td><td>3.648<br>2.9</td></tr><tr><td>a (au)</td><td></td><td>0.0553<br>0.071</td><td>0.0929<br>0.187</td><td>0.12<br>0.095</td><td>0.427<br>0.13</td><td>0.718<br>0.094</td></tr><tr><td>e</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>i<sub>inv</sub> (°)</td><td></td><td>0.6138</td><td>0.1138</td><td>0.1138</td><td>0.1662</td><td>0.08615</td></tr><tr><td>ω<sub>inv</sub> (°)</td><td></td><td>268.175</td><td>228.074</td><td>241.127</td><td>42.615</td><td>6.277</td></tr><tr><td>Ω<sub>inv</sub> (°)</td><td></td><td>270.002</td><td>270.011</td><td>270.012</td><td>89.992</td><td>89.984</td></tr><tr><td>P<sub>rot</sub> (days)</td><td></td><td></td><td></td><td></td><td></td><td>0.25–20</td></tr><tr><td>ε (°)</td><td></td><td></td><td></td><td></td><td></td><td>0–90</td></tr><tr><td>ψ + Ω (°)</td><td></td><td></td><td></td><td></td><td></td><td>0, 90,<br>180, 270</td></tr><tr><td>i<sub>sky</sub> (°)</td><td></td><td>89.2</td><td>89.7</td><td>89.7</td><td>89.98</td><td>89.9</td></tr><tr><td>ω<sub>sky</sub> (°)</td><td></td><td>178.175</td><td>138.074</td><td>151.127</td><td>312.615</td><td>276.2767</td></tr><tr><td>Ω<sub>sky</sub> (°)</td><td></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr></tbody></table>

greatly as most of our parameter space results in $J_2$ well above this value.

The *specific* polar moment of inertia of a planet, $CM^{-1}r^{-2}$ is always between 0.2 and 0.4. Here we assume the $CM^{-1}r^{-2}$ value of the Earth, 0.33 (Cook 1980). A different assumption for $CM^{-1}r^{-2}$ will change the exact values of our results. However, as the moment of inertia should be similar to this value for a terrestrial planet, the overall nature of our results should not change dramatically.

In Section 4, we employ power spectra to identify resonances between the inclination and obliquity variables. The variables used are $\zeta + \sqrt{-1}\xi$ for the obliquity variations and $q + \sqrt{-1}p$ for the inclination variations. We use the Python function *periodogram* from the package *scipy*. *signal* (Jones et al. 2001) to calculate these power spectra. This function performs fast-Fourier transforms that reduce signal leakage through Welch’s algorithm (Welch 1967) and use a window function to produce a cleaner spectrum. For our study, we use a Hanning window function.

## 2.3. Cassini State Theory

The coupling of the obliquity model and the orbital model allows for the appearance of Cassini states, and indeed we do see several of their features in our simulations. These special rotation states will have implications for obliquity evolution and climate, so it is worth reviewing the theory here.

A planet or satellite is in a Cassini state if Cassini’s laws (initially described in reference to the Moon) are satisfied. Cassini’s laws, generalized and paraphrased here, are as follows: (1) The planet or satellite is in a spin–orbit resonance, such as the synchronous rotation of the Moon or the 3:2 resonance of Mercury. (2) The planet’s or satellite’s obliquity is approximately constant with respect to an appropriate plane of reference, such as the invariable plane (the plane perpendicular to the total angular momentum) of the planetary system or the Laplace plane of the satellite (defined by the host planet’s orbit and equator). (3) Three vectors, the spin angular momentum of the planet or satellite, its orbital angular momentum, and the perpendicular to the reference plane from the second law, are approximately coplanar at all times.

Though Cassini wrote these laws in reference to Earth’s Moon, these were later generalized by Colombo (1966) and
Peale (1969) and found to apply to Mercury and other natural satellites of the solar system. The Cassini state theory developed by these later authors is elegant and mathematically rigorous and can be used to describe the trajectory of a body's spin axis as it evolves toward or around any of four unique Cassini states. The theory can be applied when the body is not actually within one of the Cassini states, provided that the obliquity evolution is dominated by the central torque and a single perturbing frequency (see Hamilton & Ward 2004; Ward & Hamilton 2004). When the theory applies, the body's spin axis will oscillate about one of the Cassini states and Cassini's third law may be satisfied, even when the first two are not. A damping force, such as tidal friction, can drive the body into its ultimate state that satisfies all three of Cassini's laws, hence the prevalence of Cassini states among the solar system satellites.

The conditions specified by the first two laws are easily identified in dynamical simulations of obliquity and orbital evolution. The third law can be identified by invoking the following relation, from Hamilton & Ward (2004):

$$ \sin \Psi = \left\| \frac{(\mathbf{k} \times \mathbf{n}) \times (\mathbf{s} \times \mathbf{n})}{|\mathbf{k} \times \mathbf{n}| |\mathbf{s} \times \mathbf{n}|} \right\|, \quad (26) $$

where **k**, **n**, and **s** are the vectors associated with the perpendicular to the appropriate reference plane (the invariable plane or Laplace plane, for example), the angular momentum of the body's orbit, and the angular momentum of the body's rotation. Alternatively, the complementary relation can be used:

$$ \cos \Psi = \frac{(\mathbf{k} \times \mathbf{n}) \cdot (\mathbf{s} \times \mathbf{n})}{|\mathbf{k} \times \mathbf{n}| |\mathbf{s} \times \mathbf{n}|}. \quad (27) $$

Close inspection indicates that $\Psi = \Lambda = \psi + \Omega$ (when $\Omega$ is measured in the invariable plane of the system), and Cassini's third law is satisfied when this angle librates about 0° or 180° instead of circulating. Since a secular spin-orbit resonance occurs when the axial precession frequency is similar to and opposite in sign to a nodal precession frequency, libration of this angle is also indicative of a strong resonance. A Cassini state is thus a form of secular spin-orbit resonance, though this type of resonance can still occur when not all of Cassini's laws are satisfied.

The locations of Cassini states can be determined from the relation (Ward & Hamilton 2004)

$$ (\alpha/s_i) \cos \theta_j \sin \theta_j + \sin(\theta_j - I_i) = 0, \quad (28) $$

where $\theta_j$ is the angle between the spin-axis vector and a normal to the planet's orbital plane for the jth Cassini state,<fnref n="6" /> $s_i$ is the ith eigenvalue associated with the Laplace–Lagrange solution for inclination, $I_i$ is the corresponding amplitude of the oscillation in inclination, and $\alpha$ is the precessional constant:

$$ \alpha = \frac{R(\varepsilon)}{\cos \varepsilon} = \frac{3\kappa^2 M_* J_2 M r^2}{a^3 \nu C} S_0. \quad (29) $$

This equation can be solved exactly using complex algebra to rewrite it as a quartic polynomial, but the approximate solutions given by Ward & Hamilton (2004) are actually very

accurate. They are

$$ \theta_{1,3} \approx \tan^{-1} \left( \frac{\sin I_i}{1 \pm \alpha / s_i} \right) \quad (30) $$

and

$$ \theta_{2,4} \approx \pm \cos^{-1} \left( \frac{-s_i \cos I_i}{\alpha} \right). \quad (31) $$

Cassini state 1 tends to occur at $|\theta_i| = \varepsilon \approx 0°$ and state 3 at $|\theta_i| = \varepsilon \approx 180°$ (retrograde spin). States 2 and 4 have the same $|\theta_i| = \varepsilon$, but the precession angles differ by $180°$ ($\psi_4 = \psi_2 \pm 180°$); that is, the spin vectors are on opposite sides of the orbit normal. State 4 is a saddle point, as shown in Figure 2 in Ward & Hamilton (2004), and is therefore unstable. Returning to the angle $\Psi = \psi + \Omega$: in the vicinity of Cassini state 2, $\psi + \Omega$ librates about zero; in the vicinity of Cassini state 4, $\psi + \Omega$ circulates or librates with large amplitude (i.e., ~360°) in a “horseshoe” trajectory.

## 2.4. Initial Conditions

The first planetary system we explore is that of Kepler-62, a K-type star with five known transiting planets (Borucki et al. 2011, 2013). The outermost planet, f, is particularly interesting for the purposes of this study because it is in the HZ and far enough from the star to potentially have any spin state (as opposed to being tidally locked and having zero obliquity; Bolmont et al. 2015). We use two sets of masses (the stable case from Bolmont et al. 2015 and the case from Kopparapu et al. 2014) and eccentricities and inclinations from Borucki et al. (2013), varying the obliquity, rotation rate, and precession angle of planet f. The longitude of pericenter, $\varpi$, was randomly chosen for each planet, while the longitude of ascending node, $\Omega$, was set to zero to keep the planets very coplanar. The orbital and spin properties are given in Table 1, where the inclination, longitude of ascending node, and longitude of periastron are listed with respect to the invariable (angular momentum) plane of the system ($i_{inv}$, $\Omega_{inv}$, and $\varpi_{inv}$) and the sky plane ($i_{sky}$, $\Omega_{sky}$, and $\varpi_{sky}$). Since Bolmont et al. (2015) showed that general relativistic corrections are important to the stability of this system, we include those here according to Equations (9) and (10).

The second system, HD 40307, we explore briefly because of its similarity (in some ways) to Kepler-62. The initial conditions are taken directly from Table 6 in Brasser et al. (2014). We vary the initial obliquity and rotation period of the outer planet, g, and examine its obliquity evolution. We run two sets of simulations: one with the initial value of $\psi + \Omega = 0°$ (in the invariable plane of the system) and one with $\psi + \Omega = 180°$.

Our final system, which we will refer to as TSYS, is one with a warm Neptune-mass planet, an Earth-mass planet in the HZ, and a $1.5M_{Jup}$ planet exterior to the HZ. This system is loosely based on HD 190360, which contains a Neptune-mass planet and a super-Jovian, and appears to have a stable HZ (assuming the minimum masses and that there are no additional planets between the two known planets). Though this is a fictitious system, we include it here in anticipation of discoveries that might be made by the PLATO (Rauer et al. 2014) and Gaia (Casertano et al. 2008; Perryman et al. 2014) missions and other future observatories. It is an extension of the work done by Atobe et al. (2004), Spiegel et al. (2010), and Armstrong

<footnote n="6">⁶ If $\theta_j$ is measured over $-180° < \theta_j < 180°$, $|\theta_j|$ is equal to the obliquity, $\varepsilon$.</footnote>
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Obliquity evolution for Earth over the next million years (upper left), from our secular model. Compare to Figure 10, top panel, in Laskar et al. (1993a). Also shown is the obliquity evolution for Earth without the Moon (lower left), from our secular model. Compare to Figure 11 top panel, in Laskar et al. (1993a). On the right, we have the obliquity evolution for Mars backwards in time with the secular orbital solution (upper right) and coupling DISTROT and HNBody (lower right).</figcaption></figure>

et al. (2014), in that it is a test of an Earth-mass planet orbiting a Sun-like star with a giant companion, which should not be ruled out as a possibility because it is undetectable by current observatories. Further, this system displays some interesting features that appear as a result of large mutual inclinations. The initial conditions and parameter space we explore for TSYS are shown in Table 2. We fix the parameters of the two known planets and vary the initial eccentricity, inclination, rotation rate, and obliquity of the Earth-mass planet. The angle $\psi$ we explore briefly in 90° intervals. For the sake of simplicity, we assume the planets 1 and 3 are roughly coplanar.

### 3. Model Validation

The orbital and obliquity models DISTORB and DISTROT are validated by comparison with N-body integrators and known results for the obliquity evolution of solar system bodies (Laskar et al. 1993a, 1993b). Figure 2 shows the resulting obliquity evolution for Earth with the Moon (upper left) and without (lower left). We simulate the effect of the Moon by forcing Earth's precession rate to equal the observed value of $50.290966 \text{ arcsec year}^{-1}$ (Laskar et al. 1993a). Both cases compare well with Laskar et al. (1993a). This figure also shows the obliquity evolution for Mars in the past 10 Myr. The fully secular model captures the primary modes of oscillation and compares well with Ward (1992). In the secular model, we do not see the transition to a higher obliquity regime at ~5 Myr ago seen in the simulations of Touma & Wisdom (1993) and

<table><caption>**Table 2**
Initial Conditions for TSYS</caption><thead><tr><th>Planet</th><th>1</th><th>2</th><th>3</th></tr></thead><tbody><tr><td><i>m</i> (<i>M</i><sub>⊕</sub>)</td><td>18.75</td><td>1</td><td>487.81</td></tr><tr><td><i>a</i> (au)</td><td>0.1292</td><td>1.0031</td><td>3.973</td></tr><tr><td><i>e</i></td><td>0.237</td><td>0.001–0.4</td><td>0.313</td></tr><tr><td><i>i</i> (°)</td><td>1.9894</td><td>0.001–35</td><td>0.02126</td></tr><tr><td><i>ω</i> (°)</td><td>353.23</td><td>100.22</td><td>181.13</td></tr><tr><td><i>Ω</i> (°)</td><td>347.70</td><td>88.22</td><td>227.95</td></tr><tr><td><i>P</i><sub>rot</sub> (days)</td><td></td><td>0.1667–10</td><td></td></tr><tr><td><i>ε</i> (°)</td><td></td><td>0–90</td><td></td></tr><tr><td><i>ψ</i></td><td></td><td>281.78</td><td></td></tr></tbody></table>

Laskar et al. (2004). Mars' obliquity evolution is irregular and, as explained in Touma & Wisdom (1993), is very sensitive to the details of the orbital model. The reason we do not reproduce this regime transition is most likely then because of the approximate nature of our orbital model. For studies of exoplanets, wherein the uncertainties in orbital parameters are much larger than for the solar system, our orbital model is still useful for broad parameter searches. We note that we can reproduce the transition at 5 Myr ago by coupling HNBody (Rauch & Hamilton 2002) to our obliquity model (lower right panel).

One comparison with N-body models, using our system TSYS from Section 4, is shown in Figure 3. The N-body
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Orbital and obliquity evolution for a case in our test system, TSYS (Section 4.3), comparing our secular model (blue) to an *N*-body model (black). There is some drift between the two solutions. Nevertheless, the secular model does an adequate job reproducing the general behavior of the system, in a small fraction of the computation time.</figcaption></figure>

<table><caption>**Table 3**
Eigenvalues for Kepler-62 (Secular Solution)</caption><thead><tr><th></th><th>Mass set A<br>g<sub>i</sub> (arcsec year<sup>-1</sup>)</th><th>s<sub>i</sub> (arcsec year<sup>-1</sup>)</th><th>Mass set B<br>g<sub>i</sub> (arcsec year<sup>-1</sup>)</th><th>s<sub>i</sub> (arcsec year<sup>-1</sup>)</th></tr></thead><tbody><tr><td>1</td><td>4923.54</td><td>-4932.41</td><td>2960.31</td><td>-2966.71</td></tr><tr><td>2</td><td>659.46</td><td>-681.49</td><td>409.23</td><td>-423.26</td></tr><tr><td>3</td><td>76.54</td><td>-1.56 × 10<sup>-5</sup></td><td>61.05</td><td>2.35 × 10<sup>-6</sup></td></tr><tr><td>4</td><td>53.24</td><td>-21.71</td><td>38.16</td><td>-14.07</td></tr><tr><td>5</td><td>16.92</td><td>-63.27</td><td>11.02</td><td>-45.17</td></tr></tbody></table>

integration was done using *Mercury* (Chambers 1999). For the *N*-body results, the obliquity evolution was done using the code used in Armstrong et al. (2014), to provide a comparison to our new code *DISTROT*.

The largest source of error between the secular models *DISTORB* and *DISTROT* comes from a slight offset in the frequencies, which we see as a drift in phase in Figure 3. This is simply because some mean-motion effects, such as minor

variations in the semimajor axes of the planets, are neglected in the secular model. Overall, the secular model captures the general behavior (amplitudes of the oscillations and approximate frequencies) very well. The advantage of the secular model is that it can run millions of years in seconds, compared to hours with the *N*-body model, allowing us to explore wide regions of parameter space. When necessary or advantageous, we also perform the integrations using the *N*-body formulation.
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 4.** Cassini state locations for Kepler-62 f, using our mass set *A*. The left panel corresponds to the fourth eigenvalue in inclination (see Table 3) and the right to the fifth eigenvalue. The fourth eigenvalue has the highest amplitude (see Table 5) and thus is the strongest forcing term on planet f's inclination and obliquity. At a critical value of $|α/s_i|$, states 1 and 4 merge and then vanish, leaving only states 2 and 3 (not shown); see Ward & Hamilton (2004). This occurs for $P_{rot} ≈ 1.2$ days for $s_4$ and $P_{rot} ≈ 0.35$ days for $s_5$.</figcaption></figure>

<table><caption>**Table 4**
Eigenvalues for Kepler-62 (N-body Solution)</caption><thead><tr><th></th><th>Mass set <i>A</i><br>g<sub>i</sub> (arcsec year<sup>-1</sup>)</th><th>s<sub>i</sub> (arcsec year<sup>-1</sup>)</th><th>Mass set <i>B</i><br>g<sub>i</sub> (arcsec year<sup>-1</sup>)</th><th>s<sub>i</sub> (arcsec year<sup>-1</sup>)</th></tr></thead><tbody><tr><td>1</td><td>4877.92</td><td>-4886.79</td><td>2950.05</td><td>-2956.55</td></tr><tr><td>2</td><td>630.38</td><td>-681.71</td><td>380.74</td><td>-423.28</td></tr><tr><td>3</td><td>75.07</td><td>-1.74 × 10<sup>-15</sup></td><td>58.94</td><td>-1.90 × 10<sup>-16</sup></td></tr><tr><td>4</td><td>53.21</td><td>-21.70</td><td>37.94</td><td>-14.06</td></tr><tr><td>5</td><td>16.90</td><td>-63.28</td><td>11.32</td><td>-45.10</td></tr></tbody></table>

# 4. Results

We apply our model to a handful of simple planetary test systems, orbiting K and G stars, which reveal a wealth of dynamical (and potentially climatic) phenomena.

## 4.1. Kepler-62 f

For our first system, Kepler-62, we vary the spin parameters of planet f with the input parameters in Table 1. Because the eigenvalues from the Laplace–Lagrange solution are useful in understanding secular resonances and Cassini states, we calculate their values for both mass sets and both integration schemes (secular and N-body). These are presented in Tables 3 and 4. The eigenvalues $g_i$ are associated with the eccentricity evolution and $s_i$ with the inclination evolution.

In Table 5, we have listed the amplitudes of the inclination oscillations associated with each eigenvalue for planet f. These are calculated from the scaled eigenvectors of the Laplace–Lagrange solution. From this, it is clear that planet f's inclination is primarily driven by the frequency $s_4$, though $s_5$ also contributes strongly. As described below, these two frequencies will be responsible for secular resonances and Cassini states.

From Equations (30) and (31), we can calculate the locations of Cassini states as a function of rotation period. These are shown in Figure 4 for the two dominant eigenvalues, $s_4$ and $s_5$. As described in Ward & Hamilton (2004), Cassini states 1 and 4 merge and vanish at a critical value of $|α/s_i|$ ($α$ depends on

<table><caption>**Table 5**
Inclination Amplitudes for Each Eigenvalue for Kepler-62 f (Secular Solution)</caption><thead><tr><th></th><th><i>I</i><sub>1</sub> (°)</th><th><i>I</i><sub>2</sub> (°)</th><th><i>I</i><sub>3</sub> (°)</th><th><i>I</i><sub>4</sub> (°)</th><th><i>I</i><sub>5</sub> (°)</th></tr></thead><tbody><tr><td><i>A</i></td><td>2.785 × 10<sup>-8</sup></td><td>1.023 × 10<sup>-4</sup></td><td>2.698 × 10<sup>-4</sup></td><td>0.1729</td><td>0.08692</td></tr><tr><td><i>B</i></td><td>4.551 × 10<sup>-8</sup></td><td>1.025 × 10<sup>-4</sup></td><td>2.838 × 10<sup>-4</sup></td><td>0.1607</td><td>0.07788</td></tr></tbody></table>

$P_{rot}$). Beyond that critical value, only states 2 and 3 exist, where $θ_2 ≈ 0°$ and $θ_3 ≈ 180°$. Cassini state 3, which is not shown, always has an obliquity of ≈180°, corresponding to a retrograde spin.

Figure 5 shows the amplitude of the obliquity oscillation for Kepler-62 f, as a function of its initial obliquity and its rotation period, with initial $ψ + Ω = 0°$. The left panel is our secular orbital solution, and the right is the N-body solution (coupling HNBody from Rauch & Hamilton 2002 to DISTROT). There are two clear secular spin–orbit resonances in this parameter space, each visible as a bright arc with shape ~cosε. The black curves show the obliquities associated with Cassini state 2 for eigenvalues $s_4$ and $s_5$. These look very similar to the secular spin–orbit resonance identified for the planet HD 40307 g in Brasser et al. (2014), which we discuss in Section 4.2. Figure 6 shows the same quantities, but for $ψ + Ω = 180°$. The black curves correspond to the Cassini state 4 obliquities. The secular resonances appear in the same
<figure><img src="image_10.png" /><figcaption>**Figure 5.** Amplitude of the obliquity oscillation for Kepler-62 f as a function of the initial obliquity and rotation period, for the stable orbital configuration from Bolmont et al. (2015) (our mass set *A*) with initial $\psi + \Omega = 0^{\circ}$. The left panel shows the solution from our secular models DISTORB and DISTROT, and the right panel shows the same solution achieved by coupling HNBody to DISTROT. Secular resonances appear at rotation periods $\lesssim 2$ days. These are minutely shifted in the *N*-body solution compared to the secular solution (compare the eigenvalues of Tables 3 and 4), because of the small orbital timescale variations in the semimajor axes of the planets. The secular resonances shown here occur in the proximity of Cassini state 2, for the two eigenvalues $s_4$ and $s_5$. The black curves show the predicted locations of the two Cassini states. The white points correspond to cases shown in Figures 7 and 8.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 6.** Same as Figure 5 but for initial $\psi + \Omega = 180^{\circ}$. Here, the secular resonances occur in the proximity of Cassini state 4, which is a saddle point.</figcaption></figure>

locations here, but the structure is somewhat different from those in Figure 5.

Figure 7 shows the orbital and obliquity evolution for planet f at point a in Figure 5. Here, the oscillation of the obliquity is very small, and the angle $\psi + \Omega$ circulates. In Figure 8, we show the evolution of the obliquity and $\psi + \Omega$ at points b, c, and d. Both points b and c are inside the secular resonance associated with frequency $s_4$, and $\psi + \Omega$ librates about $0^{\circ}$; however, the obliquity evolution is markedly different. Point c appears to be in or very near Cassini state 2; thus the $\psi + \Omega$ librates more tightly about $0^{\circ}$ and the obliquity oscillates by only $\sim 0.7^{\circ}$. At point b, the obliquity is farther from the Cassini state and oscillates about it with much larger amplitude.

At point d, we also see a large obliquity oscillation. The secular resonance here is associated with the secondary frequency $s_5$. This resonance is still able to drive the obliquity strongly, though the angle $\psi + \Omega$ does not clearly librate. This is because $s_4$ is still the stronger driver of the inclination.

Figure 9 presents another way of viewing secular resonances. The left panel shows power spectra of the obliquity variables, $\zeta + \sqrt{-1}\xi$, and the inclination variables, $q + \sqrt{-1}p$, at points b (left) and d (right). Note that we have inverted the sign of the frequencies in the inclination variables, because the resonance occurs when the frequencies in obliquity and inclination are opposite in sign. This can be seen by rewriting Equations (12) and (13) in terms of *i* and $\Omega$ instead of *p* and *q* (see Laskar 1986). The
<figure><img src="image_12.png" /><figcaption>**Figure 7.** Orbital and obliquity evolution for Kepler-62 f at point a in Figure 5, with $P_{\text{rot}} = 0.595$ days and $\epsilon_0 = 15^{\circ}$. The quantity $\psi + \Omega$ is plotted rather than $\psi$, because this angle is diagnostic of secular resonances and Cassini states. Here, the planet is well away from resonance, and $\psi + \Omega$ circulates.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 8.** Obliquity evolution at points b (left; $P_{\text{rot}} = 1.091$ days, $\epsilon_0 = 7.5^{\circ}$), c (middle; $P_{\text{rot}} = 1.091$ days, $\epsilon_0 = 15^{\circ}$), and d (right; $P_{\text{rot}} = 0.354$ days, $\epsilon_0 = 22.5^{\circ}$) in the left panel of Figure 5. The orbital evolution is identical to case a in Figure 7. At point b, the planet is in a secular spin-orbit resonance that pumps its obliquity to ≈20°, and $\psi + \Omega$ librates about 0°. It is close to Cassini state 2: the secular resonance can be described as the motion about state 2. At point c, the planet is extremely close to Cassini state 2: the argument $\psi + \Omega$ librates more closely about 0°, and the obliquity oscillation becomes much smaller, nearly satisfying two of Cassini's laws. At point d, the planet is in a spin-orbit resonance with eigenvalue $s_5$. The resonance is not clearly imprinted on the angle $\psi + \Omega$, however, because the eigenvalue $s_4$ still dominates the evolution of $\Omega$.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 9.** Power spectra for the obliquity and inclination evolution (left) for Kepler-62 f, in the secular resonances in Figure 5. The left panel corresponds to point b, the right to point d. The vertical dotted lines are the natural axial precession frequencies from Equation (14) at the minimum and maximum obliquities, and the vertical dashed lines are the inclination eigenvalues from the Laplace–Lagrange solution. In the left panel, the natural precession rate falls nearly on top of the largest inclination frequency, $s_4$. In the right panel, the precession rate lies close to a secondary cluster of frequencies associated with $s_5$. Despite the fact that the argument $\psi + \Omega$ is not seen to librate in this case (Figure 8, right panels), the power spectrum reveals the secular resonance.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 10.** Same as Figure 5 but using mass set B. The left panel has $(\psi + \Omega)_0 = 0^{\circ}$, the right $(\psi + \Omega)_0 = 180^{\circ}$. The secular resonances (and Cassini states) are shifted to longer rotation periods.</figcaption></figure>

longitude of the ascending node, $\Omega$, and the precession angle, $\psi$, appear together inside sine and cosine functions as $\Omega + \psi$.

The dotted lines in this figure, representing the minimum and maximum values of $R(\epsilon)$, lie close to peaks in the inclination. In the case at point b, the axial precession frequency is very close to the primary inclination frequency, $s_4$, while at point d, it is closer to a cluster of peaks at $s_5$. So we see that there is indeed a secular resonance at point d, even though the angle $\psi + \Omega$ does not librate.

Note the difference in the structure of the $s_4$ secular resonance between Figures 5 and 6: the initial angle $\psi + \Omega$ places the planet near Cassini state 2 in Figure 5 and near Cassini state 4 in Figure 6. Since state 2 is stable, there are “islands” of small-obliquity oscillation very close to the black curve in Figure 5. Since state 4 is a saddle point, such islands do not appear in

Figure 6; the “lumpiness” in the structure of the resonance is a result of the discrete nature of the grid of simulations and the interpolation used in making the contours.

Finally, Figure 10 shows the obliquity amplitude for mass set B. The smaller masses of the planets in this case lower all of the frequencies, as compared with mass set A (see Tables 3 and 4). Further, the $J_2$ value of planet f is increased (see Equation (25)), which increases its axial precession rate. The combination of the two effects pushes the resonances toward longer periods.

## 4.2. Comparison with HD 40307 g

In their dynamical analysis of the HD 40307 system, Brasser et al. (2014) noticed the existence of a secular resonance that can amplify the obliquity oscillation of the outermost planet,
<figure><img src="image_16.png" /><figcaption>**Figure 11.** Amplitude of the obliquity oscillation for HD 40307 g as a function of initial obliquity and rotation period. Compare to Figure 8 in Brasser et al. (2014). The left panel started with $\psi + \Omega = 0^{\circ}$, the right with $\psi + \Omega = 180^{\circ}$. There is a secular resonance at rotation rates $\lesssim 2$ days, similar to resonances for Kepler-62 f in Figure 5.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 12.** Evolution of the obliquity and the angle $\psi + \Omega$ for points a (left), b (middle), and c (right) from Figure 11. Point a approaches Cassini state 2 wherein the obliquity oscillation is smaller than the surrounding regions. Point b is librating about Cassini state 2 with a large amplitude, and its obliquity oscillation is large. Point c is near Cassini state 4.</figcaption></figure>

depending on the planet's initial obliquity and rotation rate, with a structure similar to what we find for Kepler-62 f. Note that the existence of planet HD 40307 g is currently under scrutiny (Díaz et al. 2016). Even if planet g does not exist, the system as modeled by Brasser et al. (2014) is nonetheless interesting as a test case of dynamical phenomena.

We have simulated the HD 40307 system in our model according to the initial conditions used by Brasser et al. (2014),

with a starting $\psi = 0^{\circ}$. Plotting the amplitude of the obliquity oscillation of planet g as a function of its initial rotation period and obliquity, we roughly reproduce their Figure 8 (upper left panel) in our Figure 11 (we adjusted the *x* axis to better display the resonance). The left panel of Figure 11 shows our results using our fourth-order model for $\psi + \Omega = 0^{\circ}$, the right for $\psi + \Omega = 180^{\circ}$. Including fourth-order terms reduces the amplitude of the oscillation (as compared to Brasser
<figure><img src="image_18.png" /></figure>

**Figure 13.** Amplitude of the obliquity oscillation for the Earth-mass planet in TSYS, over two slices of our parameter space. On the left, we vary $e_0$ and $i_0$ at a fixed $P_{\text{rot}} = 1$ day and $\epsilon_0 = 23.5^\circ$. A strong secular spin-orbit resonance appears as a curved arc across this space. At the minimum obliquity oscillation on the lower left (point a), Cassini's third law is satisfied. On the right, we vary the initial obliquity and rotation period at fixed $e_0 = 0.2$ and $i_0 = 20^\circ$. Large-amplitude obliquity oscillations appear as bright yellow structures throughout this space. At points c and d, the motion of the obliquity is very similar to that at point b (see text). Minima appear at points e, f, g, which are "islands" of low-amplitude obliquity oscillations. Power spectra indicate that these Cassini-like states are associated with two of the dominant inclination frequencies.

et al. 2014, Figure 8) and shifts the resonance to slightly shorter rotation periods. As Brasser et al. (2014) explain, this resonance can be easily identified in the eigenvalue solution, just as we showed for the two resonances in Kepler-62.

It is not clearly explained how those authors determined that planet g (in their simulations) is in a Cassini state inside the secular resonance, where the obliquity oscillates by the largest amount. However, much as we see for Kepler-62 f, there is a relationship between Cassini states 2 and 4 and the secular resonance. Which Cassini state the planet is near depends on the initial value of the angle $\psi + \Omega$. In the left panel of Figure 11, $\psi + \Omega = 0$, and the planet is near the stable Cassini state 2. In the right panel, $\psi + \Omega = 180^\circ$, and the planet is near the unstable state 4. Near each state, the obliquity oscillates by a large amount as the spin axis circulates about the Cassini state obliquities. However, for state 2, there is a point (as a function of rotation rate) at which the obliquity can become approximately fixed, and we see the "islands" of low-obliquity oscillation inside the resonance. Figure 12 shows the time evolution of the obliquity and the angle $\psi + \Omega$ for three cases (points a, b, and c).

### 4.3. Test System: Earth-mass Planet and Eccentric Giant

Our system, TSYS, is simulated as a test case of a longer-period Earth-mass planet perturbed by a giant planet and a test of the effects of large mutual inclinations and eccentricity. The purpose of this investigation is to demonstrate the complexity and large-amplitude obliquity variations that are possible in dynamically hot systems, which may be discovered in the coming decades. This is a fictitious case of an Earth-mass planet in the HZ of a G-dwarf, with a warm Neptune at ~0.13 au and a Jovian planet at ~4 au. Interestingly, the innermost planet turns out to have a minimal effect on the dynamical evolution of the Earth-mass planet; instead, its orbital and obliquity variations are strongly driven by the Jovian planet at ~4 au.

The panels of Figure 13 show the maximum change in obliquity of the Earth-mass planet over 2 Myr in two slices through parameter space. In general, the amplitude of the obliquity oscillation is expected to be related to the size of the inclination variation, which increases with initial inclination (see Ward 1992). We do see this amplitude correlation, particularly with an initial obliquity of 23.5°; however, there are other noteworthy regions that have substantially larger obliquity oscillations. In particular, there is a curved feature (the "arc") across the center of the left-hand plot in which the obliquity oscillates by ~70° to ~110° and a minimum that appears below the arc in the left-hand panel.

The prominent arc is a result of a secular spin-orbit resonance (Ward 1992). The minimum (point a) is possibly associated with a Cassini state; we refer to this as a "Cassini-like" state, in that Cassini's third law is satisfied, though there does not appear to be a true secular resonance, in which the natural precession rate, $R(\epsilon)$, would be commensurate with a nodal frequency. Figure 14 shows the orbital and obliquity evolution for a case with initial eccentricity $e = 0.0417$ and initial inclination $i = 10.71^\circ$ (point a in Figure 13). Here, Cassini's third law is satisfied, though the second is not: the obliquity still oscillates by ~10°.

The obliquity evolution is particularly dramatic in Figure 15, within the secular resonance. We see that over ~1.5 Myr, the planet experiences obliquities from nearly zero to 80°. Such variations, along with the somewhat large oscillations in eccentricity, would have a large influence on the climate of a planet with an atmosphere similar to Earth's. An interesting component to this evolution is that the obliquity and eccentricity evolve with nearly the same frequency, in sync with each other. It would seem then that the variation in $e$ is tightly coupled to the variation in $i$, which drives the obliquity, and hence all three evolve commensurately.

Figure 16 shows power spectra for both of the cases. The vertical dashed lines represent the natural precession frequency
<figure><img src="image_19.png" /><figcaption>**Figure 14.** Orbital and obliquity evolution for the Earth-mass planet in TSYS, near a minimum in the obliquity amplitude (point a in Figure 13). Cassini's third law is satisfied ($\psi + \Omega$ librates about 0°), though the obliquity still oscillates by $\sim 10^\circ$.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 15.** Same as Figure 14, but for a configuration within the secular resonance (point b in Figure 13). The angle $\psi + \Omega$ librates for much of the time, though not as tightly as for a point near a Cassini-like state, and periodically switches to circulation.</figcaption></figure>

as calculated from Equation (14) at the minimum and maximum obliquity over the 2 Myr simulation. The curves show the power calculated from the obliquity variables, $\zeta + \sqrt{-1}\xi$, and the inclination variables, $q + \sqrt{-1}p$. In both panels, the dominant inclination peaks appear in the obliquity power spectrum because they drive the obliquity, as explained

above. We can see in the case on the left that the natural precession rate (vertical dashed lines) falls in a valley between several strong inclination peaks, indicating that this case is not in a secular spin-orbit resonance. Interestingly, the inclination forcing of the obliquity is stronger than the natural precession rate. The libration of $\psi + \Omega$ is a natural consequence of the
<figure><img src="image_21.png" /><figcaption>**Figure 16.** Power spectra of the obliquity and inclination evolution for the Earth-mass planet in TSYS, at points a (left) and b (right) from Figure 13. The vertical dotted lines correspond to the minimum and maximum *natural* axial precession rate (Equation (14)), which varies because the obliquity varies. Outside the secular resonance, the axial precession frequency falls in between peaks in the inclination spectrum. Within the resonance, this *natural* precession frequency falls right on top of an inclination peak at ≈10 arcsec year⁻¹.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 17.** Power spectra of the obliquity and inclination evolution for the Earth-mass planet in TSYS, at points c (left) and d (right) from Figure 13, in the same format as Figure 16.</figcaption></figure>

inclination cycle dominating over the stellar torque. This case is located at a minimum in the obliquity oscillation, suggesting that Cassini’s second law could be satisfied with a slightly different set of initial conditions.

In the case shown in the right-hand panel of Figure 16, the natural precession rate falls directly on top of one of the inclination peaks and is a true secular spin-orbit resonance. Referring back to Figure 15, we can see that the obliquity response is indeed very strong as a result of the resonance, evolving from a relatively small oscillation of about ~40° to an almost 80° swing.

We also tested different starting precession angles in 90° intervals, to test this parameter’s role, still varying $e_0$ and $i_0$. The secular resonance in this parameter space is insensitive to $\psi$ and appears in all orientations of the equinox ($\psi = 119^{\circ}78, 1019^{\circ}78, 1919^{\circ}78, 2819^{\circ}78$). However, the minimum in Figure 13 (in which point a is located) disappears for $\psi \neq 2819^{\circ}78$. This minimum is purely a result of the initial conditions; nothing drives the system to this state.

Because of the large mutual inclination and eccentricities we explore in this system, the resonance and Cassini states are not easily explained in terms of the inclination eigenvalues from

the Laplace–Lagrange solution. For example, since the eigenvalues do not depend on eccentricity or inclination, they are the same at every point within the left panel of Figure 13. Yet, the secular resonance clearly depends on both orbital elements. Further, the secular resonance does not follow the $\cos\epsilon$ shape (right panel of Figure 13) that we expect from Equation (14) and that we see for Kepler-62 f and HD 40307 g. The origin of this feature is unclear, but it is probably due to the complex nonlinear dynamics of systems with large mutual inclinations.

Figure 17 shows the power spectra of the obliquity and inclination variables for points c and d from Figure 13. The two obliquity spectra are nearly identical, despite a difference in initial obliquity of 55°. In fact, the time evolution of the obliquity in both cases looks extremely similar to that of case b (Figure 15), only shifted in phase. The two cases have very different initial precession rates (~16.9 arcsec year⁻¹ at point c and ~5.1 arcsec year⁻¹ at point d), but the obliquity oscillation is so large that they end up exploring the same range of frequencies. Thus both cases cross the nodal/inclination peak at ~9 arcsec year⁻¹, passing through resonance and driving still larger obliquity oscillations.
<figure><img src="image_23.png" /><figcaption>**Figure 18.** Evolution of the obliquity and the angle $\psi + \Omega$ (left) for the Earth-mass planet at point h in Figure 13, and the power spectra of the inclination and obliquity variables (right). Vertical dashed lines represent the minimum and maximum axial precession rates, which vary with obliquity. The broad peak of the obliquity variables is indicative of chaotic motion.</figcaption></figure>

At points e, f, and g, the obliquity amplitude is reduced compared to the surrounding regions. Each of these “islands” is associated with one of the major inclination frequencies. At point e, the natural axial precession rate is close to the strongest inclination peak at $\sim$31 arcsec year$^{-1}$. At points f and g, the precession is close to the peak at $\sim$9 arcsec year$^{-1}$. This suggests, again, behavior akin to Cassini states, though involving eccentricity–inclination coupling not present in the Laplace–Lagrange eigenvalue solution.

In regions where the obliquity amplitude is $\gtrsim$100°, the motion becomes chaotic: the obliquity undergoes irregular motion, and its power spectrum displays a single, extremely broad peak. Figure 18 shows this for point h. The left-hand panels show the obliquity and $\psi + \Omega$ evolution for 10 Myr in black. The gray curves are the result of slightly different initial conditions—specifically, $\epsilon_0 = 36^{\circ}$72, a difference of 0°01. The conditions in the two simulations diverge at $\sim$2 Myr. The most probable explanation for this chaotic motion is resonance overlap (Chirikov 1979; Wisdom 1980), in which multiple resonances become active and lead to irregular, unpredictable motion. Though a rigorous calculation of the oscillation zones of the active resonances is outside the scope of this work, we can gather some information from the power spectra in Figure 18. In the region where $P_{\text{rot}} \lesssim 0.5$ days and $\epsilon_0 \lesssim 60^{\circ}$, the initial precession rate of the planet’s spin axis is between $\sim$30 arcsec year$^{-1}$ and $\sim$60 arcsec year$^{-1}$. This places the planet amidst two strong inclination frequencies at $\sim$31 arcsec year$^{-1}$ and $\sim$55 arcsec year$^{-1}$. Because the precession rate varies strongly with $\epsilon$ in this region and the obliquity oscillations are so large, the spin axis passes through both resonances easily. Referring to the right panel of Figure 18, we can see that the precession rate (vertical dashed lines) crosses both, as well as the slower frequency at $\sim$9 arcsec year$^{-1}$.

We simulated the same parameter space as in the right panel of Figure 13, but with the initial inclination reduced in 5° increments. At inclinations at and below 10°, the expected $\cos \epsilon$ shape of the secular resonances is recovered, and the obliquity oscillations are reduced. In the case we show here, where the

mutual inclination is increased to $\sim$20°, the obliquity oscillations become larger and the inclination frequency peaks broaden, so that the precession rate, $R(\epsilon)$, passes through or near several strong inclination peaks. This could explain the complex structure of Figure 13, the multiple frequencies in the obliquity evolution and intermittent circulation of $\psi + \Omega$ in Figure 15, and the chaotic motion of $\epsilon$ in some regions.

We also simulated this parameter space with $\psi_0$ rotated by 180°. It is noteworthy that in that case, the islands at points e, f, and g do *not* appear, indicating the importance of the initial orientation of the spin axis. This adds further support to the hypothesis that these regions represent Cassini-like states resulting from mutual inclinations in the fourth-order theory.

## 5. Discussion

In our simulations of Kepler-62 and TSYS, we see that the amplitude of the obliquity variations is largest in secular spin–orbit resonances, in which the natural precession frequency of the planet’s spin axis is near the frequency of inclination variations. The planet’s rotation rate and initial obliquity play an extremely important role, especially in planar systems like Kepler-62. In our TSYS, in which eccentricity and mutual inclination are allowed be somewhat large, we see the appearance of higher-order secular resonances that depend strongly on the initial $e$ and $i$.

For the coplanar system Kepler-62, secular spin–orbit resonances are the largest driver of obliquity evolution of planet f. These resonances are easily identifiable via the Laplace–Lagrange theory and depend intimately on the rotation rate and initial obliquity. While there are a number of studies that lay out possible means of determining rotation rate and obliquity for exoplanets (Pallé et al. 2008; Kawahara & Fujii 2010; Cowan et al. 2009, 2013; Fujii & Kawahara 2012; Schwartz et al. 2016), we are probably a decade or more away from doing such for small-mass planets. In the meantime, modeling studies such as this may be helpful in predicting

<figure><img src="image_24.png" /></figure>
<figure><img src="image_25.png" /><figcaption>**Figure 19.** Peak insolation curves over time for Kepler-62 f as a function of latitude, for two secular resonance cases. The left panel is point b at $P_{\text{rot}} = 1.09$ days and $\epsilon_0 = 7^{\circ}5$, the right at $P_{\text{rot}} = 0.6$ days and $\epsilon_0 = 52^{\circ}5$.</figcaption></figure>

possible climate states for a planet such as this, which may then be verified by visible phase curve observations or spectroscopic detections.

Cassini states are possible for a planet like Kepler-62 f. In the vicinity of Cassini states, the planet's obliquity oscillates with large amplitude, but if a dissipational force (like tidal friction) is at work, the planet could be driven into one of the stable states. State 2 is particularly interesting because it is stable and it can occur at high obliquity, provided the rotation rate is rapid enough. Kepler-62 f might exist in a kind of sweet spot, where the tides from the host star are not strong enough to synchronize its rotation rate, but might be strong enough to push it into a Cassini state. It may be worth testing this by coupling the orbital and obliquity models to a tidal model.

While certainly a powerful influence on climate, it is still not clear how likely terrestrial (and possibly habitable) planets are to exist in mutually inclined systems. The triple Jovian system *v* Andromedae certainly shows that strongly inclined systems exist, at least for high-mass planets (McArthur et al. 2010). Transiting surveys such as *Kepler* have shown the prevalence of compact, coplanar super-Earth systems in the galaxy, but with this detection method we may be missing many inclined systems. The only truly unbiased way of detecting mutually inclined systems is to use a combination of techniques, such as radial velocity and astrometry, or high-resolution spectrometry to detect molecular lines associated with planets. Thus far, these detection methods cannot reliably probe down to low-mass, terrestrial bodies.

On the theory end, numerous studies have shown that instabilities resulting in the loss of one or more planets from a system can result in high inclinations and high eccentricities in the remaining planets (Ford et al. 2005; Raymond et al. 2009; Barnes et al. 2011). Even our own solar system may have had such an event in its early history (Gomes et al. 2005; Batygin et al. 2012). The additional complications associated with binary stars and galactic effects (Kaib et al. 2013) add further weight to the idea that many systems will have experienced such scattering events.

Though, for simplicity, we model planets without large moons (like Earth's) here, a few words are in order. It is commonly understood that in the absence of our moon, Earth would undergo much larger obliquity oscillations than it does (Laskar et al. 1993a; Lissauer et al. 2012; Li & Batygin 2014). Before we assume that planets that have large moons are safe from large variations, we should consider how it is that Earth's moon actually stabilizes its obliquity. The moon, because of its mass and proximity, induces a torque on Earth's equatorial bulge that is stronger than the Sun's, and this torque significantly increases the planet's axial precession rate. This speed-up protects Earth from large obliquity swings, because it is far from the secular frequencies of the solar system, and hence secular spin-orbit resonances are avoided. When we consider the tidal effect of the moon, that is, the slowing of Earth's rotation rate, the matter becomes more complicated. This is because the axial precession rate is also a strong function of rotation rate; in fact, it is predicted that the moon will eventually drag our planet into a secular spin-orbit resonance (Ward 1982). Furthermore, Earth's rotation rate might be very different without the moon and its tidal tug. Thus, we argue that we should not simply assume that a large moon is beneficial (or detrimental, for that matter) to climate and habitability. Exoplanets with large moons may be driven into and out of secular spin-orbit resonances as their moons alter the planetary rotation rate.

We expect secular resonances of this nature will have a substantial impact on the climate, if the atmosphere is at all similar to Earth's. Figure 19 shows the insolation (stellar flux) received by Kepler-62 f in the secular resonance (point b in Figure 5 and a case at $P_{\text{rot}} = 0.6$ days and $\epsilon_0 = 52^{\circ}5$), for all latitudes. The insolation shown here is the peak summertime value at each latitude, calculated from the formulae in Berger (1978). For this calculation, the stellar constant is ~573.18 W m⁻², about 42% of Earth's value. The annual insolation changes dramatically in the case with low starting obliquity, especially for the poles, as the obliquity oscillates. If a planet like this has surface water, we might expect it to have severe ice-age cycles. In the high-obliquity case, the poles
receive more intense seasonal insolation than the equator. The obliquity oscillation is much smaller than the low-obliquity case, but the poles still undergo variations in peak stellar flux of ≈40 W m⁻². Our next paper will address how this impacts the surface temperature and overall habitability of the planet.

How might we determine if an exoplanet is undergoing Milankovitch cycles, given that we cannot observe planets for tens of thousands of years and are unable to travel there to collect rock records? It is worth mentioning that Milankovitch cycles have been inferred for the planet Mars from images of the polar regions taken from orbit (Laskar et al. 2002). Though we are technologically far away from the possibility of doing something similar for small HZ planets, several studies have already begun to work out techniques for creating 2D maps of a planetary surface (Pallé et al. 2008; Kawahara & Fujii 2010; Fujii & Kawahara 2012; Cowan et al. 2013; Kawahara 2016; Schwartz et al. 2016). This technique will likely require a high-contrast, large-aperture telescope, like the recently studied Large UltraViolet Optical and Infrared surveyor mission (Bolcar et al. 2016), capable of directly imaging terrestrial planets. As shown in Schwartz et al. (2016), the technique is capable of determining longitudinal brightness and color variations, obliquity, and, in some cases, latitudinal surface variations. From this, we may be able to determine some limited history of the planet's climate, if, say, the planet has large ice features at latitudes inconsistent with its current obliquity. Kepler-62 f, if it has an obliquity oscillation similar to what we find here, would be a difficult test case because the obliquity evolves rather slowly, which would allow the ice latitude to equilibrate more easily. It would be much easier to observe the inconsistency for a terrestrial planet with a giant companion, like our TSYS, because the obliquity oscillation is more rapid.

From orbital information, it is certainly reasonable to determine whether Milankovitch cycles are a possibility. If the orbital/obliquity variations are large enough to trigger a climate catastrophe (e.g., a snowball state), that might be observable through the planet's albedo, which is the subject of our follow-up study. If the rotation rate and obliquity can be constrained through surface mapping, as discussed above, the case for or against Milankovitch cycles can be further strengthened.

Milankovitch cycles may influence the climates of planets with surface oceans, but would probably have little or no effect on subsurface oceans. Hence, the habitability of planets analogous to Europa or Enceladus would not be affected by this kind of dynamical evolution, unless tidal heating of the interior is significantly modified by obliquity or eccentricity variations. However, the interaction of life with the atmosphere will likely provide the only biosignatures detectable across interstellar distances, hence our focus on surface habitability.

## 6. Conclusions

We have presented a secular model for the orbital and obliquity oscillations associated with planetary perturbations. The speed of the model allows us to explore broad regions of parameter space and to rapidly identify interesting dynamical phenomena. We apply this model to several test cases of Earth-like (or potentially Earth-like) exoplanets to understand how climate and habitability might be affected by planetary companions.

We studied three systems here, with a focus on secular spin-orbit resonances and Cassini states. Secular resonances can potentially have a dramatic impact on climate by inducing strong obliquity variations. In the next paper, we will explore whether such obliquity variations, combined with eccentricity variations and precession, make planets more or less susceptible toward climate states that might be rendered the surface inhospitable to life.

We must emphasize that planetary systems are extremely complex, and that HZs based simply on semimajor axis may be fundamentally limited. In assessing a planet's potential to host life, it is important that companion planets be looked for and considered as possible complications, especially for planets that may have not undergone strong tidal evolution.

This study demonstrates the importance of companion planets and orbital architecture when assessing the habitability of exoplanets. Used in tandem with climate models, this type of analysis will provide testable predictions of atmospheric states for upcoming exoplanet characterization missions and will aid in the interpretation of potential biosignatures. In our next paper, we will apply a simplified climate model to further explore this problem.

This work was supported by the NASA Astrobiology Institute's Virtual Planetary Laboratory under Cooperative Agreement number NNA13AA93A. This work was facilitated though the use of advanced computational, storage, and networking infrastructure provided by the Hyak supercomputer system at the University of Washington. The results reported herein benefited from the authors' affiliation with the NASA's Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA's Science Mission Directorate. We would also like to thank Ramon Brasser for extremely helpful feedback.

# Appendix

## The Disturbing Function and Its Derivatives

Here we present, for the sake of completeness, the disturbing function as used in DISTORB, in the variables *h*, *k*, *p*, and *q* (see Section 2.1). These were originally derived by Ellis & Murray (2000); we have simply applied coordinate transformations and calculated derivatives with respect to the new coordinates. This disturbing function, in its original form, can also be seen in Murray & Dermott (1999). We will not restate the semimajor axis functions, $f_1, f_2, f_3$, and so on, in this paper, as they are taken directly from Table B.3 of Murray & Dermott (1999).

The secular disturbing function, for any pair of planets, is

$$R = \frac{\mu'}{a'} R_D \qquad (32)$$

for the inner body and

$$R' = \frac{\mu}{a'} R_D \qquad (33)$$

for the outer body. Here, $a'$ is the semimajor axis of the exterior planet, and the mass factors are $\mu = \kappa^2 m$ and $\mu' = \kappa^2 m'$, where *m* is the mass of the interior planet and *m'* is the mass of the exterior planet. Finally,

$$R_D = D0.1 + D0.2 + D0.3 + \dots, \qquad (34)$$

where the terms D0.1, D0.2, and so on are given in Table 6.
<table><caption>**Table 6**
Disturbing Function</caption><thead><tr><th>Term</th><th></th></tr></thead><tbody><tr><td>D0.1</td><td>f<sub>1</sub> + (h<sup>2</sup> + k<sup>2</sup> + h'<sup>2</sup> + k'<sup>2</sup>)f<sub>2</sub> + (p<sup>2</sup> + q<sup>2</sup> + p'<sup>2</sup> + q'<sup>2</sup>)f<sub>3</sub><br> + (h<sup>2</sup> + k<sup>2</sup>)<sup>2</sup>f<sub>4</sub> + (h<sup>2</sup> + k<sup>2</sup>)(h'<sup>2</sup> + k'<sup>2</sup>)f<sub>5</sub> + (h'<sup>2</sup> + k'<sup>2</sup>)<sup>2</sup>f<sub>6</sub><br> + [(h<sup>2</sup> + k<sup>2</sup>)(p<sup>2</sup> + q<sup>2</sup>) + (h'<sup>2</sup> + k'<sup>2</sup>)(p<sup>2</sup> + q<sup>2</sup>)<br> + (h<sup>2</sup> + k<sup>2</sup>)(p'<sup>2</sup> + q'<sup>2</sup>) + (h'<sup>2</sup> + k'<sup>2</sup>)(p'<sup>2</sup> + q'<sup>2</sup>)]f<sub>7</sub><br> + [(p<sup>2</sup> + q<sup>2</sup>)<sup>2</sup> + (p'<sup>2</sup> + q'<sup>2</sup>)<sup>2</sup>]f<sub>8</sub> + (p<sup>2</sup> + q<sup>2</sup>)(p'<sup>2</sup> + q'<sup>2</sup>)f<sub>9</sub></td></tr><tr><td>D0.2</td><td>(hh' + kk')[(f<sub>10</sub> + (h<sup>2</sup> + k<sup>2</sup>)f<sub>11</sub> + (h'<sup>2</sup> + k'<sup>2</sup>)f<sub>12</sub><br> + (p<sup>2</sup> + q<sup>2</sup> + p'<sup>2</sup> + q'<sup>2</sup>)f<sub>13</sub>]</td></tr><tr><td>D0.3</td><td>(pp' + qq')[(f<sub>14</sub> + (h<sup>2</sup> + k<sup>2</sup> + h'<sup>2</sup> + k'<sup>2</sup>)f<sub>15</sub><br> + (p<sup>2</sup> + q<sup>2</sup> + p'<sup>2</sup> + q'<sup>2</sup>)f<sub>16</sub>]</td></tr><tr><td>D0.4</td><td>(h<sup>2</sup>h'<sup>2</sup> - k<sup>2</sup>h'<sup>2</sup> - h<sup>2</sup>k'<sup>2</sup> + k<sup>2</sup>k'<sup>2</sup>) + 4hh'kk')f<sub>17</sub></td></tr><tr><td>D0.5</td><td>(h<sup>2</sup>p<sup>2</sup> - h<sup>2</sup>q<sup>2</sup> - k<sup>2</sup>p<sup>2</sup> + k<sup>2</sup>q<sup>2</sup> + 4hkpq')f<sub>18</sub></td></tr><tr><td>D0.6</td><td>[hh'(p<sup>2</sup> - q<sup>2</sup>) - kk'(p<sup>2</sup> - q<sup>2</sup>) + 2pq(hk' + kh')]f<sub>19</sub></td></tr><tr><td>D0.7</td><td>(h<sup>2</sup>p<sup>2</sup> - h<sup>2</sup>q<sup>2</sup> - k<sup>2</sup>p<sup>2</sup> + k<sup>2</sup>q<sup>2</sup> + 4h'k'pq')f<sub>20</sub></td></tr><tr><td>D0.8</td><td>(h<sup>2</sup>pp' - h<sup>2</sup>qq' - k<sup>2</sup>pp' + k<sup>2</sup>qq' + 2hkp'q + 2hkpq')f<sub>21</sub></td></tr><tr><td>D0.9</td><td>[(hh' + kk')(pp' + qq') + (hk' - kh')(pq' - qp')]f<sub>22</sub></td></tr><tr><td>D0.10</td><td>[(hh' + kk')(pp' + qq') + (hk' - kh')(qp' - pq')]f<sub>23</sub></td></tr><tr><td>D0.11</td><td>[(hh' - kk')(pp' - qq') + (hk' + kh')(pq' + qp')]f<sub>24</sub></td></tr><tr><td>D0.12</td><td>(h<sup>2</sup>pp' - h<sup>2</sup>qq' - k<sup>2</sup>pp' + k<sup>2</sup>qq' + 2h'k'p'q + 2h'k'pq')f<sub>25</sub></td></tr><tr><td>D0.13</td><td>(h<sup>2</sup>p<sup>2</sup> - h<sup>2</sup>q<sup>2</sup> - k<sup>2</sup>p<sup>2</sup> + k<sup>2</sup>q<sup>2</sup> + 4hkp'q')f<sub>18</sub></td></tr><tr><td>D0.14</td><td>[hh'(p<sup>2</sup> - q<sup>2</sup>) - kk'(p<sup>2</sup> - q<sup>2</sup>) + 2p'q'(hk' + kh')]f<sub>19</sub></td></tr><tr><td>D0.15</td><td>(h<sup>2</sup>p<sup>2</sup> - h<sup>2</sup>q<sup>2</sup> - k<sup>2</sup>p<sup>2</sup> + k<sup>2</sup>q<sup>2</sup> + 4h'k'p'q')f<sub>20</sub></td></tr><tr><td>D0.16</td><td>(p<sup>2</sup>p<sup>2</sup> - p<sup>2</sup>q<sup>2</sup> - q<sup>2</sup>p<sup>2</sup> + q<sup>2</sup>q<sup>2</sup> + 4pqp'q')f<sub>26</sub></td></tr></tbody></table>

ORCID iDs

Benjamin Charnay https://orcid.org/0000-0003-0977-6545

References

Armstrong, J. C., Barnes, R., Domagal-Goldman, S., et al. 2014, AsBio, 14, 277

Armstrong, J. C., Leovy, C. B., & Quinn, T. 2004, Icar, 171, 255

Atobe, K., Ida, S., & Ito, T. 2004, Icar, 168, 223

Barker, B. M., & O'Connell, R. F. 1970, PhRvD, 2, 1428

Barnes, R., Greenberg, R., Quinn, T. R., McArthur, B. E., & Benedict, G. F. 2011, ApJ, 726, 71

Batygin, K., Brown, M. E., & Betts, H. 2012, ApJL, 744, L3

Berger, A., & Loutre, M. F. 1991, QSRv, 10, 297

Berger, A. L. 1978, JATs, 35, 2362

Bolcar, M. R., Feinberg, L., France, K., et al. 2016, Proc. SPIE, 9904, 99040J

Bolmont, E., Libert, A.-S., Leconte, J., & Selsis, F. 2016, A&A, 591, A106

Bolmont, E., Raymond, S. N., Leconte, J., Hersant, F., & Correia, A. C. M. 2015, A&A, 583, A116

Borucki, W. J., Agol, E., Fressin, F., et al. 2013, Sci, 340, 587

Borucki, W. J., Koch, D. G., Basri, G., et al. 2011, ApJ, 736, 19

Brasser, R., Ida, S., & Kokubo, E. 2014, MNRAS, 440, 3685

Bretagnon, P. 1974, A&A, 30, 141

Brouwer, D., & van Woerkom, A. 1950, Astron. Papers Amer. Ephem., Vol. 13 (Washington, D.C.: U.S. Govt. Print. Off.), 81

Casertano, S., Lattanzi, M. G., Sozzetti, A., et al. 2008, A&A, 482, 699

Chambers, J. E. 1999, MNRAS, 304, 793

Chirikov, B. V. 1979, PhR, 52, 263

Clark, P. U., & Pollard, D. 1998, PalOc, 13, 1

Colombo, G. 1966, AJ, 71, 891

Cook, A. H. 1980, Interiors of the planets (Cambridge: Cambridge Univ. Press)

Cowan, N. B., Agol, E., Meadows, V. S., et al. 2009, ApJ, 700, 915

Cowan, N. B., Fuentes, P. A., & Haggard, H. M. 2013, MNRAS, 434, 2465

Díaz, R. F., Ségransan, D., Udry, S., et al. 2016, A&A, 585, A134

Dressing, C. D., Spiegel, D. S., Scharf, C. A., Menou, K., & Raymond, S. N. 2010, ApJ, 721, 1295

Ellis, K. M., & Murray, C. D. 2000, Icar, 147, 129

Ford, E. B., Lystad, V., & Rasio, F. A. 2005, Natur, 434, 873

Forgan, D. 2016, MNRAS, 463, 2768

Forget, F., Haberle, R. M., Montmessin, F., Levrard, B., & Head, J. W. 2006, Sci, 311, 368

Fujii, Y., & Kawahara, H. 2012, ApJ, 755, 101

Goldblatt, C., Robinson, T. D., Zahnle, K. J., & Crisp, D. 2013, NatGe, 6, 661

Gomes, R., Levison, H. F., Tsiganis, K., & Morbidelli, A. 2005, Natur, 435, 466

Hamilton, D. P., & Ward, W. R. 2004, AJ, 128, 2510

Head, J. W., Marchant, D. R., Agnew, M. C., Fassett, C. I., & Kreslavsky, M. A. 2006, E&PSL, 241, 663

Huang, S. 1959, AmSci, 47, 397

Huybers, P., & Tziperman, E. 2008, PalOc, 23, PA1208

Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open source scientific tools for Python, http://www.scipy.org/

Kaib, N. A., Raymond, S. N., & Duncan, M. 2013, Natur, 493, 381

Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, Icar, 101, 108

Kawahara, H., & Fujii, Y. 2010, ApJ, 720, 1333

Kawahara, H. 2016, ApJ, 822, 112

Kinoshita, H. 1975, SAOSR, 364

Kinoshita, H. 1977, CeMec, 15, 277

Kopparapu, R. K., Ramirez, R., Kasting, J. F., et al. 2013, ApJ, 765 131

Kopparapu, R. K., Ramirez, R. M., SchottelKotte, J., et al. 2014, ApJL, 787, L29

Laskar, J. 1986, A&A, 157, 59

Laskar, J., Correia, A. C. M., Gastineau, M., et al. 2004, Icar, 170, 343

Laskar, J., Joutel, F., & Boudin, F. 1993a, A&A, 270, 522

Laskar, J., Joutel, F., & Robutel, P. 1993b, Natur, 361, 615

Laskar, J., Levrard, B., & Mustard, J. F. 2002, Natur, 419, 375

Li, G., & Batygin, K. 2014, ApJ, 790, 69

Lissauer, J. J., Barnes, J. W., & Chambers, J. E. 2012, Icar, 217, 77

Maslin, M. 2016, Natur, 540, 208

McArthur, B. E., Benedict, G. F., Barnes, R., et al. 2010, ApJ, 715, 1203

Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge: Cambridge Univ. Press)

Pallé, E., Ford, E. B., Seager, S., Montañés-Rodríguez, P., & Vazquez, M. 2008, ApJ, 676, 1319

Peale, S. J. 1969, AJ, 74, 483

Perryman, M., Hartman, J., Bakos, G. Á., & Lindegren, L. 2014, ApJ, 797, 14

Prša, A., Harmanec, P., Torres, G., et al. 2016, AJ, 152, 41

Rauch, K. P., & Hamilton, D. P. 2002, BAAS, 34, 938

Rauer, H., Catala, C., Aerts, C., et al. 2014, ExA, 38, 249

Raymond, S. N., Armitage, P. J., & Gorelick, N. 2009, ApJL, 699, L88

Roe, G. 2006, GeoRL, 33, L24703

Rogers, L. A. 2015, ApJ, 801, 41

Rose, B. E. J., Cronin, T. W., & Bitz, C. M. 2017, ApJ, 846, 28

Schwartz, J. C., Sekowski, C., Haggard, H. M., Pallé, E., & Cowan, N. B. 2016, MNRAS, 457, 926

Selsis, F., Kasting, J. F., Levrard, B., et al. 2007, A&A, 476, 1373

Sharaf, S. G., & Boudnikova, N. A. 1969, TrLen, 14, 48

Shields, A. L., Barnes, R., Agol, E., et al. 2016, AsBio, 16, 443

Shields, A. L., Meadows, V. S., Bitz, C. M., et al. 2013, AsBio, 13, 715

Spiegel, D. S., Menou, K., & Scharf, C. A. 2008, ApJ, 681, 1609

Spiegel, D. S., Menou, K., & Scharf, C. A. 2009, ApJ, 691, 596

Spiegel, D. S., Raymond, S. N., Dressing, C. D., Scharf, C. A., & Mitchell, J. L. 2010, ApJ, 721, 1308

Touma, J., & Wisdom, J. 1993, Sci, 259, 1294

Tziperman, E., Raymo, M. E., Huybers, P., & Wunsch, C. 2006, PalOc, 21, PA4206

Ward, W. R. 1974, JGR, 79, 3375
Ward, W. R. 1982, Icar, 50, 444

Ward, W. R. 1992, in Mars, Long-term orbital and spin dynamics of Mars, ed.
H. H. Kieffer et al., 298

Ward, W. R., & Hamilton, D. P. 2004, AJ, 128, 2501

Welch, P. D. 1967, in IEEE Trans. Audio & Electroacoust., Volume AU-15
(New York: IEEE), 70

Williams, D. M., & Kasting, J. F. 1997, Icar, 129, 254

Williams, D. M., & Pollard, D. 2002, IJAsB, 1, 61

Williams, D. M., & Pollard, D. 2003, IJAsB, 2, 1

Wisdom, J. 1980, AJ, 85, 1122

Wunsch, C. 2004, QSRv, 23, 1001

Yoder, C. F. 1995, Icar, 117, 250