{'fraction_non_alphanumeric': 0.03620659738957759, 'fraction_numerical': 0.024311245638049702, 'mean_word_length': 4.634368889854443, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The fates of Solar system analogues with one additional distant planet

The potential existence of a distant planet ("Planet Nine") in the Solar system has prompted a re-think about the evolution of planetary systems. As the Sun transitions from a main sequence star into a white dwarf, Jupiter, Saturn, Uranus and Neptune are currently assumed to survive in expanded but otherwise unchanged orbits. However, a sufficiently-distant and sufficiently-massive extra planet would alter this quiescent end scenario through the combined effects of Solar giant branch mass loss and Galactic tides. Here, I estimate bounds for the mass and orbit of a distant extra planet that would incite future instability in systems with a Sun-like star and giant planets with masses and orbits equivalent to those of Jupiter, Saturn, Uranus and Neptune. I find that this boundary is diffuse and strongly dependent on each of the distant planet's orbital parameters. Nevertheless, I claim that instability occurs more often than not when the planet is as massive as Jupiter and harbours a semimajor axis exceeding about 300 au, or has a mass of a super-Earth and a semimajor axis exceeding about 3000 au. These results hold for orbital pericentres ranging from 100 to at least 400 au. This instability scenario might represent a common occurrence, as potentially evidenced by the ubiquity of metal pollution in white dwarf atmospheres throughout the Galaxy.


INTRODUCTION
Thousands of planets outside of our Solar system have already been discovered despite our potential ignorance of what planets may reside in our own backyard. This arresting notion (Brown, Trujillo & Rabinowitz 2004;Gladman & Chan 2006;Iorio 2012) was given added credence with the discovery of the second Sednoid, 2012 VP113 (Trujillo & Sheppard 2014), because both that object and Sedna have arguments of pericentre clustered around −50 • (see their table 1). This clustering is not likely to be a result of observational bias (de la Fuente Marcos & de la Fuente Marcos 2014), and may be explained by one or more planets which exist beyond the orbit of Neptune but are so far invisible to us (Iorio 2014;Luhman 2014;Trujillo & Sheppard 2014;Gomes, Soares & Brasser 2015;Iorio 2015). However, the existence of these two scattered disc objects do not necessarily require the presence of additional planets (Jílková et al. 2015).
The recent finding that several distant Kuiper Belt Objects exhibit this same clustering in orbital space, as well a clustering in physical space , has prompted a febrile response amongst the public and scientific community. Constraining the location and size of a putative "Planet Nine" has been the focus of many subsequent publications, through both orbital dynamics (Beust 2016;; de la Fuente Marcos & de la Fuente Marcos 2016a,b; de la Fuente Marcos, de la Fuente Marcos & Aarseth 2016; Fienga et al. 2016;Holman & Payne 2016a,b;Lawler et al. 2016;Malhotra, Volk & Wang 2016) and intrinsic physical properties (Cowan, Holder & Kaib 2016;Fortney et al. 2016;Ginzburg, Sari & Loeb 2016;Linder & Mordasini 2016;Toth 2016).
Overall and roughly, these studies suggest that Planet Nine is more massive than the Earth and resides at hundreds or thousands of au away from the Sun (other similarly massive planets could lie further away). Other studies investigated potential origins for Planet Nine, including in-situ formation , scattering into its current (theorized) orbit ) and capture from other stars in the Sun's birth cluster (Mustill, Raymond & Davies 2016). The survivability of Planet Nine due to passing stars has also been investigated (Li & Adams 2016). Additional planets which orbit the Sun also could exist but evidence for their presence has yet to be marshalled.
Two points not emphasized in the above studies are: (1) how our knowledge of the fate of the Solar system changes with the presence of planets beyond Neptune's orbit, and (2) the deeper fundamental issue of what the implications are for similarly-constructed systems and more generally exoplanetary science. These two concepts are linked through white dwarf planetary systems, and represent the motivation for this paper.

The fate of an eight-planet Solar system
The Sun will leave the main sequence in about 6.5 Gyr, and undergo drastic changes (Fig. 1). Its radius will increase by a factor of about 230, it will lose almost half of its current mass, and its luminosity will reach a peak value which is about 4000 times its current value (see e.g Schröder & Connon Smith 2008, Veras & Wyatt 2012and fig. 3 of Veras 2016. The Sun will become so large that its radius will extend just beyond where the Earth currently sits. These major changes will occur in two phases. The red giant branch phase will last about 800 Myr. In this timespan, the Sun will gradually lose about a quarter of its mass. The second phase, when the Sun becomes an asymptotic giant branch star, is quicker: lasting just 5 Myr. Another quarter of the Sun's mass will be lost during this period. During both phases, the radius of the Sun will extend out to nearly the Earth's distance. The consequences for the inner Solar system will be profound. The terrestrial planets, which are likely to remain in stable orbits until the end of the main sequence at the approximately 99% level (Laskar & Gastineau 2009;Batygin, Morbidelli & Holman 2015;Zeebe 2015), will be in danger. Mercury and Venus will be engulfed, and the Earth will be on the edge of survivability (Rybicki & Denis 2001;Schröder & Connon Smith 2008). Mars will be roasted, but should survive, because it will escape being ensnared by the tidal reach of the Sun (Villaver & Livio 2009;Kunitomo et al. 2011;Mustill & Villaver 2012;Adams & Bloch 2013;Nordhaus & Spiegel 2013;Villaver et al. 2014;Staff et al. 2016). Asteroid belt constituents between 100 m and 10 km in radius will be spun up to breakup speed (Veras, Jacobson & Gänsicke 2014a), creating a sea of debris, some of which may be water-rich (Jura & Xu 2010, 2012Farihi, Gänsicke & Koester 2013;Raddi et al. 2015;Malamud & Perets 2016).
The consequences for the giant planets, however, will be more benign. Jupiter, Saturn, Uranus and Neptune will increase their semimajor axes by a factor of about two each, and not undergo scattering nor instability (Duncan & Lissauer 1998), even though the chemistry of at least Jupiter's atmosphere will be fundamentally altered (Villaver & Livio 2007;Spiegel & Madhusudhan 2012). The giant planet eccentricities will remain effectively fixed because they reside within the adiabatic limit, beyond which stellar mass loss changes both eccentricity and semimajor axis (Veras et al. 2011). The present-day "adiabatic" 1 limit for the Solar system specifically lies between about 10 3 and 10 4 au (Veras & Wyatt 2012). Further, even though mass loss changes stability limits (Debes & Sigurdsson 2002;Mustill, Marshall & Villaver 2013;Portegies Zwart 2013;Voyatzis et al. 2013;Mustill, Veras & Villaver 2014;Veras & Gänsicke 2015;Veras et al. 2016a), this change will not be large enough to affect the giant planets. I perform some simulations here also to back up this statement.
Moons, the Kuiper Belt, scattered disc and Oort Cloud will also be affected. Moons of planets will become more entrenched in the Hill spheres of the host planets, and would stay there in the absence of a planetary scattering event (Payne et al. 2016a,b). Although known Kuiper Belt and scattered disc objects are within the adiabatic limit, many are likely to become unstable as the stability limits between Neptune and Kuiper Belt objects change as the Sun's mass decreases (e.g. ) and its luminosity increases (Veras, Eggl & Gänscike 2015a). The Oort Cloud will be both excited and decimated (Veras & Wyatt 2012;Veras et al. 2014b), which will alter the influx of comets into the inner Solar system (Alcock, Fristrom & Siegelman 1986;Parriott & Alcock 1998;Veras, Shannon & Gänscike 2014c;Stone, Metzger & Loeb 2015), and these comets can be subsequently perturbed by radiation from the Solar white dwarf (Veras, Eggl & Gänsicke 2015b).

The effect of additional planets
The presence of additional planets (one or more) can trigger instability purely through Lagrange instability as the Sun becomes older. This phenomenon in main sequence extrasolar systems has been well-explored (e.g. Chambers, Wetherill & Boss 1996;Barnes & Greenberg 2006Chatterjee et al. 2008;Deck, Payne & Holman 2013;Veras & Mustill 2013;Davies et al. 2014;Marzari 2014;Petrovich 2015;Pu & Wu 2015) and will have greater predictive power as asteroseismology continues to pinpoint the ages of old stars (Campante, Barclay & Swift 2015;Silva Aguirre et al. 2015;Veras et al. 2015c;North et al. 2016). Mass loss from giant branch evolution will change the limits further, as outlined above. Additionally, a planet residing in the non-adiabatic regime will have its eccentricity changed during this mass loss, altering the limits even further. When combined with Galactic tides and stellar flybys, the result may be a significantly more dynamically active evolved system than previously envisaged.
Hence, although Jupiter, Saturn, Uranus and Neptune by themselves will remain stable throughout the Sun's giant branch and white dwarf phases, the presence of additional planets might trigger instability after the Sun has left the main sequence, but not before. This scenario might be common throughout the Milky Way if extrasolar planetary system architectures from about 5 to 30 au resemble the Solar system's, and contain a sufficiently distant and massive extra planet. In this paper, I quantify this scenario, and place Figure 1. The future evolution of the Sun's mass (M⋆), radius (R⋆), effective temperature (Teff⋆) and luminosity (L⋆) during the transition between the red giant branch phase (yellow background) and the asymptotic giant branch phase (pink background). These evolutionary tracks were computed with the SSE code (Hurley et al. 2000). The η parameter (Equation 1) determines the extent of the mass loss on the red giant branch phase, which has implications for evolution on the asymptotic giant branch phase. The phase transition (7.8587 Gyr) is shown for η = 0.5; other phase transitions are 7.8571 Gyr for η = 0.2 and 7.8738 Gyr for η = 0.8. The Sun eventually becomes a white dwarf with a mass of 0.511 − 0.534M ⊙ . rough bounds on the mass and orbit of one additional planet which can trigger future instability.

The polluted planet-host star
The importance of considering extrasolar systems in general arises from the outcome of the instability during the white dwarf phase. We know that planetary systems around white dwarfs are just as common as those orbiting main sequence stars, primarily through the detection of metal debris in the atmospheres of the degenerate stars. Such "polluted" white dwarfs now number in the thousands (Dufour et al. 2007;Kleinman et al. 2013;Kepler et al. 2015, 2016 andHollands et al. in preparation). High-sensitivity observational surveys reveal that between one-quarter and onehalf of the known Milky Way white dwarfs are estimated to host planetary debris (Zuckerman et al. 2003(Zuckerman et al. , 2010Koester, Gänsicke & Farihi 2014).
The specific primary mechanism which transports the rocky bodies to the white dwarf is still debatable (see section 7 of Veras 2016) but relies on interactions between minor planets and larger bodies Bonsor & Wyatt 2012;Debes, Walsh & Stark 2012;Frewen & Hansen 2014). Instabilities which arise in white dwarf planetary systems Voyatzis et al. 2013;Mustill et al. 2014;Veras & Gänsicke 2015;Bonsor & Veras 2015;Payne et al. 2016a,b;Hamers & Portegies Zwart 2016;Petrovich & Muñoz 2016) promote an architecture conducive to pollution by placing planets and liberated moons in orbits that sweep through regions of space which may access reservoirs of debris. Examples of such reservoirs are those that arise from planet-planet collisions (Shannon et al. in preparation) or YORP-induced asteroid breakup (Veras et al. 2014a). Therefore, the ubiquity of white dwarf pollution throughout the Milky Way (Zuckerman et al. 2003(Zuckerman et al. , 2010Koester et al. 2014) suggests that post-main-sequence instability is common, an idea that would be aided by the presence of extrasolar analogues of Planet Nine: i.e. potentially massive planets on wide orbits that remain stable during the main sequence phase but begin to scatter gravitationally during or after giant branch mass loss. If the presence of planets at hundreds or thousands of au is the norm, regardless of their dynamical origin, then that scenario is broadly consistent with the observed end states of planetary systems.

Outline for this paper
In this paper I characterize the mass and orbital parameters of an additional distant planet that would create gravitational instability amongst giant planet analogues orbiting a Sun-like star after the latter turns off of the main sequence. This result would also apply to the Solar system if no Planet Nine is found now: at some point later in the Sun's main sequence lifetime the Sun may capture an additional planet (Sumi et al. 2011;Perets & Kouwenhoven 2012;Varvoglis, Sgardeli & Tsiganis 2012).
In Section 2, I provide more detail about the evolution of Sun-like stars, before, in Section 3, isolating the varied effects that come into play during the giant branch and white dwarf phases. My computational method is detailed in Section 4. I explain the initial conditions and report on my simulation results in Section 5 before summarizing the conclusions of this work in Section 6.

THE SUN'S FUTURE EVOLUTION
Amongst the Sun's many physical properties (such as radius and luminosity), its mass is the most important one for my study. On the main sequence, the Sun's mass loss rate (≈ 2.4 × 10 −14 M⊙ yr −1 ) is negligible (Vial 2013). However, the Sun's future evolution beyond the main sequence is unknown. A potentially good guess is to consider the typical evolution of other 1M⊙ stars with metallicities of Z = 0.02. Although these two parameters are well-established for the Sun, an important unknown is the future rate of mass loss during the red giant branch phase.
A traditional formulation of this rate is the Reimers prescription (Reimers 1975(Reimers , 1977, which has now been improved (Schröder & Cuntz 2005), and may be parametrized as where M⋆, L⋆, R⋆, T⋆ and g⋆ are the mass, luminosity, radius, temperature and surface gravity of the star. The value of η = 0.2 reproduces the coefficient given in Schröder & Cuntz (2005), and η = 0.5 reproduces the coefficient given in the traditional formulation, which does not contain the terms dependent on T⋆ nor g⋆.
Equation (1) applies to the giant branch phase only. Mass loss evolution on the asymptotic giant branch phase is qualitatively different, and is typically characterised by the prescription of Vassiliadis & Wood (1993): where dM⋆/dt is computed in M⊙ yr −1 and such that where the value P is computed in years. The continued common use of this prescription even over two decades after it was published helps indicate its robustness in the face of new observations. This formulation also importantly includes the "superwind" (peak mass loss) which occurs at the "tip" (end-point) of the asymptotic giant branch (Lagadec & Zijlstra 2008).
After the asymptotic giant branch phase, the Sun will become a white dwarf. The Solar white dwarf will not lose mass nor change its radius, but will gradually dim. Figure  1 illustrates the Sun's evolution for (η = 0.2, 0.5, 0.8) from the SSE code, which uses the traditional Reimers formulation and will be described in more detail in Section 4. I adopted these values in order to encompass a realistic range, as outlined above, and further to conform to the range adopted in Veras & Wyatt (2012). Fig. 1 displays the Sun's mass evolution during the transition between the red giant branch and asymptotic giant branch for all three values (also included for added perspective are the Sun's radius, luminosity and effective temperature changes). The resulting Solar white dwarf mass for each of these values is 0.534M⊙ (η = 0.2), 0.519M⊙ (η = 0.5) and 0.511M⊙ (η = 0.8).

ISOLATING DIFFERENT EFFECTS
The motion of any planet orbiting a star at a distance of at least hundreds or thousands of au could be affected by a variety of effects. These include post-main-sequence mass loss, Galactic tides, and flybys from passing stars (see fig.  2 of Veras 2016). Combining these effects with the gravitational interactions from a set of giant planets like Jupiter, Saturn, Uranus and Neptune is non-trivial, especially during close encounters. I am interested in the consequences of these effects on the orbits of giant planets, and in the below subsections, isolate these physical processes.
Consider first a single planet with an arbitrary substellar mass in isolation on an arbitrary orbit whose pericentre lies beyond 10 2 au from the Sun-like star.

Stellar flybys
Stars pass by the Sun on a regular basis, and occasionally enter the Sun's gravitational sphere of influence. Sometimes these intrusions are deep. In fact, a Sun-like star can expect to encounter another star at a distance of a few hundred au over its main sequence lifetime (Zakamska & Tremaine 2004;Veras & Moeckel 2012). The consequences for a distant planet may be any outcome (orbit perturbation, collision, ejection, engulfment) depending on the velocity, direction and mass of the flyby star. On Gyr-timescales, objects with semimajor axes of the order 10 4 au have a significant chance of being stripped away by passing stars (Li & Adams 2016). Figure 2. The final semimajor axis and eccentricity of an isolated planet orbiting the Solar white dwarf, assuming no other Solar system planets exist and the only active force was and is gravity from the Sun. The left panels show the case for a main sequence, or, "initial" e pl = 0.5 value and the right panels for an initial e pl = 0.8. The adiabatic limits, given by the straight horizontal lines, are reliable predictors typically only for a pl ≪ 500 au (specifically see Eq. 15). These plots illustrate the non-monotonicity of the planet's evolution as a function of initial a pl .
These flybys will occur during all phases of stellar evolution. However, predicting where and when these flybys will occur is generally not possible because the stellar kinematic memory in the Galaxy is lost well-before 1 Gyr in the future (Nakajima, Morino & Fukagawa 2010). In fact, the significant distance errors on the closest known Solar flyby, which occurred about 70 +10 −15 × 10 3 yr ago, are 52 +23 −14 × 10 3 au (Mamajek et al. 2015). Consequently, accurately predicting the incidence of flybys on the giant branch or white dwarf phases of Sun-like stars is not possible, and hence I do not consider flybys any further.

Galactic tides
The effect of Galactic tides is more predictable. In the Solar neighbourhood 2 , the vertical component of these tides is an order of magnitude stronger than the planar components (Heisler & Tremaine 1986;Matese & Whitman 1989, 1992Matese et al. 1995;Breiter, Dybczynski & Elipe 1996;Brasser 2001;Breiter & Ratajczak 2005), greatly simplifying the equations of motion (see the contrast within table 1 of Veras et al. 2014b). For most orbits except those near the escape boundary, Galactic tides change a planet's eccentricity e, inclination i, argument of periastron ω and longitude of ascending node Ω -but not its semimajor axis a (Brasser 2001;Fouchard 2004;Fouchard et al. 2006;Veras & Evans 2013a). However, near the escape boundary (within a factor of 3-10 in orbital period), all of the orbital elements vary ( fig. 3 of Veras et al. 2014b).
The strength of these tidal changes depends primarily on the inclination of the planetary orbit with respect to the Galactic plane, with the maximum effect achieved for polar orbits. In this extreme case, the planet's eccentricity may reach unity and consequently escape the system. In the other extreme, coplanar, case the minimum eccentricity increase and decrease factors generated by Galactic tides in the Solar neighbourhood are roughly 1 ± 10 −8 (a/au) 3/2 (eq. 44 of Veras & Evans 2013a). This value corresponds to an eccentricity change of about 10 −5 for a = 10 2 au and 10 −2 for a = 10 5 au.
The Solar system, however, does not fit either extreme. The ecliptic is inclined to the Galactic plane by about 60 • , and the orbits of the outer planets therefore have a relative inclination of about 57.5 • − 62.5 • with respect to the Galactic plane. This orientation dictates that an object like Sedna (with a semimajor axis of about 500 au and eccentric-  Fig. 2, except here I did not vary the main sequence, or initial, values of a pl and e pl (left panels: a pl = 10 3 au and e pl = 0.8 and right panels: a pl = 10 4 au and e pl = 0.8) while instead varying ω pl (argument of pericentre) and Π pl (mean anomaly). The purple spades and orange triangles respectively demonstrate the effect of varying initial ω pl and Π pl . The difference is dramatic, and helps exhibit the sensitive dependence of mass loss on planet position along the orbit. This sensitivity is further emphasized by the pink stars, which result from the initial values of ω pl and Π pl being equal to each other and to the x-axes values. Because the pink stars and orange triangles become more coincident as a pl is increased (from the left panels to the right panels), non-adiabaticity directly correlates with sensitivity on Π pl , a known theoretical result (Veras et al. 2011). The larger value of a pl also allows for some systems to become unbound. Overall then, the inclination of a distant planet with respect to the Galactic plane may importantly be affected by Galactic tides, with a potential subsequent effect on the stability of the remaining planets in the system.

Stellar mass loss
Stellar mass loss will likely represent the key dynamical driver of the future Solar system, as well as for most exoplanetary systems. The motion of a distant planet will be determined by the solution of the following set of differential equations (Omarov 1962;Hadjidemetriou 1963;Veras et al. 2011) where f pl and E pl refer to the true anomaly and eccentric anomaly of the planet. I have denoted the star's (changing) mass as M⋆ (M⊙ is just the current value of the mass of a Sun-like star) and the planet's mass as M pl . The evolution of the anomalies is where the mean motion n pl = G 1/2 (M⋆ + M pl ) 1/2 a −3/2 pl , and generally,Ṁ⋆ < 0. Because the planet accretes some of the stellar ejecta (Hadjidemetriou 1963 and section 4.2.2. of Veras 2016)Ṁ pl > 0. However, for Solar system analogues, M pl /Ṁ⋆ ≪ 1.
Equations (5)- (13) are not known to be solvable analytically, and assume that the Sun-like star loses mass isotropically. Veras, Hadjidemtriou & Tout (2013b) showed that this approximation is excellent because latitudinal mass-loss variations anisotropically affect the planet's motion only if the mass loss is asymmetric about the stellar equator. Unrealistic longitudinal variations (of over 0.1%) must be maintained for long periods (over 1 Myr) in order to produce semimajor axis variations of order 0.1 au and eccentricity variations of order 0.01. See their fig. 2 for rough magnitudes, although their choice of reference frame alters the interpretation of the orbital elements (Dosopoulou & Kalogera 2016a,b).
The smaller the value of a pl , the better that the semimajor axis evolution can be approximated in the adiabatic limit as The The "adiabatic regime" is a useful characterisation for the regime which typically occurs when Ψ ≪ 0.1 − 1. Consequently, a planet located at 10 3 au should be approaching the edge of the adiabatic limit. In order to quantify this sentiment, I have run simulations of a planet at different semimajor axes and two different eccentricities, for each of the three η = 0.2, 0.5, 0.8 stellar models. Fig. 2 illustrates the results.
The figure demonstrates how non-adiabaticity increases with initial (or, main sequence values of) a pl . For η = 0.8, which has the slowest rate of mass loss of the three tracks that I sampled despite releasing the most mass overall (see Fig. 1), the adiabat is followed most closely (within a few per cent) for distances including and under 10 3 au. In all cases, the departure from adiabaticity is non-monotonic with increasing a pl , and one system in the left panels and one in the right become unbound. Such behaviour is explained further in Veras et al. (2011) and Veras & Wyatt (2012).
The systems sampled in Fig. 2 fix the planet's initial argument of pericentre and mean anomaly (≡ Π pl ). However, the value of these orbital elements is potentially of crucial importance. Hence, in Fig. 3, I fixed the initial values of a pl and e pl and instead varied orbital angles. The almost straight lines of purple spades in all of the plots reveal that (i) this choice produces a notable deviation from the adiabatic limit, and (ii) for the particular initial fixed Π pl value of 0 • , changing the initial value of ω pl has little effect on the overall outcome. Conversely, the scattered orange triangles and pink stars demonstrate a strong dependence on the planet's location along its orbit during the Solar giant branch phase. As shown by the right panels, if the initial a pl is large enough, then values of Π pl may be found which would lead to escape.
In both Figs 2 and 3 the mass of the planet is not reported. The reason is that its mass is largely insensitive to the process of mass loss. In particular, because M pl always appears in summation with M⋆ in equations (5)-(15) for all M pl ≪ M⊙, the influence of planet mass on its resulting motion should enter at the 0.1% level for Jupiter-mass planets.

Galactic tides and stellar mass loss together
The interplay between Galactic tides and mass loss can be partitioned because each acts on a different timescale. Veras et al. (2014b) concluded that in the Solar neighbourhood the tidal timescale is orders of magnitude longer than the mass loss timescale, allowing one to decouple the equations of motion in each phase. In other words, one can use the tidal equations along the main sequence and white dwarf phases, and the mass loss equations alone along the giant branch phases.
The consequences for long-term stability are highlighted in Bonsor & Veras (2015). Consider a wide-orbit (> 10 3 au) planet whose pericentre during the main sequence is high enough to not perturb an inner planetary system. After mass loss, the orbit expands enough to cause Galactic tides to excite the distant planet's eccentricity to the extent that the inner planetary system is now affected (see, in particular, their fig. 2). The eccentricity of the distant planet is likely to further change if during mass loss the planet is in the non-adiabatic regime (see Figs 2 and 3).

Planet-planet scattering
By themselves, Jupiter, Saturn, Uranus and Neptune are not assumed to scatter off one another during the Sun's evolution. 3 In order to confirm this notion, I have performed three simulations (η = 0.2, 0.5, 0.8) with the four giant planets only plus Galactic tides. I ran these simulations throughout the giant branch phases and for 10 Gyr on the white dwarf phase. In no case did the giant planets commence strong scattering.
Now consider the inclusion of a distant planet. Although unimportant for Galactic tides or stellar mass loss, the mass of the planet will crucially determine whether Jupiter, Saturn, Uranus and Neptune will be perturbed enough to cause large-scale instability in the Solar system analogue. A distant planet which is comparable in mass to any of the known gas giants will likely create such an instability during a close encounter. What happens if the planet is a super-Earth-mass planet is less clear. Before exploring these possibilities, I describe in detail the computational method used in this work and lay out the initial conditions for my simulations.

Numerical codes
In order to model the future evolution of a Sun-like star concurrently with the dynamical evolution of Jupiter, Saturn, Uranus, Neptune analogues and a distant planet, I utilized a facility which has combined the stellar evolution code SSE (Hurley et al. 2000) with a heavily modified version of the Bulirsch-Stoer integrator in the planetary dynamics code Mercury (Chambers 1999). This combination has been shown to yield converging errors as the accuracy parameter decreases (e.g. fig. 1 of , and the resulting code has been used in several previous studies (e.g. Veras et al. 2016a). I adopt the version of the code that includes effects from Galactic tides, and a realistic Hill ellipsoid, in order to correctly track escape (Veras & Evans 2013a;Veras et al. 2014b). For these computations, I assumed that the Sun-like star resides at exactly 8 kpc away from the Galactic Centre.
The Hill ellipsoid is stretched towards the centre of the Galaxy, and compressed along the other two axes. At a distance of 8 kpc, the semi-axes of this ellipsoid are (1.49, 1.92, 2.89) × 10 5 au. Any planet which the code found to reside outside of this ellipsoid at any timestep was considered to have escaped from the planetary system. None of the M pl values that I chose were large enough to alter these semi-axes.
The SSE code produces a single evolutionary track given a set of zero-age-main-sequence values of M⋆, Z and η. As previously mentioned, η is a key unknown (Fig. 1). Rather than adopting a single value, I chose three (0.2, 0.5, 0.8) for the reasons outlined in Section 2.

Simulation duration
The timespan over which I ran the simulations was limited by computational resources because of the number of orbits that needed to be accurately modelled. This value is dependent on the parent star's mass and size of the orbits. Hence, evolution during the white dwarf phase can proceed at a rate tens of times faster ( fig. 8 of ) than on the main sequence.
This computational slowness along the main sequence (together with the difficulty of self-consistently modelling multiple phases of stellar evolution) partly illustrates why many post-main-sequence exoplanetary evolution studies have restricted their simulations to a timespan of under about 1 Gyr (Debes & Sigurdsson 2002;Bonsor et al. 2011;Debes et al. 2012;Kratter & Perets 2012;Frewen & Hansen 2014). Although multi-Gyr simulations were later performed by  and Mustill et al. (2014), they had to restrict the main sequence stellar masses which they sampled to M⋆ 3M⊙ because of the correspondingly short main sequence lifetimes ( 500 Myr). Nevertheless, Veras & Gänsicke (2015) and Veras et al. (2016a) achieved full-lifetime (14 Gyr) simulations by adopting stars with main sequence masses between 1.5−2.5M⊙, and hence main sequence lifetimes of 600 Myr to 3 Gyr.
However, the main sequence lifetime for a 1.0M⊙ star is much longer, and approaches 11 Gyr, meaning that the Sun will remain on the main sequence for another 6.5 Gyr or so. This timespan was too prohibitive for my numerical integrations. Consequently, I began all my simulations just before the red giant branch phase. Nevertheless, I ran these simulations for a total of 11.6 Gyr, which included roughly 1.5 Gyr on the giant branch phases, and 10 Gyr on the white dwarf phase. In order to maximize resources, the occurrence of instability (defined as ejection or collision) terminated a simulation. I ran the simulations which generated Figs 2 and 3 for just 1.55 Gyr, because only the giant branch phases needed to be sampled.

Planet parameters
Regarding the planets, as mentioned previously, Jupiter, Saturn, Uranus and Neptune are expected to survive until the end of the main sequence in approximately their current orbits. Their exact positions 6.5 Gyr from now are not known (e.g. Zeebe 2015), even if the errors in current measurements were as small as a Planck length, because the outer Solar system itself is chaotic with a Lyapunov time of order 10 Myr. Hence, I simply adopted as initial conditions the values of a, e and f of the planets at the Julian Date 2451000.5 that are provided within the default version of Mercury.
However, I treated their values of i, Ω, and ω differently. The code considers the Galactic Plane to exist at i = 0 • , and the relative inclination of each planet to that plane may be important when coupled with the distant planet's inclination and the effect of Galactic tides. Because the ecliptic is approximately inclined at an angle of 60 • to the planets, I reproduced faithful analogues by tilting the planetary systems so that the orbital planes are approximately centred about this value. I did so by applying the rotation matrix to the Cartesian elements (positions and velocities) of the four planets, and then transforming back into orbital element space (this transformation does not affect a, e nor f ). The result was that the inclination values of the four planets ranged from 58.8 • to 60.3 • . Because the Sun-like stars were treated as point masses, they were not affected by, nor affected, this rotation.

MONTE-CARLO SIMULATIONS
Having laid out the background for my study, I am now ready to present my Monte-Carlo simulations.

Initial conditions
In order to establish initial conditions, I consider each parameter in turn. In all cases, when I use the term "initial", that refers to the starting point of my simulations (10941.0429 Myr from the zero-age-main-sequence) but also more generally to main sequence values. Time is measured from the current epoch, which was assumed to occur exactly 4.6 Gyr after the zero-age-main-sequence. My adopted ranges for the distant planet's mass, semimajor axis and inclination (below) are all larger than those assumed by recent studies attempting to constrain the properties of a (potential) Planet Nine in the Solar system. These expanded ranges allow me to help identify when instability would occur for planetary systems that generally resemble the Solar system. Naturally if Planet Nine is discovered and its mass, orbit and position are pinpointed, then a dedicated post-mainsequence study which adopts those parameters would be of interest.
• I chose to sample M pl randomly from a logarithmic distribution with range 10 1 − 10 4 M⊕. I assumed that lower masses would not likely be able to significantly perturb the giant planets even during a close encounter (the masses of Uranus and Neptune are about 14-17M⊕). Higher masses would represent stars rather than planets; in fact the typically-assumed planet-brown dwarf boundary of 13 Jupiter masses is approximately equal to 4100M⊕.
• I sampled a pl (initial) randomly from a logarithmic distribution with range 5 × 10 2 − 5 × 10 4 au, and the planet's initial orbital pericentre ≡ q pl (initial) = a pl (initial) × [1 − e pl (initial)] from a uniform distribution with a range of 100 to 400 au. These choices allowed me to sample the entire system from 10 2 au out to the edge of the Hill Ellipsoid (at ≈ 10 5 au).
• I sampled i pl (initial) randomly from a uniform distribution with range 40 • −80 • , which is centred around the approximately 60 • tilt that I have imposed for the giant planets (recall that i = 0 • corresponds to the Galactic plane). My range is somewhat arbitrary: I recognize that a distant planet may have any inclination, particularly if it is a captured object, but for computational purposes restricted the range to an order of magnitude greater than the mutual in-clinations amongst the giant planets. Sampling non-zero intermediate values of inclination is anyway particularly valuable in order to probe the effect of Galactic tides.
• The orbital angles ω pl (initial) and Ω pl (initial), and Π pl (initial) were sampled randomly from a uniform distribution over their entire ranges.

Results
I present my results in two stages. The first is an analysis of specific systems, and the second is an analysis of the ensemble.

Specific systems
I have analyzed the outcomes for each of nearly 300 simulations on an individual basis, and decided to present here a flavour of the wide variety of outcomes. In the following figures, I consistently use the same horizontal and vertical ranges, even for unstable simulations that were terminated.
• Figure 4 (stable) One of the simplest cases was evolved with a η = 0.5 Sun-like star and features a distant planet on a short-enough orbit that neither Galactic tides nor non-adiabatic mass loss ever become factors. Consequently, the semimajor axis and eccentricity changes were predictable, and all planets were far away enough from each other such that they did not suffer close encounters. The eccentricities of the analogues of the four known giant planets mirrored their main sequence values and variations, whereas the distant planet's eccentricity was slightly noticeably altered: the amplitude of e pl increased just after mass loss, and Galactic tides did cause a barely perceptible secular change over 10 Gyr of white dwarf evolution.
The Cartesian cross-sections (bottom panels) reveal how mass loss changes planetary orbits in space. For each planet, there are two distinct tori of points, corresponding to main sequence and white dwarf values. For the distant planet, a sparse intermediate ring of points can be seen, which corresponds to a transition state during the giant branch phase.
• Figure 5 (unstable) Here the distant planet's initial semimajor axis of about 2030 au was (i) small enough to achieve adiabatic mass loss, (ii) small enough to be unaffected by Galactic tides on the main sequence and giant branch phases, (iii) and large enough to be affected by Galactic tides non-negligibly on the white dwarf phase. Regarding the first point, this semimajor axis lies within the Solar system adiabaticity transition region of 10 3 − 10 4 au, but remained adiabatic because here η = 0.8 (see fig. 2 of Veras & Wyatt 2012). The modulation of e pl due to Galactic tides was high enough to be seen on the middle plot. This increase of eccentricity (by about a tenth from 0.85) decreased q pl (bottom plot) gradually, until close encounters pumped up the eccentricities of the Uranus and Neptune analogues, eventually ejecting the Neptune-like planet and stopping the simulation.
The value of q pl became much lower than the orbital pericentres of the other planets before ejection. Indicated on the right axis of the bottom plot are Rmax and R Roche , which represent the maximum radius that the star attained during the giant branch phase, and the Roche, or disruption radius, Figure 5. Tidally-induced instability: Here the distant planet (M pl ≈ 109M ⊕ ) has a large enough initial semimajor axis a pl (initial) ≈ 2030 au for Galactic tides to have a noticeable effect during the white dwarf phase. The tides create an initial increase in the already high value of e pl (initial) ≈ 0.85, triggering ejection of the analogue of Neptune as q pl approaches the location of the other four planets (bottom plot). Eventually, q pl < Rmax, causing the distant planet to sweep through any remaining debris in the inner system, which may pollute the eventual white dwarf once inside its disruption radius R Roche . Figure 6. Stability despite tidal excitation: A super-Earth-sized distant planet (M pl ≈ 12M ⊕ ) leaves the main sequence in an orbital configuration which causes non-adiabatic evolution, yielding an order of magnitude increase in semimajor axis. Consequently, Galactic tides become strong, producing oscillations in eccentricity on the order of tenths. The non-uniformity of the oscillations and jumps in a pl indicate effects from perturbations with the other four planets, whose eccentricities remain under 0.2. The repeated deep radial incursions of the distant planet near or within Rmax might provide pathways to pollute the eventual white dwarf. Figure 7. The Uranus-like planet is kidnapped, then ejected: In this unusual case, mass loss triggers scattering between the Uranus and Neptune analogues and the distant planet (M pl ≈ 806M ⊕ ), which harbours a cometary-like orbit with a pl (initial) ≈ 17000 au and q pl (initial) ≈ 222 au. Because Uranusanalogue's semimajor axis is increased by two orders of magnitude, both the Uranus-like planet and the distant planet are then strongly affected by Galactic tides (middle panel). Eventually, soon after the Uranus analogue approaches the orbits of the Neptune, Saturn and Jupiter analogues at a time of about 13 Gyr (bottom panel), the Uranus-like planet is ejected. A caveat to this figure is that stellar tidal effects on the Uranus analogue (not modelled) may have played a role during strong scattering on the star's red giant branch phase. Figure 8. The Neptune-like planet is kidnapped, but the system remains stable: For this rare case, the tip of the asymptotic giant branch coincides with a pericentre passage of the distant planet (bottom panel). During this time, the Neptune and Uranus analogues are perturbed enough by the slightly larger distant planet (M pl ≈ 30M ⊕ ) to force Neptune-analogue's semimajor axis to increase by two orders of magnitude and increase Uranus analogue's eccentricity to about 0.4. The planet's already non-adiabatic value of [a pl (initial) ≈ 27000 au] combined with this scattering event results in a negligible change in a pl . Subsequently, both the Neptune-like and the distant planet are affected by Galactic tides and sweep through the Solar system analogue, achieving q pl < Rmax. However, these actions do not disturb Uranus, Saturn nor Jupiter, and the system retains all planets until the end of the simulation. of the white dwarf. For all plots, I computed this radius assuming a Jupiter-mass planet. In this figure, q pl < Rmax. Consequently, the distant planet swept through an area of space that had been untouched since the end of the giant branch phase. Hence, extant debris from the inner system (such as e.g. planets like Mars, and remnants of asteroid belt-like breakups; Veras et al. 2014a), could be perturbed towards the white dwarf and pollute it. This same behaviour was observed in Veras & Gänsicke (2015) and Veras et al. (2016a).
• Figure 6 (stable) Here the distant planet's initial semimajor axis of 2230 au is very similar to that in Fig.  5, except in this instance the mass loss evolution was nonadiabatic. The reason partly is, for this simulation, η = 0.5 (see fig. 2 of Veras & Wyatt 2012). The resulting large jump in the value of a pl , by nearly a factor of 10, created a situation where Galactic tides played a large role during white dwarf evolution.
Whereas Galactic tides slightly decreased e pl on the main sequence and white dwarf phases, the large jump in a pl created discernable oscillations during the white dwarf phase. If Jupiter, Saturn, Uranus and Neptune analogues were not present, then these oscillations would have had a constant amplitude and frequency. However, mutual interactions triggered semimajor axis "jumps" (upper plot) plus a corresponding change in shape of the eccentricity oscillations. Further, over the first 7 Gyr of white dwarf evolution, Neptune-analogue's eccentricity increased by over an order of magnitude, to about 0.1-0.2. Nevertheless, over the duration of the simulation, the system remained stable.
The bottom plot illustrates deep radial incursions (within or close to Rmax) of the distant planet towards the white dwarf: one at 10 Gyr, and then four more from 14-18 Gyr. These illustrate how a dynamically-active inner planetary system environment can be triggered at any white dwarf age, potentially explaining white dwarf pollution at different epochs (Bonsor & Veras 2015;Hamers & Portegies Zwart 2016;Petrovich & Muñoz 2016).
• Figure 7 (unstable) A highly eccentric ( 0.9) distant planet with a pl > 10 4 au is likely to generate interesting dynamical behaviour because it is sure to evolve nonadiabatically from mass loss and always is affected significantly by Galactic tides. For this system, these attributes provided an environment which was conducive to strong scattering events, allowing the analogue of Uranus to be "captured" into a 10 4 au orbit. This orbit existed for an appreciable 5 Gyr -longer than the Solar system's current age -before the analogue of Uranus was ejected.
Some other noteworthy dynamical aspects of this η = 0.5 system are: (i) the change in frequency of the distant planet's eccentricity oscillations from the giant branch to white dwarf phases (due to semimajor axis increase), (ii) the Uranus-like planet achieving an orbit with a pericentre of nearly 1 au at the end of the red giant branch phase (but before the asymptotic giant branch phase); this fact places a significant caveat on this figure because tidal interactions with the star (not modelled) could have retarded the planet's evolution during this stage (Schröder & Connon Smith 2008), (iii) Neptune, Jupiter and Saturn analogues settling into orbits whose secular eccentricity evolution was well-behaved, compared to that from any other figures.
• Figure  8 (stable) In this figure, Figure 9. The outcomes from all simulations. Stable simulations are indicated in green symbols, and unstable ones are in red. The different shapes refer to different Solar models; these plots illustrate that the results are independent of η. Both stable and unstable simulations populate all areas of all plots, except for the blue triangle in the upper-left plot; this triangle contains all stable simulations.
Orange rectangles indicate regions where stability is more likely to occur than not. The outcomes are largely independent of q pl (initial) and i pl (initial).
a pl (initial) ≈ 27000 au, a value 59 per cent greater than the corresponding value in Fig. 7, and the initial eccentricities in both figures coincidentally differed by just 0.05 per cent. Despite these values, the amplitude of the eccentricity oscillations due to Galactic tides in this figure was not as high. The reason is because of the distant planet's other orbital parameters: i pl , ω pl , and Ω pl all played a significant, non-trivial role. Regardless, the pericentre passages were deep enough to coincide with the orbits of the four giant planet analogues. One of these passages occurred close to the tip of the asymptotic giant branch (whereas in Fig. 7 a passage occurred close to the tip of the red giant branch). The resulting scattering event helped contribute to a pl 's net negligible change at the giant branch / white dwarf boundary, and triggered a semimajor axis increase for the analogue of Neptune to almost 10 4 au after a few hundred Myr. Neptune-analogue's eccentricity then became strongly affected by Galactic tides. Subsequently, both the Neptune-like planet and the distant planet experienced repeated radial incursions -at different rates -providing a dynamically active and ever-changing environment potentially conducive to pollution. The Uranus, Saturn and Jupiter analogues were not bothered by this activity, and the system remained stable for the duration of the simulation.

System ensemble
The primary goal of this study is to determine the stability boundary: for what mass and orbital parameters of a distant planet will instability occur. Amongst the different potential combinations of variables to plot, I found that the most revealing combinations are those seen in Fig. 9. The figure reveals the following points.
• Any potentially sharp boundary is limited to the blue triangle in the upper-left plot; in this triangle every system sampled was stable. Even so, just a few tens of simulations within that triangle were performed. Otherwise, unstable and stable simulations appear on all regions of all plots, showcasing how sensitive the simulation outcome is to all orbital parameters, including Ω pl , ω pl and Π pl .
The blue triangle makes sense: Although the notion that smaller values of a pl would yield greater stability might seem counterintuitive, in fact such values guarantee that the mass loss evolution is adiabatic and that Galactic tides play no significant role in the dynamics. Also, the greater the value of M pl , the more likely that the distant planet could scatter off of the other planets and create instability.
• Mass aside, stability is more generally likely to occur than not when a pl (initial) is within a few thousand au, as indicated by the orange regions. The reason is the same as above.
• There is no discernable stability dependence on η (for values between 0.2 and 0.8).
• There is no discernable stability dependence on q pl (for values between 100 and 400 au).
• There is no discernable stability dependence on i pl (for values within 20 • below and above the ecliptic).
Other trends from the ensemble of plots include details of the instability: 97 per cent of all unstable simulations featured ejections. The remainder were engulfments (defined to occur when a planet intersected the star's Roche radius). In no case did the Jupiter-like planet become unstable. The analogues of Saturn, Uranus and Neptune and the distant planet were either ejected or engulfed 8, 53, 24 and 15 per cent of the time, respectively. This result makes sense given that in almost every simulation, the analogue of Uranus was the least massive planet. Veras et al. (2016a) showed that scattering amongst unequal-mass planets across all phases of stellar evolution will preferentially eject the least massive planet (just as along the main sequence alone).

CONCLUSIONS
I demonstrated that a distant planet with an orbital pericentre under 400 au could pose a serious danger to the stability of Solar system analogues during a Sun-like star's giant branch and white dwarf phases. This statement holds true for a distant planet which is at least as massive as Jupiter and harbours a semimajor axis beyond about 300 au, or for a super-Earth when its semimajor axis exceeds about 3000 au. The driver for the instability is a combination of Galactic tides and stellar mass loss, which together or separately may induce close encounters amongst the five planets, with the distant planet always representing the trigger.
These results have implications for both the Solar system and for extrasolar systems. The existence of a trans-Neptunian planet could eventually eliminate (likely through ejection) at least one of the giant planets (most likely Uranus, then Neptune) and rearrange the others, but only if this "extra" planet is massive enough, distant enough and its orbit is appropriately oriented with respect to the Galactic plane and the existing giants. This planet may represent the purported Planet Nine, a hidden other planet, or a planet that will be captured later during the Sun's main sequence evolution. Because the Sun will become a white dwarf and contain strewn-about debris from a destroyed asteroid belt (Veras et al. 2014a) plus a charred Mars and perhaps some liberated moons (Payne et al. 2016a,b), an ample reservoir of extant material in the inner Solar system would be available to be perturbed into the Solar white dwarf, "polluting" it.
The consequences for other planetary systems are profound. Multiple planets beyond about 5 au (such as analogues of Jupiter, Saturn, Uranus and Neptune) may be common, but are so far unfortunately effectively hidden from detection by Doppler radial velocity and transit photometry techniques, the two most successful planet-finding techniques. If more distant, trans-Neptunian-like planets are also common, then the ingredients may exist to regularly generate instability and a frequently-changing dynamical environment during white dwarf phases of evolution. Such movement could provide a natural way to perturb inner system debris into white dwarfs at a variety of white dwarf ages and rates, helping to explain current observations.