Friday 4 December 2020

Simulating the orbital dynamics of the Kepler-1649 System and the climate of Kepler-1649c.

Over the past several decades, thousands of exoplanets have been discovered enabling statistical studies of exoplanet occurrence rate. These studies include the investigation of planets that lie within the Venus Zone and Habitable Zone. A primary purpose of such studies is to place our solar system within the context of planetary system architectures. One aspect of this context is determining the frequency of systems that have terrestrial planets in both the Venus Zone and Habitable Zone. Such systems may provide potential analogs for Venus and Earth, that in turn may provide key insights into the divergence of the Venus/Earth atmospheric evolutions. Examples of such systems include TRAPPIST-1, LHS 1140, and Kepler-186. Amongst the factors that influence the climate evolution of terrestrial exoplanets are the eccentricity variations, which may have contributed to the water loss of Venus. Thus, a pathway to understanding the major factors behind the Venus/Earth divergence is exploring the dynamical impacts of eccentricity variations that may have occurred (and continue to occur) in exoplanetary systems.

A particularly strong contributor to the currently known inventory of terrestrial exoplanets was the Kepler mission, from which a catalog of Habitable Zone planets was constructed. The Kepler-1649 system was discovered in 2017, which discovery identified the observed planet as a potential Venus analog due to the planetary size and similarity of insolation flux to that of Venus. The prospect of the planet as a Venus analog was strengthened by a set of detailed climate simulations that showed the rapid rise in surface temperature from a variety of initial conditions. The system became a strong comparative planetology target when a further terrestrial planet was discovered in the Habitable Zone earlier this year. Thus, Kepler-1649 has joined a select group of systems that are known to harbour terrestrial planets in both the Venus Zone and Habitable Zone. This class of system also provide opportunities to study the role of complex dynamical evolution on the orbital variation of planets, including those within the Habitable Zone.

In a paper published on the arXiv Database at Cornell University on 18 November 2020, Stephen Kane and Zhexing Li of the Department of Earth and Planetary Sciences at the University of California, Riverside, Eric Wolf of the Laboratory for Atmospheric and Space Physics at the University of Colorado, Boulder, and Colby Ostberg, and Michelle Hill, also of the Department of Earth and Planetary Sciences at the University of California, Riverside, present a dynamical and climate study of the Kepler-1649 system, describing the potential eccentricity variations that result from the interactions between the planets, and the possible impact on the climate of the Habitable Zone planet.

The Kepler-1649 system is particularly fascinating from a comparative planetology aspect, since it contains both potential analogs to Venus and Earth within the same system. The stellar properties of Kepler-1649A include its mass (0.1977 times that of the Sun), radius (0.2317 times that of the Sun), luminosity (0.00516 times that of the Sun), and effective temperature (3240 K, compared to 5778 K for our Sun). The radii of the known planets in the system are 1.017 times that of the Earth for Kepler-1649b and 1.06 times that of the Earth for Kepler-1649c. To calculate the semi-major axes of the planets (average distance at which the planets orbit), Kane et al. used information provided in the paper which described the discovery of Kepler-1649c, resulting in semi-major axes of 0.0842 AU (i.e. 8.42% of the distance at which the Earth orbit's the Sun) for Kepler-1649b, and 0.0827 AU for Kepler-1649c. The Conservative Habitable Zone boundaries are defined by the runaway greenhouse limit at the inner edge and the maximum greenhouse limit at the outer edge. The Optimistic Habitable Zone boundaries are determined empirically based on the assumption that Venus may have had surface liquid water as recently as 1 billion years ago, and that Mars may have surface liquid water 3.8 billion years ago. Based on these system properties, Kane et al. calculated the extent of the Conservative Habitable Zone and Optimistic Habitable Zone boundaries for Kepler-1649 as 0.075–0.147 AU and 0.059–0.155 AU respectively. The outer edge of the Venus Zone corresponds to the inner edge of the Conservative Habitable Zone (the runaway greenhouse limit) at 0.075 AU.

A top-down view of the Kepler-1649 system, where the star is location at the center of the crosshairs and the orbits of the known planets are indicated by solid circles. The extent of the Conservative Habitable Zone is shown as the light green region and the Optimistic Habitable Zone extension to the Conservative Habitable Zone is shown as the dark green region. Kane et al. (2020).

Previous descriptions of the Kepler-1649 system have not constrained the eccentricities in any significant way, though smaller planets tend to have small eccentricities. It has further been suggested that there could be an additional planet in an orbit that lies between the two known planets. This suggestion was partially motivated within the context of mean motion resonances that may exist within the system. The known b and c planets are close to the 9:4 mean motion resonance, which is a relatively weak resonance in terms of planetary perturbations. However, an additional planet between the known planets at a location of 0.063 AU would create a 3:2 resonant chain with high stability potential. Furthermore, an analysis of Kepler systems has found that planetary architectures may preferentially contain similar size planets, with the vast majority of Kepler planets having a spacing greater than 10 mutual Hill radii (the Hill radius of a body is the region in which it dominates the attraction of satellites). Assuming an Earth mass for all three planets, a planet located at the 3:2 resonant chain location between planets b and c would be 12.5 mutual Hill radii away from either planet, consistent with the observed spacings in other Kepler systems.

To investigate the prospect of an additional planet as a function of eccentricity, Kane et al. conducted a dynamical simulation using the REBOUND orbital dynamics package. Such an investigation requires an estimate of the planetary masses. There are numerous planetary mass-radius relationships that have been developed from which such an estimate may be derived For their dynamical simulations, Kane et al. adopted the mass estimates provided by the Forecaster tool, which produces masses estimates of 1.06 and 1.21 times that of the Earth for planets b and c respectively.

The dynamical simulation was carried out using the MEGNO (Mean Exponential Growth of Nearby Orbits) package within REBOUND with the symplectic integrator WHFast. MEGNO is a chaos indicator that is useful in dynamical simulations to distinguish quasi-periodic or chaotic orbital time evolution of planetary systems. The final numerical value returned by a MEGNO simulation indicates whether the system would end up in a chaotic state or not, where a chaotic state is less likely to maintain long-term stability. For their simulation, Kane et al. explored the possibility of an additional planet in between the orbits of planet b and c by testing how varying eccentricity and semi-major axis values for an injected planet of one Earth mass affects the stochasticity of the system. Kane et al. conducted four sets of simulations that varied eccentricity, each providing a grid of orbital parameters for the injected planet’s eccentricity (0.0–0.3) and semi-major axis (0.0482–0.0827 AU). Since the eccentricities of the two known planets are unconstrained, Kane et al. tested planet b and c eccentricities of 0.0, 0.1, 0.2, and 0.3. Each simulation duration was 10 million orbits for the outer most planet with a time step of 0.001 years (0.365 days). The integration was set to stop and return a large MEGNO value if the planet was ejected beyond 100 AU. 

There is a large parameter space between approximately 0.052 AU and 0.075 AU with low eccentricity values where the injected Earth-mass planet could be located with long-term stability. The largest range of eccentricities for the injected planet that allow a stable configuration occur at 0.063 AU which, as stated above, is the location of the 3:2 mean motion resonance chain. The additional spikes visible in the stability distribution, most obvious in the top-left panel, also correspond to mean motion resonance locations, such as the 5:4 resonance location at 0.071 AU. However, as the eccentricities for both planet b and c increase to higher values, the stable locations for the injected planet indicated gradually diminishes and the locations shift to higher eccentricity values. When the eccentricities for both planet b and c exceed 0.3, the stable locations for the injected planet vanish. Even for an eccentricity of 0.1, the system is most stable when the pericenters of the planets are approximately aligned. It is worth noting that the described simulations assumed system coplanarity, and the introduction of mutual inclinations may present additional regions of stability, particularly for scenarios that mitigate potential overlap of orbits.

A non-zero eccentricity for the known planets can place significant restrictions on the presence of an additional terrestrial planet between them. Kane et al. next investigate the time-dependent dynamical effects of non-zero eccentricities. Such investigations have been carried out in numerous other works in a more general context. Kane et al.'s N-body integrations were performed with the Mercury Integrator Package. The simulations used a hybrid symplectic/Bulirsch-Stoer integrator with a Jacobi coordinate system to provide increased accuracy for multi-planet systems. The time resolution of the simulations were set to 0.2 days in order to comply with the recommendation that resolutions be 1/20 of the shortest orbital period (8.69 days for planet b). 

Kane et al. performed a series of dynamical simulations using the above criteria that explore a range of eccentricities. Specifically, they started each simulation with a circular orbit for planet b and eccentricities for planet c that ranged from 0.0 to 0.5 in steps of 0.05. Each simulation was allowed to run for 1 million years, or 18 million orbits of planet c, with orbital parameters for both planets out-put every 10 simulation years. The suite of simulations that explored the full range of starting eccentricities for planet c were conducted using Earth masses for both planets b and c, and also using the Forecaster masses.

The full set of dynamical simulations demonstrated that the system is able to remain stable for 1 million years up to a starting eccentricity for planet c of 0.325, beyond which the system rapidly loses dynamical integrity. The 0.325 eccentricity upper limit was true for both the Earth mass and Forecaster mass scenarios. In particular, Kane et al. found that the eccentricity of planet c induces a near-equivalent eccentricity in planet b, whereby each  planet oscillates in eccentricity in order to conserve the system angular momentum.

To investigate this further, Kane et al. conducted a fourier analysis of each of the dynamical simulation outputs described above. The periods of the eccentricity oscillations, for both Forecaster and Earth masses, as a function of starting planet c eccentricity. There are several aspects to note about this. Firstly, the period of the oscillations decreases with increasing eccentricity. Secondly, the period of the oscillations is systematically smaller for the larger Forecaster masses. Thirdly, the dynamical instability that results beyond eccentricities of 0.3 is clearly visible and is likely related to the diminishing oscillation period with increasing eccentricity. These simulation results agree well with similar calculations that may be performed using the analytical expressions derived by Sam Hadden and Yoram Lithwick. It is worth noting that the period of the eccentricity oscillations and the eccentricity at which instability occurs is likely influenced by the proximity of the planets to the 9:4 mean motion resonance. Additional factors include the mutual inclinations of the planets combined with non-zero stellar obliquity. Regardless, the high frequency of eccentricity oscillations may result in significant climate consequences for the planets. 

Kane et al. discuss the implications of the variable eccentricity for the possible climate scenarios of planet c. Given the system parameters produced the insolation flux received by planet c at its semi-major axis is about 75% of the mean Earth insolation flux. The eccentricity oscillations described will result in a correspondingly variable insolation flux for the planet, whose amplitude and period depend on the eccentricity. For example, consider the case of a starting eccentricity of 0.2 for planet c. Since, the eccentricity of planet c never reaches zero, the maximum insolation flux is always larger than that for a circular orbit. The maximum insolation in this case oscillates approximately around an Earth equivalent flux with a variation of ±20%.

The Resolving Orbital and Climate Keys of Earth and Extraterrestrial Environments with Dynamics (ROCKE-3D), is a three-dimensional climate system model developed at the NASA Goddard Institute for Space Studies specifically for studying planetary and exoplanetary atmospheres. Kane et al. used ROCKE-3D to constrain the climatic effects of oscillations in the orbital-rotational properties of Kepler-1649c. ROCKE-3D has previously been used for studying the climate and habitability of terrestrial extrasolar planets, including studies of Proxima Centauri b, moist and runaway greenhouse states, the effects of planet orbital-rotational properties, and the effects of ocean circulation. While the atmospheric composition and surface properties of Kepler-1649c remain unconstrained, here we have conducted a set of simulations using ROCKE-3D to illuminate the possible role of eccentricity evolution on climate. For a starting eccentricity of 0.3, Kepler-1649c could experience eccentricity oscillations on 1200 year timescales. To constrain the bounds of the potential climate states experienced by Kepler-1649c during its dynamical evolution, Kane et al. simulated scenarios with 0 eccentricity and synchronous rotation, and with an eccentricity of 0.3 and a 3:2 resonant rotation state. Note, for these cases Kane et al.

were not simulating the continuous evolution of Kepler-1649c over many thousands of years, but rather were simulating the equilibrium climate states for potential end-member orbital-rotational configurations. However, if the spin dynamics become chaotic, these end-member rotation states may not be fully realised. Thus, Kane et al. also conducted simulations that assume a constant, non-resonant asynchronous spin rate in order to approximate chaotic spin-dynamics. For these cases, they assumed a 16 Earth-day rotation period, which is an intermediate value between the synchronous rotation period of 19.535 days, and the 3:2 resonant rotation period of 13.0235 days. Each of Kane et al.'s simulations achieve radiative equilibrium after 100–200 years, well within the timescales of possible eccentricity variations described. This is also comparable to the timescales provided by the continuously variable eccentricity approach.

Kepler-1649c lies within the Conservative Habitable Zone, close to the inner edge. In all climate simulations conducted, Kane et al. assume a globally ocean covered planet with no emergent continents. While this represents only one of many potential surface configurations for Kepler-1649c, the study of ocean covered worlds provides a first order approximation for the study of climate evolution on potentially habitable worlds.

Kane et al. tested two different assumptions for the ocean layer. 'Slab' ocean simulations assume a 100 m deep ocean layer with no horizontal transports, but with the thermodynamic exchange of heat and moisture between the ocean and atmosphere. 'Dynamic' ocean simulations utilise a 900 m deep fully coupled ocean circulation model capable of predicting horizontal and vertical transports, in addition to thermodynamic exchanges with the atmosphere. Both ocean configurations allow for sea-ice formation. Salinity is held equal to that of the modern Earth’s ocean. Atmospheric compositions also remain unconstrained from observations and therefor constitute a large uncertainty on the climate state. Kane et al. chose two nominal atmospheric compositions which are 'Earth-like' (1 bar total with 400 ppm carbon dioxide and an nitrogen background), and a 1 bar carbon dioxide atmosphere. For all climate simulations, water vapor and clouds are variable constituents in the atmosphere, controlled by temperature, evaporation, convection, advection, and condensation. For planets residing in the middle of the habitable zone, as does Kepler-1649c, an Earth-like atmospheric composition will result in a cold climate state, dominated by ice and where ocean transport processes are of particular importance, while a 1 bar CO2 atmosphere will result in a warm climate state, dominated by the hydrological cycle, clouds, and atmospheric transport. Thus the selection of these two atmospheric compositions allow Kane et al. to capture some diversity of habitable zone climates. A wide ranging parameter study of the infinite variety of possible surface and atmospheric compositions for Kepler-1649c is beyond the scope of this work. In total, we have conducted 16 different climate simulations of Kepler-1649c

Kane et al. produced surface temperature maps from our climate simulations with ROCKE-3D. The results used for these are from an average of the last 100 orbits of Kepler-1649c (i.e. 1950 Earth days, 5.3 Earth years). The bold white line indicates where the surface temperature is 273 K, and thus serves as a close proxy for the sea ice margin. The columns show climate results from different rotation states, including synchronous, resonant, and chaotic. Synchronous rotation assumes zero eccentricity and a rotation period equaling the orbital period. Resonant rotation assumes a rotation period in a 3:2 resonance with the orbital period and an eccentricity of 0.3. Chaotic rotation, as mentioned previously, is approximated by a 16 day asynchronous rotation period. The global mean surface temperature and surface temperature maps are largely unchanged for chaotic rotators with eccentricities of zero or 0.3.

Surface temperatures maps from Kane et al.'s climate simulations of Kepler-1649c. The left column assumes synchronous rotation with zero eccentricity, the center column assumed 3:2 resonant rotation with an eccentricity of 0.3, and the right column assumes chaotic rotation approximated by a 16 day asynchronous orbit. For the chaotic states, eccentricity does not have a significant effect on the climate and thus zero and 0.3 eccentricity results are averaged together. Rows (labeled) illustrate different surface an atmospheric configurations. All simulations performed with ROCKE-3D. Kane et al. (2020).

Global mean surface temperatures range between 244.7 K and 304.8 K, and surface liquid water is maintained over at least some of the planet surface in all cases simulated. Thus, all climates simulated are at least partially habitable, with many states being globally habitable. If the orbital-rotational state of Kepler-1649c can fully oscillate between synchronous and resonant modes, significant temporal changes to the global mean surface temperatures and the total habitable area of the planet, which here Kane et al. use the open ocean fraction as a proxy for the habitable surface area, would occur. The shift from synchronous rotation to resonant rotation affects climate both by changing the longitudinal patterns of stellar flux incident on the planet, and also by changing the Coriolis parameter which feeds back on atmospheric and ocean circulations. The breaking of synchronicity and the increase in rotation rate tends to result in increased zonal transports resulting in longitudinally banded patterns, as opposed to the classic eyeball climate states found for slow and synchronous rotators. The inclusion of non-zero eccentricity with a 3:2 orbital resonance produces a characteristic double hot-spot evident in temporal mean climate statistics, due to the fact that two opposite faces of the planet are periodically irradiated at perihelion. The double hot spot pattern is clearly evident in slab ocean cases. The orbital period of Kepler-1649c is short (19.5 days) compared to thermal timescales of the ocean, thus opposing sides of the planet remain permanently warm. For dynamic ocean cases, a faint signal of the double hot spot phenomenon persists, but zonal heat transports by the ocean results in nearly zonally uniform patterns of surface temperature.

For corresponding end-member scenarios of synchronous and resonant modes, global mean surface temperatures could change by as little as 8.9 K for an Earth-like atmosphere with dynamic ocean circulation, or by as much as 28.4 K for a carbon dioxide dominated planet with a slab ocean. For both atmospheric compositions, the climate sensitivity between end-member orbital-rotational states is reduced when considering a full dynamic ocean compared to the slab ocean cases. Perhaps surprisingly, the optically thicker 1 bar carbon dioxide simulations are more sensitive to the orbital-rotational state than are the Earth-like atmosphere. While further investigation is warranted, we surmise that this is due to more vigorous water vapor and cloud feedbacks that occur for our 1 bar carbon dioxide atmospheres, due to their generally warmer climate states compared to the Earth-like atmospheres which remain cold and ice dominated.

However, as noted above, these end-member scenarios may never be fully realised. The relatively rapid changes in orbital-rotational states may instead drive Kepler-1649c into chaotic rotation, which would result in asynchronous planet rotation and not into any favoured resonant state. In such a scenario, the climate patterns of Kepler-1649c become zonally uniform for all combinations of atmospheric composition and ocean assumptions, as star light would irradiate all longitudes of the planet equally with no preferred faces receiving more insolation on time-average than any other.

In all cases simulated, the change in surface temperature caused by the dynamical evolution is relatively small on the substellar hemisphere, but is much larger for the antistellar hemisphere. Kane et al. define substellar and antistellar hemispheres based on the synchronous rotating orbital-rotational reference frame. Naturally, for non-synchronous cases (resonant and chaotic), all sides of the planet receive periodic stellar flux, meaning that the opposing hemisphere is no longer permanently cloaked in darkness and enjoys its moment in the starlight, resulting in considerable warming. The substellar-antistellar dichotomy is most strongly seen for the slab ocean cases. In both synchronous and resonant rotating instances, the inclusion of dynamic ocean circulations effectively transports heat from the day-side to the night-side of the planet, muting the substellar-antistellar dichotomy.

Kane et al.'s 3-D climate simulations reveal that the modulation of orbital-rotational states does not trigger true climate catastrophes for Kepler-1649c (i.e. runaway glaciation or runaway greenhouse states) given the atmospheric compositions explored. If Keperly-1649c oscillates between synchronous and resonant states, significant millennial scale climate changes would occur that are far in excess of anything that the Earth has experienced over such a short timescale. Given the timescales of variation, life as we know it may struggle to adapt to such wide temperature variations indicated by Kane et al.'s simulations. Extant lifeforms may remain confined to the substellar hemisphere due its relative climatic stability compared to the significantly more intense variations that occur on the antistellar hemisphere. However, if instead Kepler-1649c settles into chaotic rotation, millennial scale climate changes would be negligible.

Other factors not explored here may also affect the sensitivity of Kepler-1649c’s climate to oscillations in its orbital-rotational state. Kane et al. have modeled only selected states in equilibrium. If, the orbital-rotational state of Kepler-1649c instead evolves monotonically between our illustrated end-member cases, the climates states may instead gradually smear between synchronous and resonant. Chaotic rotation states may also experience time-evolution of the rotation period. However, differences in non-resonant asynchronous rotation periods are not likely to drive significant changes in climate, compared to toggling into, and out of, synchronous and resonant states.

The atmospheric composition also presents a significant unconstrained parameter. While the atmospheric compositions studied here all result in temperate climates, if the atmosphere of Kepler-1649c were considerably different, the climate sensitivity between states may deviate considerably. For instance, if Kepler-1649c possesses an overwhelmingly dense, cloud and haze filled atmosphere, the surface may receive little direct starlight and become decoupled from changes in the patterns of stellar flux, and surface temperatures may trend towards a globally uniform state regardless of rotation state. Alternatively, if Kepler-1649c instead has a dry, clear, optically transparent atmosphere, the change in rotation state may also cause little difference to the surface temperature as the planet’s energy balance and thus climate would be governed simply by the albedo and emissivity properties of the planet’s solid surface, along with the instantaneous stellar flux input. Finally, our results indicate that the ocean treatment makes a considerable difference in the modulation of surface temperatures between synchronous and resonant states, but makes little difference for the chaotic rotation state. Still, changes in the ocean depth, and the presence and location of continents would significantly influence the results shown here through the modulation of ocean circulations. Jérémy Leconte argued that tidally locked planets with emergent continents will tend to have their continents migrate to the substellar or antistellar points. If continents cluster near the substellar point, recent work with ROCKE-3D suggests that this would mute the effects of ocean circulations. The varied location of continents remains a significant uncertainty in modern terrestrial exoplanet climate modeling, due to their tight coupling to ocean circulations, ocean heat transport, and ultimately climate.

A major focus of exoplanetary studies lies in the exploration of the full diversity of planetary system architectures. However, some of these architectures are similar enough to our own solar system such that they may shed insight into the orbital and atmospheric evolution of our terrestrial planets. Thus, systems that contain terrestrial planets within both the Venus Zone and the Habitable Zone are of particular interest for the investigation of potential Venus/Earth analogs that exist within similar yet different contexts of age, composition, host star type, and orbital dynamics.

The Kepler-1649 system is one of the few known systems to harbour such a potential Venus/Earth analog, where the similar age and, presumably, composition of the planets allows the removal of several free parameters from which other comparative studies may be conducted. Kane et al. have presented the results of dynamical analysis in which they have demonstrated that additional terrestrial planets could exist in orbits that lie between the two known planets. This is not an unlikely scenario given the dynamical packing observed for compact planetary systems and the necessity of coplanarity for transits to be observed. Additionally, the locations of significance resonance, such as the potential 3:2 >mean motion resonance chain, lend credence to a possible planet located between the known b and c planets. Kane et al.'s dynamical simulations further demonstrate that varying eccentricity can be used as a tool to exclude either eccentricities of the known planets or the presence of the hypothetical additional planet. 

Kane et al. also explored the effect of eccentricity on the dynamical evolution of the two known planets, the results of which show that the system is stable up to eccentricities of 0.3 and that the frequency of the eccentricity oscillations can be exceptionally high (1000–3000 years). Such high frequency eccentricity oscillations are a fundamental component of Milankovitch cycles for exoplanets and may play an important role in climate evolution.

Kane et al.'s ROCKE-3D simulations for the climate of Kepler-1649c under various eccentricity and ocean circulation assumptions demonstrate that in fact the surface temperature variation remains relatively small, even for eccentricity at the high end. The region of highest climate stability consistently lies in the substellar hemisphere for all of our simulations, indicating a potential longitude dependence of metabolic processes for any life that may exist on such a planet. Thus, the overall combination of both dynamical and climate evolution studies for this planetary system show that even though their orbits may change far more rapidly than that observed within the Solar System, atmospheric characterisation may yet reveal robust environments that maintain temperate surface conditions. Such dynamically active systems should therefore not be ruled out in the continuing search for potentially habitable worlds.

See also...

Follow Sciency Thoughts on Facebook.

Follow Sciency Thoughts on Twitter.