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

HAL
open science

# Future trajectories of the Solar System: dynamical simulations of stellar encounters within 100 au

Sean Raymond, Nathan Kaib, Franck Selsis, Herve Bouy

► To cite this version:

Sean Raymond, Nathan Kaib, Franck Selsis, Herve Bouy. Future trajectories of the Solar System: dynamical simulations of stellar encounters within 100 au. Monthly Notices of the Royal Astronomical Society, 2024, 527 (3), pp.6126-6138. 10.1093/mnras/stad3604. hal-04727886

HAL Id: hal-04727886

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

Submitted on 10 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_2.png" /></figure>

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

# Future trajectories of the Solar System: dynamical simulations of stellar encounters within 100 au

Sean N. Raymond<sup>1,1★</sup> Nathan A. Kaib,<sup>2</sup> Franck Selsis<sup>1</sup> and Herve Bouy<sup>1</sup>

Accepted 2023 November 13. Received 2023 November 7; in original form 2023 September 8

## ABSTRACT

Given the inexorable increase in the Sun's luminosity, Earth will exit the habitable zone in ~1 Gyr. There is a negligible chance that Earth's orbit will change during that time through internal Solar System dynamics. However, there is a ~1 per cent chance per Gyr that a star will pass within 100 au of the Sun. Here, we use N-body simulations to evaluate the possible evolutionary pathways of the planets under the perturbation from a close stellar passage. We find a ~92 per cent chance that all eight planets will survive on orbits similar to their current ones if a star passes within 100 au of the Sun. Yet a passing star may disrupt the Solar System, by directly perturbing the planets' orbits or by triggering a dynamical instability. Mercury is the most fragile, with a destruction rate (usually via collision with the Sun) higher than that of the four giant planets combined. The most probable destructive pathways for Earth are to undergo a giant impact (with the Moon or Venus) or to collide with the Sun. Each planet may find itself on a very different orbit than its present-day one, in some cases with high eccentricities or inclinations. There is a small chance that Earth could end up on a more distant (colder) orbit, through re-shuffling of the system's orbital architecture, ejection into interstellar space (or into the Oort cloud), or capture by the passing star. We quantify plausible outcomes for the post-flyby Solar System.

**Key words:** astrobiology – methods: numerical – Moon – Oort cloud – planets and satellites: dynamical evolution and stability.

## 1 INTRODUCTION

Earth has about a billion years of habitable surface conditions remaining. Our planet is perched perilously close to the inner edge of the habitable zone (Kasting, Whitmire & Reynolds 1993). If our planet were just 5 per cent–10 per cent closer to the Sun, there would no longer be a stable balance allowing for liquid water to remain on Earth's surface (Selsis et al. 2007; Kopparapu et al. 2013; Leconte et al. 2013). At the same time, the Sun's luminosity is slowly increasing, pushing the inner edge of the habitable zone inexorably outward, toward Earth's orbit (e.g. Wolf & Toon 2014).

Yet the Solar System is not a closed system. While the orbital evolution of the planets is largely determined by secular and resonant perturbations (Laskar 2012), passing stars can have a consequential influence on the planets' orbits. Based on the local galactic density of stars and their velocity dispersion, most stellar flybys are distant and only very weakly perturb the planets' orbits (e.g. Brown & Rein 2022). Statistically speaking, flybys closer than 100 au, which would strongly affect the planets' orbits, only take place roughly once per 100 Gyr in the current Galactic neighbourhood (Zink, Batygin & Adams 2020; Brown & Rein 2022). This amounts to a ~1 per cent probability of a sub-100 au encounter in the next Gyr, which is comparable to the well-studied probability of a chaos-driven dynamical instability among the terrestrial planets before the

Sun becomes a red giant in ~5 Gyr (Laskar & Gastineau 2009; Zeebe 2015; Abbot et al. 2021).

Here, we use a suite of N-body simulations to evaluate the influence of a close stellar encounter on the orbital evolution of the Solar System. We explore the diversity of outcomes and calculate their relative probabilities. We pay particular attention to the future habitability of Earth, and show that there exist pathways by which Earth will end up on a cooler orbit (or a hotter one). The most exotic outcomes involve capture of Earth or other planets by the passing star, or the scattering of planets into the Sun's Oort cloud.

Our study was inspired in part by the pioneering work of Laughlin & Adams (2000), who simulated the dynamical response of the Solar System to flybys of binary stars. We build on their results in a number of ways, for example, by including the Moon as a separate particle in a subset of simulations (Section 4), which allows for additional outcomes (such as the Moon crashing onto the Earth).

## 2 SIMULATIONS

### 2.1 Initial conditions

Our N-body simulations started with the eight planets orbiting the Sun, and a single star set to fly past. The planets were placed at their positions as of Julian Date 2451000.5.<fnref n="1" /> The exact position

<footnote n="1"><sup>1</sup>Laboratoire d'Astrophysique de Bordeaux, CNRS and Université de Bordeaux, Allée Geoffroy St. Hilaire, F-33165 Pessac, France</footnote>

<footnote n="2"><sup>2</sup>Planetary Science Institute, 1700 E. Fort Lowell, Suite 106, Tucson, AZ 85719, USA</footnote>

<footnote n="1">¹These are the planets' positions given as an example in the MERCURY code input files (Chambers 1999).</footnote>

<footnote>*E-mail: rayray.sean@gmail.com</footnote>
at time zero has no influence on our results, because the timing and geometry of the closest stellar approach to the Sun varies from simulation to simulation between about 10 and 30 kyr due to differences in the stellar velocity. The Earth and Moon were included as a single particle orbiting at their centre of mass with the sum of their masses.

Stellar encounters were generated to match the local Galactic neighbourhood conditions. The flyby star's mass was chosen to follow the observed mass function of nearby stars determined by Reid, Gizis & Hawley (2002). The lower stellar mass limit was set to 0.05 M☉, and any stars more massive than 4 M☉ were set to 4 M☉ to account for their short lifetimes. The orbital velocity of each star was chosen in a mass-dependent fashion to match the present-day velocity distribution in the Solar neighbourhood (following Rickman et al. 2008). The resulting velocity distribution has a median of 40.2 km s⁻¹ and a standard deviation of 16.9 km s⁻¹ (averaged over all stellar types). We only include single stars. At encounter speeds faster than ~20 km s⁻¹, which are typical given the velocity dispersion in the Galactic field, the Solar System's interaction with a binary can be approximated as separate interactions with each stellar component (Li & Adams 2015). This contrasts with the approach of Laughlin & Adams (2000), who only modelled binary encounters.

The final parameter for our stellar flybys was the impact parameter *b*. The impact parameter is essentially the closest approach of the flyby star before gravity is taken into account (bringing the closest approach closer). The probability of a close encounter scales as *b*<fnref n="2" />, but we chose instead to sample following a log(*b*) distribution. This choice allowed us to explore low probability, but interesting outcomes that arise from very close flybys. We took this sampling into account when calculating probabilities (see Section 2.3).

The flyby star's orbit was placed at the edge of our numerical domain, just interior to the 100 000 au ejection radius. While this is slightly smaller than the Sun's tidal radius (which depends on the local Galactic density; Tremaine 1993), it still allows for capture of planets in the Oort cloud.

## 2.2 Integration code

Our code is based on the MERCURY integration package (Chambers 1999), to which we added two effects: general relativistic precession (using the implementation of Saha & Tremaine 1992)<fnref n="2" /> and the Galactic tide (following Levison, Dones & Duncan 2001). In our calculation of the Galactic tide, we assumed a local galactic density of 0.1 M☉ pc⁻<fnref n="3" /> (Chakrabarti et al. 2021), and assigned values for Oort's constants of A = 14 km s⁻¹ kpc⁻¹ and B = 12 km s⁻¹ kpc⁻¹ (Feast & Whitelock 1997).<fnref n="3" />

In each of our 12 000 simulations, the orbits of the planets and flyby star were integrated with a Bulirsch-Stoer integrator (Press et al. 1992), using an accuracy parameter of 10⁻¹⁵. Collisions were treated as inelastic mergers conserving linear momentum. Each simulation was integrated for 20 Myr.

As a control, we also integrated the Solar System with no stellar flyby for 2 Gyr with the same code and saw no hints of instability.

## 2.3 Probability determination

As discussed above, our 12 000 simulations sampled the impact parameters *b* of stellar flybys using a logarithmic distribution instead of a physically motivated, *b*<fnref n="2" /> distribution. We transform the probability of a simulation in our chosen sampling as follows.

Within our simulations, $dP_{sim}/d(\log(b))$ is constant, such that $dP_{sim}/db \sim b^{-1}$. In the real Universe, $dP_{Universe}/db \sim b$. To apply this to our simulations requires a linear distribution in *b*, so $dP_{Universe}/db \sim b_{sim}/(dP_{sim}/db)$, such that $dP_{Universe}/db \sim b_{sim}^2$.

To evaluate the relative probability of our simulations, we therefore weight each one by its impact parameter squared. We verified that this converges to the expected probability distribution for $b^2$ sampling. We normalize such that the sum of all of the probabilities is equal to one, keeping in mind that the total probability of a single encounter within 100 au is ~1 per cent Gyr⁻¹ given current local Galactic conditions (Zink et al. 2020; Brown & Rein 2022).

# 3 SIMULATION OUTCOMES

## 3.1 Influence of flyby parameters

The key parameter in determining the outcome of a stellar flyby is the 'impulse gradient', which measures the radial gradient in acceleration felt by planets orbiting the Sun. An impact parameter of 100 au and an encounter speed of 40 km s⁻¹ yields an encounter time-scale of 12 yr, or the orbital period of Jupiter. This implies that for such encounters, perturbations to the outer three giant planets can be approximated well within the impulsive regime. Although this is not necessarily true for Jupiter and the inner planets, we find that parametrizing our encounter results in the impulsive regime yields well-behaved empirical scalings. In particular, we find that encounter outcomes correlate well with the gradient of the stellar impulse measured at the Sun in the direction of the perpendicular to the stellar flyby path. The magnitude of this gradient is defined as:

$$ \text{Impulse gradient} = \frac{2GM_{\star}}{v_{\star}b^2}, \qquad (1) $$

where $M_{\star}$ is the mass of the flyby star, $v_{\star}$ is its speed relative to the Sun (at infinity), *b* is the impact parameter of the flyby and *G* is the gravitational constant. The units of the impulse gradient are generally s⁻¹, but we express them in km s⁻¹ au⁻¹, to give a more intuitive sense of the perturbations felt by the planets, as their orbital speeds are measured in tens of km s⁻¹.

The normalized angular momentum deficit is a measure of the orbital excitation level of a planetary system, in terms of the non-circularity and non-coplanarity of the planets' orbits (Laskar 1997). It is defined as:

$$ \mathrm{AMD} = \frac{\sum_j m_j \sqrt{a_j} (1 - \cos(i_j) \sqrt{1 - e_j^2})}{\sum_j m_j \sqrt{a_j}}, \qquad (2) $$

where $a_j$, $e_j$, $i_j$, and $m_j$ refer to planet *j*'s semimajor axis, eccentricity, inclination with respect to a fiducial plane, and mass. The AMD of the Solar System is 0.00128 (or 0.0015 when measured relative to the Earth's orbital plane).

Fig. 1 (left panel) shows that the angular momentum deficit of surviving systems is a strong function of the impulse gradient in our simulations. For weak stellar encounters (with impulse gradients below ~10⁻<fnref n="2" /> km s⁻¹ au⁻¹), the Solar System is barely perturbed by the passing star, and the planets' orbits almost always maintain their same low-eccentricity, low-inclination orbits. A planet is lost in only 3.4 per cent of simulations with IG < 10⁻<fnref n="2" /> km s⁻¹ au⁻¹ (although,

<footnote n="2">²See Brown & Rein (2023) for an evaluation of the effect on Solar System stability of different implementations of general relativity.</footnote>

<footnote n="3">³See Raymond et al. (2023) for a study of the effect of varying the Galactic density on the trapping of scattered planets in the Oort cloud.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Distribution of the outcomes of our simulations as a function of the impulse gradient of the stellar flyby. The symbol colour corresponds to the number of surviving planets in a given system, and the size is proportional to the mass of the flyby star. The angular momentum deficit AMD (Laskar 1997; Laskar & Petit 2017) is a strong function of the impulse gradient. However, the relative energy conservation of individual simulations (as measured by the fractional energy error budget $dE/E$) is independent of the impulse gradient, lending credence to our numerical methods.</figcaption></figure>

when weighted by the relative probability of each encounter, only 1.4 per cent of encounters with $IG < 10^{-2}$ km s$^{-1}$ au$^{-1}$ lose a planet. These encounters are not universally at large distances; while the median impact parameter is 47.5 au among flybys with $IG < 10^{-2}$ km s$^{-1}$ au$^{-1}$, a handful have $b < 10$ au and one quarter of stars crossed the planets' orbits, with $b < 30$ au. For stronger encounters, with $IG > 10^{-2}$ km s$^{-1}$ au$^{-1}$, the strength of the instability still scales with the impulse gradient. The level of devastation can be measured in terms of both the AMD and simply by the number of surviving planets. In the most extreme cases, one (or no) planets survive. We will explore the evolution of such cases below.

The degree of orbital energy conservation $dE/E$ is virtually independent of the impulse gradient (Fig. 1, right panel). This is reassuring, as it would be concerning if the strength of the instability governed the level of accuracy of our simulations. In addition, the overall level of energy conservation is excellent, with a median of $1.9 \times 10^{-7}$, which is orders of magnitude lower than the acceptable threshold used in simulations of giant planet scattering (e.g. Chatterjee et al. 2008; Raymond, Armitage & Gorelick 2010). This demonstrates that our integration method is reliable.

## 3.2 Number of surviving planets

If a star passes within 100 au of the Sun, there is still a very high chance that all eight Solar System planets will survive. Fig. 2 shows a probability of 95.6 per cent that no planets will be lost, at least during the 20 Myr following the stellar flyby. There is a 92 per cent probability that all eight planets will survive with each pair of neighbours on Hill-stable orbits with separations larger than 3.5 mutual Hill radii (Marchal & Bozis 1982; Gladman 1993), and the system will have an angular momentum deficit less than twice the present-day value. There is a monotonically decreasing probability of fewer planets surviving, with only a $3.7 \times 10^{-6}$ probability of zero remaining planets.

The angular momentum deficit AMD of post-flyby systems scales strongly with the number of surviving planets $N_p$ (Fig. 2). For systems with eight surviving planets, the median AMD is only 10 per cent higher than the present-day Solar System's. The orbits of systems with fewer surviving planets are progressively more and more excited as $N_p$ decreases, with higher and higher eccentricities. This same

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Probability of $N_p$ planets surviving a stellar flyby within 100 au (left axis; red). The right (blue) axis shows the angular momentum deficit distribution as a function of $N_p$. The central tick shows the median of each distribution, and the other ticks show the maxima and minima of the AMD distribution for each value of $N_p$. The present-day Solar System's AMD is ~0.00157, or -2.8 in this plot.</figcaption></figure>

trend is seen among systems of known exoplanets (Limbach & Turner 2015; Turrini, Zinzi & Belinchon 2020).

## 3.3 Examples

Fig. 3 shows the evolution of nine example simulations that span the range of outcomes. The top row shows three cases in which all of the planets survived, but with different levels of orbital excitation. The simulation at the top left underwent a weak encounter and had a final AMD that was indistinguishable from its starting value. The top centre simulation underwent a stronger instability and ended up with an AMD about three times higher than the pre-encounter Solar System. However, the planets' orbits remained stable with no orbital crossing. In the top right simulation, a still stronger instability perturbed the system to a very high AMD, with clear signs of impending instability, as Uranus and Neptune's orbits are crossing and Mercury's eccentricity has drastically increased.

The centre row of Fig. 3 shows cases in which one or more planets were lost from the Solar System. In the centre left simulation, perturbations to the giant planets were transmitted to the terrestrial
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Nine example simulations illustrating the diversity of outcomes. Each panel shows the orbital radius of each planet and the flyby star as a function of time; note that this is not the semimajor axis, so the spread in a planet's orbital radius over time is a measure of its orbital eccentricity. The top row includes three simulations in which all eight planets survived; the top left simulation is basically unperturbed, while the top right simulation is in the early stages of dynamical instability. The middle row includes systems that lost 1–2 planets. The bottom row presents simulations that were catastrophically disrupted, with 0–2 surviving planets.</figcaption></figure>

planets, leading Mercury to eventually collide with the Sun after ~10 Myr. In the centre middle simulation, a stronger stellar encounter significantly altered the ice giants' orbits and ended up destabilizing the terrestrial planets. Mercury once again collided with the Sun, and Venus and Earth underwent a giant impact. In the centre right panel, the stellar perturbations were largely confined to the outer Solar System: Uranus and Saturn were both ejected from the Solar System, but the remaining planets' orbits were stable for the rest of the simulation.

The bottom row of Fig. 3 shows a few of the most destructive outcomes. In the bottom left simulation, the stellar encounter rapidly drove Uranus into the Sun, ejected Neptune, and triggered a dynamical instability between Saturn and Jupiter. Saturn was eventually ejected by Jupiter, and the gas giant perturbations spread to the terrestrial planets, leading to Mercury, Venus, and Mars colliding with the Sun. The final Solar System contained just Jupiter and Earth, on eccentric orbits ($e_{Earth} \approx 0.4$ and $e_{Jup} \approx 0.2$). The bottom centre simulation shows a case in which only Jupiter survived. The very strong stellar encounter (with $b = 2.6$ au) immediately triggered a strong dynamical instability across the entire Solar System. Mercury, Earth, and Neptune collided with the Sun. Venus, Mars, Saturn, and Uranus were ejected after close encounters with Jupiter. Jupiter was the lone survivor, winding up with an orbital eccentricity of 0.55.

Finally, in the bottom right simulation from Fig. 3, an extremely strong close encounter with an impact parameter of 1.5 au and an impulse gradient of $2 \times 10^{-7}$ caused all eight planets to be ejected, leaving the Sun planet-less and alone.

## 3.4 Post-flyby orbits

Fig. 4 shows the final (post-flyby) orbital distribution of the planets in all of our simulations. Each planet has a distribution of final orbits that extends to indefinitely high orbital eccentricities (up to $e \sim 1$). Their semimajor axis–eccentricity distribution is reminiscent of those of the survivors of giant planet instabilities such as those thought to be responsible for the large eccentricities observed among giant exoplanets (Chatterjee et al. 2008; Raymond et al. 2010). A key difference is that, in our simulations, only the giant planets have high enough escape speeds to eject other planets (or, strictly speaking, large enough Safronov numbers $\Theta$, where $\Theta$ is the ratio of the escape speed from a planet’s surface to the escape speed from the system at the planet’s orbital radius; see Ford & Rasio 2008). Instabilities among the terrestrial planets usually lead to collisions, not ejections. Yet perturbations from the giant planets can clearly leave the terrestrial planets on orbits that are excited well beyond their self-excitation limit.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Final orbits of all surviving planets in our simulations. The tail of planets with semimajor axes of $10^4 - 5$ au are those trapped in the Oort cloud.</figcaption></figure>

Fig. 5 shows the probability-weighted eccentricity and inclination distributions for all eight planets (when they survived). The distribution is strongly weighted toward simulations in which the Solar System remained stable and with a low AMD, with the planets on orbits similar to their present-day ones, because that is the most probable outcome of a sub-100 au flyby (see Fig. 2). However, there is a tail for each planet that extends to high eccentricities and inclinations. These correspond to dynamical instabilities triggered by the flyby that are violent enough to strongly excite the planets' orbits. Of course, when instabilities are so strong that planets are ejected from the Sun, they no longer contribute to the distributions in Fig. 5. While low in overall probability, outcomes in which the planets survive on high eccentricity or highly inclined orbits are intriguing in terms of their influence on planetary climate. We will return to this issue in Section 5, with regards to the Earth.

A striking population of surviving planets in Fig. 4 are those with semimajor axes of $10^4 - 5$ au, but small enough eccentricities that they are not in dynamical contact with the other planets. These are Oort cloud planets (see Raymond, Izidoro & Kaib 2023, which is entirely dedicated to Oort cloud planets). They were scattered outward during gravitational instabilities, and spent long enough at very high orbital radius to have had their perihelion distances increased via torques from the Galactic tide, until they were no longer being scattered by the inner planets. As such, their origins are similar to those of Oort cloud comets (e.g. Brasser & Morbidelli 2013; Dones et al. 2015).

We performed additional simulations to capture the implantation of Oort cloud planets. Starting from the final orbits after 20 Myr, we selected candidate systems in which a planet's semimajor axis was larger than 1000 au and its perihelion distance was larger than 50 au. We integrated these systems out to 200 Myr. In many cases, the planet was simply ejected, but in 28 simulations a planet was captured in the Oort cloud. In most cases, it was Saturn or the ice giants that were captured in the Oort cloud, but there was one in which Mars ended up in the Oort cloud, and there were two simulations in which Earth was captured as an Oort cloud planet. The total probability of a planet being captured in the Oort cloud in our simulations was $8.9 \times 10^{-5}$, which, after taking into account the encounter likelihood, amounts to a net probability of $8.9 \times 10^{-7}$ Gyr$^{-1}$. This is likely an underestimate; at the end of our 20 Myr simulations there were an additional 182 planets with semimajor axes larger than 1000 au and perihelion distances smaller than 50 au. Assuming an Oort cloud

trapping rate of ~5 per cent based on previous results (Raymond et al. 2023), this would increase the total probability of trapping a planet in the Oort cloud by roughly 30 per cent, to $1.2 \times 10^{-6}$ Gyr$^{-1}$.

Fig. 6 shows the evolution of a simulation in which Earth was trapped in the Oort cloud. After being scattered by the giant planets to large orbital radius the Galactic tide increased the Earth's perihelion distance on a ~100 Myr time-scale. Earth finished the simulation on a stable orbit in the Oort cloud with an orbital semimajor axis of 54977 au. The long-term survival of Earth in the Oort cloud is not guaranteed. Torques from the Galactic tide continue to affect the Oort cloud planet's orbit, eventually decreasing its perihelion distance and bringing it back into dynamical contact with the inner planets (Neptune and Jupiter in this case). Earth could potentially be ejected from this simulation on a time-scale of a few hundred million years (see Raymond et al. 2023). Also, additional stellar encounters are likely to strip a 50 000 au-wide orbit after a Gyr or so (see fig 18 of Kaib, Roškar & Quinn 2011).

## 3.5 Planetary destruction pathways

We determined the most common pathways through which planets may be lost, keeping in mind that there is a ≳ 95 per cent chance that no planet will be lost if a star passes within 100 au (see Fig. 2). Table 1 shows these probabilities for each planet.

The 10 most likely destruction planetary pathways are, in order of decreasing probability:

(i) Mercury collides with the Sun (probability of 2.54 per cent).

(ii) Mars collides with the Sun (1.21 per cent).

(iii) Venus impacts another planet (1.17 per cent).

(iv) Uranus is ejected (1.06 per cent).

(v) Neptune is ejected (0.81 per cent).

(vi) Mercury impacts another planet (0.80 per cent).

(vii) Earth impacts another planet (0.48 per cent).

(viii) Saturn is ejected (0.32 per cent).

(ix) Mars impacts another planet (0.27 per cent).

(x) Earth collides with the Sun (0.24 per cent).

Examples of many of these outcomes are apparent in Fig. 3. These outcomes are not mutually exclusive; in many cases multiple planets are lost in the same simulation.

Mercury's destruction rate is higher than that of all the giant planets put together. It is the most fragile planet in the Solar System from a dynamical point of view, even through perturbations from passing stars. It is already well established that, in the absence of a stellar flyby, Mercury is the most likely planet to be lost from the Solar System, by having its eccentricity dramatically increased via a secular resonance with Jupiter, causing it to collide with the Sun (Laskar & Gastineau 2009; Batygin, Morbidelli & Holman 2015; Zeebe 2015; Abbot et al. 2021). Mercury's eccentricity increase can also result in collisions between the rocky planets, which are also common outcomes in simulations with stellar flybys. In such simulations, Mercury's destruction does not happen until long after the stellar flyby. We interpret this as evidence that its eccentricity increase is induced by indirect perturbations, through secular forcing by the giant planets — likely via the same secular resonance with Jupiter that is capable of destabilizing its orbit in the absence of flybys.

It is not surprising that the ice giants are the most susceptible to being ejected, as they are more loosely bound to the Sun than the terrestrial planets, inhabiting the realm of parameter space in which ejections are favoured over collisions (e.g. Ford & Rasio 2008), and they are particularly vulnerable to gravitational scattering by Jupiter.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Eccentricity (left) and inclination (right) distributions of surviving planets, weighted by the probability of each outcome. Each plot shows a reverse cumulative distribution; the y-axis shows the probability that a given planet has an eccentricity or inclination larger than a specified value.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** A simulation in which Earth was scattered out by the gas giants and captured in the Oort cloud.</figcaption></figure>

Saturn, on the other hand, was typically only ejected after an ice giant's orbit was destabilized (and likely already ejected), bringing it into dynamical contact with Jupiter.

A single planet was destroyed in 1130 of our 12 000 simulations. The total probability of losing just one planet if a star passes within 100 au is 2.7 per cent. Among simulations that lost a single planet, the most likely destruction pathways were collision of a planet with the Sun (55 per cent; 92 per cent of the time this was Mercury); impacts between two planets (30 per cent; Venus was involved 86 per cent of the time); and ejections (15 per cent; Uranus or Neptune were the ejected planet in 92 per cent of cases).

Could the Solar System's orbital architecture – specifically, Mercury's and the ice giants' survival, and the absence of planets interior to Mercury – constrain the properties of the Sun's birth cluster? We do not think so, because the Solar System's orbital architecture during the stellar birth cluster phase was different than the present-day one. The typical time-scale for the duration of the birth cluster is comparable to that of gas-dominated protoplanetary discs (e.g. Williams & Cieza 2011; Pfalzner, Steinhausen & Menten 2014). The terrestrial planets – perhaps including Mercury – may not have been completely formed during the cluster phase (e.g. Raymond et al. 2014). In addition, the giant planets' orbital architecture was likely different, with secular resonances in different locations (see

Raymond et al. 2009; Brasser et al. 2009), because their dynamical instability had not yet taken place (e.g. Morbidelli et al. 2018).

In 392 simulations, just one planet survived. Fig. 7 shows the distribution of planetary orbits in single-planet Solar Systems. Put together, a single planet has a total probability of $4.4 \times 10^{-4}$ of surviving if a star passes within 100 au. Jupiter was the sole surviving planet in 278 of these cases, although there are simulations in which each planet is the lone survivor, including 23 solo-Earth systems (total probability of $1.06 \times 10^{-5}$). The surviving planet is rarely found in the habitable zone – surviving terrestrial planets are almost universally found interior to the habitable zone and surviving giant planets exterior. None the less, there are a handful of habitable-zone planets, including a few systems with Jupiter on a very high eccentricity orbit – this naturally provokes questions related to the habitability of Jupiter's Galilean moons (see Williams, Kasting & Wade 1997) and whether they would survive system-wide dynamical instabilities (see Hong et al. 2018). Given the very high eccentricities characteristic of single-planet systems, the climates of any habitable zone planets or moons would certainly undergo strong orbital variations (e.g. Williams & Pollard 2002; Dressing et al. 2010; Bolmont et al. 2016).

## 3.6 Capture by the flyby star

In some simulations, a planet was captured by the passing star. We identified these cases, and re-ran the simulations with a much higher output frequency to capture the dynamics of the capture.

Fig. 8 shows a simulation in which Earth was captured during the flyby of a $0.987 M_{\oplus}$ star with an impact parameter of 2.6 au at relatively low speed ($21.4 \text{ km s}^{-1}$). Jupiter was the lone surviving planet, ending up with an orbital semimajor axis of 11.6 au and an eccentricity of 0.94. Saturn, Uranus, and Neptune were ejected, while Mercury, Venus, and Mars collided with the Sun. Earth's new captured orbit was colder than its current one, with a semimajor axis of 2.1 au. Its large orbital eccentricity of 0.73 would lead to strong climate variations over an orbit (e.g. Williams & Pollard 2002; Dressing et al. 2010).

Capture by another star is a low probability event. All together, 163 planets were captured by a flyby star in our simulations (including a handful of cases in which multiple planets were captured). The total probability of any planet being captured if a star passes within 100 au is $1.6 \times 10^{-4}$. Saturn and Uranus had the highest capture
<table><caption>**Table 1.** Probabilities of outcomes for each planet if a star passes within 100 au of the Sun – these should be multiplied by $0.01 \text{ Gyr}^{-1}$ for absolute probabilities.</caption><thead><tr><th>Planet</th><th>Collision with Sun</th><th>Impact with other planet</th><th>Ejection</th><th>Capture by flyby star</th><th>Trapped in Oort cloud</th></tr></thead><tbody><tr><td>Mercury</td><td>0.025</td><td>0.0080</td><td>6.2 × 10<sup>-4</sup></td><td>4.0 × 10<sup>-6</sup></td><td>—</td></tr><tr><td>Venus</td><td>0.0033</td><td>0.012</td><td>6.0 × 10<sup>-4</sup></td><td>1.2 × 10<sup>-5</sup></td><td>—</td></tr><tr><td>Earth</td><td>0.0024</td><td>0.0048</td><td>5.8 × 10<sup>-4</sup></td><td>1.1 × 10<sup>-5</sup></td><td>2.1 × 10<sup>-6</sup></td></tr><tr><td>Mars</td><td>0.012</td><td>0.0027</td><td>8.6 × 10<sup>-4</sup></td><td>2.0 × 10<sup>-5</sup></td><td>1.7 × 10<sup>-6</sup></td></tr><tr><td>Jupiter</td><td>9.2 × 10<sup>-6</sup></td><td>5.5 × 10<sup>-4</sup></td><td>5.5 × 10<sup>-4</sup></td><td>2.4 × 10<sup>-5</sup></td><td>—</td></tr><tr><td>Saturn</td><td>5.0 × 10<sup>-4</sup></td><td>3.5 × 10<sup>-4</sup></td><td>0.0032</td><td>4.0 × 10<sup>-5</sup></td><td>1.0 × 10<sup>-5</sup></td></tr><tr><td>Uranus</td><td>0.0017</td><td>5.4 × 10<sup>-4</sup></td><td>0.011</td><td>4.6 × 10<sup>-5</sup></td><td>1.5 × 10<sup>-5</sup></td></tr><tr><td>Neptune</td><td>7.2 × 10<sup>-4</sup></td><td>2.9 × 10<sup>-4</sup></td><td>0.0081</td><td>5.4 × 10<sup>-6</sup></td><td>6.3 × 10<sup>-5</sup></td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Final orbits of the 392 single-planet Solar Systems that were produced in our simulations. The approximate boundaries of the habitable zone for the present-day Sun are shaded (Kopparapu et al. 2013), included a scaling of the flux with orbital eccentricity (Williams & Pollard 2002; Barnes et al. 2009).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** A simulation in which Earth was captured by a ~ Solar-mass star onto a wider, eccentric orbit. The vertical axis shows the distance between each planet and the encounter star, which passed close to the planets at $t \approx 11\,000$ yr, capturing only the Earth.</figcaption></figure>

probabilities, of $4 \times 10^{-5}$ and $4.6 \times 10^{-5}$, respectively. Neptune and Mercury had the lowest capture probabilities, roughly an order of magnitude below those of Saturn and Uranus (see Table 1).

Fig. 9 shows the orbital radii of captured planets as a function of their new host star’s mass. The orbits of captured planets are

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Orbits of planets captured by flyby stars. For each planet, the horizontal error bar represents the radial excursion over its orbit due to orbital eccentricity. The habitable zone was estimated using the stellar mass–luminosity relation of Scalo et al. (2007) for low-mass stars and a simple $L^* \sim M_*^4$ scaling for more massive stars, and assuming that the planetary albedo varied linearly from 5 per cent for the lowest mass stars, to 30 per cent for Solar-mass stars, to 50 per cent for the highest mass stars in our sample.</figcaption></figure>

generally highly eccentric and follow a roughly thermal distribution with a median of ~0.7. The orbital radii of captured planets around their new host stars are generally somewhat wider than their original values. Still, many planets (generally the terrestrials) end up on orbits within their new host’s habitable zone.

Fig. 10 shows that the encounters that captured planets were typically with high-mass stars moving at low relative velocity. The median stellar mass among capture stars was $0.81 M_\odot$, almost five times higher than the median in all simulations of $0.16 M_\odot$. The median stellar velocity (at infinity) was just $11 \text{ km s}^{-1}$, almost four times lower than the median in all simulations of $40.2 \text{ km s}^{-1}$. Put together, the median impulse gradient of planet-capturing encounters was $47 \text{ km s}^{-1} \text{ au}^{-1}$, 283 times higher than the median among simulations of $0.17 \text{ km s}^{-1} \text{ au}^{-1}$. It is interesting to note that systems that ended up with no surviving planets at all followed the same trends as planetary captures, with stellar masses (median of $0.94 M_\odot$), velocities (median of $12.6 \text{ km s}^{-1}$), and impulse gradients (median of $70.7 \text{ km s}^{-1} \text{ au}^{-1}$). This is not a surprise, as captures are a subset of outcomes of planetary stripping, in which the planet’s liberated orbit happens to align with the perturbing star’s.

It is worth keeping in mind that our simulations assumed that the flyby star had no planets of its own. We know from decades of exoplanets surveys that virtually all stars host at least one planet (e.g. Cassan et al. 2012; Fressin et al. 2013; He, Ford & Ragozzine 2019), so this is statistically unlikely. One can imagine Earth being captured by a new host star on an orbit that is friendly to life, only to undergo
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Characteristics of stellar flybys that led to the capture of a Solar System planet. The size of each symbol is proportional to the log of the impulse gradient.</figcaption></figure>

a sterilizing giant impact with another planet. Alternately, the Sun could capture one or more planets from another star that end up wreaking havoc on the Solar System.

## 4 SIMULATIONS INCLUDING THE FULL EARTH-MOON SYSTEM

In our main set of simulations, the Earth and Moon were treated as a single particle with their combined mass. To explore the importance of the Earth–Moon system, we re-ran a subset of our simulations including the Earth and Moon as separate particles. We included them at their present-day orbital separation, maintaining their centre of mass at the same position as the single particle in our other simulations. The orbits of all other planets were identical.

We re-ran simulations in which, as a single Earth + Moon particle, Earth underwent the full range of outcomes, including ejection, collision with the Sun, impact with another planet, and regular, stable evolution. We re-ran each simulation 12 times, varying the orbital phase of the Moon around the Earth in 30° increments. While this approach does not strictly capture the exact configuration of the Earth–Moon system relative to the planets’ orbits, it allows us to get a statistical view of likely outcomes. As a sanity check, we ran 12 simulations of the Solar System with no stellar flyby and the full range of Earth-Moon configurations. All remained stable for 20 Myr.

Fig. 11 shows a simulation in which Earth survived without a Moon. First, a close stellar flyby (impact parameter $= 6.2$ au) triggered a dynamical instability that quickly spread through the whole Solar System. The terrestrial planets began a phase of strong interactions, and Mercury collided with the Sun. The Moon was stripped from the Earth after a close encounter with Venus, and collided with the Sun shortly thereafter. Mars and then Venus also collided with the Sun, and Earth survived on a closer-in, eccentric orbit ($a = 0.82$ au and $e \approx 0.3$).

It is not surprising that the Moon’s orbit should be destabilized during a dynamical instability. The Moon orbits Earth at roughly half the maximum stable orbital distance of ~0.5$R_H$, where $R_H$ is the Hill radius, defined as $R_H = a(m/3M_*)^{1/3}$, where *a* is the semimajor axis, and *m* and $M_*$ are the planet and star’s masses, respectively (Domingos, Winter & Yokoyama 2006). There is a decent chance of orbital disruption if a planet’s close encounter distance is comparable to a satellite’s orbital radius (e.g. Hong et al. 2018). This can lead to the satellite colliding with the planet,

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Evolution of a system in which the Earth and Moon were each included as separate particles. In this case, the Moon was stripped from the Earth after a close encounter with Venus, and later collided with the Sun. Earth survived as a Moon-less world.</figcaption></figure>

modification of the satellite’s orbit, or stripping of the satellite from the planet entirely. In some cases, it can even lead to capture of the satellite (as is thought to have been the case for Neptune’s moon Triton; Agnor & Hamilton 2006).

In total we re-ran 71 cases including the Earth and Moon as separate particles, for a total of 852 simulations (with 12 different Moon phases per case). These cases were divided between five broad categories of outcomes, based on the outcome for the single-particle Earth:

(i) Stable systems, in which all eight planets remained on well-behaved orbits with low AMD.

(ii) ‘Calm-ish’ systems, in which one or more planets were destroyed, but Earth survived (although sometimes on a perturbed orbit, or even one that might eventually be destabilized; the simulation from Fig. 11 is from the ‘calm-ish’ batch).

(iii) Simulations in which Earth underwent a collision with another planet.

(iv) Simulations in which Earth was ejected.

(v) Simulations in which Earth collided with the Sun.

Fig. 12 shows the outcomes for Earth in simulations that were re-run including the Earth and Moon as separate particles. The outcomes were strongly dependent on the category of outcome. All ‘Stable’ simulations remained stable when the Earth and Moon were included as separate particles.

The presence of the Moon greatly increased the danger to Earth among the ‘Calm-ish’ cases. In more than 40 per cent of simulations, Earth either underwent a collision (with the Moon in 2/3 of impacts), was ejected, or hit the Sun. In another 10 per cent of cases, the Moon was stripped from the Earth and usually collided with the Sun (or in one case, with Venus). These cases – in which a single-particle Earth escaped orbital disruption despite instability elsewhere in the system – are those in which including the full Earth–Moon system makes the biggest difference. This is because, compared with a single particle, the Earth–Moon system has a dramatically larger cross-section for interaction with other planets that may roam through the inner Solar System.

Among simulations that were more destructive for the Earth as a single particle, the outcomes were broadly similar. Among cases in which Earth as a single particle suffered an impact with another planet (usually Venus or Mars), it still underwent a giant impact in
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Outcomes for the Earth in simulations that were re-run including the full Earth-Moon system, labelled by the outcome in simulations with a single Earth + Moon particle. Situations in which Earth first underwent a collision and later was ejected or hit the Sun were included as collisions (because that was the first process that would have destroyed the habitable conditions on Earth's surface).</figcaption></figure>

57 per cent of simulations. However, the most common impactor was the Moon, which hit the Earth in 48 per cent of simulations (even though Earth subsequently was impacted by another planet in many cases).

When considering cases in which Earth as a single particle was ejected from the Solar System, the outcome with the full Earth-Moon system depended on the mode of ejection. When the ejection happened very early, triggered by the flyby star, the Earth and Moon were both ejected in all cases. However, among systems in which Earth as a single particle had been ejected later, after a system-wide dynamical instability, less than 10 per cent of simulations led to the ejection of the Earth or Moon. This was comparable to the ejection rate in the ‘Collision’ and ‘Hit Sun’ simulations. Rather, these simulations more often led to collisions with the Sun or with other planets. This is because ejection of the Earth requires a series of gravitational scattering events leading to an increase in orbital energy. Including the Earth and Moon as separate particles means that the exact series of events from the single-particle case cannot be reproduced.

In principle, the Earth-Moon system can survive being ejected through gravitational scattering by a giant planet (Debes & Sigurdsson 2007). However, among systems in which the Earth was ejected late, as a result of a system-wide dynamical instability, we only found a single case in which the Earth-Moon system survived (out of 11 ejections). This is likely because of the large number of scattering events leading to their ejection. Among system in which the Earth and Moon were both ejected early, due to the strong impulse from the flyby star, the Moon was retained in orbit around the Earth in almost every case.

Among cases in which Earth as a single planet hit the Sun, Earth still hit the Sun in 59 per cent of simulations when the Earth and Moon were included as separate particles. However, almost half of the time, Earth underwent a giant impact (usually with the Moon) before hitting the Sun. We classified these outcomes as ‘collisions’ in Fig. 12, since it was these collisions that first would have sterilized the Earth. In 30 cases, Earth (and the Moon) remained on stable orbits, bound to the Sun and to each other, without strong enough perturbations to destabilize the system (at least during the 20 Myr integration).

<figure><img src="image_16.png" /><figcaption>**Figure 13.** A simulation in which the Earth and Moon were both captured during the flyby of a low-mass star, but the Moon was later ejected.</figcaption></figure>

We also re-ran some simulations in which Earth was captured by the flyby star. The outcomes with the Earth and Moon as separate particles were system- and phase-dependent. In some cases, all 12 simulations resulted in Earth's capture by the star; the Moon was also captured in more than half of simulations.<fnref n="4" /> In other cases, Earth was only captured in a minority of simulations for a given flyby. This is simply because capture depends on a precise alignment of velocity vectors, and these vary depending on the Moon's exact orbital phase. However, even among systems in which both the Earth and Moon were captured by the flyby star, the Moon was always lost from the Earth and essentially became its own planet. While we did not follow the long-term evolution of the orbits around the flyby star, the Moon was in all cases on an unstable orbit crossing that of the Earth. We expect such systems to usually end in collision, or sometimes in ejection.

Fig. 13 shows a case in which the Earth and Moon were both captured by a 0.26 M☉ star. However, the Moon was no longer in orbit around Earth, and after a series of gravitational encounters, the Moon was ejected into interstellar space. Earth survived on a much colder and more eccentric orbit than its present-day one.

## 5 EARTH'S FUTURE HABITABILITY

Fig. 14 shows the distribution of post-flyby orbits of Earth, excluding cases in which Earth was captured around another star or ejected into interstellar space. Earth's post-flyby semimajor axis tends to be close to its present-day one at 1 au, although scattering events can strand it closer or further from the Sun. Of course, orbits far from 1 au tend to be very eccentric, although this is not a deal breaker for habitability (Williams & Pollard 2002; Barnes et al. 2009; Dressing et al. 2010; Bolmont et al. 2016).

The surviving Earths in Fig. 14 have a range of orbital inclinations relative to Earth's original orbital plane. While the probability of Earth surviving with a large inclination is small (see Fig. 5), this would certainly have an effect on Earth's climate, by changing its effective obliquity and, if other planets are present in the system, make it oscillate in time (Spiegel, Menou & Scharf 2009; Raymond et al. 2011; Armstrong et al. 2014).

<footnote n="4">⁴While there were many cases in which Earth was captured but the Moon was not, there were no simulations in which the Moon was captured without the Earth.</footnote>
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Post-flyby orbits of Earth around the Sun. A handful of the high-semimajor axis orbits are unlikely to be stable on long time-scales, but the two beyond $10^4$ au represent stable trapping in the Oort cloud – the outermost case is the example from Fig. 6. The colour of each planet corresponds to its inclination relative to Earth's original orbital plane.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 15.** The probability that Earth survives on a hotter or cooler orbit as a function of the number of surviving planets, in the event of a sub-100 au stellar flyby.</figcaption></figure>

What are the odds that a stellar flyby could improve the prospects for the long-term future of life on Earth? The best-case scenario is if Earth's post-flyby orbit is cooler than its present-day one. Unfortunately, a sub-100 au flyby only has a 0.28 per cent probability of making that happen. We considered an orbit to be 'cooler' if the incident Solar flux dropped to less than 90 per cent of the flux on a circular orbit at 1 au. There was a substantially higher (0.79 per cent) chance of Earth's post-flyby orbit becoming hotter by at least 10 per cent in flux. In the vast majority of cases (> 99 per cent), there was no significant change in flux.

Fig. 15 shows the probability of Earth surviving on a cooler (or hotter) orbit as a function of the number of surviving Solar System planets. There is a roughly comparable probability (of ~0.0005, or 1 part in 2000) of Earth ending up on a cooler orbit when the post-flyby Solar System contained 4, 5, 6, 7, or 8 planets. This takes into account a number of factors, including the probability of Earth surviving (including relative probabilities of different outcomes), and Earth's final orbit. There was a much higher chance of Earth ending up on a cool orbit among systems with 4–5 surviving planets than among systems with eight surviving planets, but the odds of

that happening at all was much lower, given the correlation between outcomes and the impulse gradient of the flyby star (Fig. 1).

Earth's typical eccentricity progressively increased in systems with fewer surviving planets, and was markedly higher in cases when Earth's orbit was cooler than its present-day one. In systems in which 7–8 planets survived, Earth typically was perturbed onto a slightly wider orbit, with a median semimajor axis of 1.08 au and a median eccentricity of 0.10. However, in systems in which fewer planets survived, cooler Earths tended to have significantly higher eccentricities, with a median value of 0.18, 0.21, and 0.27 for 6, 5, and 4 planets surviving, respectively. The typical eccentricities among Earths that ended up on hotter orbits were even higher, with median values of 0.19, 0.28, and 0.34 for 6, 5, and 4 planets surviving, respectively. These represent strongly perturbed systems in which Earth's orbital radius was shifted (either inward or outward) by gravitational scattering with the other planets. Such strong scattering events invariably increase the planets' eccentricities.

A plausible long-term prospect for life on Earth may be for our planet to be ejected into interstellar space. Laughlin & Adams (2000) explored the idea of the frozen Earth, and modelled its long-term thermal evolution (without any assumptions related to the atmosphere). They found that it would take roughly 1 Myr for Earth's surface to freeze over completely, although they note that life could continue to thrive in hydrothermal vents and in the deep subsurface.

The outer edge of the habitable zone may effectively extend to infinity if a planet like Earth has a thick enough thermal blanket in the form of a many-bar, hydrogen-rich atmosphere (Stevenson 1999; Pierrehumbert & Gaidos 2011) or a several km-thick layer of ice (Abbot & Switzer 2011). In such cases, the planet's internal heat flux can provide enough heat to maintain stable liquid water. If a free-floating Earth also kept its moon (as was the case in many of our simulations; see Section 4), then tides could provide an additional source of heat, which would be much higher than the present-day tidal heat flux if the Moon's orbit was modestly perturbed during ejection (Debes & Sigurdsson 2007).

Free-floating planets are known to be abundant. To date, hundreds have been detected via direct imaging (e.g. Miret-Roig et al. 2022) and gravitational microlensing (Mróz et al. 2017, 2020; Sumi et al. 2023), with masses down to almost Earth's mass. If they possess sufficient water and adequate thermal blankets, some of these (or their moons) may represent viable habitats for life.

# 6 LONG-TERM STABILITY OF THE SOLAR SYSTEM

Studies evaluating the dynamical stability of the Solar System date back hundreds of years (for a review, see Laskar 2012). The stability of a pair of planets can be evaluated analytically (Marchal & Bozis 1982) or numerically (Gladman 1993). The AMD can itself be used as a proxy for dynamical stability of a pair of planets, given that it controls that maximum eccentricity that each planet can obtain and therefore the possibility of the two planets' orbits crossing (Laskar & Petit 2017; Petit, Laskar & Boué 2018).

Brown & Rein (2022) simulated the long-term stability of the Solar System under the influence of weak, distant stellar flybys. They found that the probability of the system becoming unstable increased drastically when Neptune's semimajor axis was perturbed by more than 1 part in a thousand (roughly 0.03 au). Fig. 16 shows the fractional change in the angular momentum deficit AMD of our systems with 6, 7 or 8 surviving planets as a function of the fractional change in Neptune's semimajor axis (calculated using Jacobi coordinates, to maintain maximum stability in semimajor
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Fractional change in the angular momentum deficit of final 8-, 7-, and 6-planet Solar Systems as a function of the fractional change in Neptune's orbit.</figcaption></figure>

axis). We see the same general behaviour as in the simulations of Brown & Rein (2022), with a strong trend toward higher AMD in systems with a larger change in Neptune's semimajor axis. This same trend is seen both within sets of simulations with 6, 7, or 8 surviving planets, and between the sets themselves. This essentially boils down to the same trend that we showed in Fig 1: stronger stellar perturbations (which have a higher impulse gradient, and which lead to larger changes in Neptune's semimajor axis) lead to more strongly disrupted systems.

# 7 DISCUSSION

## 7.1 Comparison with previous work

To our knowledge, the paper that is most similar in scope to our own is Laughlin & Adams (2000), who simulated the passage of binary stars by the Solar System and evaluated the probability of different outcomes. Our paper differs from Laughlin & Adams (2000) in several ways: (i) they used only binary flyby stars, while we used only single flyby stars (see discussion in Section 2.1); (ii) our flybys followed a logarithmic distribution in impact parameter *b* (which we accounted for in calculating probabilities; see Section 2.3), while theirs followed a physically motivated *b*² distribution; (iii) we included a prescription to account for the Galactic tidal field and therefore allow for the capture of planets into the Oort cloud; (iv) we included the full eight planets in all of our simulations and ran them for 20 Myr, whereas, due to computational limitations, Laughlin & Adams (2000) only included a subset of planets in most integrations, which only lasted 1 Myr; and (v) we ran a subset of simulations that also included the Moon (see Section 5).

Our results roughly match those of Laughlin & Adams (2000) in terms of the different outcomes and general trends. Of course, there are some outcomes that were out of reach for our simulations; for instance, given that our flyby stars were single, it was impossible for them to be captured by the Sun, which happened in the simulations of Laughlin & Adams (2000). Likewise, some of our outcomes were out of reach for Laughlin & Adams (2000), such as capture in the Oort cloud and interactions that involved the Moon.

There are specific outcomes (such as those from Table 1) that we can compare, to gauge how closely our results match those of Laughlin & Adams (2000). They found a probability of 1 in 400 000 that Earth would be ejected from the Solar System in the next 3.5 Gyr, which amounts to $7.1 \times 10^{-7}$ Gyr$^{-1}$. Our simulations yielded a

significantly higher probability of $5.8 \times 10^{-6}$ Gyr$^{-1}$. They found the odds that Earth is captured by a passing star to be roughly 1 in 2 million, or $1.4 \times 10^{-7}$ Gyr$^{-1}$, which is very close to our calculated value of $1.1 \times 10^{-7}$ Gyr$^{-1}$. Laughlin & Adams (2000) found a probability of $\sim 9 \times 10^{-6}$ Gyr$^{-1}$ that Earth's post-flyby eccentricity would be larger than 0.5. We found a somewhat higher value of $2.4 \times 10^{-5}$ Gyr$^{-1}$. Despite these modest differences, we consider our results to be broadly consistent.

## 7.2 Limitations

Our numerical simulations are admittedly imperfect on several fronts. First, we did not include the gravitational influence of all Solar System bodies. Laskar et al. (2011) showed that the largest asteroids provide a non-negligible contribution to chaos among the terrestrial planets on a $\sim 60$ Myr time-scale. Given the relatively short duration of our simulations and the Gyr or longer time-scale for chaos to affect the stability of the terrestrial planets' orbits (Laskar & Gastineau 2009; Zeebe 2015; Abbot et al. 2021), it seems unlikely that neglecting asteroids has a significant impact on the stability of the planets' orbits. Yet we can imagine the possibility of asteroid or comet showers triggered by the destabilization of small body orbits during or after a flyby.

Second, as discussed above, we did not include binaries among the encounter stars. Multiple-star systems are extremely common (Duquennoy & Mayor 1991), and binary encounters can lead to certain interesting outcomes such as stellar capture (see Laughlin & Adams 2000), which we could not capture in our simulations. Li & Adams (2015) showed that the flyby of a binary can be approximated as simply the flyby of two single stars for encounter speeds faster than $\sim 20$ km s$^{-1}$. Many of the encounters that resulted in capture or complete planet stripping were at low encounter speed. The cross-section for interaction of a binary is more than twice that of a single star at low speeds (Li & Adams 2015), which implies that our results may somewhat underestimate the probability of these outcomes.

Third, we did not include the full Galactic environment, but only a single close encounter with one passing star, and a simple recipe for torques from the Galactic tidal field. Studies such as Portegies Zwart (2021) developed a more consistent treatment for the Galactic environment, which does indeed affect the outskirts of planetary systems. On a similar note, we only included a simplified approach for general relativity (taken from Saha & Tremaine 1992), although other approaches exist (see Brown & Rein 2023).

On a similar note, our simulations are limited to the conditions of the Sun's Galactic neighbourhood. The stellar number densities and velocity dispersions each vary by orders of magnitude among different Galactic locations such as globular or open (or embedded) star clusters, and the Galactic Bulge and core (see Table 1 of Brown & Rein 2022). The frequency of stellar flybys and their impact on planetary systems (via the impulse gradients) will vary accordingly. The Solar System has experienced at least two different Galactic environments: its birth cluster and the local neighbourhood. It would be interesting to take the Sun's full Galactic history into account when assessing the stability of the Solar System.

# 8 CONCLUSIONS

In this paper, we evaluated the dynamical evolution of the Solar System if a passes within 100 au of the Sun. There is a $\sim 1$ per cent chance of this happening in the next Gyr (Zink et al. 2020; Brown & Rein 2022), which is the same time frame for the Earth to leave the habitable zone due to the increasing Solar luminosity. Each of our
12 000 simulations simulated the flyby of a single star within 100 au the Solar System. We sample the impact parameter *b* in a logarithmic fashion, allowing us to capture rare outcomes and still derive realistic probabilities (see Section 2.3).

Our main results are:

(i) The most likely outcome of a sub-100 au stellar flyby is for all eight planets to survive on orbits similar to their present ones (see Fig. 2). There is a 92 per cent probability that the Solar System's angular momentum deficit (see equation 2) will also remain smaller than twice the current value with all pairs of neighbouring planets on Hill stable orbits.

(ii) The key parameter that controls the outcome is the impulse gradient, a measure of the strength of the perturbation from the stellar flyby (see Fig. 1 and equation 1).

(iii) The post-flyby orbits of surviving planets can be very different than their present-day ones, in some cases with very high eccentricities and inclinations (Fig. 4). While extremely excited orbits are a low-probability outcome (Fig. 5), they would have consequences for planetary climate.

(iv) If a star passes within 100 au of the Sun, the most likely planetary destruction pathways are: Mercury hitting the Sun (2.5 per cent probability), Mars colliding with the Sun (1.2 per cent), Venus colliding with another planet (likely Earth or Mercury; 1.2 per cent), or the ejection of Uranus (1.1 per cent) or Neptune (0.8 per cent). Table 1 lists probabilities for all major outcomes for each planet.

(v) One of the most dramatic outcomes were simulations in which all eight planets were stripped from the Solar System by a particularly strong stellar encounter (see bottom-right panel in Fig. 3). Luckily, the probability of this happening is just $8 \times 10^{-8}$ Gyr$^{-1}$.

(vi) In some simulations, a planet became trapped in the Sun's Oort cloud. This is most likely to happen to the ice giants or Saturn (see Table 1), although Earth and Mars were each trapped in the Oort cloud in certain simulations.

(vii) In some simulations, one or more planets was captured into orbit around the flyby star (see Figs 8 and 9). The eccentricities of captured planets tend to be very high. Of course, most flyby stars would likely have their own planetary systems, which means that a captured planet may well undergo a phase of collisions – or at least strong dynamical interactions – with the flyby star's native planets (or, for that matter, with other captured planets; see Fig. 13).

(viii) Including the Moon as a separate particle (rather than having a single Earth–Moon particle with the sum of their masses) makes no apparent difference in the stability of the Solar System. However, the presence of the Moon does open a window of new outcomes, the most common and destructive of which is a collision with the Earth (see Fig. 12 and Section 4).

(ix) A stellar flyby is unlikely to allow Earth to remain habitable beyond a horizon of ~1 Gyr. In the event of a sub-100 au flyby, there is only a 0.28 per cent probability that Earth will survive on an orbit that receives less than 90 per cent of the present-day stellar flux. The odds of Earth being ejected into interstellar space are even smaller, just $6 \times 10^{-6}$ in the next Gyr. If Earth is indeed ejected – or captured on a wide orbit around a low-luminosity star – the surface will freeze over completely on a time-scale of roughly 1 million years (Laughlin & Adams 2000).

(x) Regarding the long-term survival of the Solar System, our simulations are consistent with the work of Brown & Rein (2022). We see a clear correlation between the fractional change in AMD and the fractional change in Neptune's semimajor axis (Fig. 16), although we argue that the impulse gradient is really the key factor controlling the outcome.

Despite the diversity of potential evolutionary pathways, odds are high that our Solar System's current situation will not change. Earth's progressive heating from the brightening Solar luminosity (not to mention human-driven carbon emissions) will continue unabated. The Universe is statistically unlikely to help us out by providing a stellar flyby that will lead to a cooler Earth. Humanity's best solution is to help itself.

## ACKNOWLEDGEMENTS

The eventuality for the Earth to be ejected and its possible fates after such event initially came from a discussion between FS and the author luvan (https://www.luvan.org/luvan-grav/about), who is currently working on a science fiction verse novel. This research has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 682903, PI: H. Bouy), and from the French State in the framework of the 'Investments for the future' Program, IdEx Bordeaux, reference ANR-10-IDEX-03-02. SNR and FS acknowledge funding from the French Programme National de Planétologie (PNP). SNR, FS and HB acknowledge funding in the framework of the Investments for the Future programme IdEx, Université de Bordeaux/RRI ORIGINS. NAK acknowledges funding from NSF CAREER award 1846388 and NASA Exoplanets Research Program grant 80NSSC19K0445.

## DATA AVAILABILITY

All simulation data and analysis codes will be made available upon reasonable request.

## REFERENCES

Abbot D. S., Switzer E. R., 2011, ApJ, 735, L27

Abbot D. S., Webber R. J., Hadden S., Seligman D., Weare J., 2021, ApJ, 923, 236

Agnor C. B., Hamilton D. P., 2006, Nature, 441, 192

Armstrong J. C., Barnes R., Domagal-Goldman S., Breiner J., Quinn T. R., Meadows V. S., 2014, Astrobiology, 14, 277

Barnes R., Jackson B., Greenberg R., Raymond S. N., 2009, ApJ, 700, L30

Batygin K., Morbidelli A., Holman M. J., 2015, ApJ, 799, 120

Bolmont E., Libert A.-S., Leconte J., Selsis F., 2016, A&A, 591, A106

Brasser R., Morbidelli A., 2013, Icarus, 225, 40

Brasser R., Morbidelli A., Gomes R., Tsiganis K., Levison H. F., 2009, A&A, 507, 1053

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

Brown G., Rein H., 2023, MNRAS, 521, 4349

Cassan A. et al., 2012, Nature, 481, 167

Chakrabarti S., Chang P., Lam M. T., Vigeland S. J., Quillen A. C., 2021, ApJ, 907, L26

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

Chatterjee S., Ford E. B., Matsumura S., Rasio F. A., 2008, ApJ, 686, 580

Debes J. H., Sigurdsson S., 2007, ApJ, 668, L167

Domingos R. C., Winter O. C., Yokoyama T., 2006, MNRAS, 373, 1227

Dones L., Brasser R., Kaib N., Rickman H., 2015, Space Sci. Rev., 197, 191

Dressing C. D., Spiegel D. S., Scharf C. A., Menou K., Raymond S. N., 2010, ApJ, 721, 1295

Duquennoy A., Mayor M., 1991, A&A, 248, 485

Feast M., Whitelock P., 1997, MNRAS, 291, 683

Ford E. B., Rasio F. A., 2008, ApJ, 686, 621

Fressin F. et al., 2013, ApJ, 766, 81

Gladman B., 1993, Icarus, 106, 247

He M. Y., Ford E. B., Ragozzine D., 2019, MNRAS, 490, 4575

Hong Y.-C., Raymond S. N., Nicholson P. D., Lunine J. I., 2018, ApJ, 852, 85
Kaib N. A., Roškar R., Quinn T., 2011, *Icarus*, 215, 491

Kasting J. F., Whitmire D. P., Reynolds R. T., 1993, *Icarus*, 101, 108

Kopparapu R. K. et al., 2013, *ApJ*, 765, 131

Laskar J., 1997, A&A, 317, L75

Laskar J., 2012, preprint (arXiv:1209.5996)

Laskar J., Gastineau M., 2009, *Nature*, 459, 817

Laskar J., Petit A. C., 2017, A&A, 605, A72

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

Laughlin G., Adams F. C., 2000, *Icarus*, 145, 614

Leconte J., Forget F., Charnay B., Wordsworth R., Pottier A., 2013, *Nature*, 504, 268

Levison H. F., Dones L., Duncan M. J., 2001, *AJ*, 121, 2253

Li G., Adams F. C., 2015, *MNRAS*, 448, 344

Limbach M. A., Turner E. L., 2015, *Proc. Natl. Acad. Sci.*, 112, 20

Marchal C., Bozis G., 1982, *Celest. Mech.*, 26, 311

Miret-Roig N. et al., 2022, *Nat. Astron.*, 6, 89

Morbidelli A., Nesvorny D., Laurenz V., Marchi S., Rubie D. C., Elkins-Tanton L., Wieczorek M., Jacobson S., 2018, *Icarus*, 305, 262

Mróz P. et al., 2017, *Nature*, 548, 183

Mróz P. et al., 2020, *ApJ*, 903, L11

Petit A. C., Laskar J., Boué G., 2018, *A&A*, 617, A93

Pfalzner S., Steinhausen M., Menten K., 2014, *ApJ*, 793, L34

Pierrehumbert R., Gaidos E., 2011, *ApJ*, 734, L13

Portegies Zwart S., 2021, *A&A*, 647, A136

Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P., 1992, Numerical Recipes in FORTRAN. The Art of Scientific Computing. University Press, Cambridge

Raymond S. N., O'Brien D. P., Morbidelli A., Kaib N. A., 2009, *Icarus*, 203, 644

Raymond S. N., Armitage P. J., Gorelick N., 2010, *ApJ*, 711, 772

Raymond S. N. et al., 2011, *A&A*, 530, A62

Raymond S. N., Kokubo E., Morbidelli A., Morishima R., Walsh K. J., 2014, Protostars and Planets VI. University of Arizona Press, Tucson, AZ, p. 595

Raymond S. N., Izidoro A., Kaib N. A., 2023, *MNRAS*, 524, L72

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

Rickman H., Fouchard M., Froeschlé C., Valsecchi G. B., 2008, *Celest. Mech. Dyn. Astron.*, 102, 111

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

Scalo J. et al., 2007, *Astrobiology*, 7, 85

Selsis F., Kasting J. F., Levrard B., Paillet J., Ribas I., Delfosse X., 2007, *A&A*, 476, 1373

Spiegel D. S., Menou K., Scharf C. A., 2009, *ApJ*, 691, 596

Stevenson D. J., 1999, *Nature*, 400, 32

Sumi T. et al., 2023, *AJ*, 166, 108

Tremaine S., 1993, in Phillips J. A., Thorsett S. E., Kulkarni S. R., eds, ASP Conf. Ser. Vol. 36, Planets Around Pulsars. Astron. Soc. Pac., San Francisco, p. 335

Turrini D., Zinzi A., Belinchon J. A., 2020, *A&A*, 636, A53

Williams J. P., Cieza L. A., 2011, *ARA&A*, 49, 67

Williams D. M., Pollard D., 2002, *Int. J. Astrobiol.*, 1, 61

Williams D. M., Kasting J. F., Wade R. A., 1997, *Nature*, 385, 234

Wolf E. T., Toon O. B., 2014, *Geophys. Res. Lett.*, 41, 167

Zeebe R. E., 2015, *ApJ*, 811, 9

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

This paper has been typeset from a *TEX/LATEX* file prepared by the author.