Exo-Milankovitch Cycles. I. Orbits and Rotation States

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)

$$\begin{eqnarray}&&h=e\sin \varpi \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&k=e\cos \varpi \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&p=\sin \displaystyle \frac{i}{2}\sin {\rm{\Omega }}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&q=\sin \displaystyle \frac{i}{2}\cos {\rm{\Omega }},\end{eqnarray} \tag{ 4 }$$

where e is the orbital eccentricity, i is the inclination, Ω is the longitude of the ascending node, and $\varpi =\omega +{\rm{\Omega }}$ is the longitude of periastron (see Figure 1).

Zoom In Zoom Out Reset image size

Lagrange's equations for secular theory are then

$$\begin{eqnarray}\displaystyle \frac{{dh}}{{dt}} & = & \displaystyle \frac{\sqrt{1-{e}^{2}}}{{{na}}^{2}}\displaystyle \frac{\partial { \mathcal R }}{\partial k}+\displaystyle \frac{{kp}}{2{{na}}^{2}\sqrt{1-{e}^{2}}}\displaystyle \frac{\partial { \mathcal R }}{\partial p}\\ & & +\displaystyle \frac{{kq}}{2{{na}}^{2}\sqrt{1-{e}^{2}}}\displaystyle \frac{\partial { \mathcal R }}{\partial q}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}\displaystyle \frac{{dk}}{{dt}} & = & -\displaystyle \frac{\sqrt{1-{e}^{2}}}{{{na}}^{2}}\displaystyle \frac{\partial { \mathcal R }}{\partial h}-\displaystyle \frac{{hp}}{2{{na}}^{2}\sqrt{1-{e}^{2}}}\displaystyle \frac{\partial { \mathcal R }}{\partial p}\\ & & -\displaystyle \frac{{hq}}{2{{na}}^{2}\sqrt{1-{e}^{2}}}\displaystyle \frac{\partial { \mathcal R }}{\partial q}\end{eqnarray} \tag{ 6 }$$

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

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

where ${ \mathcal R }$ is the disturbing function (see Appendix), and a, n, 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, and q, 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

$$\begin{eqnarray}&&{\left.\displaystyle \frac{{dh}}{{dt}}\right|}_{\mathrm{GR}}={\delta }_{R}k\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{\left.\displaystyle \frac{{dk}}{{dt}}\right|}_{\mathrm{GR}}=-{\delta }_{R}h,\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray}&&{\delta }_{R}=\displaystyle \frac{3{n}^{3}{a}^{2}}{{c}^{2}(1-{e}^{2})},\end{eqnarray} \tag{ 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.

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{eqnarray}\displaystyle \frac{d\psi }{{dt}} & = & R(\varepsilon )-\cot (\varepsilon )\,[A(p,q)\sin \psi \\ & & +B(p,q)\cos \psi ]-2{\rm{\Gamma }}(p,q)-{p}_{g}\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d\varepsilon }{{dt}}=-B(p,q)\sin \psi +A(p,q)\cos \psi ,\end{eqnarray} \tag{ 13 }$$

where ψ is the "precession angle" (see Figure 1 and the following paragraph), ε is the obliquity, and

$$\begin{eqnarray}&&R(\varepsilon )=\displaystyle \frac{3{\kappa }^{2}{M}_{\star }}{{a}^{3}\nu }\displaystyle \frac{{J}_{2}{{Mr}}^{2}}{C}{S}_{0}\cos \varepsilon \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{S}_{0}=\displaystyle \frac{1}{2}{(1-{e}^{2})}^{-3/2}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&A(p,q)=\displaystyle \frac{2}{\sqrt{1-{p}^{2}-{q}^{2}}}[\dot{q}+p{\rm{\Gamma }}(p,q)]\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&B(p,q)=\displaystyle \frac{2}{\sqrt{1-{p}^{2}-{q}^{2}}}[\dot{p}-q{\rm{\Gamma }}(p,q)]\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{\rm{\Gamma }}(p,q)=q\dot{p}-p\dot{q}.\end{eqnarray} \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, ε 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)), κ is the Gaussian gravitational constant, M_⋆ is the mass of the host star in solar units, ν is the rotation frequency of the planet in rad day⁻¹, ${{CM}}^{-1}{r}^{-2}$ is the specific polar moment of inertia of the planet, and J₂ 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 ψ refers to the precession angle, defined as $\psi ={\rm{\Lambda }}-{\rm{\Omega }}$, where Λ is the angle between the vernal point ϒ, the position of the Sun/host star at the planet's northern spring equinox, and the location of the ascending node, Ω, measured from some reference direction ϒ₀ (often taken to be the direction of the vernal point at some reference date, hence the use of the symbol ϒ);⁵ see Figure 1. The convention of defining ϒ 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 ψ for obliquity modeling. One can equivalently refer to ϒ 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 +{\rm{\Lambda }}=\varpi +\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:

$$\begin{eqnarray}&&\xi =\sin \varepsilon \sin \psi \end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&\zeta =\sin \varepsilon \cos \psi \end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&\chi =\cos \varepsilon .\end{eqnarray} \tag{ 21 }$$

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

$$\begin{eqnarray}&&\displaystyle \frac{d\xi }{{dt}}=-B(p,q)\sqrt{1-{\xi }^{2}-{\zeta }^{2}}+\zeta [R(\varepsilon )-2{\rm{\Gamma }}(p,q)-{p}_{g}]\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d\zeta }{{dt}}=A(p,q)\sqrt{1-{\xi }^{2}-{\zeta }^{2}}-\xi [R(\varepsilon )-2{\rm{\Gamma }}(p,q)-{p}_{g}]\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d\chi }{{dt}}=\xi B(p,q)-\zeta A(p,q).\end{eqnarray} \tag{ 24 }$$

The value of J₂ 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₂, 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₂ according to

$$\begin{eqnarray}&&{J}_{2}={J}_{2\oplus }{\left(\displaystyle \frac{\nu }{{\nu }_{\oplus }}\right)}^{2}{\left(\displaystyle \frac{r}{{r}_{\oplus }}\right)}^{3}{\left(\displaystyle \frac{M}{{M}_{\oplus }}\right)}^{-1},\end{eqnarray} \tag{ 25 }$$

where ν 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⁻¹, ${r}_{\oplus }=6.3781\,\times {10}^{6}$ m, and ${M}_{\oplus }=5.972186\times {10}^{24}$ kg are used for reference (J₂ from Cook 1980; r_⊕ and M_⊕ 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₂ 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 greatly as most of our parameter space results in J₂ 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.

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):

$$\begin{eqnarray}&&\sin {\rm{\Psi }}=\parallel \displaystyle \frac{({\boldsymbol{k}}\times {\boldsymbol{n}})\times ({\boldsymbol{s}}\times {\boldsymbol{n}})}{| {\boldsymbol{k}}\times {\boldsymbol{n}}| | {\boldsymbol{s}}\times {\boldsymbol{n}}| }\parallel ,\end{eqnarray} \tag{ 26 }$$

where ${\boldsymbol{k}}$, ${\boldsymbol{n}}$, and ${\boldsymbol{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:

$$\begin{eqnarray}&&\cos {\rm{\Psi }}=\displaystyle \frac{({\boldsymbol{k}}\times {\boldsymbol{n}})\cdot ({\boldsymbol{s}}\times {\boldsymbol{n}})}{| {\boldsymbol{k}}\times {\boldsymbol{n}}| | {\boldsymbol{s}}\times {\boldsymbol{n}}| }.\end{eqnarray} \tag{ 27 }$$

Close inspection indicates that ${\rm{\Psi }}={\rm{\Lambda }}=\psi +{\rm{\Omega }}$ (when Ω 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)

$$\begin{eqnarray}&&(\alpha /{s}_{i})\cos {\theta }_{j}\sin {\theta }_{j}+\sin ({\theta }_{j}-{I}_{i})=0,\end{eqnarray} \tag{ 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,⁶ 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 α is the precessional constant:

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{R(\varepsilon )}{\cos \varepsilon }=\displaystyle \frac{3{\kappa }^{2}{M}_{\star }}{{a}^{3}\nu }\displaystyle \frac{{J}_{2}{{Mr}}^{2}}{C}{S}_{0}.\end{eqnarray} \tag{ 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

$$\begin{eqnarray}&&{\theta }_{\mathrm{1,3}}\approx {\tan }^{-1}\left(\displaystyle \frac{\sin {I}_{i}}{1\pm \alpha /{s}_{i}}\right)\end{eqnarray} \tag{ 30 }$$

and

$$\begin{eqnarray}&&{\theta }_{\mathrm{2,4}}\approx \pm {\cos }^{-1}\left(\displaystyle \frac{-{s}_{i}\cos {I}_{i}}{\alpha }\right).\end{eqnarray} \tag{ 31 }$$

Cassini state 1 tends to occur at $| {\theta }_{i}| =\varepsilon \approx 0^\circ$ and state 3 at $| {\theta }_{i}| =\varepsilon \approx 180^\circ$ (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^\circ$); 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 ${\rm{\Psi }}=\psi +{\rm{\Omega }}$: in the vicinity of Cassini state 2, $\psi +{\rm{\Omega }}$ librates about zero; in the vicinity of Cassini state 4, $\psi +{\rm{\Omega }}$ circulates or librates with large amplitude (i.e., $\sim 360^\circ$) in a "horseshoe" trajectory.

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, ϖ, was randomly chosen for each planet, while the longitude of ascending node, Ω, 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, ${{\rm{\Omega }}}_{\mathrm{inv}}$, and ${\varpi }_{\mathrm{inv}}$) and the sky plane (i_sky, ${{\rm{\Omega }}}_{\mathrm{sky}},$ and ${\varpi }_{\mathrm{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).

Table 1.  Initial Conditions for Kepler-62

Planet		b	c	d	e	f
m (M⊕)	${ \mathcal A }$	2.72	0.136	14	6.324	3.648
	${ \mathcal B }$	2.3	0.1	8.2	4.4	2.9
a (au)		0.0553	0.0929	0.12	0.427	0.718
e		0.071	0.187	0.095	0.13	0.094
iinv (°)		0.6138	0.1138	0.1138	0.1662	0.08615
${\varpi }_{\mathrm{inv}}$ (°)		268.175	228.074	241.127	42.615	6.277
${{\rm{\Omega }}}_{\mathrm{inv}}$ (°)		270.002	270.011	270.012	89.992	89.984
Prot (days)						0.25–20
ε (°)						0–90
$\psi +{\rm{\Omega }}$ (°)						0, 90, 180, 270
isky (°)		89.2	89.7	89.7	89.98	89.9
${\varpi }_{\mathrm{sky}}$ (°)		178.175	138.074	151.127	312.615	276.2767
${{\rm{\Omega }}}_{\mathrm{sky}}$ (°)		0	0	0	0	0

Download table as:  ASCIITypeset image

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 +{\rm{\Omega }}=0^\circ$ (in the invariable plane of the system) and one with $\psi +{\rm{\Omega }}=180^\circ$.

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.5{M}_{\mathrm{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 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 ψ we explore briefly in 90° intervals. For the sake of simplicity, we assume the planets 1 and 3 are roughly coplanar.

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.

Table 3.  Eigenvalues for Kepler-62 (Secular Solution)

	Mass set ${ \mathcal A }$		Mass set ${ \mathcal B }$
	gi (arcsec year−1)	si (arcsec year−1)	gi (arcsec year−1)	si (arcsec year−1)
1	4923.54	−4932.41	2960.31	−2966.71
2	659.46	−681.49	409.23	−423.26
3	76.54	−1.56 × 10−5	61.05	2.35 × 10−6
4	53.24	−21.71	38.16	−14.07
5	16.92	−63.27	11.02	−45.17

Download table as:  ASCIITypeset image

Table 4.  Eigenvalues for Kepler-62 (N-body Solution)

	Mass set ${ \mathcal A }$		Mass set ${ \mathcal B }$
	gi (arcsec year−1)	si (arcsec year−1)	gi (arcsec year−1)	si (arcsec year−1)
1	4877.92	−4886.79	2950.05	−2956.55
2	630.38	−681.71	380.74	−423.28
3	75.07	−1.74 × 10−15	58.94	−1.90 × 10−16
4	53.21	−21.70	37.94	−14.06
5	16.90	−63.28	11.32	−45.10

Download table as:  ASCIITypeset image

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₄, though s₅ also contributes strongly. As described below, these two frequencies will be responsible for secular resonances and Cassini states.

Table 5.  Inclination Amplitudes for Each Eigenvalue for Kepler-62 f (Secular Solution)

	I1 (°)	I2 (°)	I3 (°)	I4 (°)	I5 (°)
${ \mathcal A }$	$2.785\times {10}^{-8}$	$1.023\times {10}^{-4}$	$2.698\times {10}^{-4}$	0.1729	0.08692
${ \mathcal B }$	$4.551\times {10}^{-8}$	$1.025\times {10}^{-4}$	$2.838\times {10}^{-4}$	0.1607	0.07788

Download table as:  ASCIITypeset image

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₄ and s₅. As described in Ward & Hamilton (2004), Cassini states 1 and 4 merge and vanish at a critical value of $| \alpha /{s}_{i}|$ (α depends on P_rot). Beyond that critical value, only states 2 and 3 exist, where ${\theta }_{2}\approx 0^\circ$ and ${\theta }_{3}\approx 180^\circ$. Cassini state 3, which is not shown, always has an obliquity of $\approx 180^\circ$, corresponding to a retrograde spin.

Zoom In Zoom Out Reset image size

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 $\psi +{\rm{\Omega }}=0^\circ$. 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 $\sim \cos \varepsilon$. The black curves show the obliquities associated with Cassini state 2 for eigenvalues s₄ and s₅. 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 $\psi +{\rm{\Omega }}=180^\circ$. The black curves correspond to the Cassini state 4 obliquities. The secular resonances appear in the same locations here, but the structure is somewhat different from those in Figure 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 +{\rm{\Omega }}$ circulates. In Figure 8, we show the evolution of the obliquity and $\psi +{\rm{\Omega }}$ at points b, c, and d. Both points b and c are inside the secular resonance associated with frequency s₄, and $\psi +{\rm{\Omega }}$ librates about 0°; however, the obliquity evolution is markedly different. Point c appears to be in or very near Cassini state 2; thus the $\psi +{\rm{\Omega }}$ librates more tightly about 0° and the obliquity oscillates by only $\sim 0\buildrel{\circ}\over{.} 7$. At point b, the obliquity is farther from the Cassini state and oscillates about it with much larger amplitude.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At point d, we also see a large obliquity oscillation. The secular resonance here is associated with the secondary frequency s₅. This resonance is still able to drive the obliquity strongly, though the angle $\psi +{\rm{\Omega }}$ does not clearly librate. This is because s₄ 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 Ω instead of p and q (see Laskar 1986). The longitude of the ascending node, Ω, and the precession angle, ψ, appear together inside sine and cosine functions as ${\rm{\Omega }}+\psi$.

Zoom In Zoom Out Reset image size

The dotted lines in this figure, representing the minimum and maximum values of $R(\varepsilon )$, 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₄, while at point d, it is closer to a cluster of peaks at s₅. So we see that there is indeed a secular resonance at point d, even though the angle $\psi +{\rm{\Omega }}$ does not librate.

Note the difference in the structure of the s₄ secular resonance between Figures 5 and 6: the initial angle $\psi +{\rm{\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 ${ \mathcal B }$. The smaller masses of the planets in this case lower all of the frequencies, as compared with mass set ${ \mathcal A }$ (see Tables 3 and 4). Further, the J₂ 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.

Zoom In Zoom Out Reset image size

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, 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 +{\rm{\Omega }}=0^\circ$, the right for $\psi +{\rm{\Omega }}=180^\circ$. Including fourth-order terms reduces the amplitude of the oscillation (as compared to Brasser 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 +{\rm{\Omega }}$. In the left panel of Figure 11, $\psi +{\rm{\Omega }}=0$, and the planet is near the stable Cassini state 2. In the right panel, $\psi +{\rm{\Omega }}=180$, 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 ψ + Ω for three cases (points a, b, and c).

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 235; 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 $\sim 70^\circ$ to $\sim 110^\circ$ and a minimum that appears below the arc in the left-hand panel.

Zoom In Zoom Out Reset image size

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(\varepsilon )$, 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\buildrel{\circ}\over{.} 71$ (point a in Figure 13). Here, Cassini's third law is satisfied, though the second is not: the obliquity still oscillates by $\sim 10^\circ$.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

Figure 16 shows power spectra for both of the cases. The vertical dashed lines represent the natural precession frequency 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 +{\rm{\Omega }}$ is a natural consequence of the 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.

Zoom In Zoom Out Reset image size

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 $\sim 40^\circ$ to an almost 80° swing.

We also tested different starting precession angles in 90° intervals, to test this parameter's role, still varying e₀ and i₀. The secular resonance in this parameter space is insensitive to ψ and appears in all orientations of the equinox ($\psi \,=11\buildrel{\circ}\over{.} 78,101\buildrel{\circ}\over{.} 78,191\buildrel{\circ}\over{.} 78,281\buildrel{\circ}\over{.} 78$). However, the minimum in Figure 13 (in which point a is located) disappears for $\psi \ne 281\buildrel{\circ}\over{.} 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 \varepsilon$ 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.

Zoom In Zoom Out Reset image size

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 ∼31 arcsec year⁻¹. At points f and g, the precession is close to the peak at ∼9 arcsec year⁻¹. 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^\circ$, 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 +{\rm{\Omega }}$ evolution for 10 Myr in black. The gray curves are the result of slightly different initial conditions—specifically, ${\varepsilon }_{0}=36\buildrel{\circ}\over{.} 72$, a difference of 001. The conditions in the two simulations diverge at ∼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}_{\mathrm{rot}}\lesssim 0.5$ days and ${\varepsilon }_{0}\lesssim 60^\circ$, the initial precession rate of the planet's spin axis is between ∼30 arcsec year⁻¹ and ∼60 arcsec year⁻¹. This places the planet amidst two strong inclination frequencies at ∼31 arcsec year⁻¹ and ∼55 arcsec year⁻¹. Because the precession rate varies strongly with ε 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 ∼9 arcsec year⁻¹.

Zoom In Zoom Out Reset image size

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 \varepsilon$ 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^\circ$, the obliquity oscillations become larger and the inclination frequency peaks broaden, so that the precession rate, $R(\varepsilon )$, 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 +{\rm{\Omega }}$ in Figure 15, and the chaotic motion of ε 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.