Obliquity Variations of Habitable Zone Planets Kepler 62f and Kepler 186f
Abstract
Obliquity variations could play important roles in the climate and habitability of a planet. Orbital modulations caused by planetary companions and the planetary spin axis precession due to the torque from the host star may lead to resonant interactions and cause large amplitude obliquity variability. Here, we consider the spin axis dynamics of Kepler62f and Kepler186f, both of which reside in the habitable zone around their host stars. Using Nbody simulations and secular numerical integrations, we describe their obliquity evolution for particular realizations of the planetary systems. We then use a generalized analytic framework to characterize regions in parameter space where the obliquity is variable. We find that the locations of variability are finetuned over the planetary properties and system architecture in the lower obliquity regimes (). As an example, the obliquities of both Kepler62f and Kepler186f are stable below , whereas the high obliquity regions () allow moderate variabilities, assuming they are Earth analogues. However, if the rotation period of Kepler62f or Kepler186f differs from Earthlike, the lower obliquity regions could become more variable. Even small mutual inclinations () could stir up variability amplitudes up to . Extra undetected planetary companions and/or the existence of a satellite could also destabilize the low obliquity regions. In all cases, the high obliquity region allows for moderate variations.
Keywords: Exoplanets: dynamics – Exoplanets: habitability – Methods: numerical – Methods: analytical
1. Introduction
The rapidly growing arsenal of exoplanet detections has greatly improved our understanding on the occurrence, orbital and structural properties of planetary systems (e.g., Lissauer et al., 2014; Winn & Fabrycky, 2015). In particular, the NASA Kepler mission has discovered thousands of planetary candidates, and identified terrestrial planets in the Habitable Zone (HZ) of their host stars (e.g., Borucki et al., 2013; Barclay et al., 2013; Quintana et al., 2014; Torres et al., 2015). The HZ is conventionally defined as the region where liquid water may exist on the surface of the planets with atmospheres similar to that of the Earth (Kasting et al., 1993; Kopparapu et al., 2013). The Kepler Habitable Zone Working Group has provided a list of habitable zone planets, based on various HZ boundaries and planetary radii, which includes 104 planetary candidates in the optimistic HZ and 20 planets with radii less than 2 in the conservative HZ (Kane et al., 2016). Many of such systems contain multiple planets, and their dynamical interactions could play important roles in determining the habitability of such systems. Kane et al. (2016) analyzed the dynamical stability of the habitable multiplanetary systems, which also serves to validate the planetary systems.
In addition to orbital stability, planetary obliquity and its variations also play major roles in the habitability of a planet. Obliquity (or axial tilt) measures the angle between a planet’s spin and orbital axes. These values are known for planets in the solar system. However, no reliable values for exoplanets have been claimed to date. Methods to infer the spin axis direction from light curves via the effect of rotational flattening on transit depth and infrared phase curves have been proposed (e.g. Gaidos & Williams, 2004; Carter & Winn, 2010). With very high quality light curve data and sophisticated modelling, obliquity measurements may be possible for the most favourable exoplanets in the future.
The obliquity determines the latitudinal distribution of solar radiation on a planet and affects the modulation of its climate (Williams & Kasting, 1997; Chandler & Sohl, 2000; Jenkins, 2000; Spiegel et al., 2009). According to the Milankovitch theory, the ice ages on the Earth are closely associated with the variation in insolation at high latitudes, which depends on the orbital eccentricity and orientation of the spin axis (Weertman, 1976; Hays et al., 1976; Imbrie, 1982; Berger et al., 1992). At present, the obliquity variation of the Earth is regular and only undergoes small oscillations between and with a 41000 year period (Vernekar, 1972; Laskar & Robutel, 1993). This is not to say obliquity instability is wholly incompatible with life – based on a simple 1D energybalance atmospheric model, Armstrong et al. (2014) suggest that large and frequent obliquity variations could help maintain higher surface temperatures and extend the outer boundary of the traditional HZ.
Evolution in a planet’s obliquity is governed by the orbital perturbation from its companion planets, as well as the torque from the host star and any moons acting on the planetary spin axis (e.g., Laskar et al., 1993). When the perturbing frequencies from the companion planets match with the precession frequency due to the torques, the obliquity variation amplitude will increase due to resonant interactions. For instance, without the Moon, the torquing frequency of the Sun matches that from the companion planets, and the obliquity variation of the hypothetical Earth is large (though constrained between ) over billion year timescales (Laskar et al., 1993; Lissauer et al., 2012; Li & Batygin, 2014). With the Moon, the obliquity variation of the Earth is significantly suppressed (Laskar & Robutel, 1993). For an early Venus with less atmospheric tides, the obliquity variation in the low obliquity range is small (Barnes et al., 2016). On the other hand, the obliquity of Mars can vary with large amplitudes (Ward, 1973; Laskar & Robutel, 1993; Touma & Wisdom, 1993). The obliquity changes of Mars likely resulted in runaway condensation of in the atmosphere, rendering Mars uninhabitable (Toon et al., 1980; Fanale et al., 1982; Pollack & Toon, 1982; Francois et al., 1990; Nakamura & Tajika, 2003; Soto et al., 2012).
Here, we will consider the obliquity variation of Kepler62f and Kepler186f, which are Earthlike habitable planets orbiting around a K2Vtype star and a M1type star, respectively (Borucki et al., 2013; Quintana et al., 2014). Both of these two planets stay far from their host stars, where the tidal interactions of their host star are weak. In particular, Bolmont et al. (2014, 2015) and Shields et al. (2016) demonstrated, using Nbody simulations including tidal interactions based on the ‘Constant Time Lag’ model, that Kepler62f and Kepler186f have not yet evolved to be tidally locked, which allow them to keep high obliquities and short rotation periods.
In this article, we focus on the spin axis variability of the HZ planets on shorter timescales. We start by assuming all planets in the system have been detected and study the evolution of the two 5planet systems. We relax the mutual inclination assumption, and include a wide range of planetary system parameters consistent with observational constraints. We present a secular analytical framework applicable in the situation of small orbital eccentricities and inclinations. Such a framework is powerful because it allows us to visualize and make predictions for the nature of obliquity variations in a large parameter space as well as examine the sensitivity of conclusions to errors in the observed parameters. Specifically, we consider different planetary rotation rates, additional planets and satellites, and we characterize regions in this parameter space where the resonant interactions between the HZ planet and its companions may cause large obliquity variations. We briefly discuss the prospects of longterm obliquity evolution subjected to the gradual but inevitable tidal synchronization process.
The paper is organized as follows: planetary system properties used throughout the paper are explained in section §2. In section §3, we use Nbody simulations coupled with secular integration to illustrate the evolution of the obliquity. We consider the obliquity variation as a product of resonant interactions using an analytical approach and compare with our numerical results in section §4. Section §5 considers the effects of undetected planets, satellites, and the possible path to tidal synchronization. A summary is presented in section §6.
2. Planetary System Parameters
We anchor our analysis on the measured properties of the two Kepler systems with potentially habitable planets, and infer the rest of the relevant parameters under assumptions discussed below. Parameters and their representative values used in this work are given in Table 1.
Planet  Period,  Inclination,  Radius,  Mass, 

(days)  ()  ()  
Kepler62  
b  5.715  89.2 0.4  1.31 0.04  
c  12.442  89.7 0.2  0.54 0.03  
d  18.164  89.7 0.3  1.95 0.07  
e  122.387  89.98 0.02  1.61 0.05  
f  267.291  89.90 0.03  1.41 0.07  
Kepler186 

b  3.887  88.9 0.7  1.19 0.08  
c  7.267  89.3 0.4  1.38 0.09  
d  13.343  89.4 0.3  1.55 0.11  
e  22.408  89.7 0.2  1.41 0.10  
f  129.946  1.17 0.11  

Notes: , and for Kepler62 are obtained from Borucki et al. (2013). for Kepler186 are obtained from Quintana et al. (2014). and for Kepler186 are calculated based on the updated stellar mass and radius from Torres et al. (2015) and the planetary system parameters (impact parameters, orbital periods and transit depths) measured by Quintana et al. (2014). The planetary masses for both Kepler62 and Kepler186 are estimated using Chen & Kipping (2017)’s publicly available code Forecaster.
Most direct observables in exoplanet systems are combinations of properties of the planet and the host star. Therefore, any indeterminacy in the stellar properties is directly propagated into errors in the planet properties. For instance, the radius of a planet and its semimajor axis are derived from the transit depth and orbital period in concert with the stellar radius () and mass (), respectively. We set the mass and radius of Kepler62 to be and (Borucki et al., 2013), and the mass and radius of Kepler186 to be and , respectively (Torres et al., 2015). Measurements of orbital periods () from the transiting technique tend to be highly accurate. The semimajor axes () can in turn be calculated from and using Kepler’s Third Law. We adopt the semimajor axes for the Kepler62 planets from Borucki et al. (2013), and we calculated the semimajor axes for Kepler186 based on the updated stellar mass from Torres et al. (2015) and the orbital period measurements from Quintana et al. (2014). The age of the Kepler62 system is determined to be Gyr (Borucki et al., 2013), though recently Morton et al. (2016) arrived at a much younger age of Gyr. Kepler186 is estimated to be aged Gyr (Quintana et al., 2014).
The lineofsight orbital inclinations () are computed from the measured impact parameters (), as well as the stellar radius () and semimajor axes.
(1) 
We obtain the for planets in Kepler62 from Borucki et al. (2013) directly, and we calculate those for Kepler186 based on the measurements of the values in Quintana et al. (2014) and the updated stellar parameters from Torres et al. (2015). Although the line of sight inclinations may not directly translate into the inclination of the planets measured from a reference plane, the fact that the planets transit means they are most likely arranged in nearly coplanar configurations. We use as minimum orbital inclinations in our numerical study, consistent with previous works (e.g. Bolmont et al., 2014; Shields et al., 2016). While the eccentricities are not well constrained, they are all reasonably close to 0. For simplicity, we assume all orbits are initially circular in the numerical simulations.
Since the planetary masses () have not been directly measured for these systems, we use Chen & Kipping (2017)’s publicly available code Forecaster based on a probabilistic approach that predicts masses for a variety of celestial bodies from their radii. We directly use the published planetary radii () and their errors as given in Table 1 of Borucki et al. (2013) for Kepler62, and we calculate using transit depth measurements from Quintana et al. (2014) and stellar parameters from Torres et al. (2015), propagating both errors for Kepler186. For a given input radius, the output is a probability distribution of the possible masses. Due to ignorance over the planetary composition, the resultant uncertainty in the forecasted masses can be quite large and asymmetric (on order of ). The calculated masses are broadly consistent with those considered for rocky compositions in previous studies (Bolmont et al., 2014, 2015; Shields et al., 2016). We use the median radii and their uncertainty to obtain representative values for the planetary masses and their errors, shown in Table as the median and 68% confidence intervals of the mass posteriors 1. In our Monte Carlo trials, we marginalize over the full mass distributions (see Section 4.4).
3. The Evolution of Obliquity Over Time: Numerical Results
An illustration of the planetary spinorbit misalignment (obliquity) is shown in Figure 1. and denote the angular momentum of the orbit and the planet separately, and the obliquity angle () represents the mutual inclination between and . The inclination and the longitude of node characterize the orientation of the orbital plane, or the direction of . Similarly, the obliquity angle and the longitude of the spin axis determine the orientation of the planetary spin axis with respect to the orbital plane.
The phenomenon of planetary obliquity variations is induced by the torques of the the host star on the equatorial bulge of the planet, as well as by periodic forcing from other planets. Therefore, its dynamics depends on the orbital configuration of all the bodies in the system. In this section, we use a numerical approach to study the obliquity variation. We present the framework for the numerical studies in section §3.1, and numerical results on the obliquity variation of a few representative manifestations for the two 5planet systems are shown in §3.2. A qualitative comparison to existing spindynamical studies on these systems is given in §3.3.
3.1. Numerical Method Framework
To calculate the obliquity evolution of a planet, we first perform Nbody simulations to obtain its orbital evolution. Then we compute the spinaxis dynamics, which is coupled to the planetary orbit, based on the Nbody results. This approach implicitly assumed that the planets’ orbital angular momenta are dominant over that of their spin, which allows us to ignore the feedback of the spin axis to the planetary orbit.
We use the hybrid integrator in the publicly available Mercury Nbody code (Chambers, 1999) to simulate representative manifestations for each planet system for years, stepping in increments of 0.5 days. The initial values of semimajor axes, eccentricities and mutual inclinations are discussed in Section 2.
The spinorbit coupling is described by the secular Hamiltonian of obliquity variation, which is welldocumented in the literature (e.g. Goldreich & Peale, 1966; Wisdom et al., 1984; Laskar et al., 1993; Neron de Surgy & Laskar, 1997):
(2) 
Here and are the Andoyer canonical variables where and is the longitude of the planet’s spin axis shown in Figure 1 (Andoyer, 1923; Kinoshita, 1972). is the obliquity, and and are functions of and , where represents orbital inclination and is the longitude of ascending node of the planet:
(3) 
(4) 
is the precession coefficient, defined for a given planet as:
(5) 
where and denote the semimajor axis and eccentricity of the planetary orbit around the star. is the angular velocity of the planet’s spin. is the dynamical ellipticity, where , and are the moment of inertia along the three principle axes. is related to the oblateness (flattening) of the planet and generally scales with (Lambeck, 1980). For a moonless Earth (i.e. only considering the torque from the Sun), (Laskar et al., 1993).
As the internal structure and the dynamical ellipticity of the exoplanets are unknown, we assume the dynamical ellipticity is the same as the Earth if the planet rotates with the same period as the Earth.
(6) 
Note . This gives and for a planet with a rotation period of 24 hours, scaling the precession coefficient with the host stellar mass and the semimajor axis of the planet. This approach yields similar results to that of Lissauer et al. (2012) (see also Quarles et al., 2017), who assumed that the Moment of Inertia coefficient () of the planet is the same as the Earth, where is the planetary mass, and is the equatorial radius of the planet. In particular, with a rotation period of 24 hours, and based on the approach by Lissauer et al. (2012), for planets with masses and radii assumed in Table 1. We set and as our default, Earthlike values. For rotation periods deviating from that of the Earth, we scale the ’s with the inverse of the planet’s rotation period (), as in Eqn.(6). The effects of varying the rotation period as well as the presence of a satellite are explored in more detail based on an analytical approach in section §4.
3.2. Obliquity Evolution for Variable Rotation Periods, and Orbital Inclinations
We explore the dynamical evolution of the obliquity starting with a wide range of initial values, .
For particular values of the planetary rotation period (i.e. proxy for planetary oblateness), the lower region of the obliquity allows larger variabilities. Examples of such rotation periods are approximately hours () for Kepler62f and hours for Kepler186f (). The resultant obliquity evolutions are shown in the middle panels of Figure 2. We will illustrate in section §4 that the increase in the obliquity variation amplitudes are due to resonant interactions.
In the lowest panels, we consider slightly higher inclinations, requiring the mutual inclinations of the planets to be within , in accordance with studies on the mutual inclination of multiplanetary systems (Fang & Margot, 2012; Fabrycky et al., 2014; Ballard & Johnson, 2016; Moriarty & Ballard, 2016). We set the longitude of ascending nodes of the planets such that all of the planets transit. The maximum inclination of Kepler62f and Kepler186f reach over 10 Myr. We adopt the same precession coefficient (i.e. assume the same rotation periods and interior mass distribution) as those in the middle panels. With slightly higher mutual inclinations, the obliquity variation is enhanced (up to in peaktopeak amplitude), as the perturbation to the planetary orbit is stronger when the mutual inclination is larger. We obtain an analytical expression for the amplitude of variabilities in section §4.2.
3.3. Comparison to Previous Studies
The obliquity evolution for these two exoplanets were previously explored, albeit usually for either a narrow range of initial values and/or for a much longer term under tidal influence. Directly comparable with our study are the 10 Myr results, in which Bolmont et al. (2014) find that Kepler186f’s obliquity is stable if it starts with an Earthlike obliquity () and rotation period (1 day). Bolmont et al. (2015) reach the same conclusion for Kepler62f. These are very consistent with our outcome for the same set of assumptions, as shown in the top panels of Figure 2. However, more complex behaviour is possible with different , rotation periods, and mutual inclinations, as illustrated in the bottom panels in this figure.
During the final preparation of our manuscript, we noticed that Quarles et al. (2017) worked on the similar problem related to the obliquity variation of Kepler62f, using a different approach. Specifically, adopting direct Nbody simulations including spinorbit coupling, Quarles et al. (2017) found that the low obliquity region of Kepler62f is stable, assuming it is an Earthanalogue. This is consistent with our results.
4. The Theory of Obliquity Instablity & Results from an Analytical Framework
The resonant interactions between the torque from the host star acting on the planetary spin axis and the orbital perturbation from the companion planets can be quantified in a straightforward, analytical framework. This can in turn be used to approximate the obliquity evolution for given properties of the orbital architecture and that of the planet of interest. In this section, we summarize the analytical approach to characterize the locations of the resonant regions (section §4.1), and we derive an analytical expression for the size of the resonant region for Kepler62f and Kepler186f (section §4.2). Then, we interpret the numerical results based on the analytical approach (section §4.3). In the end, we characterize the parameter space of obliquity evolution considering the uncertainties in the observational orbital parameters and planetary masses based on the analytical approach (section §4.4). The basis of this analysis lies in the simplified Hamiltonian in Eqn. (2), which is analogous to that of a physical pendulum.
4.1. Location of Resonances
The expected resonant locations can be calculated analytically by identifying the obliquity values that allows resonant coupling between the precession of planetary spin axes and the planetary orbital variations, as discussed extensively in the literacture (e.g., Laskar et al., 1993; Touma & Wisdom, 1993; Lissauer et al., 2012). Specifically, the form of Eqn.(2) is analogous to the Hamiltonian describing a physical pendulum with angular position and angular momentum :
(7) 
The term here is akin to in Eqn. 2, while can be compared with , as and are expected to be on the same order.
A natural frequency for the pendulum system is given by the characteristic rate of variation for :
(8) 
and resonance occurs when a perturbing angular frequency (), coincides with this natural frequency. Comparing to Eqn. (2), we see that is analogous to the natural frequency (it is the frequency of the axial precession). Therefore, for obliquity to vary, this condition amounts to requiring . corresponds to the modal frequencies of the forcing terms and . Since and are functions of orbital inclination and the longitude of ascending node, their characteristic frequencies should follow that of the inclination vector,
(9) 
Of course, only modes with frequency leads to physical values of , and can result in resonant interactions. This equation can be inverted to read . Since, given , there is a corresponding value for every arbitrary , we conclude that for every obliquity value, each mode can induce instability at one finetuned rotation rate. This will become visually apparent in section §4.4.
4.2. Resonant Widths
Resonances are not pointlike – they have a finite extent generally centred around . The width of the resonances determines the variability amplitude of the obliquity. In this section, we calculate the resonant widths analytically and use them to characterize the amplitude of variability. Note that the actual amplitude may not equal the full width of the resonance, as it also depends on the initial longitude of the spin axis, . To derive the width of the resonant zone, we again invoke the similarity between obliquity dynamics and that of the physical pendulum. For the physical pendulum (Eqn.7), the halfwidth of a resonant zone in is (e.g. Li & Batygin, 2014).
For the obliquity variation in a resonant region, can be approximated as the following when the inclination is small:
(10) 
Thus, the halfwidth of the resonant region can be expressed as the following:
(11) 
Note that in both Equations (10) and (11), is measured in radians, and we assume to be near in the resonant region. Thus, the approximation fails when the resonant zone is large (). The width is independent of because it is canceled out in the amplitude calculation. Note that Eqn.(11) fails when , which corresponds to a completely rigid sphere with zero dynamical ellipticity. In this case, no precession is expected (hence no resonance), and the obliquity variation follows that of the orbital inclination.
Figure 3 illustrates the halfwidth of the resonant zones as a function of the resonant obliquity values. The width is much larger in the low resonant obliquity region due to the conversion from to being arccosine. The dashed, solid, and dotted line corresponds to mode amplitudes , , and , respectively. When the modal amplitude is small (e.g. ), the resonant width is also small. Therefore, the fact that a given obliquity value lies in the resonant region does not necessarily imply that its variation amplitude must be large – depending on the associated modal importance, it could very well be confined to within a few degrees. In general, planets whose orbits have a larger mutual inclination with a given planet will induce modes with greater amplitudes on that planet’s inclination. This is responsible for the enhanced variability amplitudes found in the numerical results for the more inclined cases in Figure 2.
4.3. Modal Properties based on Nbody Simulations & Predicted Obliquity Variations
As explained in section §4.1, obliquity resonance occurs where the precession rate of the spin axis coincides with a modal frequency of the orbital inclination. Therefore, the properties of the orbital inclination oscillation modes can be used to predict the behaviour of obliquity variations. Empirically, the modes can be obtained from a Fourier Transform (FT) on the time series of orbital inclination modulated by the ascending node, . Using the Nbody results described in section §3, we calculate the FTs for Kepler62f and Kepler186f in the near coplanar cases, and present them in Figure 4 in black. The locations of peaks correspond to the modal frequency values, , and determine the centres of the obliquity resonances via Eqn.(9). The peak amplitudes reflects the power associated with each mode, , which in turn determines the width of the corresponding obliquity resonance through Eqn.(11).
For Kepler62f, we identify two dominant nonzero frequency modes (amplitude ), with frequencies and , and FT peak amplitudes and for the coplanar configuration. These peak amplitudes reflect the power associated with each mode.
Following Eqn.(9), we can predict the locations of the resonant regions of for given values of the precession coefficient . In particular, for an Earthlike , we would expect Kepler62f to exhibit a resonant region at induced by the mode ( yields an unphysical for this ). Substituting the amplitude of the mode in Eqn.(11), the halfwidth of the resonance is expected to be . With a larger precession coefficient of ( hrs), corresponding to the faster rotator as presented in the middle and lower left panels in Figure 2 for Kepler62f, both modes yield physical resonant regions, at and . In the coplanar case, the resonance halfwidths are and respectively. These analytical results agree well with the numerical results shown in the top and middle left panel of Figure 2, except that the analytical results underestimate the variability at for the faster rotator. This is because the numerical FT does not resolve the modal frequencies, and tends to underestimate the modal amplitudes from the sharp peaks in the Fourier transform. Note that the variability amplitude may not span the full width of the resonance. The extent to which the variability occupies the full width of the resonance depends on the exact initial obliquity, , as well as the initial spin axis longitude, . Thus, we expect the analytical widths to be only a proxy to the variability amplitudes.
Increasing the orbital mutual inclinations tends to boost the amplitude associated with each inclination mode, which leads to enlarged resonant widths (§4.2). This phenomenon is exemplified in the Nbody results corresponding to the bottom panels in Figure 2, where the orbital inclination of Kepler62f reaches . While the modal frequencies remain almost unchanged from that of the near coplanar case ( and ), the modal amplitudes increase to and , respectively (not shown in Figure 4). Thus, according to Eqn.(11), the halfwidths of the resonances increase to and . Again, the location of the resonant regions agree very well with the numerical results shown in section §3. However, in this case the analytical approximation again underestimates the widths, in part due to the limited resolution of the discrete FT, but also because the expression for the halfwidth is no longer a good approximation when the variability is large (i.e. when deviate from in the resonant region).
Since every dominant mode can in theory cause lowobliquity instability for an appropriate value (see section §4.1 and §4.4), we should expect Kepler62f to have two rotation periods that lead to destabilized lowobliquity zones. One of the modes () is responsible for the lower obliquity instability at hr shown in Figure 2. The other mode () results in a similar instability region at hr.
Via this method, for Kepler186 we can only identify one dominant mode (with amplitude ) at yr, where the amplitude is . Then, for an Earthlike , the resonance is located at and the halfwidth of the resonance is , according to the analytical approximation in Eqn.(11). For the slower rotator with , the resonant region is at . The halfwidth in the nearly coplanar case is . For the case with higher mutual inclinations, where the inclination of Kepler186f reaches , the modal amplitude is , corresponding to a halfwidth of . Similar to Kepler62f, the analytical results of the resonant location are consistent with the numerical values (shown in Figure 2), while the analytical approximation slightly underestimates the amplitude of the obliquity variations. Another small peak occurs at and is actually an alias of a more distant mode at (see §4.4 and Table 2).
4.4. Dependence on Planetary Rotation Rates and Observational Uncertainties
In addition to Fourier transforming the Nbody results, one can alternatively estimate the relevant fundamental forcing frequencies of the planetary systems using a completely analytical approach following the LaplaceLagrange (LL) secular theory (Murray & Dermott, 1999), which allows us to estimate the characteristic oscillation frequencies as well as their approximate amplitudes in planet ’s orbital plane induced by the other planets in the system. This is done by solving for the eigenvalues and eigenvectors of a matrix constructed from the masses and orbital semimajor axes present in the planetary system. The LL approximation is accurate for systems with no mean motion resonances and containing nearly circular and coplanar orbits, conditions that are satisfied here. Such an approach is powerful because it allows us to characterize vast swaths of parameter space and visualize the results.
We include the ’s computed using the LL approach in Figure 4, shown by vertical dashed lines following a colour code for direct comparison to Figure 5. The modal properties are also summarized in Table 2. According to the LL approach, the most dominant modes for Kepler62f are /yr and /yr, with normalized eigenvector elements and respectively.
Kepler62f  

()  ()  ()  ()  
0  0.45  0.096  0.63  90 
10.0  0.32  0.071  0.29  71.8 
38.2  0.62  0.38  1.05  
278  
1857  
Kepler186f  
()  ()  ()  ()  
0  0.45  0.15  0.93  90 
11.6  0.81  0.47  1.29  85.2 
737  0.0027  0.0026  0.0056  
1520  
2407 
Column headings:
: normalized eigenvector elements
: components of the inclination eigenvector that characterize the size of ’s for the near coplanar configuration
: same as , but for the higher mutual inclination () configuration
: centre of obliquity resonance region corresponding to Earthlike
As shown in Eqn.5, the precession coefficient is highly sensitive to the planet’s spin rate, as it determines the oblateness of the planet. To investigate the dependence of obliquity variation as a function of the planetary spin rate, we present the resonant obliquity values as a function of the precession coefficient shown in Figure 5 for Kepler62 and Kepler186, calculated using Eqn.(9). We label the corresponding planetary rotation period in the top xaxis of Figure 5, assuming as discussed in Section 3.1. In general, ignoring the oblateness of the host star, a system of planets will induce nonzero frequency inclination modes. The colored solid lines represent the resonant obliquities corresponding to the five different inclination modes. For a given value of , only a subset of the modes will correspond to physical values of . The valid resonant obliquities appear in Figure 5 as intersections between an vertical and the curves. Chaotic zones, where the obliquity evolution becomes unpredictable, occur when two resonant zones overlap with each other, which becomes more likely where the resonant widths are large and curves are dense.
The values of depend on the stellar mass , as well as the masses and semimajor axes of each planet in the system (). While the orbital period is measured with high precision, the masses are often underconstrained. To investigate the sensitivity of the resonance locations to measurement errors, we conduct Monte Carlo experiments over 1000 system realizations whose stellar and planetary masses were varied within their uncertainties, and obtain a measure of the error in the determinations of the median resonant obliquity. In Figure 5, the median absolute deviation (MAD) of the resonant obliquities are overplotted in dashed lines. The MADs of the modal frequencies only vary within of theirs median value, and the modes with frequencies close to yr and yr (/yr) are also dominant for Kepler62f (Kepler186f) in the Monte Carlo simulations. It appears that the indeterminacy in these curves are relatively small in the log plots.
In Figure 5, we use green solid lines to indicate the precession frequency corresponding to Earth analogues, with rotation periods of 1 day. Both Kepler62f and Kepler186f avoid instability zones at Earthlike obliquity values (rather narrowly for Kepler62f). However, one can expect higher obliquities to undergo mild variations. Any moderate mutual orbital inclination could also enlarge the widths of each resonant zone. These predictions are consistent with the numerical results in Section 3. We also overplot in cyan dashdotted lines values corresponding to the lower panels in Figure 2, and show that they coincide with regions where lower obliquities undergo resonances. The resonance obliquity in this picture indicates that the instability in the low obliquity region is highly finetuned in rotation period.
One rotation rate of interest is that corresponding to the planet’s breakup velocity, which is 0.0586 days for Earth analogues (shown in pink dashdotted lines in Figure 5). Rotation periods shortward of this value are unphysical. If rotating close to breakup, both Kepler62f and Kepler186f intersect with multiple resonant curves at small obliquities. Nevertheless, since the resonant amplitudes are low for these modes (see Table 2), the associated variability would be quite limited. Of course, given the probable ages of the systems and the tidal interactions overtime, it is unlikely that the planets are such extreme rotators today anyway. Higher precession coefficients than the one corresponding to the breakup velocity are allowed for planets with moons (see section §5.2).
Another regime of relevance in the long term is that of synchronized rotation with the orbital periods of 267 days and 130 days for Kepler62f and Kepler186f respectively (see also Section 5.3). The synchronized rotation periods do not correspond to physical resonant obliquity values. As pointed out by Bolmont et al. (2014, 2015), and Shields et al. (2016), planetary obliquities eventually settle to 0, though the rate at which tidal synchronization occurs is dependent on many factors. Kepler62f and Kepler186f have not likely reached this state, but are inevitably marching towards it. Bolmont et al. (2014, 2015) and Shields et al. (2016) have demonstrated that many paths are possible. From Figure 5 it is visually apparent that, enroute to the synchronized states, the spin rate and therefore the precession frequency of the planet declines and allows the system to sweep across the resonant regions. This allows obliquity variations in a wide range of obliquity regions similar to the future evolution of the Earth as discussed by Neron de Surgy & Laskar (1997).
5. Discussion
Precise determination of many properties of the exoplanets remains out of reach with current techniques. As will be discussed below, a planet’s obliquity evolution can be affected by many factors that are not wellconstrained. This includes the presence of additional planets and satellites, as well as the rotation periods and the oblateness of the planets. Therefore, the parameter space governing a planet’s obliquity evolution is large, and to accurately pinpoint the location occupied by a given exoplanet is challenging.
What is possible is to map out representative regions of this parameter space and infer the general behaviours. Such characterization could provide guidance over the range of considerations towards assessing the planet’s obliquity evolution. In this section we investigate the effect of the existence of additional planets, including giant Jupiter/Saturn analogues as well as internal rocky planets. We also consider the influence of a satellite, which plays a critical role in stabilizing the Earth’s own obliquity (Laskar et al., 1993; Neron de Surgy & Laskar, 1997). Finally, we discuss the planets’ longterm obliquity dynamics as it gradually synchronizes its rotation with its orbital period due to tidal interactions.
5.1. Extra Planets
Thus far, our analysis has been anchored on the assumption that all the relevant planets that exist in the systems have been detected. In reality, the sensitivity of the transiting technique falls with increasing planet distance or inclination. The presence of additional planets, especially at high orbital inclinations, can influence the spin dynamics in a given system dramatically through the introduction of additional, potentially strong modes. In this subsection we map the resonant zones for system configurations involving external giant planets (§5.1.1) as well as an internal planet (§5.1.2).
External Giant Planets
The current estimate on the occurrence rate of Jupitermassed planets orbiting at Jupiterlike distances around Mdwarfs is (Clanton & Gaudi, 2016; Meyer et al., 2017). Therefore, while atypical, external giant planets are not wholly uncommon in these systems. Jointly motivated by the architecture of the solar system, we characterize the influence of hypothetical external giant planets on the obliquity evolution of Kepler62f and Kepler186f using the analytic techniques outlined in previous sections, and confirm them with numerical integrations.
As an example, we place a Jupiter and/or Saturn analogue in each system and calculate the resonant obliquities induced by each using the LL theory. In this context, ‘analogue’ signifies a planet with identical mass, orbital inclination, and a similar semimajor axis, in the presentday. For Jupiter and Saturn, the orbital distances are large ( to AU) and their inclinations are small (). Including only a Jupiteranalogue, there is a new dominant inclination variation mode introduced at yr for Kepler62f, and at yr for Kepler186f. The corresponding normalized eigenvector elements for this mode is for Kepler62f and for Kepler186f. These results are qualitatively similar for a range of orbital distances of Jupiter between 3 and 7 AU. Including both the Jupiter and Saturn analogues, there is another mode introduced at yr for Kepler 62f and yr for Kepler 186f. The modes attributed to the Saturnanalogue are much weaker than the modes introduced by Jupiter, with normalized eigenvector elements of for Kepler62f and for Kepler186f.
The LL resonant curves in the (or ) plane are shown in Figure 6 for systems including both the Jupiter and the Saturn analogues where the orbital distances are fixed to that of Solar System’s Jupiter (5.2 AU) and Saturn (9.6 AU). The solid lines represent the default orbital parameters as presented in Table 1, and the dashed lines represent the MAD values based on the Monte Carlo simulations taking into account the observational uncertainties of the planetary masses. An Earthanalogue in such systems would encounter more obliquity instability zones occurring at lower obliquities. For Kepler62f, the resonant regions caused by the dominant mode /yr and the less dominant mode introduced by Saturn /yr are only slightly separated in the low obliquity regime. Thus, one might expect the obliquity evolution of such a planet to be chaotic due to higher probability of overlap between the resonances. The other dominant modes for both Kepler62f and Kepler186f correspond to low frequencies, which implies greater chance for instability and chaos at higher obliquity angles for Earth analogues, as well as more severe impact on lower obliquity regions for slowly rotating planets.
For comparison, we numerically compute the obliquity evolution for Earthanalogues (i.e. hr in systems including both Jupiter and Saturnanalogues, and present them in Figure 7. Consistent with the analytical expectations, the lower obliquity regions of Kepler62f allow larger variabilities due to two closely spaced inclination oscillation modes. The higher obliquity region also exhibits larger amplitude variabilities than the case without the outer giant planets (as shown in Figure 2). This is because the giant planets induce larger mutual inclinations between the planets. For Kepler186f, the variabilities at low obliquities is still low, due to the lack of additional resonant regions in that region. Similar to Kepler62f, high obliquity variation is enlarged by the presence of the giant planets. This is consistent with the conclusions drawn by Quarles et al. (2017), who also investigated the influence of outer giant planets on the spin axis variability of Kepler62f.
Additional Internal Planet to Kepler186f
In addition to distant planetary companions, it is also possible to have internal nontransiting planetary companions. In particular, the separation between Kepler186e and Kepler186f is large, and it is likely to have an extra planet between these two planets, based on accretion disk simulations (Bolmont et al., 2014). Using dynamical simulations, Bolmont et al. (2014) characterized the mass of the extra planet. Considering a planetary mass ranging from to , Bolmont et al. (2014) found that the planet cannot be more massive than the Earth, in order to keep the mutual inclinations between the planets low so as to allow its companions to transit. Here, we adopt the orbital configuration of the extra planet assumed by Bolmont et al. (2014) to study the obliquity evolution of Kepler186f, where AU, and to avoid transiting. We assumed the extra planet to be Earthmassed ().
We adopt the nearly coplanar configuration of Kepler186, as discussed in section 3, and set the precession coefficient for planet f to be , assuming the planet is Earthlike with a rotation period of 24 hours. The obliquity evolution of Kepler186f is shown in Figure 8. The low obliquity region allows some variabilities, since the extra planet introduces a modal frequency at , which leads to a resonant region at .
Comparing with the obliquity variation without the extra planet (upper right panel of Figure 2), the variabilities in the high obliquity region is also larger. This is because the higher inclination of the extra planet also excites the mutual inclination between the planets, which leads to a stronger perturbation to the planetary spinaxis. The increase in the obliquity variation due to the extra planet is consistent with the discussions by Bolmont et al. (2014).
5.2. Presence of a Satellite
The solar system is teeming with moons. However, little is known about moons elsewhere, in part due to the challenges associated with their low expected signal. Thorough searches in transiting exoplanet data have been conducted (e.g. Kipping et al., 2013), but thus far have revealed only one possible exomoon candidate (Teachey et al., 2017). Sasaki & Barnes (2014) suggest that the tidal decay lifetimes of typical large moons around habitable planets of smaller stars tend to be shorter. However, for Kepler62f and Kepler186f, this timescale should exceed the best estimated system ages (Shields et al., 2016).
In any study of planetary obliquity evolution, it is important to consider the possible presence of moons, for they can have a large influence on a planet’s precession coefficient, described by following modified version of Eqn.(5):
(12)  
In the second term, is the mass of the moon. Similarly, the subscript on the orbital elements , , and indicates that these quantities pertain to moon’s orbit around the planet. Each additional moon contributes one such term, and their collective effect is additive.
For the Earth, the precession due to the Moon is about twice that of the Sun. The same is true for Kepler62f and Kepler186f, if we assume each harbours a moon analogous to that of the Earth. That is, a satellite that preserves the mass ratio and orbits at the same fraction of the planetary Hill radius () as that in the EarthMoon system. For a given planetstar pair, a moon for which and boosts the moonless version of by the following factor:
(13) 
assuming and . Here and refer to with and without the moon, respectively. Figure 9 shows contours of in and . The boost is enhanced at low and high values.
The EarthMoon arrangement turns out to be critical for the obliquity stability of the Earth. Without the moon, the Earth sits in a large obliquity resonant zone spanning 0 to (Lissauer et al., 2012; Li & Batygin, 2014). Our moon pushes the Earth away from this hazardous region (Laskar et al., 1993). In general, changes in the lunar parameters could be fortuitous or catastrophic, as illustrated in Figure 5 and 6, which shows that the existence of resonant obliquity regions depends sensitively on the value of . Moreover, a satellite could alter the rotation rate of the planet overtime due to tidal interactions, as our moon has done to the Earth, which would further modify the oblateness of the planet, resulting in different values (see 5.3).
Assuming that Kepler62f and Kepler186f are Earth analogues with rotation period of hours, the existence of a moon increases the torque of the planetary spin axes, and hence increases their precession coefficient. Thus, the existence of a moon could destabilize the obliquity of Kepler62f, if the precession caused by the moon leads to resonant interactions with the inclination oscillation mode of yr. However, the precession coefficient of Kepler186f is already larger than that of the frequency of the dominant inclination oscillation mode. Thus, the existence of a moon will not cause obliquity variations for Kepler186f.
5.3. Longterm Tidal Evolution
Both Kepler62f and Kepler186f are relatively far from their host stars, thereby likely have weak tidal interactions with their host stars. At their suspected system ages ( Gyr for Kepler62 and Gyr for Kepler186, they may still be evolving towards a tidallylocked state (Bolmont et al., 2014, 2015). Using different tidal models, Shields et al. (2016) showed that the rotation of Kepler62f can be synchronized within billion year timescales. For illustration purposes, assuming the planets are Earth analogues with the same love number and tidal time lag and (Lambeck, 1980), the tidal synchronization timescales are Gyr and Gyr for Kepler62f and Kepler186f respectively, assuming a constant timelag tidal model following the equation below (e.g., Hut, 1982; Ogilvie, 2014):
(14) 
where is the planetary rotational angular velocity, is the orbital angular velocity, is the orbital angular momentum, and is the spin angular momentum.
As the rotation rate of the planet evolves, the oblateness of the planet changes, which leads to a varying precession coefficient (). Thus, it is possible for the planet to move across different resonant regions during the tidalsynchronization, as illustrated in Figure 5 and 6. The fact that the tidal timescale for the alignment of the planetary spinaxis is much longer than the obliquity variation timescales in the resonant zones means that longterm tidal effects cannot suppress shortterm obliquity fluctuations. Consequently, the obliquity can still vary due to the resonant interactions. This is similar to the future obliquity evolution of the Earth as discussed in Neron de Surgy & Laskar (1997). In the case when the planets are synchronized, the spin periods are too long to allow resonant interactions for both Kepler62f and Kepler186f. Thus, the obliquities are stable for both Kepler62f and Kepler186f in the synchronized stage.
6. Conclusions
In this article, we have investigated the shortterm obliquity variability of habitable zone planets in two multiplanet transiting systems, Kepler62f and Kepler186f, over a large parameter space of possible planet properties and orbital architectures allowed by observational constraints. Using Nbody simulations coupled with secular spinorbit coupling analysis, we have shown in section §3 that lowobliquity regions of Kepler62f and Kepler186f are stable over Myrtimescales while higher obliquity regions allow small variabilities, assuming the planets are Earth analogues (i.e. same rotation rate and interior structure, obeying Eqn.(6)).
We have also presented an analytical framework to characterize the nature of obliquity instabilities from first principles (§4). The basic elements of the method are as follows:

Present the nature of obliquity instability as arising from resonant interactions between the planetary spin axis and the orbital axis (e.g., Laskar et al., 1993). Wherever there is a match between the spin axis precession and inclination oscillation frequency, obliquity variation could occur (§4.1). Given a forcing frequency, the location at which a resonance occurs can be calculated from Eqn.(9).

Deduce modal frequencies and amplitudes in the orbital inclination vector of the planet of interest (e.g., Kepler62f and Kepler186f). This could be done in two ways:
When applicable, the analytical technique has the decided advantage of being a rapid, straightforward, and transparent way to compute the regions harbouring resonant obliquities over a large parameter space. It provides a visualization for the behaviour of obliquity variablity as these parameters vary. For Kepler62f and Kepler186f, we have shown good agreement between the numerical and analytical approach.
Different planetary spin rates and orbital configurations of the Kepler systems could affect the obliquity variations. For instance, for Kepler62f, the lower obliquity region () can be unstable when the rotation period is hrs or hrs. For Kepler186f, the lower obliquity region () can be unstable when the rotation period is hours. For both planets, instability in the higher obliquity regions () occur for a wider range of rotation periods ( days). The specific values of the rotation period also depend on the properties of the assumed planetary interior structure. In general, the instability in the lower obliquity region is fine tuned, while the higher obliquity region can be unstable for a wider parameter space. The amplitude of variability is dependent on the mutual orbital inclination of the planets in the system. Higher mutual inclinations () can generate moderate () fluctuations in the low obliquity ranges.
Orbital architectures and planet properties are often difficult to measure and/or subjected to update. Our analytical approach enables us to characterize the overall obliquity variations including observational uncertainties, different planetary oblateness (which leads to different precession coefficients), extra planets and the existence of satellites. We find that the observational uncertainties in the stellar mass and in the estimates of the planetary mass do not change our conclusion qualitatively. In investigating the impact of extra planets, we find that Jupiter and Saturnanalogues can induce larger obliquity variations in the lower obliquity range for Kepler62f, assuming an Earth analogue. The obliquity variations of Kepler186f are not strongly affected by the external giant planets. However, an extra planet between Kepler186e and Kepler186f may induce stronger obliquity variations for Kepler186f in the low obliquity region. The existence of large Moons for Kepler62f could destabilize the spin axis of Kepler62f, but they cannot destabilize the spin axis of Kepler186f. Long term tidal interactions between the planet and the host star will synchronize the planetary spin axis, and reduce the oblateness of the planet and its precession coefficient to move it across resonant regions. Thus, one would expect the obliquity of the planets to vary with large amplitude during tidal synchronization, before reaching the obliquitystable regions at synchronization.
Based on a simplified energy balance model, Armstrong et al. (2014) showed that rapid and large obliquity variability can be favourable to life by keeping a planet’s global average temperature higher than it would have been otherwise, thereby systematically extending the outer edge of a host star’s HZ by . However, it is by no means clear whether large obliquity variation is necessarily beneficial to life under all circumstances. For instance, it is believed that large obliquity variation for Mars may cause the collapse of its atmosphere and render Mars inhabitable (e.g., Toon et al., 1980; Soto et al., 2012). At the very least, obliquity variability can have substantially affect transitions between multiple climate states. Recently, Kilic et al. (2017) mapped out the various equilibrium climate states reached by an Earthlike planet as a function of stellar irradiance and obliquity. They find that, in this parameter space, the state boundaries (e.g. between cryo and aquaplanets) are sharp and very sensitive to the climate history of the planet. This suggests that a variable obliquity can easily move the planet across state divisions as well as alter the boundaries themselves, which would translate into a dramatic impact on instantaneous surface conditions as well as longterm climate evolution. The detailed effects on the climate due to obliquity variations still need more investigation. While atmospheric modeling is beyond the scope of this study, our work can help provide input parameters to existing global climate models (GCMs) as another factor influencing the habilitability in a multiplanet system.
Acknowledgments
The authors would like to thank Konstantin Batygin, Jack Lissauer and Billy Quarles for helpful discussions. We are also grateful to Jennifer Yee and Matthew Holman for providing detailed feedback on our draft. This work and GL were supported in part by Harvard William F. Milton Award. YS is supported by a Doctoral Postgraduate Scholarship from the Natural Science and Engineering Research Council (NSERC) of Canada.
Footnotes
 affiliationmark:
 affiliationmark:
 affiliationmark:
 With the ‘Constant Phase Lag’ tidal model and an Earthlike tidal quality factor, Shields et al. 2016 found that Kepler62f could have become synchronized on fewGyr timescales.
 With very high precision light curves, it may become possible to measure exoplanetary oblateness from transit depth variations (Carter & Winn, 2010; Biersteker & Schlichting, 2017) and ingress/egress anomalies (Zhu et al., 2014), though it would be very difficult in general.
 In all our integrations, we started with initial longitude of spinaxis at (i.e. ). The amplitude of variability depends on (e.g., Quarles et al., 2017). Different can produce variation amplitudes that differ by a factor of up to .
 It turns out negative frequencies correspond to prograde planets and positive frequencies to retrograde ones.
 The Nyquist frequency for our Nbody sampling rate (1/1000 yr) is . Therefore, frequencies outside of this window can manifest as aliases folded upon this value.
 The eigenvector element corresponding to a given mode is an indicator for its relative dominance.
 A synchronouslyrotating body could also sustain large obliquity variations if it lies in a chaotic zone (Wisdom et al., 1984).
References
 Andoyer, H. 1923, Cours de mecanique celeste
 Armstrong, J. C., Barnes, R., DomagalGoldman, S., Breiner, J., Quinn, T. R., & Meadows, V. S. 2014, Astrobiology, 14, 277, 1404.3686
 Ballard, S., & Johnson, J. A. 2016, ApJ, 816, 66, 1410.4192
 Barclay, T. et al. 2013, ApJ, 768, 101, 1304.4941
 Barnes, J. W., Quarles, B., Lissauer, J. J., Chambers, J., & Hedman, M. M. 2016, Astrobiology, 16, 487, 1602.00053
 Berger, A., Loutre, M. F., & Laskar, J. 1992, Science, 255, 560
 Biersteker, J. B., & Schlichting, H. 2017, ArXiv eprints, 1708.08990
 Bolmont, E., Raymond, S. N., Leconte, J., Hersant, F., & Correia, A. C. M. 2015, A&A, 583, A116, 1507.04751
 Bolmont, E., Raymond, S. N., von Paris, P., Selsis, F., Hersant, F., Quintana, E. V., & Barclay, T. 2014, ApJ, 793, 3, 1404.4368
 Borucki, W. J. et al. 2013, Science, 340, 587, 1304.7387
 Carter, J. A., & Winn, J. N. 2010, ApJ, 716, 850, 1005.1663
 Chambers, J. E. 1999, MNRAS, 304, 793
 Chandler, M. A., & Sohl, L. E. 2000, J. Geophys. Res., 105, 20737
 Chen, J., & Kipping, D. 2017, ApJ, 834, 17, 1603.08614
 Clanton, C., & Gaudi, B. S. 2016, ApJ, 819, 125, 1508.04434
 Fabrycky, D. C. et al. 2014, ApJ, 790, 146, 1202.6328
 Fanale, F. P., Salvail, J. R., Banerdt, W. B., & Saunders, R. S. 1982, Icarus, 50, 381
 Fang, J., & Margot, J.L. 2012, ApJ, 761, 92, 1207.5250
 Francois, L. M., Walker, J. C. G., & Kuhn, W. R. 1990, J. Geophys. Res., 95, 14761
 Gaidos, E., & Williams, D. M. 2004, New Astroomy, 10, 67
 Goldreich, P., & Peale, S. 1966, AJ, 71, 425
 Hays, J. D., Imbrie, J., & Shackleton, N. J. 1976, Science, 194, 1121
 Hut, P. 1982, A&A, 110, 37
 Imbrie, J. 1982, Icarus, 50, 408
 Jenkins, G. S. 2000, J. Geophys. Res., 105, 7357
 Kane, S. R. et al. 2016, ApJ, 830, 1, 1608.00620
 Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, Icarus, 101, 108
 Kilic, C., Raible, C. C., & Stocker, T. F. 2017, ApJ, 844, 147
 Kinoshita, H. 1972, PASJ, 24, 423
 Kipping, D. M., Hartman, J., Buchhave, L. A., Schmitt, A. R., Bakos, G. Á., & Nesvorný, D. 2013, ApJ, 770, 101, 1301.1853
 Kopparapu, R. K. et al. 2013, ApJ, 765, 131, 1301.6674
 Lambeck, K. 1980, The earth’s variable rotation: Geophysical causes and consequences
 Laskar, J., Joutel, F., & Robutel, P. 1993, Nature, 361, 615
 Laskar, J., & Robutel, P. 1993, Nature, 361, 608
 Li, G., & Batygin, K. 2014, ApJ, 790, 69, 1404.7505
 Lissauer, J. J., Barnes, J. W., & Chambers, J. E. 2012, Icarus, 217, 77
 Lissauer, J. J., Dawson, R. I., & Tremaine, S. 2014, Nature, 513, 336, 1409.1595
 Meyer, M. R., Amara, A., Reggiani, M., & Quanz, S. P. 2017, ArXiv eprints, 1707.05256
 Moriarty, J., & Ballard, S. 2016, ApJ, 832, 34, 1512.03445
 Morton, T. D., Bryson, S. T., Coughlin, J. L., Rowe, J. F., Ravichandran, G., Petigura, E. A., Haas, M. R., & Batalha, N. M. 2016, ApJ, 822, 86, 1605.02825
 Murray, C. D., & Dermott, S. F. 1999, Solar system dynamics
 Nakamura, T., & Tajika, E. 2003, Geophys. Res. Lett., 30, 1685
 Neron de Surgy, O., & Laskar, J. 1997, A&A, 318, 975
 Ogilvie, G. I. 2014, ARA&A, 52, 171, 1406.2207
 Pollack, J. B., & Toon, O. B. 1982, Icarus, 50, 259
 Quarles, B., Barnes, J. W. Lissauer, J. J., & Chambers, J. 2017, submitted
 Quintana, E. V. et al. 2014, Science, 344, 277, 1404.5667
 Sasaki, T., & Barnes, J. W. 2014, International Journal of Astrobiology, 13, 324
 Shields, A. L., Barnes, R., Agol, E., Charnay, B., Bitz, C., & Meadows, V. S. 2016, Astrobiology, 16, 443, 1603.01272
 Soto, A., Mischna, M. A., & Richardson, M. I. 2012, in Lunar and Planetary Institute Science Conference Abstracts, Vol. 43, Lunar and Planetary Institute Science Conference Abstracts, 2783
 Spiegel, D. S., Menou, K., & Scharf, C. A. 2009, ApJ, 691, 596, 0807.4180
 Teachey, A., Kipping, D. M., & Schmitt, A. R. 2017, ArXiv eprints, 1707.08563
 Toon, O. B., Pollack, J. B., Ward, W., Burns, J. A., & Bilski, K. 1980, Icarus, 44, 552
 Torres, G. et al. 2015, ApJ, 800, 99, 1501.01101
 Touma, J., & Wisdom, J. 1993, Science, 259, 1294
 Vernekar, A. D. 1972, in Atmospheric Radiation, 228
 Ward, W. R. 1973, Science, 181, 260
 Weertman, J. 1976, Nature, 261, 17
 Williams, D. M., & Kasting, J. F. 1997, Icarus, 129, 254
 Winn, J. N., & Fabrycky, D. C. 2015, ARA&A, 53, 409, 1410.4199
 Wisdom, J., Peale, S. J., & Mignard, F. 1984, Icarus, 58, 137
 Zhu, W., Huang, C. X., Zhou, G., & Lin, D. N. C. 2014, ApJ, 796, 67, 1410.0361