As lower-mass stars often host multiple rocky planets, gravitational interactions among planets can have significant effects on climate and habitability over long timescales. Here we explore a specific case, Kepler-62f (Borucki et al. , 2013 ), a potentially habitable planet in a five-planet system with a K2V host star. N -body integrations reveal the stable range of initial eccentricities for Kepler-62f is 0.00 ≤ e ≤ 0.32, absent the effect of additional, undetected planets. We simulate the tidal evolution of Kepler-62f in this range and find that, for certain assumptions, the planet can be locked in a synchronous rotation state. Simulations using the 3-D Laboratoire de Météorologie Dynamique (LMD) Generic global climate model (GCM) indicate that the surface habitability of this planet is sensitive to orbital configuration. With 3 bar of CO 2 in its atmosphere, we find that Kepler-62f would only be warm enough for surface liquid water at the upper limit of this eccentricity range, providing it has a high planetary obliquity (between 60° and 90°). A climate similar to that of modern-day Earth is possible for the entire range of stable eccentricities if atmospheric CO 2 is increased to 5 bar levels. In a low-CO 2 case (Earth-like levels), simulations with version 4 of the Community Climate System Model (CCSM4) GCM and LMD Generic GCM indicate that increases in planetary obliquity and orbital eccentricity coupled with an orbital configuration that places the summer solstice at or near pericenter permit regions of the planet with above-freezing surface temperatures. This may melt ice sheets formed during colder seasons. If Kepler-62f is synchronously rotating and has an ocean, CO 2 levels above 3 bar would be required to distribute enough heat to the nightside of the planet to avoid atmospheric freeze-out and permit a large enough region of open water at the planet's substellar point to remain stable. Overall, we find multiple plausible combinations of orbital and atmospheric properties that permit surface liquid water on Kepler-62f. Key Words: Extrasolar planets—Habitability—Planetary environments. Astrobiology 16, 443–464.

1. Introduction

NASA's Kepler mission (Borucki et al., 2006), launched in 2009, identified more than 4700 transiting planet candidates—over 2300 of which have been confirmed as planets—in its first 5 years of observations1. Recent statistical surveys estimate ∼40% of planetary candidates to be members of multiple-planet systems (Rowe et al., 2014). Analyses also suggest a low false-positive probability of discovery, indicating that the clear majority of these multiple-planet candidates are indeed real, physically associated planets (Lissauer et al., 2012b, 2014). Kepler data also indicate that smaller stars are more likely to host a larger number of planets per star (Swift et al., 2013), and smaller planets are more abundant around smaller stars (Howard et al., 2012; Mulders et al., 2015). These statistical data suggest that systems of multiple small planets orbiting low-mass stars are a major new planetary population, and one that may be teeming with habitable worlds (Anglada-Escudé et al., 2013).

One of the systems of particular interest for habitability orbits Kepler-62 [Kepler Input Catalog (KIC) 9002278, Kepler Object of Interest (KOI) 701], a K2V star. With a mass of 0.69 M ⊙ and a radius of 0.63 R ⊙ , Kepler-62 hosts five planets—Kepler-62b–f—with orbital periods ranging from ∼5 to ∼267 days and radii from 0.54 to 1.95 R ⊕ (Borucki et al., 2013). The inclinations of all five planets are edge-on (i ∼ 89–90°). However, neither their orbital eccentricities nor their planetary obliquities are constrained. The two outermost planets, Kepler-62e (1.61 R ⊕ ) and Kepler-62f (1.41 R ⊕ ) receive 120% and 41% of the amount of flux that Earth receives from the Sun, respectively. Kepler-62e, with an equilibrium temperature of ∼270 K without an atmosphere, is likely to be too hot (with an atmosphere) to have liquid water on its surface, although a significant amount of cloud cover could reflect incident radiation and cool a planet, if it is synchronously rotating (Yang et al., 2013, 2014). In contrast, Kepler-62f sits near the outer edge of the habitable zone, with a relatively low incoming stellar insolation (hereafter “instellation”). However, the planet could avoid freezing with a sufficient greenhouse effect (Kaltenegger et al., 2013), although there are few measured planetary characteristics that could help constrain whether habitability is likely.

Based on planet size and instellation, Kepler-62f is the most likely candidate for a habitable world in this system. While Kepler-62f has a radius of 1.41 R ⊕ (Borucki et al., 2013), its mass is unknown, leaving its density also unconstrained. Recent statistical surveys of a population of Kepler planets with known masses [from companion radial velocity (RV) surveys] found that volatile inventory increases with planet radius (Marcy et al., 2014; Weiss and Marcy, 2014). However, the most recent Bayesian analysis of Kepler planets with RV follow-up data found that planets with radii ≤1.6 R ⊕ are likely sufficiently dense to be rocky (composed of iron and silicates), though this study was limited to close-in planets with orbital periods less than 50 days (Rogers, 2015).

There are many planetary properties that affect the presence of surface liquid water and perhaps life, and in the absence of sufficient observational data, many of the effects of these properties can only be explored using a global climate model (GCM). Previous work on the climate modeling of Kepler-62f found that 1.6–5 bar of CO 2 would yield surface temperatures above the freezing point of liquid water, depending on the mass and surface albedo of the planet (Kaltenegger et al., 2013). However, this work was done using a 1-D radiative-convective atmospheric code and did not include a treatment of clouds beyond a scaling of surface albedo. Additionally, the effect of the orbital architecture and evolution of this multiple-planet system on the climate of Kepler-62f was not explored. The eccentricity of a planet could be pumped to high values in the presence of additional companions in the system (Mardling, 2007; Correia et al., 2012), and the eccentricity of Kepler-62f is poorly constrained. The effect of different rotation rates on atmospheric circulation was also not examined. Changes in planetary obliquity (Ward, 1974; Williams, 1975; Williams and Kasting, 1997; Dobrovolskis, 2013) and eccentricity (Berger et al., 1993, 2006) will affect the seasonality and average instellation received by a planet throughout its orbit, respectively. Simulations have found that a large moon around Kepler-62f could be long lived over timescales necessary for biological evolution and could help stabilize the planet's obliquity and climate (Sasaki and Barnes, 2014), as the Moon has for Earth (Laskar et al., 1993). However, large moons may not be required to stabilize a planet's obliquity over long timescales, given that variations in a planet's obliquity without a moon would be constrained (Lissauer et al., 2012a) and slowly evolving (Li and Batygin, 2014). Large moons could even be detrimental at the outer edge of the habitable zone (Armstrong et al., 2014). Additionally, planetary rotation rate has been shown to affect atmospheric circulation (Joshi et al., 1997; Merlis and Schneider, 2010; Showman and Polvani, 2011; Showman et al., 2011, 2013). Quantifying the effect of orbital and rotational dynamics on the climate of Kepler-62f using a 3-D GCM is therefore crucial for a more accurate assessment of its habitability.

An understanding of the effect of tides on orbital evolution and planetary habitability is vital given the close proximity of low-mass planets orbiting in the habitable zones of their stars, and the likelihood of additional planetary companions in these systems. Tides raised between a star and a close orbiting planet introduce torques on the planet, resulting in changes in semimajor axis and eccentricity (Barnes et al., 2008, 2009). With an orbital period of ∼267 days, tidal effects are expected to be weak but could affect the rotation rate of the planet (Heller et al., 2011). Recent work has shown that planets orbiting in the habitable zones of lower-mass stars may not necessarily be synchronously rotating (Leconte et al., 2015). Bolmont et al. (2014, 2015) found it likely that the rotation period of Kepler-62f is still evolving (and therefore not synchronized with its orbital period) and that it could have a high obliquity. However, they used only one tidal model, and others exist.

Here, we combine constraints obtained with an orbital evolution model with a 3-D GCM to explore the effect of orbital configuration on the climate and habitability of Kepler-62f, for different atmospheric compositions and planetary rotation rates. We calculated the orbital locations of the inner four planets in the Kepler-62 system relative to the position of Kepler-62f, for a range of possible initial eccentricities and longitudes of pericenter for this planet. We then used an n-body model to integrate the orbits of all five planets given different planetary masses, to identify the maximum initial eccentricity possible for Kepler-62f while still maintaining stability within the planetary system. With these results we calculated the rotational and obliquity evolution possible for this planet. The potential habitability as a function of its stable eccentricity was then explored with a GCM, for atmospheres with Earth-like and high CO 2 levels, focusing on the presence of surface liquid water on the planet over an annual cycle.

2. Models

Here we describe the orbital and climate models used for this study. The existing n-body model HNBody (Section 2.1) is used to integrate the orbits of multiple planets in a system orbiting a central, dominant body (in this case a star). Our method for estimating the masses for the planets is described in Section 2.2. In Section 2.3 we describe the models used to explore the tidal evolution of Kepler-62f. The two GCMs used to run climate simulations of the potentially habitable planet Kepler-62f (with output from the n-body model as input) are described in Section 2.4.

2.1. HNBody: Modeling the orbital evolution of multiple-planet systems

The Hierarchical N-Body (HNBody) package (Rauch and Hamilton, 2002), a set of software utilities we employed as our n-body model, integrates the orbits of astronomical bodies governed by a dominant central mass. It is based on the technique of symplectic integration, an n-body mapping method, developed by Wisdom and Holman (1991), and performs standard point mass orbital integrations for a given number of planets in a system.

In the n-body mapping method, the Hamiltonian—a mathematical function used to generate the equations of motion for a dynamical system—is the sum of the Keplerian and the interacting contributions to the motion of orbiting planetary bodies (Wisdom and Holman, 1991). The former describes the Keplerian motion of the orbiting planetary bodies around a central mass (the host star), and the latter describes the gravitational interactions between the planets themselves.

The evolution of the complete Hamiltonian is determined by alternately evolving the Keplerian and interacting parts separately in a sequence of steps leading to new n-body maps of the system, which are composed of the individual Keplerian evolution of the planets and the kicks due to the perturbations the planets deliver to one another. The output from HNBody consists of a series of data files that describe the evolution of selected orbital parameters over time. In the next section we describe the inputs that HNBody requires to generate this information.

2.1.1. HNBody: Model inputs

While we know from transit timing data where each planet in the Kepler-62 system is relative to our line of sight when it transits its host star, the locations of the other four planets at each planet's individual transit time are unconstrained. An accurate integration of their orbits requires the location of all planets at the same epoch.

The Keplerian orbit of a planetary body can generally be described by a set of six parameters that characterize the orbit. The parameters we chose as input to HNBody for each of the five planets in the Kepler-62 system were the semimajor axis a, the orbital eccentricity e, the inclination of the orbit i, the longitude of the ascending node Ω, the longitude of pericenter ω, and the true anomaly f (see Fig. 1 and Table 1 for definitions of these parameters).

FIG. 1. Geometry of the elliptical orbit of a body of mass m 2 around m 1 . The ellipse has semimajor axis a, semiminor axis b, eccentricity e, and longitude of pericenter ω. The true anomaly f denotes the angle subtended by an imaginary line connecting m 1 with the location of m 2 in its orbit and one connecting m 1 with pericenter (the point of closest approach to m 1 ). This assumes that m 1 >> m 2 .

Table 1. Orbital Elements Used to Describe Kepler-62 Planetary Orbits Parameters Definition Initial Values (Kepler-62b, 62c, 62d, 62e, 62f) Semimajor axis a (AU) Half the longest diameter of an orbit 0.0553, 0.0929, 0.120, 0.427, 0.718 Eccentricity e The degree of ellipticity (oval shape) of an orbit (e = 0, circular; e = 1, parabolic) 0.0, 0.0, 0.0, 0.0, 0.0–0.9 Inclination i (°) The angle between the planet's orbital plane and the sky plane 89.2, 89.7, 89.7, 89.98, 89.90 Longitude of ascending node Ω The angle between the sky plane's 0° longitude and the point where the planet passes from in front of to behind the star 0.0, 0.0, 0.0, 0.0, 0.0 Longitude of pericenter ω The angle between a reference direction and the point of closest approach of a planet to its star (pericenter) 0.0, 0.0, 0.0, 0.0, 0.0 → 2π (Δω = π/6) True anomaly f The angle between a radius vector through pericenter and the planet's location on an orbital ellipse around its host star Calculated in the Appendix

While a and i are constrained from transit observations of the Kepler-62 system (Borucki et al., 2013), e, Ω, ω, and f are poorly constrained. Given the close proximities of planets 62b–e (0.05–0.43 AU), which are all within their tidal circularization orbital radii after 4.5 Gyr (Kasting et al., 1993), and the estimated age of Kepler-62 (∼7 billion years, Borucki et al., 2013), we assumed e = 0 for these inner four planets. However, we acknowledge that in a multiple-planet system, the tidal evolution of a close-in planet is coupled to secular interactions with other planets in the system, and these interactions can cause planet eccentricities to vary (Greenberg and Van Laerhoven, 2011; Laskar et al., 2012; Hansen and Murray, 2015).

Following from this assumption, ω is undefined, as all points along the orbit are equidistant from the star. We also assumed Ω = 0 for all planets. This is reasonable given that i ∼ 89–90° for all planets, thus constituting an edge-on orbit capable of yielding a transit observable by Earth. We determined the locations of Kepler-62b–e relative to Kepler-62f at the same point in time, assuming a range of possible eccentricities between 0.0 and 0.9 and a range of longitudes of pericenter ω between 0 and 2π for Kepler-62f. In the Appendix we outline the equations that were used to generate the locations of all planets in the system, using transit data for these planets (Borucki et al., 2013), and the aforementioned values for the other orbital elements.

2.2. Planetary masses of the Kepler-62 system

HNBody also requires a seventh parameter as input—the planet's mass relative to its host star. There are several mass-radius relations in the literature, and we chose two to explore the effect of planetary mass on the orbital evolution of the Kepler-62 system. We used the following mass-radius relation determined by Kopparapu et al. (2014) derived from the exoplanets.org database (Wright et al., 2011):

This yields a planetary mass of ∼3 M ⊕ for Kepler-62f. We also ran HNBody integrations with planet masses derived using the following mass-radius (for Kepler-62d, e, and f) and density ρ (for Kepler-62b and c) relations from Weiss and Marcy (2014):

The Weiss and Marcy (2014) relationship resulted in higher masses for Kepler-62b and f and lower masses for Kepler-62c, d, and e, compared to those using the Kopparapu et al. (2014) treatment. Finally, we ran additional orbit integrations with the maximum mass limits for all planets, as determined by Borucki et al. (2013). Although these maximum masses are likely unphysical, they should result in a similar constraint on eccentricity as more physically plausible mass estimates (such as iron or Earth-like compositions). This varied approach allows us to address the spread in mass-radius relationships that can arise given different planetary compositions (Wolfgang and Lopez, 2015). The masses for all five Kepler-62 planets used as input to HNBody integrations are given in Table 2.

Table 2. Masses (M ⊕ ) Used as Inputs to HNBody for Different Sets of Orbital Integrations of the Planets Kepler-62b–f Kepler-62b Kepler-62c Kepler-62d Kepler-62e Kepler-62f Relation/Source 2.3 0.1 8.2 4.4 2.9 Kopparapu et al. (2014) 2.8 0.1 5.0 4.2 3.7 Weiss and Marcy (2014) 9.0 4.0 14 36 35 Borucki et al. (2013)

Additional model specifications include the preferred coordinate system, the class of particles (based on the scale of the masses and how their interaction is to be taken into account), and the time step and total length of integration. We specified a bodycentric coordinate system (ex. heliocentric reference frame in the case of the Solar System), which treats the system as dominated by the mass of the central star. We ran our HNBody integrations for interactions between “heavy weight particles” (which includes the star and the planets) for 106 years, with a time step equal to 1/20th and 1/100th of the orbital period of the innermost planet in the Kepler-62 system. As Kepler-62b has an orbital period of 5.7 days, we used time steps of 0.29 and 0.06 days, respectively. See Deitrick et al. (2015) for more details on HNBody. Barnes and Quinn (2004) showed that in simulations that last one million orbits, only ∼1% of ejections occurred in the timescale between 105 and 106 years of simulation, with the vast majority of unstable configurations occurring within 104 years. While orbital instabilities can arise at any timescale, one million orbital periods has emerged as a practical reference.

We defined orbital stability as a successful integration in which stable, periodic amplitude oscillation was present for all planets throughout the entire million years, the energy was conserved to better than one part in 104, and no planets were ejected from the system. Eccentricities spanning the full range of stability for Kepler-62f were then used as input to our GCMs. We then ran climate simulations of this planet with a variety of atmospheric compositions, orbital configurations, and rotation rates to explore and assess its possible climate states and to determine the best possible combination of these parameters for surface habitability.

2.3. Tidal model

In this section we consider the tidal evolution of Kepler-62f. Bolmont et al. (2014, 2015) examined the rotational evolution of planets e and f and found that the rotation period of f is not synchronized with the orbital period and that a wide range of obliquities is possible. However, they only considered circular orbits for planet f. Furthermore, they used only one equilibrium tide model in which the time interval between the passage of the perturber and the passage of the tidal bulge is constant, a model we call the constant time lag (CTL) model (Hut, 1981; Ferraz-Mello et al., 2008; Leconte et al., 2010). In this section we relax some of these choices and find that synchronous rotation of Kepler-62f is possible.

In addition to the CTL model, we also employed a model in which the angle between the perturber and the tidal bulge, as measured from the center of the planet, is constant, a model we call the constant phase lag (CPL) model (Goldreich and Soter, 1966; Ferraz-Mello et al., 2008; Heller et al., 2011). These two equilibrium tide models are mathematical treatments of tidal deformation, angular momentum transfer, and energy dissipation due to tidal friction, and both have known flaws (Greenberg, 2009; Efroimsky and Makarov, 2013). However, they provide a qualitative picture of tidal evolution and can be used to infer possible rotation states of Kepler-62f.

Both models are one-dimensional with the tidal properties encapsulated in either the time lag τ for CTL or the tidal quality factor Q for CPL. For both models we assumed that Kepler-62f has the same tidal properties as Earth: τ ⊕ = 640 s (Lambeck, 1977) and Q ⊕ = 12 (Yoder, 1995). The values for Kepler-62f are assuredly different, but these choices are likely not off by more than an order of magnitude, assuming Kepler-62f has oceans and continents. A Q of 100 would result in a tidal locking timescale that is less than 10 times longer.

Rather than present the full set of equations, the reader is referred to Appendix E of Barnes et al. (2013). The salient features are that each model contains six coupled differential equations that track the orbital semimajor axis and eccentricity, as well as the rotation rate and obliquity of both bodies. The effects of the other planets in the system are ignored. Therefore, we consider limiting cases to explore the range of these effects. The equations conserve angular momentum, and energy is dissipated by the tidal deformation of the bodies. In this case, tidal effects on the star are negligible.

We used calculated minimum and maximum initial eccentricities for Kepler-62f (Section 3); initial spin periods of 8 h, 1 day, or 10 days; and initial obliquities of 5°, 23.5°, or 80°. We chose this range of initial conditions for Kepler-62f with the goal of setting boundary conditions for the climate simulations.

2.4. Climate modeling of Kepler-62f

Our primary goal in the climate modeling of Kepler-62f was to identify the most favorable combination of planetary parameters that would result in areas of the planet with warm enough surface temperatures for liquid water. To determine the scenario for which the largest habitable surface area is possible, we varied the atmospheric composition, orbital eccentricity, planetary obliquity, the angle of the vernal equinox (the point in the planet's orbit where both hemispheres receive equal instellation) relative to pericenter (closest approach to the star), and the rotation rate of the planet in our GCM simulations.

2.4.1. Model inputs to CCSM4

We used version 4 of the Community Climate System Model (CCSM42), a fully-coupled, global climate model (Gent et al., 2011). We ran CCSM4 with a 50 m deep slab ocean (see, e.g., Bitz et al., 2012), with the ocean heat transport set to zero, as done in experiments by Poulsen et al. (2001) and Pierrehumbert et al. (2011). The ocean is treated as static but fully mixed with depth. The horizontal resolution is 2°. There is no land, hence we refer to it as an “aqua planet.” The sea ice component to CCSM4 is the Los Alamos sea ice model CICE version 4 (Hunke and Lipscomb, 2008). We made the ice thermodynamic only (no sea-ice dynamics) and use the more easily manipulated sea-ice albedo parameterization from CCSM3, with the surface albedo divided into two bands, visible (λ ≤ 0.7 μm) and near-IR (λ > 0.7 μm). We used the default near-IR and visible band albedos (0.3 and 0.67 for cold bare ice and 0.68 and 0.8 for cold dry snow, respectively). For more details, see Shields et al. (2013).

Because CCSM does not easily allow a 267-day orbital period (the actual orbital period of Kepler-62f), the orbital period was set to 365 days, so the model would still simulate a full annual cycle. As atmospheric radiative and convective adjustment timescales are short compared to either orbital period, we do not expect this to make much difference in the overall climate. We used CCSM4 for continuity with our previous work, to evaluate the general climate behavior given changes in orbital parameters, and used the Laboratoire de Météorologie Dynamique (LMD) Generic GCM to simulate the climate of a planet with physical characteristics more closely like those of Kepler-62f. The details of LMD Generic GCM are given in the following section.

Kepler-62f receives 41% of the modern solar constant from its star, therefore significant amounts of CO 2 or other greenhouse gases may be required to keep temperatures above the freezing point of water on the surface, as is widely assumed for planets near the outer edge of their host stars' habitable zones (Kasting et al., 1993). An active carbon cycle capable of generating increased amounts of atmospheric CO 2 in response to decreasing surface temperatures (Walker et al., 1981) would be a relatively straightforward means of maintaining habitable surface temperatures on the planet. More than ∼2 bar of CO 2 could accumulate in the atmosphere of a planet with an active carbon cycle before the maximum greenhouse limit for CO 2 is reached (Pierrehumbert, 2010). Large increases in atmospheric CO 2 concentration begin to have significant effects on convection and the manner in which it adjusts the temperature lapse rate (the rate at which atmospheric temperature decreases with increasing altitude) on a planet. Additionally, CO 2 condensation becomes likely at levels of 1–2 bar, and collisional line broadening becomes important, increasing the infrared opacity of the atmosphere (Pierrehumbert, 2005). These effects are neglected in Earth-oriented GCMs such as CCSM4. We therefore used the CCSM4 model to simulate only scenarios with an Earth-like atmospheric CO 2 concentration (400 ppmv) and allowed water vapor to vary throughout each simulation according to evaporation and precipitation processes. The rest of the atmospheric composition is preindustrial. We then used an additional GCM—LMD Generic GCM—to simulate the climate of Kepler-62f with Earth-like as well as 1–12 bar of CO 2 , as LMD Generic GCM contains parameterizations for addressing atmospheres with high CO 2 content.

It was important to consider the possibility that Kepler-62f may not have sufficient atmospheric CO 2 to keep surface temperatures above freezing. We therefore explored alternate means of creating habitable areas of the planet with lower, Earth-like CO 2 levels. Given the effects of planetary obliquity (Ward, 1974; Williams, 1975) and eccentricity (Berger et al., 1993, 2006) on seasonality and annual global instellation, the best possible scenario for habitability in the low-CO 2 case may be one in which Kepler-62f has a high obliquity and a high eccentricity. Additionally, an orbital configuration in which the hotter, summer months in a given hemisphere coincide with the pericenter of the planet's orbit could amplify the effects of high obliquity and eccentricity. To test this prediction, we ran simulations with CCSM4 (and LMD Generic GCM), assuming an aqua planet, with a range of different values for these parameters.

We ran 30-year GCM simulations using CCSM4 with the input maximum initial eccentricity for stable HNBody integrations for Kepler-62f (see Section 3). We also ran additional simulations with lower eccentricity values to explore the effect of eccentricity on instellation and planetary surface temperature. We used a synthetic stellar spectrum from the Pickles Stellar Atlas (Pickles, 1998), with flux in the range of 1150–25000 Å, and an effective photospheric temperature of 4887 K, which is close to the estimated effective temperature of Kepler-62 (4925 K, Borucki et al., 2013). The percentage of the total flux from the star was specified in each of the 12 incident wavelength bands in CAM4 (see Table 3), as done in Shields et al. (2013, 2014).

Table 3. CAM4 Spectral Wavelength Bands Specifying Shortwave (Stellar) Incoming Flux into the Atmosphere, and the Percentage of Flux within Each Waveband for a Synthetic Spectrum of a K Dwarf Star with a Similar Photospheric Temperature to Kepler-62, from the Pickles Stellar Atlas (Pickles,1998) Band λ min λ max K2V star % flux 1 0.200 0.245 0.128 2 0.245 0.265 0.075 3 0.265 0.275 0.054 4 0.275 0.285 0.056 5 0.285 0.295 0.070 6 0.295 0.305 0.091 7 0.305 0.350 1.076 8 0.350 0.640 27.33 9 0.640 0.700 6.831 10 0.700 5.000 64.35 11 2.630 2.860 0.000 12 4.160 4.550 0.000

We also ran 40-year CCSM4 simulations with an Earth-like obliquity of 23°, and with an obliquity of 60°. To explore the influence of the location of summer solstice relative to pericenter, we varied the angle of the vernal equinox relative to the longitude of pericenter (VEP), which governs the difference in instellation between southern hemisphere summer and northern hemisphere summer (Fig. 13). Because CCSM4 is parameterized for Earth-like conditions, we kept the radius, mass, and surface gravity of the planet equal to those of Earth in these simulations.

2.4.2. LMD generic GCM

We used the Laboratoire de Météorologie Dynamique (LMD) Generic GCM (Hourdin et al., 2006), developed to simulate a wide range of planetary atmospheres and climates. It has been used in studies of the early climates of Solar System planets (Charnay et al., 2013; Forget et al., 2013; Wordsworth et al., 2013) and in previous studies of the climates of extrasolar planets (Wordsworth et al., 2011; Leconte et al., 2013; Wordsworth, 2015).

Like CCSM4, LMD Generic GCM solves the primitive equations of fluid dynamics using a finite difference method and a 3-D dynamical core. The radiative transfer scheme is based on a correlated-k method, with absorption coefficients calculated from high-resolution spectra generated using the HITRAN 2008 database (Rothman et al., 2009). A smaller database of correlated-k coefficients was then generated from this spectra using 12 × 9 × 8 temperature (T = 100, 150, …, 600, 650 K, in 50 K steps), log-pressure (p = 0.1, 1, 10, …, 107 Pa), and water vapor mixing ratio (q H 2 O = 10−7, 10−6, …, 1) grids. These correlated-k coefficients ensure expedient radiative transfer calculations in the GCM. The spectral intervals included 38 shortwave (incident stellar radiation) bands and 36 longwave (outgoing radiation) bands, with a 16-point cumulative distribution function for integration of absorption data within each band. The radiative transfer equation is then solved in each atmospheric layer using a two-stream approximation (Toon et al., 1989). Parameterizations for convective adjustment, and CO 2 collision-induced absorption, are included based on work by Wordsworth et al. (2010). Simulations were run with a two-layer ocean, including a 50 m deep top layer and an underlying 150 m deep second layer. Both layers are assumed to be well mixed, with horizontal diffusion used to approximate ocean heat transport by large-scale eddies. Adiabatic adjustment is also included, whereby the ocean lapse rate is adjusted to maintain a warmer top ocean layer at all times.

We ran our LMD Generic GCM simulations for 60–200 years (depending on the amount of CO 2 and required model equilibration time) with a horizontal spatial resolution of 64 × 48 (corresponding to 5.625° longitude × 3.75° latitude), with 20 vertical levels. A blackbody spectrum with the stellar effective temperature of Kepler-62 (4925 K, Borucki et al., 2013) was used as the host star spectrum. The albedo of snow was set to 0.55 for the high-CO 2 simulations and 0.43 for the simulations with Earth-like CO 2 (to match the two-band albedo parameterization in CCSM4 for accurate comparison). The albedo of sea ice was allowed to vary between a minimum of 0.20 and a maximum of 0.65, depending on ice thickness. We assumed an aqua planet with different amounts of atmospheric CO 2 , and water vapor that varied throughout each simulation as with CCSM4. We used the radius of Kepler-62f (1.41 R ⊕ , Borucki et al., 2013), a surface gravity of 14.3 m/s2, based on a ∼3 M ⊕ planet (see Section 2.2), a 267-day year, Earth-like (24 h) and synchronous (267-day) rotation rates, and the minimum and maximum stable initial eccentricities possible for Kepler-62f (see Section 3). Additional input parameters for simulations using CCSM4 and LMD Generic GCM for our work on Kepler-62f are given in Tables 4 and 5. The results of these simulations and the implications for the habitability of Kepler-62f are presented and discussed in the following section.

Table 4. Additional Input Parameters Used in CCSM4 Climate Simulations of Kepler-62f Parameter Run1 Run2 Run3 Run4 Run5 Run6 Run7 Orbital eccentricity 0.0 0.1 0.2 0.32 0.32 0.32 0.0 Obliquity (°) 23 23 23 23 60 23 0 Rotation rate 24 h 24 h 24 h 24 h 24 h 24 h 365 days VEP (°) 90 90 90 90 90 0 90

Table 5. Additional Input Parameters Used in High-CO 2 LMD Generic GCM Climate Simulations of Kepler-62f Parameter Run1 Run2 Run3 Run4 Run5 Run6 Run7 Run8 Run9 Run10 Run11 Run12–27 Ecc 0.32 0.00 0.32 0.32 0.32 0.32 0.32 0.00 0.00 0.00 0.00 0–0.32 Obl (°) 23.5 23.5 23.5 23.5 23.5 23.5 23.5 0.00 0.00 0.00 0.00 0–90 CO 2 (bar) 1 3 3 5 8 10 12 1 1 3 3 3 Rot 24 h 24 h 24 h 24 h 24 h 24 h 24 h 24 h 267 days 24 h 267 days 24 h VEP 0 0 0 0 0 0 0 0 0 0 0 0–90

3. Results

3.1. N-body simulations

Figure 2 shows the fraction of stable orbital integrations performed using HNBody, assuming different initial eccentricities for Kepler-62f and zero eccentricity for Kepler-62b–e (see e.g., Barnes and Quinn, 2001). The maximum initial eccentricity for which stable integrations were possible for greater than 90% of the simulated longitudes of pericenter was e = 0.32, assuming the Kopparapu et al. (2014) mass-radius relation. A higher upper eccentricity limit is 0.36 assuming a smaller mass (equal to that of Earth) for Kepler-62f and a larger mass for Kepler-62e. Simulations using a shorter time step of 1/100th of the orbital period of the innermost planet reduced the percentage of longitudes of pericenter with stable integrations at e = 0.32 by ∼30%, though we still found over half of our simulated configurations at this eccentricity to be stable. Stable integrations were possible for 23% of the simulated longitudes of pericenter at e = 0.33. However, we did not consider this a large enough fraction of stable configurations, and consider e = 0.32 to be the conservative maximum. The maximum stable initial eccentricity decreased by ∼3%, to e = 0.31, when the Weiss and Marcy (2014) mass-radius and density relations were used, and decreased by 28%, to e = 0.23, using the maximum mass limits for all planets (Borucki et al., 2013). The evolution of the eccentricities of all five planets for a stable integration at e = 0.32 for Kepler-62f is shown in Fig. 3. All planets exhibited eccentricity oscillations—a characteristic of multiple-planet systems—with a regular period that remained constant over the entire 106-year integration. Integrations assuming higher initial eccentricities for Kepler-62f yielded eccentricity evolution that exhibited irregular oscillatory motion indicative of the occurrence of significant orbital shifts. Consequently, we considered 0.00 ≤ e ≤ 0.32 to be the maximum allowable range of initial eccentricities for Kepler-62f.

FIG. 2. Fraction of stable configurations after a 106-year HNBody integration for initial eccentricities between 0.0 and 0.9 for Kepler-62f, assuming the Kopparapu et al. (2014) mass-radius relationship. The eccentricities of all other planets in the Kepler-62 system were set to zero.

FIG. 3. Evolution of the eccentricities of Kepler-62b–f as a function of time calculated with HNBody. Initial eccentricities and longitudes of pericenter for Kepler-62b–e were set to zero. The initial eccentricity of Kepler-62f was set to 0.32 with a longitude of pericenter equal to π.

3.2. Tidal evolution of Kepler-62f

We began our simulations with e = 0.00 or 0.32 and spin periods and obliquities as given in Section 2.3. In general, tidal evolution is faster for larger eccentricity, larger spin period, and smaller obliquity. Thus the e = 0.32, 10-day spin period and 5° obliquity cases should evolve most rapidly toward an equilibrium spin state.

In Fig. 4 we show the evolution of spin period and obliquity for the two tidal models. The difference between the two models naturally produces different torques on the planetary rotation and thus different timescales to reach equilibrium. On the left, the CTL model shows results very similar to those in Bolmont et al. (2014, 2015), and it is unlikely that the system has reached a tidally locked state, though tidal de-spinning is significant over 10 Gyr. The right panels show that tidal locking is possible, which is easily shown by the spin period evolution flattening out at either 267 or 178 days. The latter corresponds to a 2:3 spin-orbit resonance, which is closer to the equilibrium spin period for an eccentricity of 0.32. Only the extremely fast-rotating and high-obliquity cases do not tidally lock within 5 Gyr in the CPL model.

FIG. 4. Rotational evolution of Kepler-62f due to tidal processes. Left: Evolution of the spin period (top) and obliquity (bottom) for the CTL model. The gray curves assume a circular orbit, and the black curves assume an eccentricity of 0.32. Solid lines represent planets that begin with an 8 h rotation period and an obliquity of 80°. Dashed lines assume the planet began with modern Earth's rotational state. Dotted lines assume an initial spin period of 10 days and an obliquity of 5°. Right: Same as left panels, but for the CPL model.

As the estimated age of Kepler-62 (7 Gyr, Borucki et al., 2013) is uncertain, as well as the initial eccentricity and rotation state, we conclude that the rotational period of Kepler-62f can lie anywhere from less than 1 day all the way to 267 days. Similarly, the obliquity can have a range of values from 0 to at least 90°. For our assumptions, the rotation period and obliquity should lie between the solid gray line and the dotted black line. Therefore, habitability and climate studies of this planet should account for this range of rotational states, as well as both small and large obliquities. We present the results of such studies in the following section.

3.3. Climate simulations

With the present atmospheric level (PAL) of CO 2 on Earth (400 ppmv), all CCSM4 simulations with eccentricity e = 0.00–0.32 resulted in completely ice-covered conditions, with global mean surface temperatures below 190 K. While our CCSM4 simulations were run with a 365-day orbital period rather than the actual orbital period of Kepler-62f (267 days), a comparison of LMD Generic GCM sensitivity tests run with 267- and 365-day orbital periods showed equivalent results. Figure 5 shows the annual mean instellation as a function of latitude for these different eccentricities, assuming an obliquity of 23°. The average instellation over an annual cycle increases with eccentricity (Berger et al., 1993, 2006; Williams and Pollard, 2002), in accordance with the following relation:

where S is the average instellation at the mean star-planet distance; S a is the instellation at a given distance a from the star during a planet's orbit; and e is the eccentricity of the planet's orbit (Berger et al., 2006).

FIG. 5. Instellation as a function of latitude for Kepler-62f (using CCSM4), assuming e = 0.0 (black), e = 0.1 (blue), e = 0.2 (green), and e = 0.32 (red). The plots show 12-month averages. The obliquity of the planet was set to 23°. The angle of the vernal equinox relative to pericenter was set to 90°, similar to that of Earth (102.7°). The larger the eccentricity, the larger the annually averaged instellation received at a given latitude.

If CO 2 is efficiently outgassed, and the silicate weathering rate is weak, higher levels of CO 2 may be expected to accumulate in the planet's atmosphere prior to reaching the maximum greenhouse limit (Kopparapu et al., 2013a, 2013b). Since PAL CO 2 was clearly insufficient to yield habitable surface temperatures for Kepler-62f at its value of incoming stellar flux, we ran simulations using LMD Generic GCM with CO 2 concentrations of up to 12 bar, at both limits of the stable eccentricity range for Kepler-62f. We let the amount of water vapor vary (through transport, evaporation, and circulation) during the course of each simulation. Global mean, minimum, and maximum surface temperatures for each simulation are plotted in Fig. 6. With 5 bar of CO 2 in the atmosphere, Kepler-62f exhibits a global mean surface temperature similar to modern Earth, at ∼282 and ∼290 K for the lower and upper stable eccentricity limits, respectively. Surface temperatures increase with CO 2 concentration. However, while the surface temperature on the planet is 46–50 K warmer (depending on eccentricity) with 5 bar of atmospheric CO 2 compared to the 3 bar CO 2 case, it is only an additional 20–25 K warmer (∼308 K) with 8 bar of CO 2 in the atmosphere. Further increases in atmospheric CO 2 yield successively smaller increases in surface temperature. We also found that the planetary albedo reached a minimum at 8 bar and then began to slowly increase (Fig. 7). This indicates that we are likely approaching the point where the effects of Rayleigh scattering begin to dominate over the greenhouse effect at this level of CO 2 . Additionally, the difference between minimum and maximum surface temperatures decreases with CO 2 concentration, highlighting the role of a denser atmosphere in evening out temperature contrasts (Pierrehumbert, 2011; Wordsworth et al., 2011).

FIG. 6. Mean (plus symbols), minimum (triangles), and maximum (diamonds) surface temperature for Kepler-62f after 100- to 200-year LMD Generic GCM simulations, assuming e = 0.00 (blue) and e = 0.32 (red), an Earth-like rotation rate, and different levels of atmospheric CO 2 . The mean values take into account the location of measured surface temperature values relative to the total surface area of the planet. An obliquity of 23.5° and VEP = 0° are assumed.

FIG. 7. Planetary albedo of Kepler-62f as a function of simulated atmospheric CO 2 concentration after 100- to 200-year LMD Generic GCM simulations. The 3 bar CO 2 simulation, which was ice-covered, with a planetary albedo of 0.491, is not shown, to enlarge the turning point in planetary albedo at 8 bar CO 2 . An obliquity of 23.5°, an eccentricity of zero, and VEP = 0° are assumed.

While the majority of our simulations run with 3 bar of CO 2 , including that run at Earth's present obliquity, yielded globally ice-covered conditions (Fig. 8), we did find that at an extremely high planetary obliquity (90°), at the upper eccentricity limit (0.32), surface temperatures exceeded the freezing point of water over close to 20% of the planet, though the global mean surface temperature was well below freezing, at 246 K. As shown in Fig. 8, the regions with zero ice cover occur in the polar regions, which receive more instellation than the equator at high obliquities. We also found that, while high-eccentricity simulations with an obliquity of 60° ultimately yielded ice-covered conditions, maximum surface temperatures reached 273 K during the planet's orbit at small orbital distances from its star, which could result in melting of ice formed when the planet is at more distant points in its orbit.

FIG. 8. Sea ice cover fraction as a function of latitude for Kepler-62f after 160-year LMD Generic GCM simulations, with 3 bar of CO 2 in the atmosphere, the maximum stable initial eccentricity (e = 0.32), Earth-like (blue) and 90° (green) obliquities, and an Earth-like rotation rate. A VEP of 0° is assumed.

Figure 9 shows the surface temperature as a function of latitude and longitude for an LMD Generic GCM simulation with 5 bar of CO 2 , variable water vapor, and e = 0.00. With a global mean surface temperature of 282 K, only ∼17% of the planet is ice-covered. Simulations run at e = 0.32 resulted in a higher global mean surface temperature by ∼8 K, and ∼10% less ice cover on the planet.

FIG. 9. Surface temperature as a function of latitude for Kepler-62f after a 120-year LMD Generic GCM simulation, assuming e = 0, an obliquity of 23.5°, 5 bar of CO 2 , and an Earth-like rotation rate. A VEP of 0° is assumed.

Given that Kepler-62f could have a wide range of rotation periods, including a synchronous rotation rate, using CCSM4 and LMD Generic GCM we ran simulations of the planet assuming an Earth-like rotation rate, and a synchronous rotation period, for 400 ppmv (Earth-like), 1 bar, and 3 bar CO 2 atmospheres, with both the eccentricity and obliquity set to zero. Our CCSM4 simulations with Earth-like CO 2 levels were completely ice covered. Figure 10 shows the surface temperature as a function of latitude and longitude for the LMD Generic GCM 1 bar and 3 bar CO 2 cases. The synchronous case with 1 bar of CO 2 has a global mean surface temperature of ∼206 K, with 99.7% of the planet covered in ice. There is a small circular region of open water at the substellar point on the planet (Fig. 10). However, as we have not included sea ice transport in our GCM, and glacial flow of thick sea ice could cause a planet to become fully glaciated (Abbot et al., 2011), this region of open water is likely unstable. Furthermore, temperatures on the nightside of the planet reach below the limit for CO 2 to condense at 1 bar surface pressure, indicating the likelihood for atmospheric freeze-out at this atmospheric concentration. The nonsynchronous case with the same CO 2 level has a global mean surface temperature that is ∼11 degrees warmer (∼217 K), although it is completely covered in ice. The 3 bar CO 2 cases, while ∼17–18 K warmer, exhibit a similar behavior, with the nonsynchronous case yielding a global mean surface temperature (∼235 K) that is ∼12 K warmer than the synchronous case (∼223 K). Surface temperatures on the nightside of the planet in the 3 bar CO 2 case are right at the boundary for CO 2 condensation at this surface pressure and CO 2 concentration.

FIG. 10. Surface temperature for a synchronous (left) and Earth-like (24 h) rotation rate for Kepler-62f, after a 65- to 100-year LMD Generic GCM simulation with 1 bar (top) and 3 bar (bottom) of CO 2 in the atmosphere. We assumed e = 0, VEP = 0°, and an obliquity of 0° for all four simulations.

High amounts of CO 2 would be a relatively straightforward means of generating habitable surface temperatures on a planet receiving low instellation from its star. However, the efficiency of the carbonate-silicate cycle—which has been shown to be sensitive to a variety of factors, including the mantle degassing rate of a planet (Driscoll and Bercovici, 2013)—is unknown for Kepler-62f. If the planet has a low amount of CO 2 in its atmosphere, and lacks an active carbon cycle to adjust the silicate weathering rate with temperature (Walker et al., 1981), the fraction of habitable surface area may decrease significantly compared to our simulated cases with 3 bar and higher CO 2 . We therefore explored different orbital configurations that could improve conditions for habitability on a planet that does not have an effective means of increasing its concentration of greenhouse gases.

Since we do not know the location of pericenter for the orbit of Kepler-62f, we explored the effect on instellation of changes in the VEP of a planet's orbit, as this can affect the hemispherical annually averaged instellation on a planet (Berger, 1978; Berger et al., 1993). Figure 11 shows the results of CCSM4 simulations assuming a VEP of 0° and 90°, with the obliquity and eccentricity held fixed at 23° and 0.32, respectively. At VEP = 0°, the point where both hemispheres receive equal amounts of instellation coincides with the planet's closest approach to its star. The difference in monthly instellation is relatively small in the southern hemisphere, as it receives equal instellation to that received by the northern hemisphere at pericenter. At a VEP of 90°, the southern hemisphere receives a significantly higher percentage of instellation compared to the northern hemisphere during its summer months when the planet is at or near pericenter. Because of the planet's obliquity, this is when the southern hemisphere is angled toward the star.

FIG. 11. Annual mean instellation as a function of latitude for Kepler-62f as a function of the month of the year after 40-year CCSM4 simulations, assuming a 12-month annual cycle and a VEP of 0° (left) and 90° (right). The obliquity and eccentricity of the planet were set to 23° and 0.32, respectively.

The obliquity of Kepler-62f is observationally unconstrained, so we also explored how different obliquities might affect the planet's climate in a low-CO 2 scenario. As shown in Fig. 12, CCSM4 at an obliquity of 60° reveals more instellation received in the high-latitude regions of the planet than in the tropics. The global mean surface temperature is still significantly below freezing at both of the simulated obliquities (23° and 60°), given the Earth-like CO 2 levels and low stellar flux. However, surface temperatures do get above freezing in the southern hemisphere during its summer months in the high-obliquity case, with VEP = 90°, in both our CCSM4 and LMD Generic GCM simulations. This means that summer in the southern hemisphere occurs near the planet's closest approach to its star, as shown in the schematic diagram in Fig. 13. This results in higher annual mean surface temperatures in this hemisphere during its summer months, compared to the northern hemisphere. The less extreme cold temperatures in the LMD Generic GCM simulations are due to a 10 m maximum thickness limit of sea ice in LMD Generic GCM, while the CCSM4 sea ice thickness near the poles is ∼30 m. This causes a greater temperature difference between the models at the poles during winter months. The reduced ice thickness in LMD Generic GCM results in a larger conductive heat flux through the ice, though this has a lower impact for the warmer climates (e.g., with 3–5 bar and higher CO 2 ), where there is less ice. We calculated the temperature difference between the models by equating the residual surface energy flux with the conductive heat flux through the ice assuming 10 and 30 m ice thicknesses, and found the difference to equal that between the coldest temperatures in LMD Generic GCM and CCSM4 (∼60 K). Both models indicate that an orbital configuration that places the summer solstice near pericenter may amplify the effects of high obliquity and eccentricity, and this may cause surface melting to occur during an annual cycle in a low-CO 2 scenario. This effect is also seen clearly in a comparison between LMD Generic GCM simulations with 3 bar of atmospheric CO 2 and VEP of 0° and 90° (Fig. 14), where the latter simulation shows warmer surface temperatures in the southern hemisphere of the planet. We noticed a decrease in planetary ice cover of ∼0.6% per 90-degree increase in VEP in our 3 bar CO 2 simulations run at 90° obliquity and e = 0.32.

FIG. 12. Top: Annual mean instellation as a function of latitude for Kepler-62f after 40-year CCSM4 simulations, assuming the PAL of CO 2 on Earth (400 ppmv), an obliquity of 23° (blue) and 60° (green). Middle: Surface temperature as a function of the month of the year, assuming a 12-month annual cycle, for an obliquity of 60° (left) and 23° (right), after 40-year CCSM4 simulations. Bottom: Surface temperature as a function of month of year for an obliquity of 60° (left) and 23° (right), after 60-year LMD Generic GCM simulations. In both models VEP was set to 90°, similar to Earth (102.7°). The eccentricity was set to 0.32. The less extreme cold temperatures in the LMD Generic GCM simulations are due to a 10 m maximum thickness limit of sea ice in LMD Generic GCM, while the CCSM4 sea ice thickness near the poles is ∼30 m.

FIG. 13. Schematic diagram of assumed orbital configuration for CCSM4 (and select LMD Generic GCM) simulations of Kepler-62f. The angle of the vernal equinox with respect to pericenter was set to 90°, similar to Earth (102.7°).

FIG. 14. Surface temperature as a function of latitude for Kepler-62f, after 160-year LMD Generic GCM simulations with 3 bar of atmospheric CO 2 and VEP 0° (dashed line) and 90° (solid line). An obliquity of 90° and an eccentricity of 0.32 are assumed.

4. Discussion

We explored the plausible range of orbital, rotational, and atmospheric states of Kepler-62f and found that some permit habitability, but some do not. A global mean surface temperature similar to modern-day Earth is possible throughout the planet's orbit with 5 bar of CO 2 in its atmosphere, for both lower and upper limits of the range of stable initial eccentricities (e = 0.00 and e = 0.32), assuming present Earth obliquity. Our simulations with e = 0 yielded lower global mean surface temperatures, consistent with average instellation decreasing with decreasing eccentricity. These results indicate that if Kepler-62f has an active carbon cycle, where CO 2 is allowed to build up in the atmosphere as silicate weathering decreases at lower surface temperatures (Walker et al., 1981), the best possible scenario for habitable surface conditions is one that combines moderate to high eccentricity with high atmospheric CO 2 . The maximum CO 2 greenhouse limit for stars with the effective temperature of Kepler-62 occurs at a stellar flux3 that is well below that received by Kepler-62f (41% of the modern solar constant, Borucki et al., 2013). Given that 3–5 bar of CO 2 are significantly below this limit (∼7–8 bar, Kopparapu et al., 2013a, 2013b), these levels of CO 2 are plausible for this planet.

Earlier work exploring the possible climates of Kepler-62f proposed habitable conditions with 1.6 bar of CO 2 or more in the planet's atmosphere (Kaltenegger et al., 2013). However, this work was not done with a 3-D GCM and therefore lacked a full treatment of clouds and atmosphere-ocean interactions. Using a 3-D GCM, we find that the only scenario that permits Kepler-62f to exhibit clement temperatures for surface liquid water with 3 bar of CO 2 in the planet's atmosphere requires stringent orbital configuration requirements (high eccentricity and obliquity). With an obliquity of 90°, open water was present over just ∼20% of the planet, at polar latitudes. This indicates that it is likely that more than 1.6 bar of CO 2 is required for surface habitability on this planet.

We assumed a five-planet system for Kepler-62 in our n-body simulations. While there is currently no observational evidence for additional planets, the existence of other planets in this system is a clear possibility, and their presence would certainly affect the eccentricity limit for dynamical stability that we have calculated here. Future discovery of additional planets in this system with transit timing variations or RV may provide additional dynamical constraints on the eccentricity of Kepler-62f.

We combined climate simulations using CCSM4 and LMD Generic GCM to provide a comprehensive exploration of the possible climates for Kepler-62f given its n-body model constraints. Previous work using LMD Generic GCM at a higher resolution (e.g., 128 × 96) yielded similar results to those at the resolution we employed here (64 × 48), and as we confirmed that the global mean surface temperature had not changed by more than 1 K in the last 20 years of simulation, we are confident that running our simulations for longer timescales would not change our results significantly. We verified that simulations assuming Earth-like CO 2 levels in both GCMs exhibited similar (ice-covered) conditions, confirming the robustness of these simulations to various assumptions.

Our simulations with higher levels of CO 2 resulted in increasingly higher surface temperatures on the planet. As CO 2 was increased, the changes in global mean surface temperatures became progressively smaller. This logarithmic relationship between CO 2 concentration and radiative forcing is long established (Wigley, 1987). As regions of the spectrum become opaque, additional CO 2 molecules become far less effective at increasing temperatures (Shine et al., 1990). Thus, the greenhouse effect starts to become less efficient as a warming mechanism. Additionally, CO 2 is 2.5× more effective as a Rayleigh scatterer than Earth's air (Kasting, 1991; Forget and Pierrehumbert, 1997), and this behavior likely contributes to the loss of warming at higher CO 2 concentrations (Kasting, 1991; Kasting et al., 1993; Selsis et al., 2007).

We have not included the effect of CO 2 condensation in our simulations. As discussed in Section 2.4, at levels of 1–2 bar, CO 2 condensation is likely to occur in the upper atmosphere (Pierrehumbert, 2005). Depending on the particle size of CO 2 ice grains, this could result in cooling of the planet due to the albedo effect of CO 2 ice clouds (Kasting, 1991) or warming by scattering outgoing thermal radiation back toward the surface of the planet (Forget and Pierrehumbert, 1997).

Kopparapu et al. (2013a, 2013b) found that the maximum CO 2 greenhouse limit is ∼7–8 bar for a star with a similar effective temperature to that of Kepler-62 (Kopparapu et al., 2013a, 2013b). We did find that the planetary albedo, which had decreased with increasing CO 2 concentration, started to increase, albeit slowly, at a CO 2 level of 8 bar. However, our simulations with 8–12 bar of CO 2 resulted in global mean surface temperatures that were still higher than those with lower CO 2 , so we conclude that we had not yet reached the maximum CO 2 limit in our simulations and that it may be higher than originally proposed. Kopparapu et al. (2013a, 2013b) used a 1-D radiative-convective model in their work and did not include the effect of water clouds or CO 2 clouds in their calculations. While water clouds could increase the planetary albedo, thereby cooling the planet further, they may also contribute to the greenhouse effect, as both H 2 O and CO 2 have strong absorption coefficients in the near-IR, which increase the amount of radiation absorbed by planets with lower-mass host stars that emit strongly in the near-IR (Kasting et al., 1993; Selsis et al., 2007; Joshi and Haberle, 2012; Kopparapu et al., 2013a, 2013b; Shields et al., 2013, 2014). Our results with a 3-D GCM do include water clouds, though not CO 2 clouds. A comprehensive study of the effect of CO 2 condensation as a function of particle size would be an important step toward identifying its ultimate effect on planetary climate. Regardless, dense CO 2 atmospheres can be expected to cause more even distributions of heat across a planet, and reduce the contrast between maximum and minimum surface temperatures. This is particularly important on synchronously rotating planets (Joshi et al., 1997; Edson et al., 2011), where the difference in instellation on the dayside and nightside of the planet is large. We also see this result in our simulations, which show smaller maximum/minimum surface temperature contrasts with larger amounts of CO 2 .

All our simulations assumed an aqua planet configuration, with no land. Previous work exploring the habitability of planets composed almost entirely of land and orbiting G dwarf stars suggests that due to their lower thermal inertia and drier atmospheres, land planets are less susceptible to snowball episodes than aqua planets, requiring 13% less instellation to freeze over entirely (Abe et al., 2011). The presence of land could certainly affect the silicate weathering rate and atmospheric concentration of CO 2 on a planet. However, Abbot et al. (2012) found that climate weathering feedback does not have a strong dependence on land fraction, as long as the land fraction is at least 0.01. Edson et al. (2012) found that the amount of CO 2 that accumulates in the atmosphere of a synchronously rotating planet could be much greater if the substellar point is located over an ocean-covered area of the planet, where continental weathering is minimal, although atmospheric CO 2 concentration could still be limited by seafloor weathering processes (Edson et al., 2012). Including land in future simulations of a synchronously rotating Kepler-62f, using a GCM with a carbonate-silicate cycle included (rather than assigning a prescribed atmospheric CO 2 concentration as we have done here) would be a valuable step toward assessing the role of surface type in regulating the atmospheric CO 2 inventory on synchronously rotating planets.

We have concentrated on CO 2 as the primary greenhouse gas in our study of the influence of atmospheric composition on the habitability of Kepler-62f. Previous work has highlighted the role of molecular hydrogen as an incondensable greenhouse gas that, due to collision-induced absorption, could allow planets to maintain clement temperatures for surface liquid water far beyond the traditional outer edge of a star's habitable zone (Pierrehumbert and Gaidos, 2011). Given that the orbital distance of Kepler-62f is interior to the limit for which an entire primordial H 2 envelope could be lost due to extreme UV-driven atmospheric escape (assuming a 3–4 M ⊕ planet), even at apocenter (at an eccentricity of 0.32), an H 2 greenhouse may not be an effective mechanism for warming this planet. Exploring this mechanism in detail in the context of the Kepler-62 system would be an interesting topic for future study.

4.1. Additional orbital influences on the climate of Kepler-62f

High amounts of CO 2 in the atmosphere of Kepler-62f would provide for consistently habitable surface conditions throughout the planet's orbit at its level of instellation, thereby offering the best chances for sustained surface liquid water and life. However, given the uncertainty of an active carbonate-silicate cycle operating on this planet, we explored orbital parameters that may periodically improve surface conditions during the course of an orbit in the event of low-CO 2 atmospheric conditions similar to Earth, though these conditions would be short-lived. Planetary habitability throughout the course of an eccentric orbit has been shown to be most strongly affected by the time-averaged global instellation—provided there is an ocean to contribute to the planet's heat capacity—which is greater at higher eccentricities (Williams and Pollard, 2002). Therefore, although high-eccentricity planets may spend significant fractions of an orbit outside of their stars' habitable zones, high eccentricity may help these planets maintain habitable surface conditions over an annual cycle (Kopparapu et al., 2013a, 2013b), though habitability could be affected by tidal forces and resultant heating at close distances from the star during a planet's eccentric orbit (Driscoll and Barnes, 2015). Planets on eccentric orbits could undergo freeze/thaw cycles depending on the orbital distance of apocenter and pericenter, respectively, and this could also have important implications for planetary habitability. Limbach and Turner (2015) found that catalogued RV systems of higher multiplicity exhibit lower eccentricities. However, as these are statistical results, they do not mean that all such systems have low eccentricities. Furthermore, they used data from RV systems of more massive planets than those in the Kepler-62 system. There are currently no observational constraints on the eccentricities of low-mass, long-period extrasolar planetary systems. Our calculated eccentricity of 0.32 is entirely plausible for Kepler-62f and consistent with observations.

Given that the obliquity, rotation rate, and the VEP for Kepler-62f are unknown, we explored how these factors might also influence surface habitability on this planet. Our results indicate that high obliquity, which has been shown to increase seasonality and stability against snowball episodes (Williams and Kasting, 1997; Williams and Pollard, 2003; Spiegel et al., 2009), results in even higher seasonality at high eccentricity, due to the larger difference between the orbital distance at pericenter and apocenter (Williams and Pollard, 2002). This effect is more pronounced in the southern hemisphere during its summer months, when the angle of the vernal equinox relative to pericenter is 90°, so that the planet's high obliquity significantly increases the instellation received by this hemisphere. If we had simulated VEP = 270°, we would expect the northern hemisphere to exhibit a similar effect, as this is a symmetric problem. Though we did not include ocean heat transport in our high-obliquity CCSM4 simulations that yielded global ice cover, previous work has found that snowball collapse is possible at high obliquity regardless of the presence of a dynamical ocean with ocean heat transport (Ferreira et al., 2014).

We identified several orbital configurations that, though rare, may cause temporarily habitable surface conditions during the course of an orbit given Earth-like levels of CO 2 . Even on an ice-covered high-obliquity planet, our simulations using both models showed that surface temperatures reached above the freezing point of water during southern hemisphere summer months, which could cause surface melting. Our LMD Generic GCM simulations with 3 bar of CO 2 also resulted in higher surface temperatures in the summer hemisphere with a higher VEP of 90° compared to 0°. The VEP can therefore amplify the warming effects of high obliquity and eccentricity, and may keep ice from forming at the poles, or reducing an ice sheet formed during the planet's orbit. If the planet experiences large oscillations in obliquity, polar ice could be prevented from forming on both hemispheres over an annual cycle, allowing habitable surface conditions to be maintained on planets with large eccentricities (Armstrong et al., 2014).

The weaker Coriolis force on synchronously rotating planets permits rapid heat transport by advection compared to radiative heat transfer for surface pressures significantly greater than ∼0.2 bar (Showman et al., 2013). Combined with sufficient greenhouse gas concentration levels (Joshi et al., 1997), this can prevent atmospheric freeze-out on the nightside of the planet, through the transport of high amounts of heat from the dayside to the nightside (Pierrehumbert, 2011). Our results have shown that a synchronously rotating Kepler-62f exhibits a lower global mean surface temperature than one that is nonsynchronous. This may underscore a lack of sufficient heat distribution between the dayside and nightside of the planet at the levels of CO 2 that we have simulated at this planet's instellation. Horizontal heat transport is less effective for smaller amounts of CO 2 (Wordsworth et al., 2011), which we confirmed in the 1 bar CO 2 case, where the ratio between the outgoing thermal fluxes on the nightside and dayside (a proxy for the redistribution efficiency) is lower. This caused larger differences between maximum and minimum surface temperatures on the dayside and nightside of the 1 bar CO 2 planet compared to the 3 bar CO 2 case.

Additionally, though we did not measure significant differences in planetary albedo between the synchronous and nonsynchronous climate simulations at a given CO 2 concentration, Kepler-62f lies near the outer edge of its star's habitable zone, receiving 41% instellation. On more temperate synchronously rotating planets than we have simulated here, the increased cloud cover that will likely result on the dayside of a synchronously rotating planet with an ocean could increase the overall planetary albedo and reduce surface temperatures (Yang et al., 2013, 2014). The strength of this effect on distant worlds would depend on the greenhouse gas concentration and the behavior of the hydrological cycle. The results presented here imply that synchronously rotating planets may require more CO 2 in their atmospheres than their nonsynchronous counterparts to generate equivalent global mean surface temperatures farther out in their star's habitable zones, and this could affect planetary habitability.

The n-body model used in this work does not include the effect of tides. With an orbital period of 267 days, tides are likely to be weak on Kepler-62f, but they could affect its rotation period (see Section 3.2), and would certainly be an important consideration for potentially habitable planets orbiting even closer to their stars. Tidal effects can lead to changes in orbital parameters, and may induce capture into resonances in spin-orbit period, depending on the planet's eccentricity and its equatorial ellipticity (the equilibrium shape attained as a result of the gravitational interaction between the planet and the host star, Rodríguez et al., 2012). Our work here explored the extreme case of spin-orbit resonance—synchronous rotation. Nonsynchronous resonance configurations, such as the 3:2 spin-orbit resonance observed on the planet Mercury (Goldreich and Peale, 1966; Correia and Laskar, 2004) are also possible and may be sustained on a planet in a noncircular orbit over long timescales (Rodríguez et al., 2012). Such configurations become more likely at larger orbital eccentricities (Malhotra, 1998; Correia and Laskar, 2004) and are worth exploring.

The habitability requirements we have determined here did not account for any uncertainties in measured parameters such as stellar luminosity and semimajor axis. Including such uncertainties might permit a larger habitable surface area for an even wider range of parameters.

Using constraints from n-body models as input to GCM simulations permits the exploration of the impact on climate of the gravitational interactions inherent to a growing population of exoplanets—potentially habitable planets orbiting in multiple-planet systems around low-mass stars—and the resulting prospects for the habitability of these worlds. The techniques presented here can be applied to planets orbiting stars of any spectral type, with a range of possible atmospheric and surface compositions and dynamical architectures. They can be used to help assess the potential habitability of newly discovered planets for which observational measurements are still limited, and can be easily modified to incorporate new observational data that are acquired for these planets in the future. While future missions such as the Transiting Exoplanet Survey Satellite (TESS, Ricker et al., 2009, 2014) and the James Webb Space Telescope (Gardner et al., 2006) will not be capable of characterizing distant planets like Kepler-62f, the procedure we have carried out here serves as excellent preparation for closer planets that could be observed and characterized with these missions. These methods will help identify which among the habitable-zone planets discovered to date could exhibit conditions conducive to the presence of surface liquid water for the widest range of atmospheric and orbital conditions, making them priorities for these and other future characterization missions.

5. Conclusions

We carried out a comprehensive exploration of the orbital evolution of Kepler-62f using an n-body model and found that the range of eccentricities that Kepler-62f could have while maintaining dynamical stability within the system was 0.00 ≤ e ≤ 0.32. The upper limit of 0.32 is consistent with the analytic Hill Stability criterion (cf. Gladman, 1993; Barnes and Greenberg, 2006). A higher upper eccentricity limit is 0.36 assuming a smaller mass (equal to that of the Earth) for Kepler-62f and a larger mass for Kepler-62e. The constraints from the n-body model were used as input to 3-D climate simulations to explore the possible climates and habitability of Kepler-62f.

At 41% of the modern solar constant, this planet will likely require an active carbonate-silicate cycle (or some other means by which to produce high greenhouse gas concentrations) to maintain clement conditions for surface liquid water. With 3 bar of CO 2 in its atmosphere and an Earth-like rotation rate, 3-D climate simulations of Kepler-62f yielded open water across ∼20% of the planetary surface at the upper limit of the stable eccentricity range possible for the planet, provided that it has an extreme obliquity (90°). With 5 bar of CO 2 in its atmosphere, a global mean surface temperature similar to modern-day Earth is possible for the full range of stable eccentricities and at the present obliquity of Earth. This higher CO 2 level is therefore optimal, as it is below the maximum CO 2 greenhouse limit, and generates habitable surface conditions for a wide range of orbital configurations throughout the entire orbital period of the planet. If Kepler-62f is synchronously rotating, CO 2 concentrations above 3 bar would be required to distribute sufficient heat to the nightside of the planet to avoid atmospheric freeze-out.

We have also shown that surface temperatures above the freezing point of water during an annual cycle are possible on the planet if it has a low (Earth-like) level of CO 2 , provided that the obliquity is high (60° or greater) compared to an Earth-like obliquity (23°), and the summer solstice at a given hemisphere occurs at or near the planet's closest approach to its star. This is a rare but possible orbital configuration that could cause surface melting of an ice sheet formed during a planet's orbit and amplify the effects of moderate to high eccentricity and obliquity. While less optimal than the high-CO 2 case, this may result in periodically habitable surface conditions in a low-CO 2 scenario for Kepler-62f.

Acknowledgments

This material is based upon work supported by the National Science Foundation under Award No. 1401554, and Grant Nos. DGE-0718124 and DGE-1256082, and by a University of California President's Postdoctoral Fellowship. This work was facilitated through the use of advanced computational, storage, and networking infrastructure provided by the Hyak supercomputer system at the University of Washington. We would also like to acknowledge high-performance computing support from Yellowstone (ark:/85065/d7wd3xhc) provided by NCAR's Computational and Information Systems Laboratory, sponsored by the National Science Foundation. This work was performed as part of the NASA Astrobiology Institute's Virtual Planetary Laboratory Lead Team, supported by the National Aeronautics and Space Administration through the NASA Astrobiology Institute under solicitation NNH12ZDA002C and Cooperative Agreement Number NNA13AA93A. The authors wish to thank Dorian Abbot and an anonymous referee for their comments and suggestions, which improved the paper. B.C. acknowledges support from an appointment to the NASA Postdoctoral Program at NAI Virtual Planetary Laboratory, administered by Oak Ridge Affiliated Universities. R.B. acknowledges support from NSF grant AST-1108882. A.S. thanks Brad Hansen, Jonathan Mitchell, John Johnson, Russell Deitrick, and Tom Quinn for helpful discussions regarding this work.

Author Disclosure Statement

No competing financial interests exist.

Appendix: Determining Initial Planet Locations The location of a planet can be expressed in terms of the planet's true anomaly, the longitude measured from the direction of the pericenter of the planet's orbit. We used the transit times for each of the five Kepler-62 planets and calculated the true anomaly values for all five planets at the same epoch, using one planet's location as a reference point. As stated in Section 2.1.1, we do not know the eccentricity of Kepler-62f (we assumed e = 0 for Kepler-62b–e). The angular velocity and location of a planet in its orbit depend on its eccentricity, following Kepler's second law; therefore, we do not know the pericenter of the orbit for Kepler-62f. We therefore calculated the true anomaly for a range of possible eccentricities for Kepler-62f. Here we outline the equations used in our model to generate the true anomaly values, given assumed values for other key orbital parameters. Time of pericenter passage The time t(f i ) that it takes for a planet to go from some initial reference point (to which it arrives at t o ) to the point of closest approach (where it arrives at t peri ), can be calculated using the following formula given by Sudarsky et al. (2005): Taking the relation that t(f i ) = t o − t peri , we rearranged A1 to yield an expression for t peri : To find each planet's true anomaly, it is necessary to first determine the location of each planet when it transits its star, f i . Using Fig. 1 for reference, with Θ = f + ω = 0 when the planet m 2 passes through the sky plane (the plane perpendicular to the reference direction), at the time the planet passes behind the star, Θ = π/2. By extension, the planet is in front of the star and transiting at Θ = 3π/2. Therefore, for any value of ω, f i = 3π/2 − ω. Since we have assumed ω = 0 for all planets except Kepler-62f, f i = 3π/2 for these planets. Values of f i for Kepler-62f varied depending on the value of ω. Values used in Eq. A2 for the orbital period P, the transit time for each planet t o , e, and f i are given in Table A1. With these values, t peri was calculated for each planet. Table A1. Key Parameters Used as Inputs toEq. A2for Kepler-62b–f Kepler-62b Kepler-62c Kepler-62d Kepler-62e Kepler-62f P (days) 5.714932 12.4417 18.16406 122.3874 267.291 t o (BJD-2454900) 103.9189 67.651 113.8117 83.404 522.710 e 0.0 0.0 0.0 0.0 0.0–0.9 f i 3π/2 3π/2 3π/2 3π/2 -π/2 → 3π/2 From mean anomaly to true anomaly The location of a planet in its orbit around a central star is generally described by one of three angular parameters: the mean anomaly, the eccentric anomaly, or the true anomaly. The mean anomaly is a function of the average angular velocity of the planet and therefore does not provide a precise location. Rather, the mean anomaly is an approximate location that is easy to calculate: where n = 2π/P. We used the first transit time for Kepler-62f as the reference point t o (see Table A1). The eccentric anomaly E is an angle measured from a line through the focus of the planet's elliptical orbit to the radius of a circle that passes through the pericenter of the ellipse. It is related to the mean anomaly and the orbital eccentricity e through Kepler's equation: As this equation cannot be solved analytically, we solved it numerically following the method provided by Danby and Burkardt (1983). With calculated values of E for our assigned eccentricities, the true anomaly was calculated using the following relation: Equation A5 allows us to calculate the relative positions of all the planets at any given time so that our n-body integrations will realistically reproduce the system's orbital evolution.

1 http://kepler.nasa.gov as of May 10, 2016

2 We used three components of CCSM4 in this work: CICE4, CAM4, and a slab ocean. For simplicity, we refer collectively to this suite of model components as CCSM4 throughout this work, as done in Shields et al. (2013, 2014). Bitz et al. (2012)—the first study to use the CCSM4 slab configuration—also used this convention.

3 ∼30% of the modern solar constant (Kopparapu et al., 2013a, 2013b).