Quantifying the Influence of Jupiter on the Earth's Orbital Cycles

Through the coming decade, the next generation of astronomical observatories should yield a wealth of planets that will, to a greater or lesser extent, seem to resemble the Earth (e.g., Ricker et al. 2015; Arya et al. 2017; France et al. 2017). At this point, there will be a significant investment of observational resources attempting to search for evidence of biosignatures on those alien worlds (e.g., Des Marais et al. 2002; Segura et al. 2005; Kaltenegger et al. 2010; Rauer et al. 2011; O'Malley-James et al. 2014). But where should we look? By the time that we are capable of looking for the evidence of life on planets beyond the solar system, we will likely have a vast catalog of potential targets, from which the most promising must be chosen for that search (e.g., Turnbull & Tarter 2003; Lammer et al. 2009; Horner & Jones 2010; Kopparapu et al. 2014; Barnes et al. 2015; Cuntz & Guinan 2016; Agnew et al. 2017, 2018a, 2018b, 2019; Lingam & Loeb 2018).

The Exoplanet Era began in the latter stages of the last millennium, with the discovery of the first planets orbiting other stars (e.g., Campbell et al. 1988; Latham et al. 1989; Wolszczan & Frail 1992; Mayor & Queloz 1995). The first planets discovered revealed that the diversity of planetary systems was far greater than we had previously imagined. Giant planets were found orbiting perilously close to their host stars—a population of planets that became known as "hot Jupiters" (e.g., Mayor & Queloz 1995; Masset & Papaloizou 2003; Bouchy et al. 2005; Hellier et al. 2011; Albrecht et al. 2012; Wright et al. 2012). Planets were found orbiting pulsars (e.g., Wolszczan & Frail 1992; Thorsett et al. 1993; Bailes et al. 2011), and others were found moving on highly elongated orbits, dramatically different to anything seen in the solar system (e.g., Wittenmyer et al. 2007, 2017b; Tamuz et al. 2008; Harakawa et al. 2015).

As time has passed, every new generation of exoplanet discoveries has once again highlighted the great diversity of exoplanetary systems, driving home the concept that the planet formation process can yield an incredible variety of outcomes (e.g., Burgasser et al. 2010; Marcy et al. 2014; Masuda 2014; Johns et al. 2018). This result has been strikingly driven home by the results from the Kepler spacecraft, which carried out the first true exoplanetary census (e.g., Borucki et al. 2010; Batalha et al. 2013; Mullally et al. 2015). Among its many other noteworthy discoveries, Kepler revealed that "super-Earths" and "sub-Neptunes," classes of planet that are not represented in the solar system (e.g., Charbonneau et al. 2009; Vogt et al. 2010; Winn et al. 2011; Howard et al. 2012; Sinukoff et al. 2016), are common. Kepler also revealed that "dynamically packed" planetary systems are common, with the planets therein packed so tightly together that it would be impossible to have any other planets orbit between them (e.g., MacDonald et al. 2016; Mills et al. 2016; Gillon et al. 2017).

At the same time, we have seen a revolution in our understanding of the formation and evolution of the solar system. These developments include a new understanding of planet formation mechanisms (Adams 2010), the dynamical history and evolution of the orbits of solar system bodies, (Duncan & Quinn 1993; Tsiganis et al. 2005; Nesvorný 2018), and the prevalence of water and water delivery mechanisms in the early solar system (Encrenaz 2008). Such an understanding of our own planetary system is essential to placing exoplanetary systems in context, particularly given the diversity of orbital architectures that have been discovered in other planetary systems (Ford 2014; Kane & Raymond 2014; Batygin 2015; Hatzes 2016; Raymond et al. 2018). Many of the compact planetary system discoveries resulted from Kepler observations (Fang & Margot 2012), and it is expected that further insights into orbital architectures will result from the discoveries that will be made by the Transiting Exoplanet Survey Satellite (TESS). The parallel developments in both solar and exoplanet system science are therefore gradually converging upon a consistent picture of planetary system architectures.

The next generation of exoplanet observatories will yield a vast fresh catch of new discoveries. Given the evidence that terrestrial planets are common in the cosmos (e.g., Wittenmyer et al. 2011; Dressing & Charbonneau 2013), it is likely that new discoveries, such as those from TESS, will include many planets in the super-Earth regime or smaller (e.g., Ricker et al. 2015; Sullivan et al. 2015). At the same time, new instruments such as the James Webb Space Telescope (Gardner et al. 2006; Kempton et al. 2018) and the next generation of ultra-large ground-based telescopes (such as the European Extremely Large Telescope and the Giant Magellan Telescope; Gilmozzi & Spyromilio 2007; Johns et al. 2012) are expected to be able to deliver the first measurements that could properly characterize such planets, and potentially detect any evidence of life upon them (e.g., Beichman et al. 2014; Barstow & Irwin 2016; Schwieterman et al. 2016; Meadows et al. 2018).

There is, however, a problem. Simply put, by the time we are ready to search for evidence of life on suspected "Earth-like" worlds, we will have far more potential targets to study than available resources to study them. For this reason, it is imperative to prepare the playing field to examine the various facets that come together to make one planet more (or less) habitable than another (e.g., Horner & Jones 2010). By doing so, we hope to inform the target selection for the search for life beyond the solar system, helping researchers to target the most promising planets, and maximize the likelihood of a positive outcome.

In this light, a number of studies have begun to investigate the myriad factors that influence planetary habitability. Such studies range from those of the impact of stellar variability and binarity on the climates of terrestrial planets with potentially temperate surface conditions (e.g., Eggl et al. 2013; Haghighippour & Kaltenegger 2013; Kane & Hinkel 2013; Forgan 2014) to investigations of the role of giant planets in determining the impact rates on habitable zone (HZ) planets (e.g., Horner & Jones 2008, 2009, 2012; Horner et al. 2010; Lewis et al. 2013; Grazier 2016; Grazier et al. 2018) to studies of the impact of planetary orbital and spin dynamics on the potential climate of a planet (Armstrong et al. 2014; Linsenmeier et al. 2015; Deitrick et al. 2018b; Georgakarakos 2018). The planetary architecture plays a fundamental role in every one of those studies as gravitational forces induced mainly by planet–planet interactions alter the orbits of planets on astronomical timescales and thereby change the stellar radiation a planet receives—one among many factors that could affect the suitability of a given planet as a host for life.

Throughout most of the solar system's history, Earth's climate has remained within the range that allows liquid water to exist on the surface, and thereby accommodated the development of life (Mojzsis et al. 2001). That is not to say that climatic conditions on our planet remained unchanged. The long-term climate of Earth is primarily driven by geological processes such as plate tectonics and volcanism that alter the composition of the atmosphere. On shorter timescales, periodic variations in climate (collectively known as Milankovitch cycles) are superimposed on the long-term trend. The gravitational interactions with objects in our solar system induce systematic flexing and tilting of Earth's orbit over time, and affect the planetary axial spin dynamics (Milankovitch 1930; Berger & Loutre 1994; Laskar et al. 2004). Subsequently, the distribution and to a lesser extent the total amount of solar flux received at the top of the atmosphere varies over timescales greater than 10⁴ yr, causing cyclic modifications to the global and regional environment (Berger 1978; Berger & Loutre 1994).

Perhaps the most famous example of Milankovitch cycles on Earth are the sequences of glacial and interglacial periods that have been particularly pronounced during the Late Pleistocene (the past million years; e.g., Hays et al. 1976; Imbrie et al. 1992; Lisiecki & Raymo 2005). It should be noted that Earth's astronomical cycles have a relatively small amplitude and low frequency, and yet produce significant climate oscillations. Given that architectures of alien planetary systems are certain to be greatly different to that of the solar system, some "exo-Earths" are likely to experience astronomical forcing of greater variability (Deitrick et al. 2018a, 2018b). Such oscillations could have important implications for the suitability of a planet as a host for life, and may also impact the long-term survival of any existing life. Ample evidence exists that climate fluctuations on Earth have impacted the biodiversity, evolution, and migration of species (e.g., Bennett 1990; Jansson & Dynesius 2002; Van Dam et al. 2006). More extreme climatic variations (such as Snowball Earth transitions) may even have facilitated the explosion of new life forms during the Neoproterozoic (e.g., Hoffman et al. 1998; Kirschvink et al. 2000).

It is therefore clearly important to examine in detail the role that the architecture of planetary systems could have on the climates of potentially habitable worlds. Such studies could help to identify systems in which planets that would otherwise be considered eminently habitable could be ruled out as targets for the initial search for life beyond the solar system on the basis of extreme climate variability on astronomical timescales, driven by interactions between the planets in that system.

In this work, we consider the influence of the architecture of our own planetary system on the astronomical cycles that Earth experiences. Essentially, we treat the Earth as though it were a candidate exoplanet, and examine the variability of its orbital elements on the basis of the architecture of its host planetary system. To sample potential architectures for the system, we choose a methodical approach, considering the question: "how would Earth's Milankovitch cycles change if Jupiter moved on a different orbit?" We expect that variations in Jupiter's orbit will be especially important in altering Earth's orbital cycles because the two leading terms in Earth's eccentricity solution are related to Jupiter's orbital geometry (Laskar et al. 2004). In Section 2, we offer a short refresher on the various Milankovitch cycles experienced by the Earth, before describing our methodology in Section 3. We present our results in Sections 4 and 5, before discussing the implications of our results in Section 6. Finally, we consider the implications of our work for the Rare Earth hypothesis in Section 7, before drawing our conclusions in Section 8, with a discussion of the direction in which we intend to take our future work.

The long-term (>10⁴ yr) variability of Earth's climate is driven by intricate interactions between our planet's orbital evolution, spin dynamics, the spatial distribution of continents, oceans, variations in biogeochemistry, the occurrence of ice sheets, and many other factors (see, e.g., Cronin 2009). Being able to fully model such long-term climate behavior for newly discovered exoplanets is currently beyond us, as most of the required information will not be available. To put this complexity in context, we describe the current understanding of astronomically forced climate cycles experienced by the Earth.

Three astronomical cycles are typically discussed in the context of Earth's climate, namely eccentricity (e, the ellipticity of Earth's orbit), obliquity (, the angle of planetary axial tilt relative to the orbital plane), and climatic precession (e sinϖ, where ϖ is the longitude of perihelion⁵ ).

Climatic precession describes which hemisphere faces toward the Sun during a particular season and thus controls the spatial distribution of incident solar flux that is responsible for the seasonal contrast between hemispheres. When the Northern Hemisphere faces toward the Sun at perihelion, summers in this hemisphere are particularly warm, while winters are extremely cool. At the same time, the Southern Hemisphere experiences reduced seasonal disparity. A precession cycle is approximately 23 kyr long, meaning that the hemisphere experiencing maximal seasonal contrast alternates every ∼11 kyr. Eccentricity and precession are tightly connected in the sense that a perfectly circular orbit (e = 0) results in an equal amount of insolation received by both hemispheres, which minimizes interhemispheric differences. Even though the eccentricity of Earth may reach values close to zero every ∼100 kyr, hemispheric contrast never completely disappears. The oceans and continents have a different albedo, heat capacity, and thermal inertia; the asymmetric continental distribution therefore maintains the hemispheric contrast to some extent, even at times of near-zero orbital eccentricity.

The ∼100 kyr (short) and ∼400 kyr (long) eccentricity cycles not only modulate climatic precession effects, they also slightly alter the total amount of insolation Earth receives. The annual mean solar flux (F) at the top of the atmosphere scales with orbital eccentricity such that F ∝ (1−e ${}^{2}{)}^{-0.5}$. Earth's eccentricity varies between ∼0 and 0.06, which is equivalent to a difference in radiative forcing of ∼0.5 W m⁻² on astronomical timescales (Berger 1978; Berger & Loutre 1994; Laskar et al. 2004).

Despite these small variations in eccentricity, power spectra of paleoclimate data often reflect a strong imprint of eccentricity frequencies. The glacial–interglacial cycles of the late Pleistocene appear to be eccentricity-paced, while ∼100 and ∼400 kyr cycles are also prominently present in mid- and early Cenozoic and Mesozoic geological records (e.g., Herbert & Fischer 1986; Lourens et al. 2005; Kirtland Turner 2014). The discrepancy between the small changes in annual mean insolation and the relatively large climate consequences that result could imply that processes intrinsic to Earth may lead to a highly nonlinear response between insolation and climate (Hays et al. 1976; Clemens & Tiedemann 1997). Positive feedbacks, for example the ice-albedo feedback (Curry et al. 1995; McGehee & Lehman 2012), may drive the climate state of a planet into its extremes, while negative feedbacks, such as CO₂ consumption through weathering (Walker et al. 1981) can dampen the effects of cyclic insolation variation on timescales greater than 10⁴ yr.

We should note that there are still many unknowns in the theory of Milankovitch forcing. One of the main ideas is that the ∼100 kyr glacial cycles of the late Pleistocene are not directly related to eccentricity, but rather indirectly through the modulation of precession effects (Raymo 1997; Ridgwell et al. 1999; Maslin & Ridgwell 2005). The long (∼10⁵ yr) residence time of carbon in the oceans and atmosphere may also transfer power from precession frequencies to the longer-eccentricity frequencies (Zeebe et al. 2017). Alternative explanations for the ∼100 kyr cycles exist that invoke intrinsic feedback mechanisms that could explain a natural climate oscillation on ∼100 kyr timescales due to inertia in parts of the Earth system that are completely unrelated to astronomical forcing (Saltzman & Sutera 1987; Wunsch 2003; Nie et al. 2008). At the same time, it has been suggested that obliquity forcing, rather than eccentricity variation, is the main driver for climate oscillations (e.g., Huybers & Wunsch 2005; Raymo et al. 2006), on the basis that obliquity has a more direct impact on the mean insolation at high latitudes than does eccentricity. It should be noted, however, that variations in eccentricity would serve to modulate the impact of such obliquity variations.

The Earth's obliquity, or axial tilt, oscillates between 215 and 245 with a period of ∼ 41 kyr. This expands and shrinks the polar circles, alters high-latitude insolation patterns, and subsequently leads to the waxing and waning of ice caps; 41 kyr cycles in paleoclimate records are particularly pronounced in the Pliocene and early Pleistocene (Gildor & Tziperman 2000; Naish et al. 2009; Lourens et al. 2010).

Most likely, all of the above suggestions play a role in the periodic climate oscillations on Earth, but the relative importance of those elements shifts depending on the state of the planet during a given period in geological history. For instance, the presence and distribution of ice sheets (Raymo et al. 2006), or the background climate state (Berger et al. 1999) could shift the spectral signal of the oscillations.

It is clear that astronomical forcing contributed to the oscillations in Earth's climate not only in the more recent time when polar ice caps could amplify the effects of insolation changes, but also in the more distant geological history when the Earth surface was mostly devoid of ice, despite the minor amplitudes of the eccentricity and obliquity cycles of Earth (e.g., Zachos et al. 2001).

Given the expected wide diversity in exoplanetary system architectures (Winn & Fabrycky 2015), it is reasonable to assume that some planets orbiting other stars undergo more extreme orbital variations and therefore experience dramatic spatiotemporal insolation changes, potentially having important implications for the climatic environment of an exoplanet and its habitability. Since the obliquity and precession dynamics of terrestrial exoplanets cannot currently be observationally constrained, we focus here on the orbital parameters that are measurable, namely the planet's orbital eccentricity and inclination (the tilt of the orbit relative to the reference plane). We would like to stress, however, that the orbital and spin dynamics are intimately connected. Any changes in the orbital inclination or eccentricity due to planet–planet interactions will directly translate, albeit nonlinearly, into perturbations to the spin dynamics of that planet (e.g., Kinoshita 1975; Laskar et al. 2004).

The calculations to estimate the evolution of the planetary spin dynamics are extremely complex and output is highly sensitive to the multitude of input variables. Not only is the angular momentum altered by variations in the orbital motion, it is also affected by tidal forces (Laskar et al. 1993) and internal planetary processes such as core–mantle friction (Néron de Surgy & Laskar 1997), atmospheric tides (Barnes et al. 1983; Volland 1996), mass displacement from plate tectonics or mantle convection (Ward et al. 1979; Forte & Mitrovica 1997), and climatic friction (Dehant et al. 1990; Rubincam 1990; Ito et al. 1995). Such information will likely never be obtained for exoplanets.

For that reason, we explore the manner in which the orbital evolution of Earth responds to architectural changes to the solar system. Specifically, in this work we consider the impact of the orbit of the solar system's most massive planet, Jupiter, on the Earth's orbital cycles.

To examine how Earth's orbital Milankovitch cycles would be altered were the solar system's architecture markedly different, we carried out an extensive suite of n-body simulations, using the Hybrid integrator within the dynamical integration package Mercury (Chambers 1999). Each individual simulation followed the dynamical evolution of the eight solar system planets for a period of 10 Myr, with the instantaneous orbital elements of the planets written to file every 1000 yr. An integration time step of one day was used to ensure that the integrations were as accurate as possible.

In addition, Mercury was modified by the addition of a user-defined force in order to take account of first-order post-Newtonian relativistic corrections (Gilmore & Ross 2008). This correction allows the orbital behavior of the innermost planets to be modeled accurately, ensuring that our results fairly reflect the physical reality of the orbital evolution of the planets in question.

A total of 159,201 individual simulations were carried out, in which the initial orbits of the planets Mercury, Venus, Earth, Mars, Saturn, Uranus, and Neptune were the same. The only elements that were changed from one simulation to the next were the initial semimajor axis, a, and eccentricity, e, of Jupiter's orbit. The process by which the simulations were created followed that established for studies of the stability of exoplanetary systems and solar system small bodies (e.g., Horner & Lykawka 2010; Horner et al. 2012; Wittenmyer et al. 2012, 2017a), such that the semimajor axis and eccentricity of Jupiter were each varied in regular steps, creating a rectangular grid of solutions in a–e space for Jupiter's orbit.

The orbital elements for the planets were obtained from the Horizons DE431 ephemeris, and converted from Cartesian to Keplerian coordinates through Mercury, by running a short simulation and outputting the Keplerian elements at t = 0. These Keplerian elements were then used as the basis for our grid of Jovian orbital solutions. In total, we tested 399 unique values of Jupiter's semimajor axis, covering a 4.0 au range centered on the initial semimajor axis of Jupiter in "our" solar system (a ∼ 5.203 au). For each of these 399 semimajor axes, we tested 399 unique values of Jovian orbital eccentricity, ranging from circular orbits (e = 0.0) to ones with moderate eccentricity (e = 0.4). While such an upper bound might seem unusually high, we note that there is a growing body of exoplanets that have been found on orbits far more eccentric than this (e.g., Wittenmyer et al. 2007, 2017b; Tamuz et al. 2008; Harakawa et al. 2015). Even in the solar system, it has been suggested that the eccentricity of Mercury's orbit can exceed this value, as part of its own long-term periodic variability (e.g., Strom & Sprague 2003; Correia & Laskar 2009).

Individual simulations were halted early if any planet was so perturbed that it collided with the Sun, another planet, or reached a heliocentric distance of 40 au. We flagged those simulations containing architectures that proved dynamically unfeasible, and recorded the time within the integration at which the simulations were halted. Once the remaining (stable) simulations were complete, we extracted the evolution of Earth's orbital elements at 1000 yr intervals to determine the frequency and amplitude of the variations that occurred in Earth's orbit.

In Figure 1, we show the evolution of the Earth's orbital elements over the last million years of our simulations for three exemplar cases. The three scenarios feature "Jupiters" initially located at 3.203 au (left), 5.203 au (center), and 7.203 au (right), while Jupiter's initial orbital eccentricity was set to 0.0 and initial values of the four other orbital elements were set to their canonical values. This figure highlights the degree to which changes in Jupiter's orbit can impact both the amplitude and frequency of the Earth's orbital cycles.

Zoom In Zoom Out Reset image size

Of 159,201 unique realizations of the solar system simulated in this work, the majority (∼74%) proved unstable. The stability of the solar system in our simulations is shown in Figure 2. It is clear, particularly when the data are plotted on a linear scale, that the stability of the solar system is a strong function of Jupiter's orbital eccentricity, but even at low eccentricities there are regions where no stable solutions are found. Equally, two narrow regions can be seen where stable solutions exist across the full range of eccentricity tested in this work—the result of the stabilizing influence of resonant interactions between Jupiter and Saturn.

Zoom In Zoom Out Reset image size

From Figure 2, it is apparent that instability occurs on a variety of timescales, and the further we move from the "stable" solutions in our work, the more rapidly things become unstable. As a result, it seems highly likely that scenarios on the fringe of stability, which in our simulations proved stable on 10 Myr timescales, would likely become unstable on timescales of tens, hundreds, or thousands of millions of years. While such instability would be of great interest, the computational challenges involved in integrating our data for one or two orders of magnitude longer make a detailed study of such edge cases impractical. Nonetheless, we caution readers that the true extent of the "unstable" region is likely slightly larger than that seen in our data, as a result of the dynamical timescale over which our integrations were performed.

For the rest of our analysis, we ignore the unstable regions of the plot, and solely focus on those regions where the system remains dynamically stable. For each of those stable solutions, we can determine the degree to which the Earth's orbital parameters vary with time, in both their amplitude and frequency. The results are visualized in Figures 3–10.

Zoom In Zoom Out Reset image size

Spectral analysis was performed on each of the stable dynamical simulations to extract the primary and secondary periods and amplitudes for Earth's orbital elements. The R package "Astrochron" was used to perform a multitaper method spectral analysis of our data using a time-bandwidth product of 2.0 (Meyers et al. 2014). To confirm the significance of spectral peaks, autoregressive (AR1) red noise models were generated for each individual time series. Spectral peaks were considered to be statistically significant when they rose above the 99% confidence level (Figures 3 and 4). The maximum amplitudes of the two most pronounced and statistically significant cycles were estimated by applying a bandpass filter to the original time series data. The width of the filter is 10% of the period that is associated with the cycle. For instance, a cycle with a period of 100 kyr will have a bandpass width of 10 kyr, covering 95–105 kyr periods. All significant cycles that have periods within this range are considered to have a similar origin. In this manner, we accounted for the temporal variability of a cycle over the simulated 10 Myr.

Zoom In Zoom Out Reset image size

For validation, we applied these methods to the simulation in which Jupiter's location and eccentricity most closely resemble its "real" values (Figure 3). The most pronounced cycles that result for Earth's eccentricity have periods of 400, 123 and 94 kyr, consistent with earlier calculations (e.g., Laskar et al. 2004). Likewise, the associated simulated amplitudes for Earth's eccentricity are also in good agreement with the expected values.

Spectral analysis reveals that the periods in Earth's orbital cycles are mainly affected by the semimajor axis of Jupiter a_J (as can be seen in Figures 4, 8, and 9 moving from left to right within a subplot). The influence of Jupiter's eccentricity e_J is minor (e.g., Figures 8 and 9, moving from bottom to top within a subplot). To examine this in more detail, we focus on the variability in Earth's orbital parameters over a range of a_J, while keeping e_J constant at zero (Figure 4).

When Jupiter is initially placed close to the Sun at a_J < 4 au, Earth's eccentricity cycles predominantly have short periods of 50–100 kyr. In scenarios where Jupiter is more distant, the dominant eccentricity cycles tend to have periods of 100–150 kyr and 1–2 Myr, as can be seen in panels (a) and (b) of Figure 4 and in Figure 8. These cycles appear to be disrupted when Jupiter is moved to 4.1 au, but this does not result in system instability within the simulated 10 Myr. However, the cycles in Earth's eccentricity at this location exhibit aperiodic and relatively large-amplitude variations that suggest the system may become truly dynamically unstable after the 10 Myr of our integrations.

Rapid changes in the long-period oscillations in orbital eccentricity occur around 4.25, 4.6, 4.8, 5.05, and 5.7 au, demonstrating that the periods of the longest-eccentricity cycles are sensitive to minor changes in Jupiter's semimajor axis just before and after an unstable region. This does not apply for the shorter cycles that are relatively stable throughout. For a_J ∼ 5.2–7.2 au, the periodicity of the short-eccentricity cycles remains at ∼140 kyr. The long-eccentricity cycles vary from 500 kyr to 2.5 Myr, depending on Jupiter's location.

Cycles in Earth's orbital inclination are generally shorter than eccentricity cycles (Figure 4, panels (c) and (d)). The period of the dominant cycle increases from 30 to 85 kyr when Jupiter is moved outward. The disruption around 4.1 au is clearly visible. Multiple short cycles converge at this point, while a strong peak appears in the longest cycle that rapidly rises from 500 kyr to 2.5 Myr, before falling back to ∼500 kyr. At this location, the Earth's inclination cycles exhibit large oscillations (sometimes greater than 10°). However, in contrast to the aperiodic variability in Earth's eccentricity, these cycles appear to be steady over time throughout the 10 Myr simulation.

Aside from the most dominant 30–85 kyr cycle, two other orbital inclination cycles occur that arise in the vicinity of the 4.1 au inconsistency. We observe a 100 kyr cycle that remains approximately constant, and a cycle that starts with a period of ∼100 kyr and gradually increases to 500 kyr as the Jupiter–Sun distance increases to 7.2 au.

Multiple significant cycles coexist in the variability of Earth's orbital eccentricity and inclination at any given Jupiter–Sun distance. This is because Jupiter is not the only planet that affects Earth's orbit—the other planets in the solar system also interact gravitationally with the Earth, as well as with Jupiter. In scenarios where Jupiter is initially placed close to the Sun, the periods of Earth's orbital eccentricity and inclination cycles are short and remarkably similar. Two orbital parameters with similar periods can potentially lead to interesting climatic behavior as external forcing and subsequent climate feedbacks may either cancel out or reinforce each other, depending on the phase of the variation of the parameters. As the parameters move in and out of phase, their combined effects could "beat," causing significantly larger variability than would otherwise be expected.

The degree to which Jupiter influences Earth's orbital cycles is clearly apparent upon examination of the three exemplar simulations shown in Figure 1. Both the amplitude and frequency of the cyclic variations in Earth's orbital elements are affected by Jupiter's orbit, which would change the spatiotemporal distribution of incident solar flux and thereby alter the seasonal patterns experienced by Earth.

In Figure 5, we show the impact of Jupiter's orbit on the time evolution of Earth's orbital eccentricity. The regions in black in those plots are those for which the solar system proved unstable in our simulations. The data plotted in the top left panel of that figure reveal that, through most of the phase-space studied, Earth's maximum orbital eccentricity remains low, less than 0.1. However, there are small regions in which the solar system is stable on the timescale of our integration, but Earth's eccentricity can be forced to larger values. In the main, these regions are those on the edge of stability, and it is quite possible that those systems would have become unstable with longer simulation time. Of perhaps more interest are the panels in Figure 5 showing the rate of change of Earth's orbital eccentricity. As Jupiter is moved closer to the Sun, both the mean and maximum rate of change of our planet's orbital eccentricity increase, a finding mirrored in the left-hand column of Figure 1.

Zoom In Zoom Out Reset image size

Given information on how Earth's orbital eccentricity changes with time, it is possible to quantify the degree to which the annual mean insolation on our planet changes. The higher Earth's orbital eccentricity, the higher the insolation averaged over the course of the year, as described in Equation (1) (where ${\bar{F}}_{\mathrm{annual}}$ is the insolation averaged over the course of the year, a is Earth's semimajor axis, e the Earth's eccentricity, and L_⊙ is the luminosity of the Sun). For that reason, in Figure 6, we show the degree to which variations in Earth's orbital eccentricity would drive variations in annual mean insolation. We plot the difference between the highest and lowest values for annual mean insolation that would be experienced across the course of our integrations at a given location. While the structure shown in these plots clearly follows the maximum eccentricity obtained by the Earth in our simulations, the nonlinear response of insolation to changes in eccentricity can be clearly seen. As a point of reference, the current situation on Earth is that the difference in insolation between eccentricity maxima and minima is less than 0.5 W m⁻². The most extreme scenarios tested in this work, for which Earth exhibits the greatest orbital eccentricities, would deliver changes in annual mean insolation that exceed those seen on Earth by more than an order of magnitude—a result that would have interesting implications for the planet's climate.

$$\begin{eqnarray}&&{\bar{F}}_{\mathrm{annual}}=\displaystyle \frac{{L}_{\odot }}{(16\pi {a}^{2}\sqrt{(1-{e}^{2})})}.\end{eqnarray} \tag{ 1 }$$

Zoom In Zoom Out Reset image size

The evolution of Earth's orbital inclination as a function of Jupiter's initial orbit is shown in Figure 7. In general, the maximum inclination to which Earth's orbit is excited remains low throughout the plot. The three exceptions to this are the regions just beyond 4 au, at around 4.75 and 6.3 au. At each of these locations, the Earth's orbit can experience significant inclination variability. The innermost of these regions also exhibits enhanced rates of inclination variability, as can be seen in the three other panels of Figure 7. The outermost band of increased i_max values, however, shows no such feature; here, the rate of change of inclination is no greater than in the regions surrounding it.

Zoom In Zoom Out Reset image size

These regions offer a cautionary tale for the assessment of potential exoplanet habitability: just because the orbit of an Earth-like planet is stable, that does not mean that the oscillations in its orbit cannot be relatively large. Such large excursions in the orbit of a dynamically stable planet could well have an impact on the evolution of its climate. To illustrate this, consider the impact of the planet's orbital inclination.

Given that the seasons are driven by the inclination of a planet's axis with respect to the plane of its orbit, rather than the absolute orientation of that axis in space, changing the inclination of a planet's orbit while the absolute orientation of the axial tilt remains fixed would result in an equivalent change in the tilt of the planet's axis with respect to the plane of its orbit. As a result, the extent of the polar circles would vary dramatically on a planet experiencing changes in orbital inclination greater than 10°, which would result in severe changes in seasonal insolation. Would such a planet still be habitable? Perhaps, but it is probably fair to argue that such a world would be less promising as a target for the search for life than another that did not exhibit such extreme variability.

The rates at which the Milankovitch-like oscillations occur may be key here. It has been suggested that rapid obliquity and eccentricity cycles could result in extreme seasonal variations that could trigger intense ice age cycles (e.g., Deitrick et al. 2018b). On the other hand, it has been suggested that such rapid oscillations may suppress the ice-albedo feedback, and act to expand the outer edge of the HZ (e.g., Armstrong et al. 2014). It seems plausible that each of those scenarios could occur. Which one is more likely to happen depends on the inertia of the system and the rate of change. Planets with a thicker atmosphere or larger ocean volume may, for example, respond more slowly to periodic alterations in the incoming stellar radiation and might therefore prove to be more resistant to oscillations in astronomical forcing than planets with thinner atmospheres and shallower oceans (e.g., Cowan et al. 2012).

Figure 8 shows the period (left) and amplitude (right) of the two dominant oscillations in Earth's orbital eccentricity, with Figure 9 showing the same information for the evolution of Earth's orbital inclination. In broad terms, Figures 8 and 9 show the same features that can be seen in Figure 4. Of particular interest is the "phase-change" around 4 au, from a regime where both dominant periods occur on timescales <100 kyr to one that features both long- and short-period components.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Variations in the incident flux received by a planet can affect its climate evolution and that of its atmosphere, which both influence overall planetary habitability. For example, the notable effect of eccentricity on climate has been previously studied for a variety of scenarios and specific exoplanets (Williams & Pollard 2002; Kane & Gelino 2012; Way & Georgakarakos 2017; Georgakarakos et al. 2018). The effects of obliquity and eccentricity on the incident flux for a planet have been quantified with application to exoplanets by Kane & Torres (2017). An exploration of how Milankovitch cycles affect exoplanet climates was recently undertaken by Deitrick et al. (2018a, 2018b), showing the significant impact on obliquity variations and subsequent effect on planetary climate.

This work is particularly important for the evaluation of the variation in flux for planets that lie within the HZ of their host stars (Kopparapu et al. 2013, 2014; Kane et al. 2016). Given the diversity of exoplanet systems and their corresponding dynamics, it is possible and likely that combinations of orbital parameters exist which would produce stellar flux variations that render a HZ planet uninhabitable over its lifetime. It is also possible for the habitability of a planet to be time-dependent, and to be driven by periodic variations like Milankovitch cycles. Studies such as ours therefore play an important role: they allow the current dynamical state of the system to be assessed, but can also be used to examine the likelihood that the planet has remained "habitable" on astronomically long timescales.

In addition to studying the short-term habitability of a system, simulations such as ours can also reveal the degree to which a given Milankovitch regime is robust against the migration of the planets in that system. For example, our simulations suggest that, were Jupiter to migrate just a short distance to ∼5 au from the Sun, the solar system could become catastrophically unstable. Such migration-driven instability is not a new concept in the narrative of solar system evolution, having been invoked in the past to explain the Late Heavy Bombardment of the inner solar system (e.g., Gomes et al. 2005; Morbidelli et al. 2010; Walsh et al. 2012).

Even on smaller scales, the fine structure visible in the plots of Earth's orbital evolution suggests that even relatively small-scale migration could cause marked shifts in the potential habitability of a planet. Given that Jupiter is continually ejecting cometary and asteroidal material from the solar system (e.g., Horner et al. 2003; Horner & Jones 2008, 2009; Grazier et al. 2018, 2019), it must still be undergoing a very gradual inward migration. The same will no doubt be true of planets orbiting other stars; nothing is truly static on astronomical timescales.

By performing simulations like those presented herein for potentially habitable exo-Earths, it may even be possible to identify those that might have experienced catastrophic or chaotic Milankovitch cycles in the past. Whether such periods of instability would be deleterious or beneficial to the development of detectable life is still open to debate. Compare, for example, the impact of the ancient episodes of "Snowball Earth," which have been suggested as potentially contributing to the explosion of new species during the Neoproterozoic (e.g., Hoffman et al. 1998; Kirschvink et al. 2000), to the effect that such an event would have on life at the modern era (which would likely be catastrophic, given that even the relatively minor recent glaciations have had a marked effect on the Earth's biodiversity; e.g., Williams et al. 1993; Hewitt 2000, 2004; Sniderman et al. 2013; Nadachowska-Brzyska et al. 2015). Despite the complexity of such situations, the ability to identify such "edge cases" may help focus the choice of the most promising targets in the search for life, and so should definitely be considered in future work.

In this work, we have solely considered the impact of the orbital parameters (eccentricity, inclination, apsidal precession) that contribute to variations in the mean annual flux received at the top of a planet's atmosphere. Elements in the axial group (obliquity and axial precession) clearly play an important role in determining the spatial distribution of incident flux on a planet, but these variables cannot be measured for exoplanets at the present time. Observations to determine those angles will be as challenging as those required to search for evidence of life on the planets considered. As our work is intended to help guide the selection of the targets for such observations, our focus lies on those parameters that we might reasonably measure in the near future.

Nonetheless, it is interesting to briefly view our results in the context of the real solar system. In the past, it has been argued that the Earth is unusually favorable for the development of life—a core tenet of the "Rare Earth" hypothesis (e.g., Ward & Brownlee 2000; Waltham 2019). In that light, it is interesting to consider the degree to which the Earth is unusual in the context of our simulations. Are the orbital variations experienced by our planet unusual or typical, when compared to our ensemble of "alternative Earths"? We stress that such a study is purely illustrative, since our results take no account of the impact on our planet's climate from elements in the axial group.

Zoom In Zoom Out Reset image size

To consider what fraction of the stable simulations would be equally as or more clement than our solar system (with low-amplitude and low-frequency oscillations), we calculated the distribution of frequencies and amplitudes for each of the orbital parameters that impact the Milankovitch cycles. The results of this analysis are shown in Figure 11. The two uppermost panels show the distribution of the periods and amplitudes of the two strongest periodic oscillations in Earth's orbital eccentricity, across our runs. It becomes apparent that the periods of oscillation span a broad range. In total, the fraction of Earths whose dominant eccentricity cycle had a period longer than Earth's 400 kyr cycle was 31%, while 51% of simulations featured a secondary eccentricity cycle with period greater than Earth's 100 kyr periodicity. The amplitude of the dominant eccentricity cycle was smaller than that seen in our solar system in 52% of cases, while the secondary cycle amplitudes were smaller than those we experience in 67% of cases. In other words, when it comes to the variability of our orbital eccentricity, the Earth seems unremarkable.

Zoom In Zoom Out Reset image size

In a similar fashion, we can examine the evolution of Earth's orbital inclination and the longitude of its perihelion in the context of the ensemble of stable runs. Once again, the Earth's orbital evolution is relatively typical of the ensemble. A total of 58% of systems featured primary inclination variability on timescales longer than that experienced by the Earth, with 29% showing secondary inclination periodicity on longer timescales; 47% of systems feature primary inclination oscillations smaller than those for our planet, with 50% of systems having smaller secondary oscillations.

The rate at which Earth's longitude of perihelion precesses is at the more sedate end of those seen in our simulations. In total, 80% of stable scenarios featured "Earths" whose dominant perihelion precession rate was more rapid than seen for our planet, with 67% exhibiting more rapid precession for the second strongest periodicity.

Finally, we note for posterity the extremes in the variability observed across our dynamical stable ensemble of simulations. Across those 41,652 runs, the largest eccentricity obtained by the Earth was 0.415 in the scenario where Jupiter was initially located at a semimajor axis of 4.04 au and had an eccentricity of 0.168. Such an extreme eccentricity value for the Earth would result in an increase in the annual mean solar flux of around 10%, or 24 W m⁻² when also accounting for our planet's modern albedo—a significant change compared to the 0.2% variation received by Earth over the eccentricity timescales associated with the glacial–interglacial fluctuations.

To put this into the context of the Earth's modern climate sensitivity, we note that the estimated radiative forcing that would be caused by a doubling of atmospheric CO₂ is 3–4 W m⁻², a value which takes into account various climate processes such as atmospheric water vapour, cloud, and lapse rate feedbacks (Andrews et al. 2012; Huang & Bani Shahabadi 2014). An increase of 24 W m⁻² would be comparable to the rise in global temperature that would result from eight doublings of atmospheric CO₂.

It is conventionally assumed that global temperatures increase by 0.8 K per W m⁻². An increase of 24 W m⁻² would therefore equate to a global warming of approximately 20 K. Similar estimates have been obtained in previous climate modeling studies that assessed temperature variability for Earth-like planets on orbits of various eccentricities (Williams & Pollard 2002; Dressing et al. 2010). Such major global temperature changes are unprecedented in the Phanerozoic era of Earth (past 500 Myr), based on estimates from paleoclimate records (e.g., Royer et al. 2004; Hansen et al. 2013).

In the coming years, exoplanet experiments will produce an improved sensitivity to Earth-sized planets orbiting within the HZ of solar-type stars. Such planets will become the targets in our efforts to search for evidence of life beyond the solar system. However, such observations will be immensely challenging and, in the early stages of that search, we will only be able to study a small subset of the discovered planets. It is therefore vital that we prioritize which of those planets are the most promising targets for intensive follow-up observations. Such prioritization will consider many factors that come together to render one planet more or less habitable than another (e.g., Horner & Jones 2010). An important example of such a factor is the orbital dynamics of the planet and the subsequent effects on the planetary climate.

In this work, we detail the results of a large suite of n-body simulations designed to examine the influence of Jupiter on the orbital cycles experienced by the Earth. We systematically varied the initial orbit of Jupiter across a region spanning ±2 au in semimajor axis around the orbit of Jupiter in our solar system. At each unique semimajor axis we tested, we varied the initial orbital eccentricity of the giant planet in the range 0.0–0.4. This yielded a grid of almost 160,000 unique variants of our solar system, each of which featured the other seven planets (Mercury, Venus, Earth, Mars, Saturn, Uranus, and Neptune) moving on identical orbits—those that they occupy in the solar system. We simulated the evolution of these modified solar systems over a period of 10 Myr, and examined how Earth's orbit subsequently varied over time.

Our results reveal the sensitivity of the solar system's stability to Jupiter's orbit, with some ∼74% of the variant systems proving catastrophically unstable within the 10 Myr of our integrations. For the subset that proved stable, we found that both the periods and amplitudes of the oscillations in Earth's orbital elements varied markedly as a function of Jupiter's initial orbital elements. When Jupiter began on an orbit closer to the Sun, the periodicity of Earth's orbital element variation was typically shorter than when Jupiter was more distant. Simultaneously, the amplitude of the Earth's orbital cycles varied as the giant planet was moved through the solar system, with some stable solar system variants featuring oscillations in Earth's orbital inclination that approached, or even exceeded, ten degrees.

Our work highlights the degree to which small changes in the architecture of a planetary system can drive large variations in the Milankovitch cycles that would be experienced by any potentially habitable planets therein. It will therefore be critically important to accurately determine the orbital elements for all planets within the system, in order to assess the scale and frequency of their orbital cycles so that the impact of Milankovitch forcing on terrestrial planets can be considered.

When the amplitudes and periods for the dominant cycles in Earth's orbital eccentricity, inclination, and the precession of Earth's perihelion are compared with the results across our ensemble of stable solutions, we find that our planet's orbital behavior is remarkably unremarkable. The variations in the Earth's orbital parameters are neither unusually fast nor unusually slow, neither unusually large nor small. As such, at least when it comes to the orbital Milankovitch cycles, it seems that the central tenet of the "Rare Earth" hypothesis does not hold true. The Earth is not unusual—and so it may be that planets with similar Milankovitch cycles are common in the cosmos.

The transit and radial velocity follow-up of TESS planets will likely reveal a number of other planets in the same systems, which will enable studies such as the one described here to be carried out for those new systems. Beyond simply allowing us to compare the potential habitability of those planets, such studies will form a fascinating complement to the atmospheric observations carried out by the James Webb Space Telescope and other facilities. The potential combination of our being able to study both the orbital architectures and atmospheric properties for these TESS planets will offer an unprecedented opportunity to directly examine the influence of orbital dynamics on planetary atmospheres outside of our solar system.

This research has made use of the Habitable Zone Gallery at hzgallery.org. This research has also made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. The results reported herein benefited from collaborations and/or information exchange within NASA's Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA's Science Mission Directorate. P.V. and S.K.T. were supported by a Heising-Simons Foundation award while this research was ongoing. The authors wish to express their gratitude to the anonymous referee, whose suggestions helped to greatly improve the flow and clarity of our work.

Software: Mercury (Chambers 1999), Astrochron R package.