Obliquity Evolution of the Potentially Habitable Exoplanet Kepler-62f

Abstract

Variations in the axial tilt, or obliquity, of terrestrial planets can affect their climates and therefore their habitability. Kepler-62f is a 1.4 R_⊕ planet orbiting within the habitable zone of its K2 dwarf host star. We perform N-body simulations that monitor the evolution of obliquity of Kepler-62f for 10-million-year timescales to explore the effects on model assumptions, such as the masses of the Kepler-62 planets and the possibility of outer bodies. Significant obliquity variation occurs when the rotational precession frequency overlaps with one or more of the secular orbital frequencies, but most variations are limited to ≲10°. Moderate variations (∼10–20°) can occur over a broader range of initial obliquities when the relative nodal longitude (ΔΩ) overlaps with the frequency and phase of a given secular mode. However, we find that adding outer gas giants on long-period orbits (∼1000 days) can produce large (∼60°) variations in obliquity if Kepler-62f has a very rapid (4 h) rotation period. The possibility of giant planets on long-period orbits impacts the climate and habitability of Kepler-62f through variations in the latitudinal surface flux, where large variations can occur on million year timescales.

1. Introduction

The Kepler mission has discovered >1000 extrasolar planets (Lissauer et al., 2014; Rowe et al., 2014) that represent a broad spectrum of possible worlds. The Kepler planets represent a reservoir of outcomes from planet formation. From this, researchers can explore and test our models of how planets behave. In this study, we focus on the evolution of a planet's axial tilt, or obliquity, which can be modified by the neighboring planets within a system through the gravitational torque that they exert on its equatorial bulge.

Borucki et al. (2013) announced the discovery of a five-planet system orbiting a K2 dwarf, which is a star slightly smaller in mass and radius than our Sun. Interestingly, the outer two planets in this system (Kepler-62e and Kepler-62f) may reside within a region of space where liquid water could exist on their surfaces given a rocky surface and an atmosphere that permits a reasonable greenhouse effect, otherwise known as the habitable zone. This star, Kepler-62, and its planets, Kepler-62b—Kepler-62f, occupy a region of parameter space that is exciting to explore and can have astrobiological implications.

Among the parameters that govern a planet's astrobiological potential is the obliquity, or axis tilt, ψ. Our Earth with a middling ψ _⊕ = 23.44° enjoys moderate seasonal weather variability (Williams and Kasting, 1997; Spiegel et al., 2009). If a planet's obliquity gets too low, then a lack of illumination at high latitudes can lead to polar glaciations. Planets with obliquities above ψ ≈ 56° would experience severe seasons, alternately baking and freezing their poles while their equators receive the least annual illumination on the planet. That high-obliquity regime applies to both present-day Pluto (White et al., 2017; Howard et al., 2017) and Mars at times in the geologic past (Mège and Bourgeois, 2011; Head and Weiss, 2014).

The study of obliquity evolution of the Solar System (SS) planets has a long history, where Ward (1973) showed that the obliquity of Mars undergoes large variations due to the perturbations of the other planets. Later, Touma and Wisdom (1993) showed that these variations are in fact chaotic, and Laskar and Robutel (1993) performed a frequency analysis to the remaining terrestrial planets of the SS. Laskar et al. (1993) showed that the obliquity of Earth would be also chaotic in the absence of the Moon over a range of initial rotation periods and obliquities. Correia et al. (2003) and Correia and Laskar (2003) showed that the past obliquity of Venus was chaotic and may have aided the planet to transition between a prograde and retrograde rotation. Ward and Hamilton (2004) and Hamilton and Ward (2004) used the planetary perturbations to explain the present high obliquity of Saturn (through a resonance with Neptune). The planetary perturbations were also proposed as the mechanism to tilt Uranus during the initial migration stages in the SS (Brunini, 2006), although this scenario was later discarded (Lee et al., 2007). Brasser and Walsh (2011) showed that the martian obliquity during the Noachian era could be outside the chaotic regime.

Lissauer et al. (2012) showed that the obliquity of a hypothetical Earth in the absence of a large Moon typically varies by approximately ±10°, in contrast to (although not in conflict with) prior calculations showing an allowed range of 0° < ψ _⊕ < 85° (Laskar and Robutel, 1993). The difference between these studies is twofold: (1) Lissauer et al. (2012) used a full n-body method where the prior calculations were based on a secular solution, and (2) Lissauer et al. (2012) evolved the SS for ±2 Gyr where the timescale required to explore the full region is longer. Recently, Saillenfest et al. (2019) updated the procedure for using the secular solutions and have highlighted possible applications to exoplanets.

Barnes et al. (2016) analyzed the possible obliquity variations of early Venus, in part, as a possible analogue for habitable exoplanets, but direct studies of exoplanet obliquity variations have been lacking. Recently, Deitrick et al. (2018b) demonstrated semianalytical methods for evaluating obliquity evolution and the connection to exo-Milankovitch cycles (Spiegel et al., 2010). Shields et al. (2016) performed a limited dynamical study of Kepler-62f starting the interior planets from circular orbits while allowing Kepler-62f to begin with a moderate eccentricity. Shields et al. (2016) used their dynamical model in three-dimensional (3D) climate simulations to estimate potential climates for the planet. However, these results depend on the model assumptions assumed in the dynamical simulations, where the orbital architecture of the system is largely incomplete.

In this study, we investigate the potential obliquity variability of Kepler-62f under a range of model assumptions, including the possible existence of outer bodies and two different mass/radius relationships to determine the planetary masses. The outer bodies we consider are analogs of the gas giants in the SS so that our results are comparable with previous studies in terms of dynamics and habitability. Our numerical method builds upon the work of Lissauer et al. (2012) and is described in Section 2, where the initial orbital architectures of our realizations are presented. Section 3 details the results of our simulations and discusses interesting facets from within the possible outcomes and the possible impact on habitability due to flux variations. Section 4 presents the general conclusions that we may draw from this study.

2. Methodology

2.1. Estimating the masses of the observed planets

Many parameters within Kepler-62 are well characterized by Borucki et al. (2013), including the orbital periods and planetary radii, but they were not able to determine masses for any of the five planets. They performed an analysis focused on transit timing measurements within the data that yielded upper limits on the masses, but these upper limits do not exclude any physically plausible compositions. Bolmont et al. (2015) performed calculations for Kepler-62 where they assume that all five planets are rocky to prescribe masses. They give caution to their study especially considering the effects of tides and rotational flattening that are sensitive to the assumed composition. In this work, we assume the smallest 3 or 4 of the planets to be rocky and the largest planet(s) to be more massive but less dense than the rocky planets.

Besides the deficit of knowledge of mass in the known planetary bodies, there is also an incomplete knowledge of the full architecture of this system. The Kepler mission provided nearly 4 years of continuous monitoring, but other planetary bodies may exist in the Kepler-62 system with larger orbital periods. Moreover, planets can also escape detection through a small misalignment of inclination relative to our line of sight. We thus perform many possible realizations in the rotational state of Kepler-62f with and without giant planets on long-period orbits to identify the sensitivity of the obliquity evolution on our model assumptions.

Following Lissauer and de Pater (2013), we estimate the masses of the planets using a mass/radius relationship in the form of a piecewise power law, $M_{p} ∕ M_{\oplus} = ξ {(R_{p} ∕ R_{\oplus})}^{\oplus}$ , where M_p is the mass of the planet, R_p is the radius of the planet, ξ is a scale factor, and $\oplus$ is the power. This power law is divided into three regimes as follows: $\frac{M_{p}}{M_{\oplus}} = ξ {(\frac{R_{p}}{R_{\oplus}})}^{\oplus} \{\begin{matrix} ξ = 1, \oplus = \frac{1}{0.3}; R_{p} < 1 R_{\oplus} \\ ξ = 1, \oplus = \frac{1}{0.27}; 1 R_{\oplus} \leq R_{p} \leq β R_{\oplus} \\ ξ = β^{\frac{1}{0.27} - 1}, \oplus = 1; R_{p} > β R_{\oplus} \end{matrix},$ (1)

where the last regime uses a parameter β that defines the assumed transition point between “rocky” and volatile-rich planets and an assumed slope in the volatile-rich regime that is within the range of those derived from fits to measured values for exoplanets (Weiss and Marcy, 2014; Wolfgang et al., 2016; Chen and Kipping, 2017).

The obliquity evolution of a planet is only of astrobiological interest if it is rocky. Planets in the habitable zone with $R_{p} = 1.41 R_{\oplus}$ may well be rocky, while those with $R_{p} = 1.61 R_{\oplus}$ are probably not (Rogers, 2015). For most of our simulations, we assume $β = 1.41$ (Case A) where we consider Kepler-62e to be volatile rich and $β = 1.61$ (Case B) that allows Kepler-62e to be rocky. Table 1 shows the nominal values in the planet radius from Borucki et al. (2013) and the masses for each planet for each assumed composition of Kepler-62e. Kepler-62d is considered to be volatile rich in both cases, but using a different transition approximately doubles its mass. However, this will likely only modify slightly the perturbations on Kepler-62f because its semimajor axis is ∼5 × greater.

Table 1.

Assumed Radii and Masses for the Kepler-62 Planets

		62b	62c	62d	62e	62f
R_p	(R _⊕)	1.31	0.54	1.95	1.61	1.41
M_p	(M _⊕), Case A	2.71	0.128	4.94	4.05	3.57
M_p	(M _⊕), Case B	2.71	0.128	7.06	5.83	3.57

Note: Assumed masses for the planets considering Kepler-62e to be either volatile-rich (Case A) or “rocky” (Case B).

2.2. Orbital solution from observations

Borucki et al. (2013) provided orbital parameters derived from a statistical analysis using a transit model to compare with the observational data. Using the results of Borucki et al. (2013), we define or derive all the orbital elements to uniquely prescribe a starting orbit and list those parameters in Table 2. Using the observed orbital period ratios, we prescribe for each planet a semimajor axis (a) where the mass of the host star is 0.69 $M_{⊙}$ . We also use the prescribed masses in Table 1 for this calculation, but their inclusion is largely negligible and augments the semimajor axes of the outer two planets by ∼10⁻⁴ AU, where AU denotes an astronomical unit or 149,597,870,700 meters.

Table 2.

Initial Orbital Parameters Used for the Kepler-62 Planets

Planet	Period (days)	a (AU)	e	i (deg.)	ω (deg.)	Ω (deg.)	M (deg.)	T_eq (K)
Kepler-62b	5.714932	0.0552802	0.07071068	0.883	278.130100	90.373788	74.589661	750
Kepler-62c	12.4417	0.0928576	0.08602325	0.560	54.462322	90.473841	327.803830	578
Kepler-62d	18.16406	0.1195013	0.09486833	0.300	288.434940	90.002197	243.096020	510
Kepler-62e	122.3874	0.4263096	0.13000000	0.303	292.619860	90.303058	63.220181	270
Kepler-62f	267.291	0.7176171	0.09433981	0.488	57.994616	90.478155	176.546390	208

Note: Initial values of the orbital elements used in our simulations for the Kepler-62 system. The radiative equilibrium temperature, T _eq, is also provided to delineate each planet relative to its potential habitability. See Borucki et al. (2013) for the true uncertainties in these parameters.

AU, Astronomical Unit.

The results of Borucki et al. (2013) provide the components of the eccentricity vector (ecosω, esinω), which we convert to the dynamicist convention through a sign change in ω. Borucki et al. (2013) also acknowledge that the nominal values provided result in an unstable configuration leading to the ejection of Kepler-62c, the Mars-sized planet. Most of our planets are assigned an eccentricity using the nominal values of the eccentricity vectors given in Borucki et al. (2013). However, we start our simulations choosing the 1σ upper bound value of esinω component for Kepler-62c (−0.07 rather than the nominal value of −0.18) and use the nominal values for ecosω. By choosing the eccentricity of Kepler-62c in this way, we found the orbits in the Kepler-62 system to be dynamically stable up to 100 Myr.

Also, in the dynamicist convention, we measure orbital inclination relative to the line of sight (90° from the sky plane) and define the ascending node (Ω) to be 90° during the time of transit. The observations tell us that all five planets will have nearly identical values of ascending node by the virtue that they all transit. However, we do not know precisely what the values are, and we assume them to be within a degree of the node defining the line of sight. By analyzing the ratios of transit durations for planets in Kepler's multiplanet systems, statistical studies show that the typical mutual inclination of planets within Kepler multiplanet systems is <2.2° (Fang and Margot, 2012; Fabrycky et al., 2014; Ballard and Johnson, 2016; Moriarty and Ballard, 2016). To define this numerically, we add randomly generated values between 0° and 0.5° to the line of sight value (90°). After defining the ascending node for each planet, we determine the orbital inclination relative to the line of sight through the spherical law of cosines. From the parameters defined thus far, we determine the mean anomaly (M) of each planet at the epoch of central transit as given by the data and then find the mean anomalies relative to a common epoch.

2.3. Possible outer bodies

Our best present methods for indirectly detecting exoplanets are inherently biased toward bodies with relatively short orbital periods, but more planets may exist at longer periods (Mills et al., 2019). The existence of these bodies could have an astrobiological impact on exoplanets within the habitable zone without rendering the overall system unstable. To determine the impact on obliquity, we include two sets of analogs to the giant planets, where one is similar to those in the SS and another is drawn from a large number of stability simulations of two planet pairs of Jupiter and Saturn-mass planets. These analogues are identical in mass to their SS giant namesakes. The set that follows the orbital elements of the SS giants (see appendix A in Murray and Dermott, 2000) is scaled in semimajor axis by their orbital period (Table 3) due to the less massive host star.

Table 3.

Initial Orbital Parameters Used for Our Gas Giants

Planet	Period (days)	a (AU)	e	i (deg.)	ω (deg.)	Ω (deg.)	M (deg.)
SS_J	4336.1086	4.598637392	0.04839266	1.3053	−85.8023	100.55615	19.65053
SS_S	10757.994	8.42784611	0.0541506	2.48446	−21.2831	113.71504	−42.48762
SS_U	30707.225	16.95857856	0.04716771	0.76986	96.73436	74.22988	142.26794
SS_N	60223.18	26.57081264	0.00858587	1.76917	−86.75034	131.72169	259.90868
RG_J	1083.6600	1.82455841	0.03328494	0.8819402	51.9209548	93.8755285	210.832984
RG_S	2781.5082	3.42044515	0.03193422	2.286905	234.538802	89.5462986	351.061437

Note: Initial values of the orbital elements used in our simulations using the gas giants of the SS and a selected stable configuration of randomly drawn RGs. The subscripts denote the mass of each planet by the respective analogue within the SS (Jupiter–Neptune).

RG, gas giant pair; SS, Solar System.

The orbital period of our Jupiter-analog (∼4300 days) is quite large and its perturbations on the inner system do not affect stability. Also, a set of giant planets at such distances may not produce the largest changes to the spin evolution of Kepler-62f. To identify the most extreme conditions for the obliquity evolution of Kepler-62f, we need to identify plausible orbital elements for gas giants where the invoked bodies could have escaped detection and the inner system remains stable. Thus, we perform a suite of ∼15,000 simulations using a pair of gas giants (Jupiter- and Saturn-mass) over a range of orbital periods and evaluate whether the system can remain stable for 10 Myr, given the random initial conditions for the pair of gas giants. In these simulations, we use the masses of the inner planets from Case A.

The initial conditions for the Jupiter mass planet are drawn from a uniform distribution in period ranging from 300 to 1600 days, Rayleigh distributions in the eccentricity (σ_e = 0.05) and inclination (σ_i = 1°), a uniform distribution in the ascending node ranging from 85° to 95°, and uniform distributions for the argument of periastron and mean anomaly ranging from 0° to 360°. The initial conditions for the Saturn mass planet are chosen in a similar manner for most of the orbital elements. Instead of drawing from a distribution in period, we use a uniform distribution ranging from 10 to 20 R_H, where R_H = a_J (M_J/(3 $M_{*}$ ))^1/3 represents the Hill Radius, to ensure the stability of the pair of gas giants (Gladman, 1993; Chambers et al., 1996).

From these simulations, we find that systems are stable when we choose a Jupiter analog with an orbital period ∼1000 days with a corresponding Saturn analog outside of mean motion resonance. Thus, we choose 1 Jupiter–Saturn pair to include in our exploration of obliquity with a Jupiter analog (orbital period ∼1084 days) and a Saturn analog separated by ∼11 R_H, which is wide of the 5:2 mean motion resonance. By choosing our setup in this manner, we will have a giant planet architecture that will substantially perturb the inner system without causing a global instability and provide a much larger perturbation to the possible obliquity of Kepler-62f. Using a hypothetical pair of giant planets is important because it provides a broader context to the study of habitability of Kepler-62f through an investigation of very extreme conditions. The initial orbital elements are given in Table 3 for the gas giants drawn from the SS and our randomly drawn gas giant pair (RG).

2.4. Numerical setup for obliquity evolution

To evaluate the obliquity evolution, we use a modified version of the smercury integration package (Lissauer et al., 2012) that has been optimized for determining the extrema in obliquity evolution up to a given integration step and uses the formalism developed in Touma and Wisdom (1993). We define the obliquity as the mutual inclination, through the spherical law of cosines, between the spin axis (i_s , Ω_s) and orbital axis (i, Ω). The nodal difference, ΔΩ = Ω_s – Ω, is set to 0° in a majority of our simulations, where a subset of our five planet systems using the Case B masses begins with ΔΩ = 90°.

The short orbital period of Kepler-62b (∼5.715 days) poses a numerical challenge for evaluating a broad and deep range of parameters. As a result, we limit our simulations to 10 Myr using a timestep (0.286 days) that is 5% of the orbital period of Kepler-62b. One avenue that could be used to reach longer times is to remove the inner two planets (Kepler-62b and Kepler-62c), where they are added to the mass of the host star effectively increasing the J ₂ of the host star. This would allow a larger timestep to be chosen relative to the orbital period of Kepler-62d. However, much of our discussion on the variation of obliquity depends on the secular frequencies of the system and removing the inner planets would shift the relevant frequencies. Farago et al. (2009) used the averaged Hamiltonian of an inner planet to evaluate the orbital evolution of more distant planets on longer timescales, but this method is beyond the scope of our present study. Our numerical code, smercury, does not include possible tidal interactions and thus keeps the rotation period constant throughout the simulation. This is justified because our simulations do not reach the timescales (∼1 Gyr) necessary for tides to be important for Kepler-62f.

Previous works (Laskar and Robutel, 1993; Li and Batygin, 2014a, 2014b; Shan and Li, 2018; Deitrick et al., 2018a) have used the secular solution for obliquity, while more recent studies have used n-body methods that include spin/orbit interactions (Lissauer et al., 2012; Bolmont et al., 2015; Barnes et al., 2016). However, we can make comparisons with historical formalisms that use secular solutions through relevant precession frequencies. Laskar and Robutel (1993) and Shan and Li (2018) use the “precession” constant $α$ measured in arcseconds per year, which is defined as follows: $α = \frac{3 n^{2}}{2 v} \frac{C - A}{C} \approx \frac{3 n^{2}}{2 v} \frac{J_{2}}{\bar{C}},$ (2)

where n represents the mean motion, ν denotes the rotational frequency, $\bar{C}$ relates to the moment of Inertia, and J ₂ is the zonal harmonic related to the flattening due to the rotation. The moment of Inertia and J ₂ are presently unknown, so we use the values assumed in Lissauer et al. (2012), where J ₂ is derived using the Darwin–Radau relationship (see appendix A in Lissauer et al., 2012).

Figure 1 shows how our assumptions on the planet masses in Case A relate to the equatorial radius (R_eq), the derived zonal harmonic (J ₂), and the approximate value of the precession constant, α, and Table 4 provides specific values for select rotation periods. We note that Fig. 1 (bottom panel) shows the precession constant for Kepler-62e (red dashed line) to be very high (>60''/year), even for slow rotation periods (P_rot > 40 h). As a result, the obliquity variations of the planet will likely be small for the rotational parameters that we consider and instead focus on the possible obliquity variations of Kepler-62f.

FIG. 1.

Rotational parameters (equatorial radius, J ₂, α) using the masses in Case A for the Kepler-62 planets as a function of an assumed rotation period in hours. The inset panel shows a zoomed view of Kepler-62f over a range of rotation periods similar to the Earth.

Table 4.

Initial Rotational Parameters Used in Our Simulations

P _rot (hours)	R _eq (km)	J₂	α (''/year)
4	9316.35	0.0331066	166.81803
6	9116.36	0.0137867	104.20265
8	9055.86	0.0076016	76.60628
10	9029.19	0.0048222	60.74512
12	9015.03	0.0033330	50.38312
14	9006.60	0.0024419	43.06447
16	9001.17	0.0018662	37.61329
18	8997.46	0.0014727	33.39277
20	8994.82	0.0011918	30.02705
22	8992.87	0.0009843	27.27958
24	8991.39	0.0008267	24.99394
26	8990.24	0.0007041	23.06249
28	8989.33	0.0006070	21.40866
30	8988.60	0.0005286	19.97652
32	8988.00	0.0004645	18.72423
34	8987.50	0.0004114	17.61988
36	8987.08	0.0003669	16.63868
38	8986.73	0.0003293	15.76111
40	8986.43	0.0002971	14.97155
42	8986.17	0.0002695	14.25739
44	8985.95	0.0002455	13.60831
46	8985.75	0.0002246	13.01579
48	8985.58	0.0002063	12.47275

Note: Initial values for rotation period (P_rot), equatorial radius (R_eq), zonal harmonic (J ₂), and precession constant (α) for Kepler-62f determined from Lissauer et al. (2012). These values are the same for both Cases A and B as the mass of Kepler-62f is unchanged between the two cases.

Significant variation of retrograde obliquities (ψ_o > 90°) takes a long time to develop computationally and dynamically, even under less computationally demanding conditions (Barnes et al., 2016), so this work focuses mainly on the obliquity evolution of prograde (ψ_o ≤ 90°) rotators. Another unknown parameter is the rotation period of the planets, where we explore a wide range (4–24 h) for the runs considering the five-planet system, Kepler-62b—Kepler-62f. We extend this range to 48 h for the cases, including the outer bodies, because longer orbital periods can introduce lower frequencies that can potentially overlap. By taking steps in the rotation period rather than precession frequency (Laskar and Robutel, 1993; Laskar et al., 1993), we seek to identify the larger structures within the parameter space.

2.5. Calculating the surface flux

The potential habitability of an exoplanet is hard to define and often depends on a range of assumed parameters that influence the exchange of energy between the subsurface, atmospheric, and local space environment (Kasting et al., 1993; Kopparapu et al., 2014, 2016; Ramirez and Levi, 2018). Therefore, our study is limited in terms of the energy received at the top layer of the atmosphere, or surface, of Kepler-62f. The atmospheric composition and albedo of Kepler-62f, which are unknown, are necessary to provide realistic estimates of the temperature variations on the surface of the planet. We consider a proxy for potential habitability the surface flux, S_p , as a function of the stellar luminosity ( $L_{*}$ in $L_{⊙}$ ), the instantaneous stellar distance (r in AU), and the solar constant ( $S_{\oplus}$ ) to get the following: $S_{p} = \frac{L_{*}}{L_{⊙}} {(\frac{1 A U}{r})}^{2} S_{\oplus},$ (3)

The instantaneous stellar distance can be obtained through numerical integration of an orbit, but the variation of orbital parameters (semimajor axis and eccentricity) is small on the timescale of a single orbit, so that it can be computed over all values of the true anomaly, f, by $r = \frac{a (1 - e^{2})}{1 + e cos f}$ (4)

However, the surface flux at any point on the planet will vary as a function of latitude on a sphere. To incorporate this effect, we define the daily mean top-of-atmosphere insolation, I_d in W/m ², at any point as $I_{d} = \frac{S_{p}}{π} [η sin δ_{*} sin δ + sin η cos δ_{*} cos δ_{*}],$ (5)

where η is the half-angle of daylight (i.e., a measure of the day length in radians) at a given latitude, δ, and the substellar latitude, $δ_{*}$ = ψ cos(f + Δf), that is determined by the obliquity (ψ), the true anomaly (f), and an orbital-phase offset (Δf), or the offset in orbital phase between periastron and the highest solar declination in the northern hemisphere (Armstrong et al., 2004, 2014; Shields et al., 2016; Kane and Torres, 2017). The orbital-phase offset for any exoplanet is unknown, where we use Δf = 0.25 throughout our work. An orbital phase of zero corresponds to the planet's periastron passage (Kane and Torres, 2017). The half-angle of daylight is computed through the following conditions: $cos η = \{\begin{matrix} - tan δ tan δ_{*}; |δ| = \frac{π}{2} - |δ_{*}| \\ - 1; δ - δ_{*} \leq - \frac{π}{2} o r δ - δ_{*} \leq \frac{π}{2} \\ 1; δ + δ_{*} \geq \frac{π}{2} o r δ + δ_{*} \leq \frac{π}{2} \end{matrix} .$ (6)

following the work of Armstrong et al. (2014). Using Eqs. 3–6, we calculate the latitudinal flux incident on Kepler-62f and average over the orbital phase to identify how the flux changes annually as a function of latitude.

3. Results

Our simulations investigate the obliquity evolution of Kepler-62f for a 10 Myr timescale. We examine two different assumptions on the mass, and thereby, compositions of the five known (transiting) planets, Cases A and B. Planets on long-period orbits have a low transit probability. Therefore, hypothetical outer planets are also included in many of our simulations. In some cases, we add to the transiting planets a set of giant planet analogues drawn from the SS and scaled by period. In others, randomly determined giant planets (RG) are also used with initial conditions drawn from the results of our stability simulations (Section 2.3).

The inclusion of giant planets alters the eccentricity and inclination of Kepler-62f over time. Figure 2 shows these effects using the five-planet system (black), the nine-planet system with the scaled SS giants (blue), and the seven-planet system with the randomly drawn giant planet pair (red). Much of our analysis relates to overlap between axial precession frequencies and the secular orbital frequencies. Sections 3.1 and 3.2 identify where certain frequencies are important in relation to our initial rotational states, and Section 3.3 focuses more on why those frequencies (and other factors) are important. Given our broad range of rotation periods, the range in the respective precession constant will also be large and our figures account for this by using a base-10 logarithmic scale for the precession constant. We highlight the areas where select frequencies may overlap by black curves in Figs. 4, 5, 6, and 9 using the f_j values given in Table 5.

FIG. 2.

Evolution of the eccentricity (top) and inclination (bottom) of Kepler-62f using the planet masses from Case A (black). The evolution of these parameters, including the scaled SS giants (blue) and our random RG (red), is also shown. RG, gas giant pair; SS, Solar System.

Table 5.

Top 10 Frequencies Determined Through Fast Modal Fourier Transform Analysis

	Case A			Case B			Case A + SS			Case A + RG
j	f_j (''/year)	B (deg.)	$γ_{j}$ (deg.)	f_j (''/year)	B (deg.)	$γ_{j}$ (deg.)	f_j (''/year)	B (deg.)	$γ_{j}$ (deg.)	f_j (''/year)	B (deg.)	$γ_{j}$ (deg.)
1	0	0.41910	90.3	0	0.39692	90.3	0	1.57025	107.7	0	1.27783	91.6
2	−59.73	0.08776	93.0	−73.96	0.10279	90.6	−3.57	0.72488	283.4	−87.25	0.38749	260.5
3	−11.73	0.01468	271.7	−89.19	0.01024	251.1	−14.64	0.34227	293.0	−87.67	0.08507	30.1
4	−79.3	0.00331	255.4	−58.74	0.00981	110.0	−4.35	0.16848	316.5	−93.89	0.0337	231.5
5	−40.15	0.00326	110.7	−17.33	0.0053	284.2	−1.00	0.08441	22.6	−14.67	0.01218	270.2
6	−95.08	0.00248	352.4	−15.52	0.00165	232.0	−65.51	0.07116	82.8	−86.82	0.01176	129.3
7	−24.37	0.00085	13.6	−30.75	0.00105	32.5	−39.30	0.04123	126.8	10.62	0.00746	211.1
8	−47.08	0.00034	351.1	−43.51	0.00052	129.4	−64.47	0.01153	170.7	−40.39	0.00516	38.9
9	−4.80	0.00029	211.3	−60.54	0.0002	342.4	−4.04	0.00606	232.9	−91.4	0.00469	225.1
10	−59.04	0.00016	284.7	−75.77	0.00017	323.5	−50.96	0.00467	107.9	−90.97	0.00377	275.0

Note: Top 10 values of the frequency (f_j ), amplitude (B), and the phase ( $γ_{j}$ ) found using the fast modal Fourier transform analysis for our five-planet systems (Cases A and B), as well as for our simulations that include the scaled SS giants (Case A + SS) or a Jupiter–Saturn pair of gas giants (Case A + RG).

3.1. Variations due to the transiting planets

The Kepler-62 system consists of five transiting planets, all whose orbital periods are shorter than the Earth's, and three of the planets have periods shorter than Mercury's. As a result, the orbits are relatively close together and can produce perturbations on neighboring planets that in turn influence the evolution of obliquity. Bolmont et al. (2015), using a tidal model and incorporating general relativity effects over Gyr timescales, showed that the rotation periods of the planets, which are interior to Kepler-62f, can be substantially slowed leading to a state near ψ = 0°. However, the changes to Kepler-62f under the same model experience much smaller changes to its obliquity and rotation period. Thus, we examine a broad range of initial rotation periods (4–24 h) for a prograde Kepler-62f using both of our assumptions for the masses of the planets (Cases A and B) on a 10 Myr timescale.

Figure 4 illustrates the results of these simulations for nominal values of the planetary masses (Case A). Using the range of obliquity variation, Δψ ≡ |ψ_max – ψ_min |, we estimate the most likely values of obliquity variation and how they depend on the initial rotation state of Kepler-62f. The most common Δψ values in Fig. 4 are less than ∼3°, which is roughly similar in the amount of variation for the present-day Earth with the Moon (∼2.4°). There are distinct regions where Δψ ∼ 1° (light gray) is more prevalent than those >1°, but <3° (red).

The largest Δψ values occur at ψ_o = 0° with a 10-h rotation period (purple) and correspond to a precession constant ∼60''/year (Table 4). This amounts to an instantaneous precession period of ∼22,000 years, which is less than the present-day Earth's. For the 10-h rotation period, Δψ decreases as the initial obliquity, ψ_o , increases until ∼16°. The formula for the expected axial precession, Ω_s, is as follows: ${. Ω}_{s} = - f_{j} = α cos ψ,$ (7)

where f_j denotes the modal frequency (Section 3.3). Using Eq. 7, the precession frequency is ∼62.5''/year for ψ = 16° and α ≈ 60''/year. Thus, the decline in obliquity variation depends on the proximity to the expected precession frequency. We show the evolution of obliquity in a representation similar to a phase portrait in Fig. 3, where a fixed point appears ψ ∼ 16°. When Kepler-62f begins with a 9-h rotation period (α ∼ 67.7''/year), the highest variation appears at ψ_o = 28°, where using Eq. 7 we find ${. Ω}_{s}$ ≈ 60''/year and the width around this peak in Δψ is smaller. This trend continues for faster (decreasing) rotation periods over the range that we simulated.

FIG. 3.

Phase-like portrait of the obliquity (ψ in deg.) and the numerical derivative (Δψ/Δt in ''/year) for three initial values of obliquity (ψ_o = 0° – 10°) assuming the masses from Case A and a rotation period of 10 h.

FIG. 4.

Obliquity variations of Kepler-62f using the five-planet system assuming a mass-radius transition at 1.41 R _⊕ (Case A). The color scale denotes the range in obliquity variation, |ψ_max – ψ_min |, obtained for each simulation over a 10 Myr timescale. The left vertical axis marks the value of the precession constant, α, on a logarithmic scale, where the right vertical axis provides the corresponding rotation period (Eq. 2). The black curves (solid and dashed) represent where the respective orbital precession frequencies (f ₂ and f ₃) in Table 5 overlap with the precession constant.

FIG. 5.

Similar to Fig. 4, but assuming a mass-radius transition at 1.61 R _⊕ (Case B). The black curves (solid, dashed, dotted, and dash-dot) represent where the respective orbital precession frequencies (f ₂, f ₃, f ₄, and f ₅) in Table 5 overlap with the precession constant.

FIG. 6.

Similar to Fig. 4 with the addition of four giant planets similar to the SS in mass and orbital architecture but scaled by period (Case A + SS). Even rotation periods beyond 24 h are included to demonstrate the additional variation due to the giant planets at smaller frequencies. As a result, the maximum value of the color scale is changed. The black curves (solid and dashed) represent where the respective orbital precession frequencies (f ₆ and f ₇) in Table 5 overlap with the precession constant.

Figure 5 demonstrates a similar exploration but considers a different set of masses for the two largest transiting planets, Kepler-62d and Kepler-62e (Case B). Since the mass of Kepler-62f remains unchanged, the values for the zonal harmonic, J ₂, and the precession constant, α, also remain unchanged. Naively, we would expect to see the same result as in Fig. 4. The overall features in Fig. 5 are similar to Fig. 4, but the largest variation appears at a shorter rotation period (8 h), where we would expect −73''/year to be the significant frequency due to the larger masses. This is likely caused by an induced precession from the increased mass for Kepler-62e (Section 3.3). Kepler-62d nearly doubles in mass between the two cases, but it is much farther away from Kepler-62f than is Kepler-62e (∼40% mass increase in Case B). Both Figs. 4 and 5 differ from the SS due to the compactness of the system. There is more variation in the similar plots of the SS (Laskar and Robutel, 1993) due to the slower orbital precession of the outer giant overlapping with plausible spin precessions of the inner planets.

3.2. Effects of outer giant planets

The full architecture of the Kepler-62 systems (i.e., number of planets and masses) is largely unknown, so outer perturbers could introduce other precession frequencies. We explore three scenarios where outer giant planets may exist: (1) a scaled version of the SS giant planets, (2) the scaled SS giant planets with their orbital inclinations doubled, and (3) a Jupiter–Saturn pair of planets drawn from stability simulations (Section 2.3). All three scenarios are performed using the masses from Case A, but only the first scenario is performed using the masses from Case B. Also, we include rotation periods beyond 24 h to show the expected variations at lower frequencies (∼20''/year). The precession frequencies when adding the SS giant planets differ from those in the SS because we did not scale masses of the giant planets and the relative semimajor axis between Kepler-62f and the giant planets is substantially different (Murray and Dermott, 2000, Chapter 7).

Including a scaled version of the SS giant planets (using the Case A masses) introduces more variation, Δψ, at higher initial obliquities and longer rotation periods. Figure 6 demonstrates this result; note that the color scale has been adjusted in response. The most common values of Δψ are ∼3–5° (red), which is a little higher than in Fig. 4. The location of the largest variation also changes in response to the perturbations of the outer planets on the orbit of Kepler-62f. The regions of large variation (Δψ > 5°) typically occupy regions of the parameter space with a short rotation period or high (>45°) initial obliquity. However, there are regions with larger variations for rotation periods longer than 24 h. Figure 7 (using the Case B masses) shows similar differences compared with Fig. 5, with larger variations in the same regions of parameter space. Similar to Fig. 6, the most common variation in Fig. 7 increases to ∼5°.

FIG. 7.

Similar to Fig. 6, where the alternate masses are used (Case B + SS). Even rotation periods beyond 24 h are included to demonstrate the large-scale variations of obliquity.

In the second and third scenarios, the differences in large-scale structure are the most interesting, where we use the masses from Case A. Thus, we evaluate only the even rotation periods for the full range from 4 to 48 h. Doubling the inclination of the SS giant planets increases the overall variation, as shown in Fig. 8, where the largest variation increases to ∼42°. This occurs in the low initial obliquity range, but for slower rotation periods (>36 h). Part of this increase can be attributed to the larger range of values possible when the orbital angular momentum vector decouples from the spin vector (Δψ ∼ 2i). One may expect that doubling the inclinations would also increase the potential for large variations in retrograde (ψ_o > 90°). We perform a set of runs exploring initially retrograde obliquities for this scenario and find that for ψ_o ≥ 95°, the obliquity varies up to 5°, which is similar to our prograde (ψ_o ≤ 90°) results (i.e., red points in Fig. 8). When ψ_o starts in the 91–94° range, more substantial variations (Δψ ∼ 26°) can occur, but they quickly decrease with increasing initial obliquity. In the third scenario, the lowest variation in obliquity is ∼3° and the most common value is ∼10°, as shown in Fig. 9. For a very fast rotator (P_rot = 4 h), the obliquity variation can be quite large, up to ∼65° even on the relatively short timescale (10 Myr) of our simulations.

FIG. 8.

Similar to Fig. 6, where the orbital inclinations of the SS giant planets are doubled. The highest obliquity variations (ψ_o < 5° and P_rot = 38 h) exceed the color scale with values of Δψ up to 42°. The range of initial obliquity is expanded to 105° to show the transition from prograde to retrograde rotators, where the retrograde obliquities not shown (ψ_o > 105°) have variations <5°. Even rotation periods are included to demonstrate the large-scale variations of obliquity. As a result, the maximum value of the color scale is changed.

FIG. 9.

Similar to Fig. 6, but considering our randomly determined Jupiter–Saturn pair (RG, Section 2.3) instead (Case A + RG). Even rotation periods are included to demonstrate the large-scale variations of obliquity. As a result, the minimum and maximum values of the color scale are changed, where obliquity variations below 5° are all colored dark gray. The black curves (solid, dashed, dotted, and dash-dot) represent where the respective orbital precession frequencies (f ₂, f ₄, f ₅, and f ₈) in Table 5 overlap with the precession constant.

3.3. Effects from overlapping frequencies

Our simulations (Figs. 4–9) show that particular regions of parameter space are more likely to exhibit larger variations of obliquity, Δψ, relative to other regions. Figure 3 indicates that 60''/year is a particularly important frequency, given that the largest variation occurred for ψ ₀ = 0°, which corresponds to the case ${\dot{Ω}}_{s}$ = α (Eq. 7). To identify the frequencies more precisely, we apply a fast Fourier transform on the inclination vector of Kepler-62f using a 10 Myr data set with samples every 1000 years. Figure 10 illustrates the resulting Fourier spectra considering four of our assumed architectures: (1) the five-planet system with masses from Case A, (2) the five-planet system with masses from Case B, (3) the nine-planet system with the SS giants (Case A + SS), and (4) the seven-planet system with our Jupiter–Saturn pair (Case A + RG). There are several active peaks (power larger than −11, log scale) in each spectrum. We note that the highest peak at 0''/year has been removed from each spectrum.

FIG. 10.

Fourier spectra illustrating the relevant orbital frequencies (''/year) using the inclination vector (icos Ω, isin Ω) of Kepler-62f while using the planet masses from Case A (black, solid) and Case B (black, dashed). The Fourier spectra of Kepler-62f, including the scaled SS giants (blue) and our random RG (red), are also shown. The vertical axis is on a logarithmic scale.

Table 5 shows the top 10 frequencies (f_j ) identified in the time series using Frequency Modified Fourier Transform* (Šidlichovský and Nesvorný, 1996). The rows of Table 5 are ordered by the amplitude of each mode, where the index j refers to the counting of the modes and does not correspond to a particular body. We list for each scenario the frequency (f_j in ''/year), the amplitude (B in degrees), and the phase of the mode (γ_j in degrees). The highest amplitude frequency (j ₁) occurs at 0''/year, which arises from a degeneracy in inclination vectors (Murray and Dermott, 2000). Shan and Li (2018) performed an analytical analysis of Kepler-62f using the Lagrange/Laplace method and found similar frequencies present. Also, the instantaneous precession period can be determined using these frequencies, where the high-frequency terms produce precession periods much faster than the present-day Earth.

Figure 3 shows large variations near 60''/year, and from Table 5 we expect this to occur most strongly at 59.73''/year. Another region of obliquity variation, although much smaller, occurs at α ≈25''/year and ψ_o ≈ 62°. At this location, we find that Ω_s = 11.7''/year, which is approximately equal to the f ₃ frequency for Case A. In Section 3.1, we found that the region of largest variation changed for Case B (Fig. 5) shifting to higher frequencies, which is because Kepler-62e has a larger assumed mass and induces a higher precession frequency. When looking at the Fourier spectra (Fig. 10) and the associated frequencies (Table 5), it is apparent that a shift to higher frequencies has occurred and explains the large-scale differences between our results in Figs. 3 and 5.

Adding the SS giants to Case A (Case A + SS) produces much larger variations in obliquity (Fig. 6) in the low-frequency regime. The region of largest variation (long rotation period, low initial obliquity) can be associated with the f ₃ frequency, where the high-obliquity regions are associated with the f ₂ frequency. The associated frequencies for the scenario where we double the SS giants' inclination are very similar, where the amplitudes (B in deg.) are approximately double. However, f ₁₀ is different with values of −2.57''/year, 0.00824°, and 70.0° for f ₁₀, B, and γ ₁₀, respectively. This distinction is important because Fig. 8 shows two regions for P_rot = 38 h with Δψ ≈ 42°, corresponding to the f ₃ – f ₅ (ψ_o ≈ 0°) and f ₃ – f ₁₀ (ψ_o ≈ 40°) frequency combinations. Our Jupiter–Saturn pair (Case A + RG) is dominated by much higher frequencies and even includes a positive frequency (10.62''/year) that would make a backward rotating (retrograde) Kepler-62f interesting.

The obliquity evolution over the first 1 Myr for an Earth-like (ψ_o = 23.44° and P_rot = 23.934 h) Kepler-62f is shown in Fig. 11 for both prograde (top) and retrograde (bottom) rotators. The evolution of the five-planet system (Case A, black) displays very small variations, where those with the SS giants added (Case A + SS, blue) are more substantial. The evolution, when including our Jupiter–Saturn pair, is much larger (∼5°) in both prograde and retrograde. We note that the frequency of variation between the prograde and retrograde rotators (Case A + RG, red) is slightly different, where this is due to the overlap of slightly different positive and negative orbital frequencies (f ₅ and f ₇; Table 5).

FIG. 11.

Evolution of a prograde (top) and retrograde (bottom) Kepler-62f with Earth-like values in initial obliquity and rotation period while using the planet masses from Case A (black). The evolution of these parameters, including the scaled SS giants (blue) and our random RG (red), is also shown.

The orbital obliquity can change depending on the assumed masses (Case A or Case B), the outer bodies (possible giant planets), and the assumed longitude of the spin node (Ω_s), where ΔΩ = Ω_s – Ω, the difference between the spin node and the orbital node, is the more important quantity. In almost all our simulations, we have assumed that ΔΩ = 0° based upon our previous results in the work of Barnes et al. (2016). Our previous finding showed little change because the phase of the secular frequencies (γ_j ) in the SS was not near the values of ΔΩ that we simulated. From Table 5, we can see that the strongest nonzero frequency (J ₂) has a phase angle near 90°.

Figure 12 shows the variation in obliquity (Δψ) considering a prograde Kepler-62f with an 8-h rotation period. The top and middle panels demonstrate the difference in variation between the five-planet systems (Cases A and B) and their respective nine-planet systems where the SS giant planets are added. Similar peaks are present that correspond to color changes in Figs. 4–7 for an 8-h rotation period. The bottom panel of Fig. 12 shows the changes due to our assumption on ΔΩ. The region of significant obliquity variations occurs near the same initial obliquity (ψ ∼ 16°), but the range of initial obliquity values, ψ_o , is broadened when we consider ΔΩ = 90°. In addition, the variation near ψ_o = 40° is larger.

FIG. 12.

Variation in obliquity, Δψ, for Cases A and B (black) as a function of the initial obliquity for a rotation period of 8 h. Variations when adding the SS giants are also shown in blue. The bottom panel shows the same results from the middle panel (Case B) but adds results of simulations that modify the relative nodal angle by 90° (orange).

3.4. Effects on the potential habitability

The average flux received by the Earth is larger than what Kepler-62f (∼60% less) receives due to differences in orbital distance relative to the difference between the host stars (Eq. 3). Thus, if Kepler-62f is habitable, irrespective of its changes in obliquity, then its atmosphere must be different such that a more significant greenhouse effect is present. Shields et al. (2016) showed, using 3D Global Circulation Models, that a carbon dioxide (CO₂)-dominated atmosphere with five bars of atmospheric pressure would allow Kepler-62f to be considered habitable by present standards. Several studies of the habitability of a planet include many such assumptions that vary between models (Kasting et al., 1993; Kopparapu et al., 2013, 2014, 2016), where we present results that survey the effect of coupled orbital and obliquity variation on the latitudinal surface flux (Williams and Pollard, 2002, 2003; Armstrong et al., 2004, 2014; Shields et al., 2016; Kane and Torres, 2017; Kilic et al., 2017). We focus on the effects of obliquity variation relative to what the modern Earth experiences using our numerical simulations (Case A, Case A + SS, Case A + RG), including some of the resonant cases that induce large obliquity variations. Some of these cases occur for faster rotation rates than Earth, which could be important because others have suggested that faster rotation rates can increase the prospects of habitability for Kepler-62f (Ramirez and Levi, 2018).

First, we examine the flux variations of the modern Earth so that our later results can be contrasted and placed into context. The mean annual flux F_avg , as shown in Fig. 13 (top row), appears largely stratified where the equator receives the bulk (∼400 W/m²) of the radiation and the poles receive substantially less radiation from the Sun (∼175 W/m²). Although the mean annual flux appears roughly constant, the latitudinal flux changes over a yearly cycle (Fig. 13; middle row), where the polar regions can experience the most dramatic effects with differences of ∼3.5 × the mean flux between the summer and winter extremes. Our simulations reproduce the expected obliquity variation of the modern Earth (±1.3° over 41,000 years), which causes regular climatic shifts (Fig. 13, bottom row). In addition, the nonperiodic shifts in the obliquity are seen in the fractional change of the flux (ΔF/F_avg ).

FIG. 13.

Latitudinal surface flux variations of the Earth as a result of obliquity variations over short (50 Kyr) and long (1 Myr) timescales. The top row illustrates the mean annual flux (F_avg ) as a function of time, the middle row identifies the relative change in flux (ΔF = F_max – F_min ) over an orbit, and the bottom row shows the evolution of obliquity for the respective timescales.

3.4.1. Earth-like and resonant spins of Kepler-62f

As noted before, Kepler-62f resides in a relatively more distant orbit than that of the Earth, and we expect the magnitude of the mean annual flux to differ. Figure 14 illustrates Earth-like conditions in terms of the spin state (ψ_o = 23.4°, P_rot = 24 h), where the mean annual flux at the equator (∼170 W/m²) more closely resembles Earth's polar regions (top row). Apart from the difference in magnitude, Figs. 13 and 14 (top rows) appear quite similar in structure. Differences appear to arise when we consider the fractional change (Figs. 13 and 14, middle rows) of the flux in the southern polar region (∼4.25 × compared with ∼3.5 × ), but the absolute differences ΔF (∼300 W/m² compared with ∼590 W/m²) show that the changes between summer and winter can be milder than those experienced at Earth's south pole due to the lower mean annual flux F_avg . However, the minimum flux at the poles for Kepler-62f is much lower than what Earth experiences. For the obliquity evolution (Fig. 14, bottom row), we find the precession period to be similar to Earth, but with more periodic variations. The similarity of precession period of Earth is purely coincidental with the spin precession period for this test case. The more periodic nature of the obliquity evolution, on the contrary, comes from the relatively weak perturbations of the neighboring planets in Kepler-62.

FIG. 14.

Similar to Fig. 13 but considering Kepler-62f (using Case A masses) with an Earth-like initial spin state (ψ_o = 23.4°, P_rot = 24 h). We note that the scale of the color code changes in this figure and subsequent figures.

Figure 3 illustrates the locations where spin/orbit coupling can play a significant role and thereby induce large obliquity variations. We examine, in Fig. 15, the flux variations when Kepler-62f begins with a shorter rotation period (10 h) and near a commensurability with f ₂ (Table 5). The relatively short-term (<50 Kyr) evolution shows the expected result for ψ ≈ 0° where the mean annual flux F_avg is extremely stratified and the poles receive negligible amounts of radiation even when accounting for the yearly variation ΔF. The obliquity slowly increases on an ∼1 Myr timescale to ∼20° from the resonant effect of the spin/orbit interactions, which causes the mean annual flux to increase at the poles up to a maximum (∼70 W/m²) and the difference between summer and winter extremes is ∼210–280 W/m². Figure 15 (middle rows) shows a ringing effect, which illustrates the effects on flux the variations due to eccentricity variations (inset panel). Flux variations due eccentricity are also apparent in the previously discussed Earth-like rotator case (Fig. 14, middle row).

FIG. 15.

Similar to Fig. 14 but considering Kepler-62f (using Case A masses) with an initial spin state near the spin orbit resonance (ψ_o = 0.4°, P_rot = 10 h).

3.4.2. Effects of giant planets on Kepler-62f

From Section 2.3, we demonstrate that the addition of giant planets on longer period orbits increases obliquity variation across the parameter space and introduces new regions at longer rotation periods where large obliquity variation is possible. Here, we examine where the spin/orbit interactions produce significant obliquity variation, including (1) the SS giants (scaled by period to Kepler-62) and (2) a random Jupiter–Saturn pair.

Figure 16 considers a slowly rotating planet (P_rot = 40 h) with a nearly Earth-like initial obliquity (ψ_o = 28°). When the scaled SS giants are included, the mean annual flux (top rows) is initially similar in structure to the Earth-like case discussed in Section 3.4.1, but the poles receive slightly more average radiation (∼100 W/m²) due to the increased obliquity. Over much longer timescales (1.5 Myr), the obliquity decreases to near 0° and the mean annual flux changes dramatically. This becomes important to habitability. In our SS, for instance, long periods at low obliquity contributed to the collapse of Mars' atmosphere, assuming it began with a more substantial CO₂ atmosphere (Forget et al., 2013). There are variations in the seasonal extremes in flux (∼3.5–4.25 × F_avg ) at the poles (Fig. 16, middle rows) that come from the eccentricity variation. Although the fractional change remains high, the magnitude of the change decreases as the obliquity decreases. The obliquity changes by ∼26° over a 1.5 Myr timescale, where the effects of this change (in terms of the flux) can be quite dramatic when including those due to the eccentricity.

FIG. 16.

Similar to Fig. 15 but considering Kepler-62f (using Case A + SS masses) with an initial spin state near the spin orbit resonance (ψ_o = 15°, P_rot = 15 h).

We examine another giant planet case but including the Jupiter–Saturn pair (Case A + RG) to see how the obliquity of a rapid rotator (P_rot = 4 h) evolves and the impact on the flux variations. Figure 17 shows that over a 10 Myr timescale, the obliquity changes (bottom rows) greatly affect the mean annual flux (top rows) and the fractional change in flux (middle rows). Initially, the mean annual flux is similar to the previous cases, but as the obliquity increases (over the first 50 Kyr), the polar regions receive more radiation per year. After ∼3.5 Myr, there is a shift toward even higher obliquity, but the year-to-year variation is less. Another shift occurs at ∼7.5 Myr that pushes the planet into a high-obliquity regime, where the mean annual flux at the poles is relatively high and there are large seasonal variations (ΔF and 400 W/m²). The evolution of these states occurs over millions of years, but there are times of stark transition that may be detrimental in the present view of habitability, if such conditions exist.

FIG. 17.

Similar to Fig. 15 but considering Kepler-62f (using Case A + RG masses) with an initial spin state near the spin orbit resonance (ψ_o = 12°, P_rot = 4 h).

4. Conclusions

We explore possible dynamical states of the Kepler-62 system focusing on the variation of obliquity for the outermost planet, Kepler-62f, due to its high astrobiological interest (Borucki et al., 2013; Ramirez and Levi, 2018). The possible obliquity variations of Kepler-62f depend on many unknowns, including the masses of the five known planets (Case A or Case B), the possible presence of outer bodies (Case A/B + SS or Case A + RG), the magnitude and direction of the spin, and the relative nodal angle, ΔΩ. Each of these assumptions can introduce variations larger than the present-day Earth (including the stabilizing effects of the Moon) on timescales of a few million years, where the largest contributor is the presence of gas giants with larger orbital periods and inclinations that could have escaped detection. We determine the magnitude of these variations using N-body simulations and their relation to the assumed rotation parameters using frequency analysis. The flux received at the top-of-the atmosphere is measured on both short (50 Kyr) and long (1 Myr) timescales for a range of representative cases, where the effects of planetary eccentricity and obliquity evolution are present.

Regions of significant variation in obliquity differ based upon our initial assumptions of the masses (and compositions) of the Kepler-62 planets. Considering a 0° initial obliquity and a 10-h rotation period for Kepler-62f (Case A) produces an ∼20° difference between the highest and lowest attained obliquity, while most other choices for initial obliquity and rotation period are limited to variations <3°. We show the latitudinal surface flux to vary in response to the orbital and obliquity evolution of the system (Fig. 14), resulting in a variation of ∼280 W/m² in the mean annual surface flux at the poles. A similar range of variation in Case B occurs at a higher initial obliquity (ψ_o ∼ 17°) when considering a different mass/radius relationship for the inner five planets, but requires a faster rotation (8-h period) due to an induced precession from Kepler-62e. The range in initial obliquity that produces moderate variations (>10°) can broaden due to overlap between a secular mode and the relative nodal longitude, ΔΩ.

Including a set of giant planets similar to the SS to either case (Case A or Case B) increases the overall obliquity variation, but not substantially for most cases. In order for strong variations (>25°) to occur, the rotation period of Kepler-62f needs to be longer than 24 h and this also depends on the initial obliquity. SS-like giant planets with double their orbital inclinations produce much broader regions of moderate obliquity variation, where strong variations (Δψ ≈ 42°) can occur for specific initial parameters. Obliquity variations (ψ_o > 90°) for the double inclination scenario (which exhibits the largest variations for large prograde obliquities) are large for initial obliquity in the range 95° > ψ_o > 90°, but typically produce relatively low-obliquity variations (Δψ < 5°) for obliquity above 95°.

Our seven-planet systems that include a pair of gas giants on close-in (1.8–3.5 AU) orbits could induce high variations of obliquity (∼66°), but this requires Kepler-62f to be a relatively rapid (<8 h) rotator due to the higher orbital precession frequency of the giant planets. The most common obliquity variation is larger than the present-day Earth–Moon system, but not extremely high (Δψ < 10°). We show the latitudinal surface flux to vary in response to the orbital and obliquity evolution of system with a 4-h rotation period (Fig. 17) resulting in a variation of ∼50–400 W/m² in the mean annual surface flux at the poles. During epochs of high obliquity, the polar regions can receive a substantial surface flux (∼400 W/m²) at the poles for nearly a third of an orbit. The obliquity can transition into different ranges on a 10 Myr timescale and thus dramatically affect the prospects of habitability.

Obliquity variation can have an impact on the potential climates of exoplanets (Williams and Kasting, 1997; Spiegel et al., 2009), where some climate model calculations have been performed specifically for the Kepler-62 system (Bolmont et al., 2015; Shields et al., 2016). Although we do not include a full climate model in our analysis, we find that the obliquity evolution can differ substantially when additional planets on long-period orbits are considered and thereby alter the amount of latitudinal surface flux that a planet receives at various epochs. The amount of obliquity variation, Δψ, can increase substantially and potentially affect the broader conclusions drawn about climates on potentially habitable worlds. Recently, Ramirez and Levi (2018) found that Kepler-62f is one of the three confirmed exoplanets that lie within a zone of habitability for water worlds called the ice cap zone and a fast rotation rate (∼8 h) would be necessary to allow for habitability by most definitions.

Our study probes the Kepler-62 system using the best estimates for the planetary masses and best-known observationally derived orbital elements. However, there remains significant uncertainty in these values when compared with those of the SS planets, which highlights the need for additional observations that could better constrain the system architecture for more robust studies. Upcoming planet surveys (Transiting Exoplanet Survey Satellite, Ricker et al., 2014) have prioritized searching for Earth-mass planets orbiting M dwarfs and it may be some time before another system dynamically similar to Kepler-62 (with a habitable zone planet) is discovered. Kane et al. (2016) produced a categorized catalog of potentially habitable planets, where Kepler-62f is included in all four categories and indicates that it would be an ideal candidate for any observational follow-up program that targets habitable zone exoplanets.

Footnotes

Acknowledgments

The authors thank the anonymous reviewers for constructive insights that greatly enhanced the quality and clarity of the article. In addition, the authors thank Gongjie Li for helpful discussions. They are also grateful to Yutong Shan and Gongjie Li for providing detailed feedback on their draft.

Author Disclosure Statement

No competing financial interests exist.

Funding Information

The authors acknowledge support from the NASA Exobiology Program (Grant No. NNX14AK31G).

Abbreviations Used

References

Armstrong

J.C.

, Barnes

, Domagal-Goldman

, Breiner

, Quinn

T.R.

, and Meadows

V.S.

(2014) Effects of extreme obliquity variations on the habitability of exoplanets. Astrobiology, 14:277–291.

Armstrong

J.C.

, Leovy

C.B.

, and Quinn

(2004) A 1 Gyr climate model for Mars: new orbital statistics and the importance of seasonally resolved polar processes. Icarus, 171:255–271.

Ballard

and Johnson

J.A.

(2016) The Kepler dichotomy among the M dwarfs: half of systems contain five or more coplanar planets. Astrophys J, 816:66.

Barnes

J.W.

, Quarles

, Lissauer

J.J.

, Chambers

, and Hedman

M.M.

(2016) Obliquity variability of a potentially habitable early Venus. Astrobiology, 16:487–499.

Bolmont

, Raymond

S.N.

, Leconte

, Hersant

, and Correia

A.C.M.

(2015) Mercury-T: a new code to study tidally evolving multi-planet systems. Applications to Kepler-62. Astron Astrophys, 583:A116.

Borucki

W.J.

, Agol

, Fressin

, Kaltenegger

, Rowe

, Isaacson

, Fischer

, Batalha

, Lissauer

J.J.

, Marcy

G.W.

, Fabrycky

, Désert

J.M.

, Bryson

S.T.

, Barclay

, Bastien

, Boss

, Brugamyer

, Buchhave

L.A.

, Burke

, Caldwell

D.A.

, Carter

, Charbonneau

, Crepp

J.R.

, Christensen-Dalsgaard

, Christiansen

J.L.

, Ciardi

, Cochran

W.D.

, DeVore

, Doyle

, Dupree

A.K.

, Endl

, Everett

M.E.

, Ford

E.B.

, Fortney

, Gautier

T.N.

3rd, Geary

J.C.

, Gould

, Haas

, Henze

, Howard

A.W.

, Howell

S.B.

, Huber

, Jenkins

J.M.

, Kjeldsen

, Kolbl

, Kolodziejczak

, Latham

D.W.

, Lee

B.L.

, Lopez

, Mullally

, Orosz

J.A.

, Prsa

, Quintana

E.V.

, Sanchis-Ojeda

, Sasselov

, Seader

, Shporer

, Steffen

J.H.

, Still

, Tenenbaum

, Thompson

S.E.

, Torres

, Twicken

J.D.

, Welsh

W.F.

, and Winn

J.N.

(2013) Kepler-62: a five-planet system with planets of 1.4 and 1.6 Earth radii in the habitable zone. Science, 340:587–590.

Brasser

and Walsh

K.J.

(2011) Stability analysis of the martian obliquity during the Noachian era. Icarus, 213:423–427.

Brunini

(2006) Origin of the obliquities of the giant planets in mutual interactions in the early Solar System. Nature, 440:1163–1165.

Chambers

J.E.

, Wetherill

G.W.

, and Boss

A.P.

(1996) The stability of multi-planet systems. Icarus, 119:261–268.

10.

Chen

and Kipping

(2017) Probabilistic forecasting of the masses and radii of other worlds. Astrophys J, 834:17.

11.

Correia

A.C.M

. and Laskar

(2003) Long-term evolution of the spin of Venus. II. Numerical simulations. Icarus, 163:24–45.

12.

Correia

A.C.M.

, Laskar

, and de Surgy

O.N.

(2003) Long-term evolution of the spin of Venus. I. Theory. Icarus, 163:1–23.

13.

Deitrick

, Barnes

, Quinn

T.R.

, Armstrong

, Charnay

, and Wilhelm

(2018a) Exo-Milankovitch Cycles. I. Orbits and rotation states. Astron J, 155:60.

14.

Deitrick

, Barnes

, Bitz

, Fleming

, Charnay

, Meadows

, Wilhelm

, Armstrong

, and Quinn

T.R.

(2018b) Exo-Milankovitch cycles II: climates of G-dwarf planets in dynamically hot systems. Astron J, 155:266.

15.

Fabrycky

, Lissauer

J.J.

, Ragozzine

, Rowe

J.F.

, Steffen

J.H.

, Agol

, Barclay

, Batalha

, Borucki

, Ciardi

D.R.

, Ford

E.B.

, Gautier

T.N.

, Geary

, Holman

M.J.

, Jenkins

J.M.

, Li

, Morehead

R.C.

, Morris

R.L.

, Shporer

, Smith

J.C.

, Still

, and Van Cleve

(2014) Architecture of Kepler's Multi-transiting Systems. II. New Investigations with Twice as Many Candidates. Astrophys J, 790:146.

16.

Fang

and Margot

J.L.

(2012) Architecture of planetary systems based on Kepler data: number of planets and coplanarity. Astrophys J, 761:92.

17.

Farago

, Laskar

, and Couetdic

(2009) Averaging on the motion of a fast revolving body. Application to the stability study of a planetary system. Celest Mech Dyn Astr, 104:291–306.

18.

Forget

, Wordsworth

, Millour

, Madeleine

J.B.

, Kerber

, Leconte

, Marcq

, and Haberle

R.M.

(2013) 3D modelling of the early martian climate under a denser CO₂ atmosphere: temperatures and CO₂ ice clouds. Icarus, 222:81–99.

19.

Gladman

. (1993) Dynamics of systems of two close planets. Icarus, 106:247.

20.

Hamilton

D.P.

and Ward

W.R.

(2004) Tilting Saturn. II. Numerical model. Astrophys J, 128:2510–2517.

21.

Head

J.W.

and Weiss

D.K.

(2014) Preservation of ancient ice at Pavonis and Arsia Mons: tropical mountain glacier deposits on Mars. Planet Space Sci, 103:331–338.

22.

Howard

A.D.

, Moore

J.M.

, Umurhan

O.M.

, White

O.L.

, Anderson

R.S.

, McKinnon

W.B.

, Spencer

J.R.

, Schenk

P.M.

, Beyer

R.A.

, Stern

S.A.

, Weaver

, Young

, Smith

K.E.

, Olkin

, and New Horizons Science

Team

(2017) Present and past glaciation on Pluto. Icarus, 287:287–300.

23.

Kane

S.R.

and Torres

S.M.

(2017) Obliquity and eccentricity constraints for terrestrial exoplanets. Astron J, 154:204.

24.

Kane

S.R.

, Hill

M.L.

, Kasting

J.F.

, Kopparapu

R.K.

, Quintana

E.V.

, Barclay

, Batalha

N.M.

, Borucki

W.J.

, Ciardi

D.R.

, Haghighipour

, Hinkel

N.R.

, Kaltenegger

, Selsis

, and Torres

(2016) A catalog of Kepler habitable zone exoplanet candidates. Astrophys J, 830:1.

25.

Kasting

J.F.

, Whitmire

D.P.

, and Reynolds

R.T.

(1993) Habitable zones around main sequence stars. Icarus, 101:108–128.

26.

Kilic

, Raible

C.C.

, and Stocker

T.F.

(2017) Multiple climate states of habitable exoplanets: the role of obliquity and irradiance. Astrophys J, 844:147.

27.

Kopparapu

R.K.

, Ramirez

, Kasting

J.F.

, Eymet

, Robinson

T.D.

, Mahadevan

, Terrien

R.C.

, Domagal-Goldman

, Meadows

, and Deshpande

(2013) Habitable zones around main-sequence stars: new estimates. Astrophys J, 765:131.

28.

Kopparapu

R.K.

, Ramirez

R.M.

, SchottelKotte

, Kasting

J.F.

, Domagal-Goldman

, and Eymet

(2014) Habitable zones around main-sequence stars: dependence on planetary mass. Astrophys J, 787:L29.

29.

Kopparapu

R.K.

, Wolf

E.T.

, Haqq-Misra

, Yang

, Kasting

J.F.

, Meadows

, Terrien

, and Mahadevan

(2016) The inner edge of the habitable zone for synchronously rotating planets around low-mass stars using general circulation models. Astrophys J, 819:84.

30.

Laskar

, Joutel

, and Robutel

(1993) Stabilization of the earth's obliquity by the moon. Nature, 361:615–617.

31.

Laskar

and Robutel

(1993) The chaotic obliquity of the planets. Nature, 361:608–612.

32.

Lee

M.H.

, Peale

S.J.

, Pfahl

, and Ward

W.R.

(2007) Evolution of the obliquities of the giant planets in encounters during migration. Icarus, 190:103–109.

33.

and Batygin

(2014a) On the spin-axis dynamics of a Moonless Earth. Astrophys J, 790:69.

34.

and Batygin

(2014b) Pre-late heavy bombardment evolution of the Earth's obliquity. Astrophys J, 795:67.

35.

Lissauer

J.J.

and de Pater

(2013) Fundamental Planetary Science: Physics, Chemistry and Habitability. Cambridge University Press, Cambridge, UK.

36.

Lissauer

J.J.

, Barnes

J.W.

, and Chambers

J.E.

(2012) Obliquity variations of a moonless Earth. Icarus, 217:77–87.

37.

Lissauer

J.J.

, Marcy

G.W.

, Bryson

S.T.

, Rowe

J.F.

, Jontof-Hutter

, Agol

, Borucki

W.J.

, Carter

J.A.

, Ford

E.B.

, Gilliland

R.L.

, Kolbl

, Star

K.M.

, Steffen

J.H.

, Torres

(2014) Validation of Kepler's multiple planet candidates. II. Refined statistical framework and descriptions of systems of special interest. Astrophys J, 784:44.

38.

Mège

and Bourgeois

(2011) Equatorial glaciations on Mars revealed by gravitational collapse of Valles Marineris wallslopes. Earth Planet Sci Lett, 310:182–191.

39.

Mills

S.M.

, Howard

A.W.

, Weiss

L.M.

, Steffen

J.H.

, Isaacson

, Fulton

B.J.

, Petigura

E.A.

, Kosiarek

M.R.

, Hirsch

L.A.

, and Boisvert

J.H.

(2019) Long-period giant companions to three compact, multiplanet systems. Astron J, 157:145.

40.

Moriarty

and Ballard

(2016) The Kepler dichotomy in planetary disks: linking Kepler observables to simulations of late-stage planet formation. Astrophys J, 832:34.

41.

Murray

C.D.

and Dermott

S.F.

(2000) Solar System Dynamics. Cambridge University Press, Cambridge, UK.

42.

Ramirez

R.M.

and Levi

(2018) The ice cap zone: a unique habitable zone for ocean worlds. Mon Not R Astron Soc, 477:4627–4640.

43.

Ricker

G.R.

, Winn

J.N.

, Vanderspek

, Latham

D.W.

, Bakos

.Á., Bean

J.L.

, Berta-Thompson

Z.K.

, Brown

T.M.

, Buchhave

, Butler

N.R.

, Butler

R.P.

, Chaplin

W.J.

, Charbonneau

, Christensen-Dalsgaard

, Clampin

, Deming

, Doty

, De Lee

, Dressing

, Dunham

E.W.

, Endl

, Fressin

, Ge

, Henning

, Holman

M.J.

, Howard

A.W.

, Ida

, Jenkins

, Jernigan

, Johnson

J.A.

, Kaltenegger

, Kawai

, Kjeldsen

, Laughlin

, Levine

A.M.

, Lin

, Lissauer

J.J.

, MacQueen

, Marcy

, McCullough

P.R.

, Morton

T.D.

, Narita

, Paegert

, Palle

, Pepe

, Pepper

, Quirrenbach

, Rinehart

S.A.

, Sasselov

, Sato

, Seager

, Sozzetti

, Stassun

K.G.

, Sullivan

, Szentgyorgyi

, Torres

, Udry

, and Villasenor

(2014) Transiting Exoplanet Survey Satellite (TESS). Space Telescopes and Instrumentation 2014:, 9143:914320.

44.

Rogers

L.A

. (2015) Most 1.6 Earth-radius planets are not rocky. Astrophys J, 801:41.

45.

Rowe

J.F.

, Bryson

S.T.

, Marcy

G.W.

, Lissauer

J.J.

, Jontof-Hutter

, Mullally

, Gilliland

R.L.

, Issacson

, Ford

, Howell

S.B.

, Borucki

W.J.

, Haas

, Huber

, Steffen

J.H.

, Thompson

S.E.

, Quintana

, Barclay

, Still

, Fortney

, Gautier III

T.N.

, Hunter

, Caldwell

D.A.

, Ciardi Edna Devore

D.R.

, Cochran

, Jenkins

, Agol

, Carter

J.A

, and Geary

(2014) Validation of Kepler's multiple planet candidates. III. Light curve analysis and announcement of hundreds of new multi-planet systems. Astrophys J, 784:45.

46.

Saillenfest

, Laskar

, and Boué

(2019) Secular spin-axis dynamics of exoplanets. Astron Astrophys, 624:C2.

47.

Shan

and Li

(2018) Obliquity variations of habitable zone planets Kepler-62f and Kepler-186f. Astron J, 155:237.

48.

Shields

A.L.

, Barnes

, Agol

, Charnay

, Bitz

, and Meadows

V.S.

(2016) The Effect of orbital configuration on the possible climates and habitability of Kepler-62f. Astrobiology, 16:443–464.

49.

Šidlichovský

. and Nesvorný

. (1996) Frequency modified Fourier transform and its applications to asteroids. Celest Mech Dyn Astr, 65:137–148.

50.

Spiegel

D.S.

, Menou

, and Scharf

C.A.

(2009) Habitable climates: the influence of obliquity. Astrophys J, 691:596–610.

51.

Spiegel

D.S.

, Raymond

S.N.

, Dressing

C.D.

, Scharf

C.A.

, and Mitchell

J.L.

(2010) Generalized Milankovitch cycles and long-term climatic habitability. Astrophys J, 721:1308–1318.

52.

Touma

and Wisdom

(1993) The chaotic obliquity of Mars. Science, 259:1294–1297.

53.

Ward

W.R

. (1973) Large-scale variations in the obliquity of Mars. Science, 181:260–262.

54.

Ward

W.R.

and Hamilton

D.P.

(2004) Tilting Saturn. I. Analytic model. Astron J, 128:2501.

55.

Weiss

L.M.

and Marcy

G.W.

(2014) The Mass-radius relation for 65 exoplanets smaller than 4 Earth radii. Astrophys J, 783:L6.

56.

White

O.L.

, Moore

J.M.

, McKinnon

W.B.

, Spencer

J.R.

, Howard

A.D.

, Schenk

P.M.

, Beyer

R.A.

, Nimmo

, Singer

K.N.

, Umurhan

O.M.

, Stern

S.A.

, Ennico

, Olkin

C.B.

, Weaver

H.A.

, Young

L.A.

, Cheng

A.F.

, Bertrand

, Binzel

R.P.

, Earle

A.M.

, Grundy

W.M.

, Lauer

T.R.

, Protopapa

, Robbins

S.J.

, and Schmitt

(2017) Geological mapping of Sputnik Planitia on Pluto. Icarus, 287:261–286.

57.

Williams

D.M.

and Kasting

J.F.

(1997) Habitable planets with high obliquities. Icarus, 129:254–267.

58.

Williams

D.M.

and Pollard

(2002) International Journal of Astrobiology, 1:61.

59.

Williams

D.M.

and Pollard

(2003) International Journal of Astrobiology, 2:1.

60.

Wolfgang

, Rogers

L.A.

, and Ford

E.B.

(2016) Probabilistic Mass-radius relationship for sub-Neptune-sized planets. Astron J, 825:19.