<figure><img src="image_1.png" /></figure>

HAL
open science

# Passing Stars as an Important Driver of Paleoclimate and the Solar System's Orbital Evolution

Nathan Kaib, Sean N. Raymond

<figure><img src="image_2.png" /></figure>

## To cite this version:

Nathan Kaib, Sean N. Raymond. Passing Stars as an Important Driver of Paleoclimate and the Solar System's Orbital Evolution. The Astrophysical Journal Letters, 2024, 962 (2), pp.L28. 10.3847/2041-8213/ad24fb. hal-04727891

HAL Id: hal-04727891

https://hal.science/hal-04727891v1

Submitted on 10 Oct 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_4.png" /></figure>

# Passing Stars as an Important Driver of Paleoclimate and the Solar System's Orbital Evolution

Nathan A. Kaib¹,² ID and Sean N. Raymond³ ID

¹Planetary Science Institute, 1700 E. Fort Lowell, Suite 106, Tucson, AZ 85719, USA

²Homer L. Dodge Department of Physics & Astronomy, University of Oklahoma, 440 W. Brooks Street, Norman, OK 73019, USA

³Laboratoire d'Astrophysique de Bordeaux, CNRS and Université de Bordeaux, Allée Geoffroy St. Hilaire, 33165 Pessac, France

Received 2023 December 18; revised 2024 January 10; accepted 2024 January 24; published 2024 February 14

## Abstract

Reconstructions of the paleoclimate indicate that ancient climatic fluctuations on Earth are often correlated with variations in its orbital elements. However, the chaos inherent in the solar system's orbital evolution prevents numerical simulations from confidently predicting Earth's past orbital evolution beyond 50–100 Myr. Gravitational interactions among the Sun's planets and asteroids are believed to set this limiting time horizon, but most prior works approximate the solar system as an isolated system and neglect our surrounding Galaxy. Here we present simulations that include the Sun's nearby stellar population, and we find that close-passing field stars alter our entire planetary system's orbital evolution via their gravitational perturbations on the giant planets. This shortens the timespan over which Earth's orbital evolution can be definitively known by a further ~10%. In particular, in simulations that include an exceptionally close passage of the Sun-like star HD 7977 2.8 Myr ago, new sequences of Earth's orbital evolution become possible in epochs before ~50 Myr ago, which includes the Paleocene–Eocene Thermal Maximum. Thus, simulations predicting Earth's past orbital evolution before ~50 Myr ago must consider the additional uncertainty from passing stars, which can open new regimes of past orbital evolution not seen in previous modeling efforts.

*Unified Astronomy Thesaurus concepts: Solar system terrestrial planets (797); Solar system evolution (2293); Solar system planets (1260); Planetary climates (2184); Earth (planet) (439); Stellar dynamics (1596); Dynamical evolution (421); Orbits (1184)*

## 1. Introduction

Over the past 3–4 decades, computational dynamical simulations have enabled direct modeling of the orbital evolution of the solar system over Myr–Gyr timescales (Laskar 1988; Quinn et al. 1991; Wisdom & Holman 1991). This work has demonstrated that the orbital evolution of the Sun's inner planets are chaotic, with Lyapunov timescales of ~5 Myr, and this chaos is thought to be driven by the proximity of multiple secular planetary frequencies as well as asteroid–asteroid encounters (Sussman & Wisdom 1988; Laskar 1989, 1990; Laskar et al. 2011b). Because we do not know the modern solar system's properties to infinite precision, chaotic divergence prevents us from knowing the planets' detailed orbital evolution to arbitrarily long times, and it is only characterized in a statistical sense beyond ~100 Myr timescales, as a broad range of behavior is possible (Laskar 1994). In fact, even the long-term stability of the solar system is not guaranteed over the Sun's lifetime, as simulations indicate a ~1% chance that Mercury will be lost via collision with the Sun or Venus (Laskar 2008).

Consequently, Earth's past or future orbital evolution can only be confidently predicted inside a time horizon much shorter than Earth's age, as small uncertainties in current planetary orbits eventually lead to dramatically divergent behavior. The time horizon set by the internal chaos among the Sun's eight planets is ~70 Myr, but additional strong chaos resulting from encounters between large asteroids shortens the

horizon by another ~10 Myr (Laskar et al. 2011b). However, inside this time horizon, backward integration of the Sun's planets has been used to predict the detailed past orbital evolution of the Earth (Laskar et al. 2004). Because variations in the Earth's paleoclimate are known to correlate with Earth's eccentricity fluctuations (Milankovitch 1941; Hays et al. 1976), this has allowed us to reconstruct the chronology of sedimentation records and assess the orbital context of notable climate events (Zeebe & Lourens 2019).

Simulations of the long-term orbital evolution of the Sun's planets have nearly always modeled the solar system as an isolated system. For many purposes, this is a very good approximation, but the solar system is of course part of the Milky Way Galaxy. Consequently, it occasionally suffers close encounters with other field stars (Öpik 1932; Oort 1950), and solar neighborhood kinematic studies predict an average of ~20 stellar passages within 1 pc of the Sun each Myr (Bailer-Jones 2018). Because the solar system cross section scales with the square of heliocentric distance, the large majority of these encounters will be distant and inconsequential to the planets' dynamics, but this is not guaranteed. In fact, there is a ~0.5% chance that a field star passage will trigger the loss of one or more planets over the next 5 Gyr (Brown & Rein 2022; Raymond et al. 2024), and such passages may actually guarantee the disruption of the planets' orbits many Gyr after the Sun becomes a white dwarf (Zink et al. 2020). Yet, encounters need not trigger an instability for them to have dynamical consequences for the planets. For instance, it has been suggested that ~one-third of Neptune's modern eccentricity has been generated through past stellar encounters (Brunini 1993), but many of the long-term dynamical effects of stellar passages remain unknown. Here we use numerical

<figure><img src="image_5.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further

distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
simulations to characterize how stellar passages resulting from
the Sun's local Galactic environment affect the chaotic
evolution of our planets. Our work and its results are described
in the subsequent sections.

## 2. Methods

We use a modified version of the MERCURY hybrid integrator to perform the simulations presented here (Chambers 1999). This modified version allows for the inclusion of an arbitrary number of stellar mass bodies (Kaib et al. 2018) and also features a simplified central force modification to induce general relativistic orbital precession (Equation (30) of Saha & Tremaine 1992). Unless specified otherwise, the simulations here are integrated for 150 Myr with a step size of 1.5 days. Orbital elements of the planets and Pluto are set to the heliocentric osculating 2000 January 1 elements specified within the JPL Horizons system. The same is done for asteroids Ceres, Vesta, Pallas, Iris, and Bamberga when asteroids are included (which is noted when this is the case). To generate unique sets of initial conditions, the mean anomaly of each body is altered by a random amount that generates an orbital drift between ±2 cm, which is obviously much smaller than the actual uncertainty of the bodies' positions and orbits.

Some of our simulations include populations of passing field stars. In these cases, stellar passages are initiated at random positions 1 pc away from the Sun with an isotropic velocity distribution (Henon 1972). Stellar masses are assigned with the present-day stellar mass function (Reid et al. 2002), and the passage frequency and velocity distributions of each stellar spectral class are specified based on characterizations of solar neighborhood kinematics (García-Sánchez et al. 2001; Rickman et al. 2008). Our stellar passage routine generates 18 stellar encounters within 1 pc per Myr, which compares favorably with recent observational estimates (Bailer-Jones 2018). However, only passages with impact parameters below 0.1 pc are actually integrated, as more distant passages require extraordinarily high masses or low velocities to generate significant perturbations (see the results). Passing stars are integrated until their heliocentric distance again exceeds 1 pc, at which point the body is removed from the simulation.

In simulations that “neglect” stars, we actually still introduce the same sets of stellar passages from our nonisolated runs. The only difference is that the stellar masses are scaled down by a factor of 1000. This is done to try to keep the numerical algorithm sequences between batches of runs as similar as possible, since it is known that the dynamical behavior of a simulated solar system can display variance with the numerical algorithm’s properties (Zeebe 2015).

Earth's Moon and the Sun's oblateness are excluded from our simulations, as well as any tidal or rigid body effects. Our simulations are not designed to generate the most accurate history of Earth's orbital evolution, and far more sophisticated models exist (Laskar et al. 2004, 2011a; Zeebe & Lourens 2019). Instead, they are designed to highlight and characterize the important effects that stellar passages have on the orbital evolution of the planets and compare these effects to those of asteroids as well as the internal chaos that planet-planet interactions drive.

## 2.1. *Secular Mode Analysis*

To decompose a planet's orbital eccentricity evolution into its eigenfrequencies and solve for their phases and amplitudes, we use the Modified Fourier Transform<fnref n="4" /> (Šidlichovský & Nesvorný 1996; Laskar 1999). For the frequencies, amplitudes, and phases measured in Figure 2, orbital elements are recorded every 100 yr for 3.2768 Myr. In all other cases discussed here, we record orbital elements every 100 yr for 6.5536 Myr.

## 3. The Effects of Stellar Passages

Given that even asteroid-asteroid interactions are known to accelerate the chaotic diffusion of Earth's orbit, it may be that perturbations from other stars can also accelerate the divergence of orbital evolution predictions. In Figure 1, we demonstrate that stellar passages do indeed significantly accelerate Earth's orbital evolution divergence. Here we back-integrate 100 realizations (or clones) of the modern solar system for 150 Myr and measure the standard deviation of Earth's orbital eccentricity among these systems over time. (Eccentricity divergences beyond 10% of Earth's mean eccentricity, or 0.0028, have been previously used to mark the time beyond which an orbital solution is unreliable; Zeebe 2017.) If we only consider the eight planets and Pluto, we find that this standard deviation reaches 10% of Earth's mean orbital eccentricity after 77 Myr of back-integration. While our simulations neglect some effects considered in other works such as the Sun's oblateness, tidal dissipation, and gravity from Earth's Moon, this result is roughly consistent with previous results (Laskar et al. 2011b). When we next repeat the integrations, including the asteroids Ceres, Vesta, Pallas, Iris, and Bamberga, we find that this time horizon shortens by 10 Myr down to 67 Myr, and this decrease in divergence time is in approximate agreement with prior results (Laskar et al. 2011b). However, when we then repeat these back-integrations in the presence of a population of passing stars based on the solar neighborhood's properties (García-Sánchez et al. 2001; Rickman et al. 2008), we find that the time horizon again shortens by another 5-7 Myr (62 Myr in the case with stellar passages and the asteroids, and 60 Myr in the case with stellar passages but without asteroids). Because we see a significantly shorter time horizon whether we include asteroids or not, we conclude that perturbations from stellar passages degrade our ability to predict Earth's past or future orbital evolution even more strongly than asteroid perturbations.

The mechanism by which stars alter Earth's orbital evolution is illustrated in Figure 2. In the top panel, we plot the difference in Earth's orbital eccentricity for pairs of solar system integrations with minutely different initial conditions (±2 cm). Among a pair that is integrated in isolation, Earth attains an eccentricity difference of ~10⁻⁷ after 20 Myr of evolution. However, when one of the pair members is subjected to a 25 km s⁻¹ passage of a 1.4 M☉ star within 7200 au, the eccentricity difference jumps up by over 1 order of magnitude to ~3 × 10⁻⁶ after 20 Myr of integration. The rapid increase in eccentricity divergence after the stellar encounter at 10 Myr is obvious. In a final pair integration, we still expose one of the members to the stellar encounter, but this time the giant planets are not included in either simulation. In this last case, the effect of the stellar encounter disappears from the plot, and the pair

<footnote n="4">⁴ http://www.boulder.swri.edu/~davidn/fmft/fmft.html</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 1.** The standard deviation of Earth's orbital eccentricity is plotted against time for 100 back-integrations of solar system clones. Four batches of solar system simulations are run. Two only consider the Sun, its eight planets, and Pluto (blue lines), and the other two (orange lines) also include five large asteroids: Ceres, Vesta, Pallas, Iris, and Bamberga. Half of our batches also perturb the solar system with passing field stars (dashed lines), while the others do not include this perturbation (solid lines).</figcaption></figure>

finishes with an eccentricity difference very near that of the original isolated pair. This strongly suggests that the giant planets play an important role in connecting Earth's orbital evolution to stellar passages.

Even in isolation, the giant planets play a critical role in Earth's orbital evolution (as well as the other inner planets' orbits). Any planet's eccentricity evolution with time can be linearly approximated within Lagrange-Laplace theory as shown in Equations (1) and (2). In this expression, $e_i$ and $\bar{\omega}_i$ are the eccentricity and longitude of perihelion of the ith planet from the Sun, while $g_j$ and $\beta_j$ are the eigenfrequency and phase at which planet j's eccentricity vector precesses, and $M_{ij}$ can be thought of as the magnitude of the eccentricity forcing planet j exerts on planet i. (An analogous decomposition can also be performed for a planet's inclination nodal regression frequencies, $s_{j,i}$.) In Earth's case ($i=3$), the magnitude of Jupiter's forcing $M_{3,5}$ is ~0.019, or more than half of Earth's mean orbital eccentricity. Thus, the instantaneous value of Earth's orbital eccentricity is highly dependent on Jupiter's orbital evolution (which, in turn, can be influenced by the other giant planets),

$$e_i \cos \bar{\omega}_i = \sum_{j=1}^{8} M_{ij} \cos(g_j t + \beta_j) \quad (1)$$

$$e_i \sin \bar{\omega}_i = \sum_{j=1}^{8} M_{ij} \sin(g_j t + \beta_j). \quad (2)$$

Figures 2(B), (C), and (D) illustrate the impact that a stellar passage can have on the eigenfrequencies, phases, and magnitudes of the giant planets' eccentricity forcing. In Figure 2(B) we plot the change in $g_j$ over time between panel (A)'s clone pair, in which one member is subjected to a stellar passage and the other is not. We also do this for the clone pair

integrated in isolation. As can be seen, for the first 10 Myr, the two pairs of systems have identical drifts in eigenfrequency for all four gas giants. However, the divergence in eigenfrequency becomes 2–3 orders greater in the stellar-perturbed pair after the stellar encounter occurs.

Figure 2(C) shows that the stellar passage has a similar effect on the change in $M_{ii}$. Within the context of Lagrange-Laplace theory, $M_{ii}$ is effectively the eccentricity that the ith planet would have in the absence of all other planetary perturbations. Here we see that the stellar encounter drives an $M_{ii}$ divergence among the giant planets that is 2–4 orders of magnitude greater when one system is perturbed by the stellar passage compared to a pair of completely isolated systems. Finally, Figure 2(D) shows that the divergence in the phases of the giant planets' eccentricity forcing is 1–3 orders of magnitudes greater when a stellar encounter perturbs a pair member than when both are isolated. Thus, the frequencies, magnitudes, and phases at which the giant planets perturb the inner planets change much more quickly in the presence of stellar passages than in isolation. In this way, the giant planets serve as a dynamical link that ultimately allows the Milky Way's stars to influence the long-term evolution of Earth's orbit.

## 4. Parameterizing Stellar Passage Importance

The degree to which the giant planets' orbits (and by extension the Earth's orbit) are perturbed of course depends on the strength of the stellar perturbation. One way to measure stellar passage strength is by the velocity impulse it imparts upon the Sun, or $\frac{2GM_*}{v_*b_\odot}$, where $b_\odot$ is the stellar impact parameter with respect to the Sun, $v_*$ is the passage velocity, and $M_*$ is the stellar mass (Rickman 1976). If a planet's dynamical timescale is much longer than the passage timescale (known as
<figure><img src="image_7.png" /><figcaption>**Figure 2.** (A) The absolute difference in Earth's orbital eccentricity between pairs of solar system clone integrations is plotted against time. One pair is integrated in complete isolation (blue line). In the other pairs (orange and black), one pair member includes the 25 km s⁻¹ passage of a 1.4 M☉ star within 7200 au of the Sun at t = 10 Myr, but the other does not include this passage. One pair (black) also excludes the giant planets. (B) For the clone pairs including the giant planets, the absolute pair difference in the proper secular eccentricity precession frequency for each giant planet (*g*<sub>j</sub> in Equations (1) and (2)) is shown vs. time. Planets from the isolated clone pair (smaller, lighter symbols) and the stellar-perturbed pair (larger, darker symbols) are shown. (C) For the clone pairs including the giant planets, the absolute pair difference in the amplitude of the proper secular eccentricity precession frequency for each giant planet is shown vs. time. Planets from the isolated clone pair (smaller, lighter symbols) and the stellar-perturbed pair (larger, darker symbols) are shown. (D) For the clone pairs including the giant planets, the absolute pair difference in the phase of the proper secular eccentricity precession frequency for each giant planet is shown vs. time. Planets from the isolated clone pair (smaller, lighter symbols) and the stellar-perturbed pair (larger, darker symbols) are shown.</figcaption></figure>
the impulsive regime), the planet's impulse relative to the Sun's
is then $\frac{2GM_*}{v_*b_*^2}\Delta x$, where $\Delta x$ is the relative separation of the Sun
and planet projected along the perpendicular to the star's path
at its nearest approach to the Sun.

In fact, most of our stellar passage timescales are far from the impulsive regime. The opposite limit, in which the encounter timescale is much longer than the planetary dynamical timescale, is known as the adiabatic regime. Typically, the planetary dynamical timescale is taken to be the orbital period. However, since we are primarily studying passing stars' perturbations on the planets' secular architecture (see Figure 2), the most relevant planetary timescale is likely the planetary precession periods, which are ~$10^{5-6}$ yr. Thus, with respect to perturbations on the solar system's secular architecture, many stellar passages will be nonadiabatic.

Although most of our stellar passages are far from the formal impulsive regime, we nonetheless find that the perturbation delivered to the solar system scales with $\frac{2GM_*}{v_*b_*^2}$, which is effectively the impulse gradient along the perpendicular to the passage path at nearest approach (Raymond et al. 2024). This is illustrated in Figure 3(A). Here we have generated 100 different 50 Myr sequences of stellar encounters. In each sequence, we select only the encounter with the largest impulse gradient and then subject a solar system clone to this encounter. The perturbation delivered to $M_{8.8}$ (approximately Neptune's eccentricity in the absence of other planetary forcing) is then plotted against the impulse gradient of these 100 different encounters. The two quantities are well-correlated across 4–5 orders of magnitude of perturbation. Perturbations to other giant planet orbital parameters display a level of correlation with stellar impulse gradient that is similar to $M_{8.8}$.

In Figure 3(B), we see that these giant planet perturbations ultimately lead to perturbations on Earth's orbit. Measuring Earth's eccentricity 1 Myr after the stellar passages, we find that the deviation of Earth's eccentricity from the ensemble's median eccentricity scales with the stellar impulse gradient in a manner similar to the giant planets. To once again highlight the importance of the giant planets' role, we repeat these stellar passage experiments without the giant planets. In these repeat runs, we see that the perturbation to Earth's orbit falls by ~2 orders of magnitude. This again confirms that Earth's sensitivity to stellar passages primarily relies on the passages' influence over the more distant giant planets. In the event that the solar system were to have even more distant planets, the dynamical impact of stellar passages could be still stronger than what is described here (Trujillo & Sheppard 2014; Batygin & Brown 2016).

5. HD 7977 and the Paleoclimate

It is clear that the stellar passages expected within the solar neighborhood significantly influence the orbital evolution of the Sun's planets, and we now assess the effects of a specific encounter known to have occurred. Among past stellar encounters inferred from Gaia Data Release 3, HD 7977 stands out as potentially the closest recent known encounter. This 1.1 $M_\odot$ star passed near the solar system ~2.8 Myr ago at ~27 km s⁻¹ (Bailer-Jones 2022; Gaia Collaboration et al. 2023). Although this encounter's median inferred impact parameter is ~13,200 au, there is a large amount of uncertainty, with a 5% probability of passage within ~3900 au. This range of impact parameter corresponds to over 1 order of magnitude

variation in impulse gradient (which governs the level of planetary perturbation).

In Figure 4(A), we back-integrate 100 solar system realizations analogous to the isolated ensemble including asteroids from Figure 1. However, this time we also subject the ensemble to an HD 7977-like passage at $t = -2.8$ Myr. This is done twice. In the first instance, we hold the impact parameter at the median value of 13,200 au, while our individual clones randomly vary the encounter velocity (±0.4 km s⁻¹) and time of passage (±0.004 Myr) within the observational uncertainty range (Bailer-Jones 2022). The second back-integration ensemble employs the same variations in stellar velocity and passage time, but the impact parameter is decreased to 3900 au. As seen in Figure 4(A), the median impact parameter for HD 7977 has no discernible effect on the divergence of Earth's eccentricity during back-integrations. Given that the real encounter's impact parameter was probably near this value or larger, HD 7977 has most likely not had a significant effect on the solar system's long-term chaotic evolution. However, this is not for certain. We find that the stellar passage's effect becomes very significant if we assume the 3900 au impact parameter, which has a 5% likelihood. In this case, the standard deviation of Earth's eccentricity exceeds 0.0028 ~ 10 Myr faster ($t = -58$ Myr) than our models that just consider planets and asteroids.

In this last back-integration ensemble that considers the smaller impact parameter with 5% likelihood, even faster divergence is also possible. If we instead plot the difference between the 98th and 2nd percentile of Earth's eccentricity (roughly the $2\sigma$ confidence bounds), we see that the difference exceeds 0.0028 in just 50 Myr, nearly 20 Myr faster than the nominal chaos time horizon set by asteroid interactions in our models. We should note that employing a 3900 au impact parameter for HD 7977 would likely make it one of the 10 most powerful encounters experienced over the solar system's history as measured by impulse gradient (owing to the star's above-average mass and below-average encounter velocity). Nevertheless, such an impact parameter is within the bounds of observational uncertainty, and better constraining of this encounter is essential to understanding the recent dynamical history of Earth and the rest of the Sun's planets.

Incredibly precise dynamical models of Earth's past orbital evolution have been created over the last three decades to assist in interpreting the past 50–70 Myr of Earth's sedimentation record and paleoclimate (Laskar et al. 2004, 2011a; Zeebe & Lourens 2019). However, beyond 50 Myr, these model outcomes are known to be sensitive to slight changes in solar system properties such as the Sun's oblateness and the number of asteroids considered (Zeebe 2017), and models that maximize agreement between Earth's eccentricity evolution and sedimentation patterns 50–60 Myr ago must invoke modern solar system parameters that are somewhat in tension with observational constraints (Zeebe & Lourens 2019). These maximal agreement models also contain a specific resonance transition 50 Myr ago (Laskar et al. 2004; Ma et al. 2017; Zeebe & Lourens 2019). Currently, the eigenfrequency difference between Earth's and Mars' nodal regression ($s_4 - s_3$) is twice that of their perihelion precession difference ($g_4 - g_3$), and entry into this resonance occurred when the perihelion precession converged to a ratio of 2 from some other value at some point in the past.
<figure><img src="image_8.png" /><figcaption>**Figure 3.** (A) The absolute difference between the observed amplitude of Neptune's proper secular eccentricity precession frequency and that measured in 100 integrations of solar system clones is plotted against the impulse gradient of the stellar encounter to which each clone is exposed. (B) The absolute difference between the Earth's eccentricity in 100 solar system clone integrations and the median eccentricity is plotted against the impulse gradient of the stellar encounter to which each clone is exposed. Filled blue data points are for clones that contain all eight planets, while open red data points do not contain the giant planets. Earth eccentricities are measured ~1 Myr after the stellar encounter in each integration.</figcaption></figure>

In Figure 5(A), we calculate the ratio ($s_4 - s_3$) to ($g_4 - g_3$) among our solar system realizations from Figure 4 at $t = -50$ Myr. We see that none of our isolated systems are far from the 2:1 resonance at this past time. Meanwhile, the systems subjected to HD 7977's 3900 au flyby are also concentrated near 2:1 but display larger variance. When we

reexamine these ratios at the earlier epoch of $t = -60$ Myr, we find that our isolated systems are systematically nearer to 2:1 than the HD 7977-perturbed systems, many of which are clearly well outside the 2:1 resonance and have yet to enter it. Thus, the close encounter with HD 7977 at $t = -2.8$ Myr can unlock a much broader possible spectrum of orbital behavior
<figure><img src="image_9.png" /><figcaption>**Figure 4.** (A) The standard deviation of Earth's orbital eccentricity is plotted against time for 100 back-integrations of solar system clones. Three batches of solar system simulations include the Sun, its eight planets, Pluto, Ceres, Vesta, Pallas, Iris, and Bamberga. One batch includes no stellar perturbations (orange). Another also includes a passage of HD 7977 at *t* = −2.8 Myr at its median impact parameter of 13,200 au, as indicated by Gaia uncertainties (solid blue; Bailer-Jones 2022). The final batch includes a passage of HD 7977 at *t* = −2.8 Myr at 3900 au, the 5th percentile of possible impact parameters according to Gaia uncertainties (dashed blue; Bailer-Jones 2022). (B) The difference between the 98th and 2nd percentile of Earth's orbital eccentricity vs. time among 100 back-integrations of solar system clones. The three batches of clones are the same as in the top panel.</figcaption></figure>

50–60 Myr ago with respect to entry into this secular
resonance. The consequences of this can be seen in
Figure 5(B), where we examine the maximum eccentricity
Earth attains 55–56 Myr ago. This highlighted period
encompasses the Paleocene–Eocene Thermal Maximum,
which has been suggested to coincide with a period of

heightened orbital eccentricity for Earth (Zeebe & Lourens 2019). Among our 100 isolated solar system realizations, Earth's maximum eccentricity during this time varies between 0.0546 and 0.0555 (±0.00045). Among our 100 HD 7977-perturbed systems, the range in maximum eccentricity is nearly 1 order of magnitude larger at 0.0500–0.0569
<figure><img src="image_10.png" /><figcaption>**Figure 5.** (A) The distribution of $(s_4 - s_3)/(g_4 - g_3)$ is calculated for our isolated back-integrated systems (blue) and our back-integrated systems perturbed by a 3900 au flyby of HD 7977 at $t = -2.8$ Myr (orange). These ratios are calculated at $t = -50$ Myr (solid) and $t = -60$ Myr (dashed). (B) The distribution of maximum eccentricity values Earth attains between −55 Myr and −56 Myr is plotted for our isolated solar system realizations (blue) and our systems perturbed by a 3900 au flyby of HD 7977 at $t = -2.8$ Myr (orange).</figcaption></figure>

($\pm0.0035$). Thus, our interpretation of Earth's orbit near the Paleocene–Eocene Thermal Maximum is dependent on the exact parameters of stellar encounters that have occurred in the past several Myr. Due to the broader range of Earth's possible orbital behavior, orbital solutions may exist that provide good matches to the geological record before 50 Myr ago but also do not necessarily have a period of notably

heightened orbital eccentricity during the Paleocene–Eocene Thermal Maximum. Conversely, the divergent behaviors illustrated in Figure 5(A) raise the possibility that geological constraints on the timing of this secular resonance entry could imply additional constraints on the magnitude of recent stellar perturbations (Ma et al. 2017; Westerhold et al. 2017).
# 6. Summary and Conclusions

Much more sophisticated simulations of Earth's past orbital evolution exist compared to those presented here (Laskar et al. 2004, 2011a; Zeebe & Lourens 2019), and our present work is not meant to make the most accurate predictions of Earth's exact orbital history. Instead, it compares the importance of stellar encounters relative to the much more intensely studied internal dynamics of the solar system. We show that stellar encounters play an important role in our solar system's long-term dynamical evolution, and our results highlight a number of notable conclusions. First, stellar encounters significantly accelerate the chaotic diffusion of Earth's orbit and the time back to which numerical simulations can confidently predict Earth's orbital evolution is ~10% shorter than previously thought. Second, this chaotic divergence that stellar passages impart on Earth's orbit results from their perturbations to the giant planets' orbits, and these perturbations roughly scale with the velocity impulse gradients of stellar encounters. Third, the known encounter with HD 7977 2.8 Myr ago has the potential to unlock new sequences of Earth's past orbital evolution beyond 50 Myr ago that have not been considered or generated in previous modeling efforts. Although it takes tens of Myr for the effects of stellar passages to significantly manifest themselves, the long-term orbital evolution of the Earth and the rest of the planets is linked to these stars.

# Acknowledgments

This work was performed with support from NSF CAREER Award 2405121. S.N.R. acknowledges support from the French Programme Nationale de Planetologie (PNP) and in the framework of the Investments for the Future program IdEx, Université de Bordeaux/RRI ORIGINS. Computing for this project was performed at the OU Supercomputing Center for Education & Research (OSCER) at the University of Oklahoma (OU).

# ORCID iDs

Nathan A. Kaib https://orcid.org/0000-0001-5272-5888

Sean N. Raymond https://orcid.org/0000-0001-8974-0758

# References

Bailer-Jones, C. A. L. 2018, A&A, 609, A8

Bailer-Jones, C. A. L. 2022, ApJL, 935, L9

Batygin, K., & Brown, M. E. 2016, AJ, 151, 22

Brown, G., & Rein, H. 2022, MNRAS, 515, 5942

Brunini, A. 1993, A&A, 276, 261

Chambers, J. E. 1999, MNRAS, 304, 793

Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023, A&A, 674, A1

García-Sánchez, J., Weissman, P. R., Preston, R. A., et al. 2001, A&A, 379, 634

Hays, J. D., Imbrie, J., & Shackleton, N. J. 1976, Sci, 194, 1121

Henon, M. 1972, A&A, 19, 488

Kaib, N. A., White, E. B., & Izidoro, A. 2018, MNRAS, 473, 470

Laskar, J. 1988, A&A, 198, 341

Laskar, J. 1989, Natur, 338, 237

Laskar, J. 1990, Icar, 88, 266

Laskar, J. 1994, A&A, 287, L9

Laskar, J. 1999, RSPTA, 357, 1735

Laskar, J. 2008, Icar, 196, 1

Laskar, J., Fienga, A., Gastineau, M., & Manche, H. 2011a, A&A, 532, A89

Laskar, J., Gastineau, M., Delisle, J. B., Farrés, A., & Fienga, A. 2011b, A&A, 532, L4

Laskar, J., Robutel, P., Joutel, F., et al. 2004, A&A, 428, 261

Ma, C., Meyers, S. R., & Sageman, B. B. 2017, Natur, 542, 468

Milankovitch, M. M. 1941, Kanon der Erdbestrahlung und seine Anwendung auf das Eiszeitenproblem, Vol. 132 (Belgrade: Royal Serbian Academy), 32

Oort, J. H. 1950, BAN, 11, 91

Öpik, E. 1932, PAAAS, 67, 169

Quinn, T. R., Tremaine, S., & Duncan, M. 1991, AJ, 101, 2287

Raymond, S. N., Kaib, N. A., Selsis, F., & Bouy, H. 2024, MNRAS, 527, 6126

Reid, I. N., Gizis, J. E., & Hawley, S. L. 2002, AJ, 124, 2721

Rickman, H. 1976, BAICz, 27, 92

Rickman, H., Fouchard, M., Froeschlé, C., & Valsecchi, G. B. 2008, CeMDA, 102, 111

Saha, P., & Tremaine, S. 1992, AJ, 104, 1633

Sussman, G. J., & Wisdom, J. 1988, Sci, 241, 433

Trujillo, C. A., & Sheppard, S. S. 2014, Natur, 507, 471

Šidlichovský, M., & Nesvorný, D. 1996, CeMDA, 65, 137

Westerhold, T., Röhl, U., Frederichs, T., et al. 2017, CliPa, 13, 1129

Wisdom, J., & Holman, M. 1991, AJ, 102, 1528

Zeebe, R. E. 2015, ApJ, 798, 8

Zeebe, R. E. 2017, AJ, 154, 193

Zeebe, R. E., & Lourens, L. J. 2019, Sci, 365, 926

Zink, J. K., Batygin, K., & Adams, F. C. 2020, AJ, 160, 232