# Simulations of two-planet systems through all phases of stellar evolution: implications for the instability boundary and white dwarf pollution

Dimitri Veras, Alexander J. Mustill, Amy Bonsor, Mark C. Wyatt

► To cite this version:

Dimitri Veras, Alexander J. Mustill, Amy Bonsor, Mark C. Wyatt. Simulations of two-planet systems through all phases of stellar evolution: implications for the instability boundary and white dwarf pollution. Monthly Notices of the Royal Astronomical Society, 2013, 431, pp.1686-1708. 10.1093/mn-ras/stt289 . insu-03617074

HAL Id: insu-03617074

https://insu.hal.science/insu-03617074

Submitted on 23 Mar 2022

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Simulations of two-planet systems through all phases of stellar evolution: implications for the instability boundary and white dwarf pollution

Dimitri Veras,<fnref n="1" />★ Alexander J. Mustill,<fnref n="2" /> Amy Bonsor<fnref n="3" /> and Mark C. Wyatt<fnref n="1" />

Accepted 2013 February 13. Received 2013 January 15; in original form 2012 November 20

## ABSTRACT

Exoplanets have been observed at many stages of their host star's life, including the main-sequence (MS), subgiant and red giant branch stages. Also, polluted white dwarfs (WDs) likely represent dynamically active systems at late times. Here, we perform three-body simulations which include realistic post-MS stellar mass-loss and span the entire lifetime of exosystems with two massive planets, from the endpoint of formation to several Gyr into the WD phase of the host star. We find that both MS and WD systems experience ejections and star-planet collisions (Lagrange instability) even if the planet-planet separation well-exceeds the analytical orbit-crossing (Hill instability) boundary. Consequently, MS-stable planets do not need to be closely packed to experience instability during the WD phase. This instability may pollute the WD directly through collisions, or, more likely, indirectly through increased scattering of smaller bodies such as asteroids or comets. Our simulations show that this instability occurs predominately between tens of Myr to a few Gyr of WD cooling.

**Key words:** planets and satellites: dynamical evolution and stability – planet-star interactions – stars: AGB and post-AGB – stars: evolution – white dwarfs.

## 1 INTRODUCTION

A planet's life may be split into four distinct stages: (1) formation and concurrent dynamical excitation, (2) main-sequence (MS) evolution, (3) evolution during post-MS stellar phase changes and (4) white dwarf (WD) evolution. The first stage generally lasts no longer than 0.1 per cent of the entire MS lifetime. The second stage is relatively dynamically quiescent, with only occasional but often important scattering interactions. In the third stage, the planet is subject to dynamical changes due to the star's violent actions as it becomes a giant. In the final stage, the star has become a WD, and the planet again enters and remains in a phase of relative dynamical quiescence occasionally punctuated by scattering interactions or external forcing. This general picture, which does not include possibilities such as the capture of free-floating planets, planetary destruction due to supernovae, or multiple host stars, describes the life cycle of the vast majority of known exoplanets.

The volume of planetary literature investigating the first two stages dwarfs the literature describing the final two stages, despite the fact that the Universe is already over 13.5 Gyr old (Jarosik et al. 2011) and that the Milky Way contains about 10⁹ WDs (Binney & Tremaine 2008, p. 2–3 and Holberg et al. 2008). Further, these final

two stages are becoming increasingly relevant given the suggestions or discoveries of exoplanets in post-MS systems (Wolszczan & Frail 1992; Wolszczan 1994; Sigurdsson et al. 2003; Silvotti et al. 2007; Mullally et al. 2008, 2009; Geier et al. 2009; Lee et al. 2009, 2012; Setiawan et al. 2010; Wickramasinghe et al. 2010; Charpinet et al. 2011; Adamów et al. 2012; Farihi et al. 2012c; Lee, Han & Park 2013; Sato et al. 2013).

Explorations of exosystem evolution in the third stage include one-planet studies (Villaver & Livio 2007, 2009; Veras et al. 2011; Kratter & Perets 2012; Mustill & Villaver 2012; Spiegel & Madhusudhan 2012; Veras & Tout 2012; Adams, Anderson & Bloch, in preparation), just a few of dedicated multiple-planet studies (Debes & Sigurdsson 2002; Portegies Zwart 2013; Voyatzis et al. 2013) and studies focusing on the evolution of comets (Alcock, Fristrom & Siegelman 1986; Parrott & Alcock 1998). Furthermore, Bonsor & Wyatt (2010) consider the effect of post-MS evolution on debris discs. Motivated by observations of metal-polluted WDs, Bonsor, Mustill & Wyatt (2011), Bonsor, Augereau & Thebault (2012) and Debes, Walsh & Stark (2012) model the interplay between a planet and a belt of smaller material amidst stellar mass-loss.

Here, we self-consistently simulate the second, third and fourth stages together. We combine stellar evolution with planetary gravitational scattering amongst multiple massive planets, and extend the work of Debes & Sigurdsson (2002) by considering full-lifetime simulations with realistic mass-loss prescriptions at each post-MS

<footnote n="1">¹Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="2">²Universidad Autónoma de Madrid, Departamento de Física Teórica C-XI, 28049 Madrid, Spain</footnote>

<footnote n="3">³UJF-Grenoble 1 / CNRS-INSU, Institut de Planétologie et d'Astrophysique de Grenoble (IPAG), UMR 5274, BP 53, F-38041 Grenoble cedex 9, France</footnote>

<footnote>*E-mail: veras@ast.cam.ac.uk</footnote>
phase. Following the evolution over the whole stellar lifetime means that we can be sure that systems whose stability is investigated on the giant and WD stages will have survived the long MS evolution. Through these integrations, we can determine what types of planetary architectures might be expected in exoplanet-hosting WD systems, and could allow us to extrapolate backwards in time from observed WD systems. We restrict our explorations to two-planet systems in this initial study given the vast phase space to explore; three-planet simulations will be presented in a follow-up paper. First, we briefly summarize our knowledge of planetary instability for one- and two-planet systems during the MS (Section 1.1) and post-MS (Section 1.2).

## 1.1 Instability in MS planetary systems

Dynamical instability in planetary systems is often said to occur when a planet suffers a close encounter with the star or another planet, or is ejected from the system. Occasionally, investigators use stricter definitions of instability, such as when the semimajor axis or eccentricity variation of a planet exceeds a certain per cent of its nominal value. Additionally, a wide body of literature has arisen characterizing chaotic orbits as a precursor to instability; Darriba et al. (2012) and references therein summarize many of these techniques.

### 1.1.1 One-planet instability

One planet orbiting an MS star will typically remain stable throughout the star's MS lifetime in the absence of external forces. Exceptions may include planets which are close enough to their parent stars to be tidally disrupted (e.g. Gu, Lin & Bodenheimer 2003) and possibly evaporated (e.g. Guillot et al. 1996). In the opposite extreme, a planet which is far enough away from its parent star may be ejected due to external forces such as passing stars (Zakamska & Tremaine 2004; Veras & Moeckel 2012) or achieve a high enough eccentricity through Galactic tides to cause a collision with the star (Veras & Evans 2013a,b).

### 1.1.2 Two-planet instability

In addition to tidal interactions and external forces, the mutual perturbations between two planets may also create instability. Partially motivated by tractable analytical solutions to the general three-body problem, the source of this instability has been studied extensively. If the orbits of two planets are guaranteed to never overlap (precluding a collision between both planets), then they are said to be ‘Hill stable’. Gladman (1993) pioneered the analytic use of Hill stability for planetary systems in specific cases and has motivated many subsequent analyses, as recently summarized by, e.g., Donnison (2009, 2010a,b, 2011). Hill stability does not guarantee that the outer planet remains bound to the system, nor does it prevent the inner planet from colliding with the star. If both planets remain bound and retain their ordering, and no collision with the star occurs, then the system is ‘Lagrange stable’.<fnref n="1" /> Unlike Hill stability, Lagrange stability does not benefit from a known analytical formulation, but rather empirical estimates based on numerical simulations.

The analytical Hill stability boundary is conservative. Two planets whose initial separation is less than the Hill stable distance may in fact remain stable. If the initial separation is greater than the Hill stable distance, then the planets are guaranteed to retain their ordering. In simulations of the HD 12661 and 47 Uma systems, Barnes & Greenberg (2006) found that pairs of planets close to the Hill stability boundary are not Lagrange stable, and hence are not generally stable. They tentatively suggest that the Lagrange stability boundary exceeds the Hill stable boundary by at least 21 per cent as measured by the semimajor axis ratio. Subsequent work (Barnes & Greenberg 2007) revealed how mean-motion commensurabilities can broaden the divide between the Hill and Lagrange stable boundaries; Kopparapu & Barnes (2010) demonstrated how the boundary between stable and unstable systems is not sharp. Hence, numerical validation of analytical stability estimates is crucial.

## 1.2 Instability in post-MS planetary systems

Mass-loss from a dying star can trigger planetary instability in different ways, which are outlined below. A common assumption amongst the studies which have considered instabilities in post-MS systems is isotropic stellar mass-loss. We also adopt this assumption here, as modelling non-isotropic mass-loss would significantly complicate both numerical and analytical modelling and is best left to separate, dedicated studies. One such dedicated post-MS study (Parrott & Alcock 1998) importantly observes that the speed of (effectively massless) comets near the boundary of a planetary system may be comparable to the recoil velocity of the parent star due to asymmetric mass-loss. That study suggests that anisotropic mass-loss will affect the details of planets being ejected after scattering but is unlikely to have a significant effect on the prior dynamics. Other studies modelling planetary dynamics due to non-isotropic mass-loss instead focus on jet accelerations present at the birth sites of planets (Namouni 2005, 2007, 2012).

Furthermore, in all cases we assume that the planets are orbiting a single star. Extensions to the multiple-star case (Kratter & Perets 2012; Veras & Tout 2012; Portegies Zwart 2013) are likely to be nontrivial.

### 1.2.1 One-planet instability

For decades, binary star investigations revealed that stellar mass-loss causes orbital semimajor axis expansion. Less well known is that when the mass-loss is rapid enough, the eccentricity of the companion's orbit can change as well (Omarov 1962; Hadjidemetriou 1963, 1966; Veras et al. 2011). If the eccentricity is great enough, a planetary companion may escape from the system.

In the more commonly used adiabatic limit, $de/dt = 0$ and $da/dt = -(a/\mu)(d\mu/dt)$. Here, $\mu \equiv G(M_* + M_p)/M_*$, $M_*$ and $M_p$ are the masses of the star and planet, $a$ is the planet's semimajor axis and $e$ is the planet's eccentricity. This limit holds when the mass-loss time-scale is much longer than the orbital time-scale. Adiabaticity will be broken, even briefly, if at any point a sudden burst of mass-loss causes the time-scales to become comparable (Veras & Tout 2012; Veras & Wyatt 2012). Hence, characterizing whether or not planetary evolution is adiabatic amidst mass-loss will be important for any post-MS scattering study.

Another source of instability for one-planet systems could come from tidal orbital decay and potentially direct engulfment by the rapidly expanding stellar envelope. Villaver & Livio (2007, 2009) and Villaver (2011) treat this effect in detail with additional physics such as frictional drag, planet accretion, planet evaporation and

<footnote n="1">¹ This type of stability has also been referred to as ‘Laplace stability’ (e.g. Kubala et al. 1993) and featured but remained unnamed in many papers published before the discovery of exoplanets and high-speed computing.</footnote>
spin. Mustill & Villaver (2012) model individual thermal pulses and demonstrate how they affect planetary stability. In this study, we only consider planets that are too distant to be affected by the stellar envelope expansion and so are not affected by tides, accretion or evaporation (see Section 3).

### 1.2.2 Two-planet instability

Debes & Sigurdsson (2002) considered adiabatic evolution of two-planet systems while exposed to a $1 M_\odot$ star losing half of its mass over 1000 planetary orbits. This foundational study considered circular and coplanar equal mass planets, with planet/star mass ratios ranging from $10^{-3}$ to $10^{-7}$. They discovered importantly that although adiabatic evolution causes both planets to move outward and maintain their initial semimajor axis ratio, their critical Hill separation changes.

The rate of change of the separation measured in units of Hill radii is equal to $\mu^{-2/3} d\mu/dt$. This dependence causes previously Hill stable planetary systems to become unstable, and incite gravitational scattering which could not occur on the MS. Their simulation results suggest that scattering instabilities may be more widespread during post-MS evolution than during MS evolution. Here, we investigate this claim in significant detail. For discussion on the high mass-loss non-adiabatic multi-planet case recently presented by Voyatzis et al. (2013), please see Section 6.4.

## 1.3 Paper outline

We begin in Section 2 with a description of the challenges of using N-body numerical simulations for gravitational scattering amidst mass-loss. In Section 3, we determine the regimes where engulfment and tides from stellar envelope expansion can be neglected for this study. Section 4 presents a general formulation of the Hill stability limit and shows how it changes due to stellar mass-loss. We use the results of Sections 2–4 to motivate the setup for our numerical scattering simulations. In Section 5, we perform these simulations, and report the results. We discuss the consequences in Section 6 and conclude in Section 7.

# 2 NUMERICAL CONVERGENCE

In this section, we highlight the difficulty in achieving accurate N-body simulations that model both central star mass-loss and gravitational scattering amongst multiple massive planets, and implement a solution.

## 2.1 Stellar evolution code

We utilize the SSE stellar evolution code (Hurley, Pols & Tout 2000), which adopts empirically derived algebraic formulations in order to quickly generate a stellar evolutionary track solely from a given progenitor mass and metallicity, and stellar model parameters such as the Reimers mass-loss coefficient. We use the same mass-loss prescriptions as described in section 7.1 of Hurley et al. (2000) with their Reimers mass-loss coefficient default value of 0.5. Their choice is observationally motivated by the horizontal branch morphology in Galactic globular clusters (Iben & Renzini 1983), and lies in the centre of the range recently considered by Veras & Wyatt (2012), who discuss this choice in light of an updated version of the Reimers law (Schröder & Cuntz 2005). The SSE code allows us to sample many different evolutionary tracks easily, and outputs the important parameters, $M_*(t)$ and $R_*(t)$, where $R_*$ is the radius of the star.

## 2.2 Planetary evolution code

We also use the MERCURY integration package (Chambers 1999), which specializes in modelling planetary dynamical evolution. In order to accurately model close encounters between planets and the parent star – a necessity for this study – we use the Bulirsch-Stoer (BS) integrator. This integrator features an adaptive time-step, which is determined by a tolerance parameter given at the start of the simulation. A tolerance of $10^{-12}$ is considered to be highly accurate (Jurić & Tremaine 2008). Smaller tolerance values should roughly converge to the same result; fig. 7 b of Smith & Lissauer (2009) demonstrates that in crowded five-planet systems separated by several Hill radii, tolerances of $10^{-12}, 10^{-13}, 10^{-14}, 10^{-15}$ and $10^{-16}$ will yield instability time-scales which are all within the same order of magnitude. Tolerances below $10^{-16}$ generally cannot be achieved because in those cases the accuracy requested is greater than machine precision.

## 2.3 Merging both codes

Veras et al. (2011) found that linearly interpolating SSE stellar mass output at each MERCURY time-step adequately models the dynamical evolution of a single planet amidst stellar mass-loss because the numerical simulations reproduced the analytical results. For multi-planet systems, this technique alone is inadequate. The interaction between both planets coupled with stellar mass-loss causes a failure of convergence of orbital parameters as the tolerance is decreased.

In order to improve the accuracy, we have performed an additional interpolation of the SSE stellar mass output in between each MERCURY time-step at each BS substep. The resulting finer gradation makes a crucial difference, as demonstrated by Fig. 1. The figure plots the final semimajor axis values for the inner $10^{-3} M_\odot$ planet and the outer $10^{-3} M_\odot$ planet in a system with initial semimajor axes of 10 and 30 au, and initial eccentricities of 0.0, 0.5, respectively. All initial orbital angles were set to 0°. The simulations were run for the entire evolution of the thermally pulsing asymptotic giant

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

**Figure 1.** Difference between interpolating SSE-outputted mass at every MERCURY time-step (blue circles) versus interpolating this mass within MERCURY time-steps (orange squares). 'Spliced' indicates the latter and 'Not Spliced' indicates the former. Shown are the final values of the inner planet's semimajor axis for a pair of 0.001 $M_\odot$-planets with initial semimajor axes of 10 and 30 au, and initial eccentricities of 0.0 and 0.5, respectively. The parent Solar-metallicity star was modelled to lose $\approx 6.22 M_\odot$ of its initial $\approx 7.66 M_\odot$ at the start of an AGB phase lasting almost $5 \times 10^5$ yr. The 'Analytical' line refers to the final semimajor axis of the inner planet predicted by adiabatic mass-loss. The convergence properties of the spliced BS integrator for systems with both gravitational scattering and mass-loss is a significant improvement.
branch (TPAGB) phase of a $Z = Z_⊙ = 0.02$ (Solar metallicity) star. We chose a progenitor mass of $8 M_⊙$ to model particularly violent mass-loss. Our simulations ran during the TPAGB phase only, when $M_*$ was reduced from 7.659 to $1.438 M_⊙$ in about 492 744 yr.

The plot contains two curves from the simulation output, representing final values of the semimajor axis due to SSE-outputted mass interpolation at each MERCURY time-step (blue circles; 'non-spliced'), and with an additional interpolation in-between time-steps (orange squares; 'spliced'). The third, green, curve is the analytic prediction for the final semimajor axis of the inner planet. This value can be determined because the mass-loss is adiabatic and the planets are not near a strong mean-motion commensurability (hence, their semimajor axes remain secularly unaffected).

Without the additional interpolation, the results do not appear to converge until perhaps at the machine precision limit for the BS tolerance.<fnref n="2" /> Further, the extent of the variance in the non-spliced curves may fundamentally change the endstate of the system if any more close encounters occur. Therefore, we use the spliced BS integrator throughout the rest of this work. Convergence with the spliced integrator is achieved in this case at an accuracy of ~$10^{-12}$; we are conservative and adopt the value of $10^{-13}$ for our integrations.

Another consideration is the ejecta-crossing lag time. Stellar ejecta will cause the inner planet's orbit to shift before the outer planet's orbit. In some cases, this 'lag time' between orbital shifts may produce a noticeable change in the dynamics that is missed by assuming that both planets simultaneously change their orbits. The weakness of this assumption is accentuated for widely spaced orbits and for systems which are not in the adiabatic regime. For the (adiabatic) systems studied here, however, this assumption likely produces a negligible effect,<fnref n="3" /> and hence is neglected for the remainder of this study.

## 2.4 Further adaptations

All orbital elements in this work are reported in Jacobi coordinates. Therefore, as MERCURY receives input in astrocentric coordinates, we performed the conversion. Furthermore, we had to modify the default version of MERCURY to account for a changing stellar mass in the output file xv.out so that the conversion from Cartesian output to Jacobi elements was performed correctly. Consequently, the size of xv.out nearly triples in size. Although this increase might be prohibitive for high-resolution studies of individual systems, here we are interested primarily in the final stability state of each system. Therefore, in our case, outputs at a Myr resolution are all that is required. Independent of the paucity of outputs, MERCURY does record the times of collisions to within a time-step.

## 3 TREATING THE STELLAR RADIUS

Additionally, we modified MERCURY to incorporate the stellar radius evolution profiles from SSE. Over its lifetime, a star's radius evolves

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Maximum stellar radius as a function of progenitor mass and metallicity. Colours denote different metallicities; the blue curve ($Z = 0.02$) represents Solar metallicity. Curves end when supernovae occur. This plot demonstrates a rule of thumb: generally, the number of au at the maximum radius is approximately equal to the number of initial $M_⊙$.</figcaption></figure>

non-linearly and non-monotonically. Because these variations are modest and all occur within 0.05 au during the MS, most previous planet scattering studies treated the radius as static and/or negligible. However, during the post-MS, the radius variations can be violent and extend beyond several au.

### 3.1 Expansion

Because of the potential for planetary collisions, evaporation and/or envelopment due to the expanding stellar envelope, the variations in stellar radius must be taken into account during post-MS scattering simulations. An important question is whether or not a planet, expanding its orbit due to mass-loss, can outrun an expanding stellar envelope. The answer is complicated by the fact that the time-scales for and amplitudes of mass-loss and radial increase are not completely in sync, although they often are similar. Additionally, a star's radius may decrease. WD radii are even more compact than MS radii.

In order to characterize these variations, we generated Figs 2 and 3, with SSE data. The figures characterize the maximum stellar radius for different metallicities and stellar phases, respectively. In particular, Fig. 2 suggests that the maximum stellar radius is largely independent of metallicity, and that roughly the number of au at the maximum radius is equal to the number of initial $M_⊙$. Fig. 3 illustrates that the stellar radius generally increases during post-MS phases, although for the $1 M_⊙$ case, there is an order-of-magnitude decrease after the red giant branch (RGB) stage. This decrease becomes progressively smaller as the progenitor stellar mass is increased until vanishing at about $3 M_⊙$.

These trends are model dependent. Other stellar evolution codes (e.g. Vassiliadis & Wood 1993) demonstrate that as a result of thermal pulses on the AGB the expanding stellar radius can expand up to 1 au more than the rule-of-thumb maximum radius from the last paragraph. Furthermore, the maximum radius at each phase might differ depending on the model used; see, for example, Villaver & Livio (2009).

### 3.2 Tides

The maximum stellar radius is just a physical boundary; stellar tides can extend beyond the reach of the star. Because modelling tides

<footnote n="2">² In the one-planet case, when the outer planet is removed from these particular simulations, then all three curves are visually indistinguishable from one another on this plot. This result reinforces the finding of Veras et al. (2011) that splicing within time-steps is generally not necessary in one-planet simulations.</footnote>

<footnote n="3">³ By using the observed mass ejecta speed in the post-MS system R Sculptoris of ≈$14.3 \text{ km s}^{-1}$ (Maercker et al. 2012), one can estimate that the ejecta will take 121 d to travel 1 au. Thus, for an inner planet at 10 au and an outer planet at 12 or 13 au, the travel time is less than a year, a small fraction of the inner planet's orbital period.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Maximum stellar radius as a function of progenitor mass and stellar phase. The y-axis is logarithmic, and the x-axis is monotonically increasing in time. Symbols denote different progenitor masses, and the stellar phases are MS = main sequence, HG = Hertzsprung Gap, RGB = red giant branch, CoHe = core helium burning, EAGB = early asymptotic giant branch and TPAGB = thermally pulsing asymptotic giant branch. For Solar metallicity, stellar radii generally increase monotonically throughout post-MS phases except for progenitor masses approximately equal to 1 M⊙.</figcaption></figure>

is both beyond the scope of this study and remains the subject of debate, we choose initial conditions for our numerical simulations where we can safely neglect tidal effects. Planets are not necessarily destroyed by tides nor by being engulfed in the stellar envelope. The remarkable sub-10 h periods of the two planets in the hot B subdwarf star KIC 05807616 system (Charpinet et al. 2011) provide strong evidence that planets can survive deep immersion into the stellar envelope.

Both Villaver & Livio (2009) and Kunitomo et al. (2011) have analysed planet engulfment by RGB stars by including a number of physical factors, and use more detailed stellar evolution models than SSE. Their results indicate that, for lower mass stars, tides can significantly affect planets with semimajor axes that are about 2.5 times as high as the maximum stellar RGB radius, and that engulfment is sensitively dependent on progenitor stellar mass. However, the RGB radius does not greatly exceed 1 au. Although AGB radii are significantly larger (Figs 2 and 3), Mustill & Villaver (2012) found that for higher mass stars, the most distant Jovian planet that becomes engulfed is initially at approximately the maximum stellar radius; tides can slightly shrink the orbits of planets for about an au beyond this maximum. Hence, we adopt 10 au as the initial orbit of the innermost planet in our simulations. This planet will certainly be safe from engulfment in the envelope, and will experience no significant tidal decay except possibly from the most massive stars.

### 3.3 Radius-based code modifications

Our simulation initial conditions are chosen such that stable planetary systems are well beyond the influence of tidal effects. However, instabilities which arise during the simulations may cause planets to approach or collide with the expanding stellar envelope. If a collision occurs, the system is flagged as unstable and is stopped. In reality, if a star engulfs a planet, the star's mass would increase slightly and as a result its radius might change as well. Remaining planets in the system would then be affected because angular momentum must be conserved.

Additionally, the MERCURY collision detection algorithm mce_cent checks to see if the pericentre of a planet's orbit lies

within the star's radius. If so, a collision is flagged. However, if the star's reflex velocity is sufficiently high, then a collision might not occur. Therefore, we modified mce_cent to account for the stellar reflex velocity. For planetary mass companions, however, the reflex velocities are low, and are not likely to factor into collision statistics.

### 4 HILL STABILITY

We now make some analytical stability estimates to identify the systems of interest for our N-body runs. We particularly seek systems that are likely to be stable on the MS and subsequently destabilized during the giant or WD stages. Donnison (2011) provides a formulation of Hill stability in Jacobi coordinates which allows for arbitrarily inclined and eccentric orbits. His treatment is fully general with one exception: the expression for the system energy is the two-body approximation. This approximation is necessary in order to obtain an analytically tractable (but not strictly correct) solution; the intractable terms appear, for example, in equation 2.27 of Veras (2007), which provides the complete expression for the energy of a three-body system in terms of Jacobi orbital elements. See Section 6.5 of this paper for further discussion on this point.

In the following, the subscripts '1' and '2' refer to the inner and outer planets, the subscript 'in' refers to the star/innermost-planet binary, and the subscript 'out' refers to the outer planet properties measured with respect to the inner binary. Let $i_m$ represent the mutual inclination of the inner and outer binaries. Then the Hill stability curve is given by (Donnison 2011)

$$ (1 + y^2) (y^2 \beta^2 + 2y\beta \cos i_m + 1) = \frac{2S_{cr}(M_* + M_1 + M_2)y^2}{M_2^3 (M_* + M_1)^3 (1 - e_{out}^2)} \quad (1) $$

where

$$ y = \sqrt{\frac{a_{in} M_2 (M_* + M_1)}{a_{out} M_* M_1}} \quad (2) $$

$$ \beta = \left( \frac{M_* M_1}{M_2} \right)^{3/2} \sqrt{ \frac{M_* + M_1 + M_2}{(M_* + M_1)^4} \frac{(1 - e_{in}^2)}{(1 - e_{out}^2)} } \quad (3) $$

and with (Donnison 2006)

$$ S_{cr} = \frac{1}{2(M_* + M_1 + M_2)} \left( M_* M_1 + \frac{M_* M_2}{1 + x_{crit}} + \frac{M_1 M_2}{x_{crit}} \right)^2 \times \left( M_* M_1 + M_* M_2 (1 + x_{crit})^2 + M_1 M_2 x_{crit}^2 \right) \quad (4) $$

such that $x = x_{crit}$ is the unique real solution to the following quintic equation:

$$ (M_* + M_1)x^5 + (3M_* + 2M_1)x^4 + (3M_* + M_1)x^3 - (3M_2 + M_1)x^2 - (3M_2 + 2M_1)x = (M_2 + M_1) \quad (5) $$

Care should be taken when choosing the root of the quartic equation in equation (1) when solving for y.

The literature is replete with special-case solutions to equations (1)–(5) (see Georgakarakos 2008 for a review)<fnref n="4" /> but typically treat the masses as static and define a 'separation' as a modulated

<footnote n="4">⁴ We have discovered two typographic errors in the previous literature: The last quantity of the LHS of equation 2.13 in Veras & Armitage (2004) should not be squared, and the sign in front of A in equation 9 of Donnison (2011) is incorrect.</footnote>
ratio of the planetary semimajor axes. Equation (2) demonstrates why. In order to model how the Hill stability curves change with mass-loss, one need only to evaluate equations (1)–(5) at different times during a star’s evolution.

In the circular, coplanar, equal-planetary mass case presented in Debes & Sigurdsson (2002), the critical Hill separation is shown to vary by a few tenths due to significant mass-loss. We have undertaken a wider parameter exploration of phase space, and have discovered that this result generally holds true for orbits of any eccentricity, inclination and stellar mass as long as the planetary masses are at most about 1 *M*<sub>J</sub> each. Our results are presented in Fig. 4, where each panel showcases a different parameter dependence.

The variation in the Hill stability limit due to stellar mass-loss is often equivalent to several au for planets which reside beyond about 10 au. Consequently, if planetary packing (e.g. Raymond et al. 2009) produces planets near the Hill stability limit during the MS, then post-MS evolution may trigger instability.

Despite these considerations, one must remember that the Hill stability criterion is a sufficient but not necessary condition for the planets to remain ordered. Hill stable planets may be Lagrange unstable, and planets failing to satisfy the Hill stability condition may nevertheless be stable. Regarding the latter case, one outstanding feature of Fig. 4 is that moderately eccentric or inclined orbits yield critical semimajor axes ratios that are high – much higher, for example, than the semimajor axis ratios of any adjacent pair of planets in the Solar system.

Large critical semimajor axis ratios may strongly influence the location where two planets become Lagrange stable. Veras & Armitage (2004) show that as the mutual inclination of the circular orbits of two equal-mass planets increases, the critical Hill stability limit becomes a progressively worse indicator of the separations at which planets may actually become Lagrange stable. Their fig. 5 illustrates that for *i*<sub>m</sub> = 35°, instability occurs effectively for values under half of the critical separation. However, their numerical simulations were run for just 2 Myr, almost certainly missing instances of longer term instability.

Therefore, determining how mass-loss affects the stability prospects of the orbits of two planets is perhaps more complex than just considering the analytic effect on the Hill stability limit. Hence, we now turn to *N*-body simulations.

# 5 *N*-BODY SIMULATIONS

Ideally, we could sample the entire two-planet/single-star phase space with numerical simulations. Realistically, we must take measures to restrict our studies to computationally feasible and insightful simulations. To better understand how to restrict the phase space, we consider typical MS lifetimes, the mass lost at each evolutionary phase and the planetary period enhancement during the WD phase, in Sections 5.1–5.3. We motivate and state our initial conditions in Section 5.4. Sections 5.5 and 5.6 present the simulation results.

## 5.1 MS time-scales

MS ages are given in Fig. 5. These ages may vary by Gyr depending on stellar metallicity, and are at least 1 Gyr long for any progenitor mass less than 2 M<sub>☉</sub>. With an MS lifetime of over 10 Gyr, a Solar-mass, Solar-metallicity star with orbiting planets is particularly prohibitive to integrate. This long time-scale helps explain the uncertainties in long-term evolution of the Solar system planets (Kholshevnikov & Kuznetsov 2007; Laskar et al. 2011). Fig. 6

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Critical Hill semimajor axis ratios as a function of stellar mass for different eccentricities (upper panel), different mutual inclinations (middle panel) and different planetary masses (bottom panel). The purple stars on the bottom of each plot, from right to left, represent the eventual WD mass (in brackets in the following) for MS progenitor masses of 8 M<sub>☉</sub>[1.44 M<sub>☉</sub>], 7 M<sub>☉</sub>[1.29 M<sub>☉</sub>], 6 M<sub>☉</sub>[1.14 M<sub>☉</sub>], 5 M<sub>☉</sub>[1.00 M<sub>☉</sub>], 4 M<sub>☉</sub>[0.87 M<sub>☉</sub>], 3 M<sub>☉</sub>[0.75 M<sub>☉</sub>], 2 M<sub>☉</sub>[0.64 M<sub>☉</sub>] and 1 M<sub>☉</sub>[0.52 M<sub>☉</sub>]. The figure illustrates that the Hill radius is sensitively dependent on *e*<sub>out</sub>, *e*<sub>in</sub>, *i*<sub>m</sub>, *M*<sub>1</sub> and *M*<sub>2</sub>, but weakly dependent on *M*<sub>*</sub>. For Jovian and terrestrial-mass planets, stellar mass-loss can change the Hill stability limit by at most a few tenths in *a*<sub>out</sub>/*a*<sub>in</sub>.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** The MS age of stars with Solar metallicity (blue, top curve) and with a very low ($Z = 0.0001$) metallicity (red, bottom curve). Higher mass progenitors significantly reduce the CPU time needed to integrate planets over a star's entire MS lifetime.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** The number of orbits taken by an isolated planet around a star throughout its MS phase. In each pair of equivalently coloured curves, the top curve is for $Z = 0.02$ and the bottom curves is for $Z = 0.0001$. The number of orbits decrease with higher stellar mass because the decreased MS time-scale dominates the shortened orbital period. These values are important both physically – to determine instability – and computationally – to assess the feasibility of integrating over the entire MS phase.</figcaption></figure>

provides estimates for the number of planetary orbits that would be achieved during the MS as a function of stellar progenitor mass, for a variety of planetary semimajor axes. The curves result from the competition between the decreased orbital time-scale with increased stellar mass, and the decreased MS lifetime with increased stellar mass; the latter overwhelmingly wins. Note that the number of orbits tails off significantly as $M_*(0)$ is increased, for all semimajor axes. A one order-of-magnitude change in semimajor axis corresponds to 1.5-order-of-magnitude change in the number of MS orbits.

## 5.2 Post-MS phase properties

The evolution time-scales of the intermediate post-MS, pre-WD phases are short compared to the MS time-scale. Fig. 7 provides time-scales for each stellar phase, and relates the phase to the percentage of the star's original mass lost, for Solar-metallicity stars. The plot demonstrates that except for the $1 M_{\odot}$ case, most mass is lost during the TPAGB and negligible percentages of mass are lost in the other phases.

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Correlating mass-loss fractions, phase durations and progenitor masses for stars of Solar metallicity and a Reimers mass-loss coefficient of 0.5. Each curve contains eight symbols representing stellar progenitor masses of 8, 7, 6, 5, 4, 3, 2 and $1 M_{\odot}$, ordered monotonically. The green squares represent the thermally pulsing asymptotic giant branch (TPAGB), the brown upward-pointing triangles the early asymptotic giant branch (EAGB), the grey downward-pointing triangles the core helium burning phase (CoreHe), the red diamonds the red giant branch (RGB), the yellow circles the Hertzsprung gap (HG) and the blue open squares the MS. Most MS mass-loss is too small for this plot. Except for the $1 M_{\odot}$ case, the most important mass-loss phases are the TPAGB, which all last of the order of 1 Myr.</figcaption></figure>

## 5.3 WD period enhancement

After the star has become a WD, the star stops losing mass and gradually cools down. Compared to its MS mass, the WD mass is greatly reduced, and cannot exceed the Chandrasekhar Limit of ≈1.4 $M_{\odot}$. The result is that the planetary period may be drastically increased. Assume that the planet's evolution is entirely adiabatic. A reduction of the star's mass by a factor of *k* will cause the planet's semimajor axis to be increased by a factor of *k*. Hence, the planetary period around the WD is $k^2$ times the period around the MS star. Fig. 8 plots this enhancement factor as a function of progenitor stellar mass, and demonstrates both that WD planets perform fewer orbits than MS counterparts in the same amount of time (with possible implications for scattering) and that WD numerical simulations may proceed much more quickly than MS simulations.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** The planetary period enhancement factor due to post-MS evolution, supposing that the planet has evolved entirely adiabatically and is isolated from any other perturbations. This enhancement factor is independent of semimajor axis. The top curve is for $Z = 0.02$ and the bottom curve is for $Z = 0.0001$.</figcaption></figure>

## 5.4 Initial conditions

### 5.4.1 Fiducial choices

The above considerations lead us to choose an integration duration of 5 Gyr and Solar-metallicity progenitor masses of 8, 7, 6, 5, 4 and 3 M⊙. This mass range extends down to the upper mass-end of the observed range of exoplanet host stars (Sato et al. 2012b). This combination allows us to sample an ensemble of multi-planet systems over every phase of evolution, including a substantial sampling (over 4.5 Gyr) of evolution in the WD phase. Simulation output occurs at a frequency of 1 Myr. Performing a statistically significant simulation ensemble for 1 M⊙ stars is well beyond our available resources; for more details on the planetary consequences of the possible evolutionary tracks of 1 M⊙ stars, see Veras & Wyatt (2012).

We adopt 1 $M_1$ for the mass of each planet ($M_1 = M_2 = M_1$), assume that the planets are on coplanar orbits ($i_m = 0°$) and assume that they have small but non-negligible MS eccentricities [$e_1(t = 0) = e_2(t = 0) = 0.1$]. These eccentricities are low compared to the observed MS values for massive exoplanets beyond the tidal circularization limit, but higher than the near-circular orbits predicted by core accretion theory. The inner planet is initially set at $a_1 = 10$ au to avoid tides with the expanding stellar envelope; over 25 known planets are thought to harbour $a > 10$ au.<fnref n="5" /> Also, this semimajor axis guarantees that orbital evolution due to mass-loss will be adiabatic unless in the presence of an event such as a supernova, which is not modelled here.

We perform 632 simulations per ensemble, where each ensemble features a different progenitor mass. In each individual simulation, the orbital angles (mean anomalies and longitudes of pericentre) of both planets are selected from a uniform random distribution. We adopt 79 values of $a_2$ so that we sample 8 different sets of orbital angles for each $a_2$ value. The range of $a_2$ values sampled encompasses both the ‘chaos limit’ and the Hill stability limit in order for us to sample many different types of dynamical behaviour.

The chaos limit refers to a maximum semimajor axis ratio separation at which mean-motion resonances do not necessarily overlap. This limit is smaller than the Hill and Lagrange stability limits, and

represents a fuzzy boundary within which instability occurs readily and quickly. Wisdom (1980) found the chaos limit to be

$$ \frac{a_2 - a_1}{a_1} \approx 1.3 \left( \frac{M_p}{M_*} \right)^{2/7} \qquad (6) $$

for two equal-mass circular-orbit planets. Recently, Quillen & Faber (2006) and Mustill & Wyatt (2012) expanded on this result by considering bodies’ eccentricity. Mustill & Wyatt (2012) discovered that for $e_m > 0.21(M_p/M_*)^{3/7}$, where $e_m$ is the eccentricity of the least massive planet, the chaos limit is (for $M_p = M_1 = M_2$)

$$ \frac{a_2 - a_1}{a_1} = 1.8 e_m^{1/5} \left( \frac{M_p}{M_*} \right)^{1/5} . \qquad (7) $$

We set our minimum value of $a_2$ to be less than the limit from the more conservative definition (equation 6) for each progenitor mass to help ensure that we have a tail of unstable simulations. Similarly, we wish to have a tail of stable simulations for large separations. Hence, we select a maximum value of $a_2$ that exceeds the MS Hill stability limit in each case.

### 5.4.2 Additional simulations

Motivated by the outcome of our fiducial simulations, we performed two additional ensembles of simulations (632 simulations per ensemble). In both, we adopted a stellar progenitor mass of 5 M⊙. The first case assumed different planetary masses; we adopted 1 M⊕ for each planet ($M_1 = M_2 = M_⊕$). In the second, we adopted $e_2 = 0.5$ to model a moderately eccentric outer planet. Doing so yields a much wider Hill stability separation (see Fig. 4) than in the fiducial case.

## 5.5 Simulation results: fiducial cases

### 5.5.1 Overview

Our goal is to identify instability and when it occurs. We define instability as Lagrange instability: if the planets at any point are found to achieve a hyperbolic orbit, cross orbits, or collide with each other or the star. Hill instability includes just a few of these possibilities: planet–planet collisions and crossing orbits. Therefore, Hill stable systems may eventually be Lagrange unstable. Those that do will feature ejection of the outer planet and/or collision of the inner planet with the star.

We plot instability time versus initial semimajor axis ratio in Fig. 9, which represents our main result. Dots indicate unstable systems. No dot appears for a system that has remained stable over the entire 5 Gyr integration. If all eight simulations for a given semimajor axis ratio are stable over 5 Gyr, then we place a blue star at $10^{10}$ yr in the appropriate horizontal position, even though the vertical position of the star has no physical meaning and is selected for visual impact.

The figure includes the six ensembles of simulations with different progenitor stellar masses. Each post-MS phase change occurs at different times on each plot. See Hurley et al. (2000) for detailed physical descriptions of each phase. Although the horizontal lines are close together in Fig. 9, they are clearly distinguishable on the zoomed-in Fig. 10 plots. Different progenitor masses also cause differences in the location of the Hill stability limit. On each plot is a black vertical dashed line, representing the Hill stability limit computed from the star’s MS mass. The black dotted lines represent the Hill stability limits computed with the star’s WD mass (see the caption to Fig. 4 for these values).

<footnote n="5">⁵ See the Extrasolar Planet Encyclopedia at http://exoplanet.eu/</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Instability times versus initial semimajor axes ratios for full-lifetime two-planet simulations. The dots indicate individual unstable systems, and the blue stars indicate that all eight systems at that separation ratio were stable over 5 Gyr. The blue stars not shown extend out to at least $a_2/a_1 = 1.747$ in an unbroken string in each plot. Each coloured horizontal line represents a stellar phase change, and is at a different position on each plot. The two vertical lines represent the Hill stability limit for the MS (dashed) and WD phase (dotted) for each progenitor mass. Any unstable systems on the MS beyond the MS Hill stability limit (such as for the 5, 4 and 3 $M_\odot$ cases) or during the WD phase beyond the WD Hill stability limit are Lagrange unstable. The plot demonstrates that instability during the WD phase can be achieved at separations that well-exceed both the MS and WD Hill stability limits.</figcaption></figure>
<figure><img src="image_11.png" /></figure>

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

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

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

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

<figure><img src="image_16.png" /><figcaption>**Figure 10.** Zoomed-in versions of Fig. 9 to show detail at times of stellar phase changes. In ascending vertical order, the phases are ‘Hertzsprung gap’ = blue, ‘red giant branch’ = red, ‘core helium burning’ = green, ‘early AGB’ = orange, ‘TPAGB’ = purple, ‘WD’ = grey. The longest pre-WD post-MS phase is the red giant branch; the most violent phase (with the greatest mass-loss, and causing the greatest amount of instability) is the TPAGB. The WD Hill stability limit acts as an effective empirical boundary for pre-WD post-MS instability.</figcaption></figure>

5.5.2 Physical description of Fig. 9

The Hill stability limits and the onset of post-MS evolution provide
boundaries within which one can understand the following different
regions in the Fig. 9 plots.

(1) During the MS, represented by the region under all the coloured horizontal lines, dots appear predominantly inside of the MS Hill stability limit and predominately at times under 10^7 yr. Hence, the limit is useful for identifying short-term instability.
(2) Some dots appear on the MS but outside of the MS Hill stability limit in the plots with progenitor stellar masses of 5, 4 and $3M_\odot$. All these dots indicate long-term instability (occurring after ~$10^7$ yr). The long-term MS unstable simulations with initial separations exceeding the Hill stability limit must be (and indeed are) Lagrange unstable: where the outer planet is ejected and/or the inner planet collides with the star. As the progenitor mass is decreased, the number of unstable MS systems beyond the MS Hill stability boundary increases. One possible reason is because the longer MS lifetimes (see Fig. 5) translate into more orbits for the planets (see Fig. 6), and hence a greater opportunity for instability to occur. Another potential reason is that smaller planet-star mass ratios broaden the boundaries between Hill and Lagrange instability.

(3) MS instability beyond the MS Hill stability limit appears to extend only as far as the WD Hill stability limit. However, this must be coincidence – due to our choice of initial conditions – as the planetary system has no knowledge of the post-MS mass-loss that will occur.

(4) The WD Hill stability limits ensure that any WD instability occurring beyond this limit must be Lagrange instability. Our simulations corroborate this statement.

(5) Each plot in Fig. 9 demonstrates that all systems become Lagrange stable for $a_2/a_1 \gtrsim 1.55$, just inside the 2:1 commensurability. Reinforcing this estimate are blue stars which were excluded from the plot (for visual clarity) that extend all the way out in an unbroken chain to $a_2/a_1 \approx [1.779, 1.775, 1.770, 1.765, 1.757, 1.747]$ for $M_\odot(0) = [8, 7, 6, 5, 4, 3 M_\odot]$. This sampled range helps to establish the robustness of the Lagrange stability boundary for our chosen 5 Gyr integration duration. This boundary lies at a distance corresponding to approximately [178, 176, 176, 172, 167, 163 per cent] of the MS Hill stability limit and [138, 134, 133, 130, 126, 123 per cent] of the WD Hill stability limit.

In all our cases, planets with initial separations that exceed the 2:1 commensurability are stable throughout the 5 Gyr integration.<fnref n="6" /> However, the Lagrange stability limit may be higher than the values reported here for progenitor masses lower than those sampled here. Especially for $1 M_\odot$ stars, the outcome is uncertain, given the long MS lifetime and strong mass-loss over both the RGB and AGB (see Fig. 7). The distribution of unstable systems in Fig. 9 shows a wide variation of instability times, and suggests that longer simulations, perhaps out to the age of the Universe, could yield additional instability.

### 5.5.3 Description of Figs 10–12

Post-MS pre-WD phase changes can prompt instability, which can be seen more clearly in Fig. 10. The phases which cause the greatest mass-loss (see Fig. 7) also tend to be the phases which are most likely to trigger instability. This tendency is indicated by the number of dots between the TPAGB and WD lines versus the number of dots appearing below the TPAGB line. Note, however, that the relatively long length of the core helium burning phase helps to increase the number of unstable systems during that phase. The abundance of dots just above the WD line indicates that the TPAGB can unsettle stable systems enough to cause slightly delayed instability. In Fig. 10, the WD Hill stability limit appears to provide an

effective boundary beyond which post-MS pre-WD instability does not occur. However, this apparent boundary again must be coincidence because the systems are unaware of forthcoming post-MS mass-loss.

In support of the above claims, we quantify the types of instability in Fig. 11 for each progenitor mass. The figure shows six normalized bar plots. The blue, or topmost, bars represent the fraction of systems (out of eight) that feature a collision between both planets; orange bars represent a collision with the central star; yellow bars represent any other type of instability (including ejection or periodic instances of a planet attaining a hyperbolic orbit); green, or bottommost, bars indicate systems which remained stable for 5 Gyr. Therefore, Lagrange unstable systems are represented by orange and yellow bars. The predominance of the yellow bars towards the right sides of the plots indicate that any type of collisions become less likely as the initial planet separations are increased. Any blue bars beyond the MS Hill stability limit indicate planet-planet collisions during the post-MS, importantly demonstrating that after leaving the MS, planets are not restricted to (but still predominately experience) Lagrange instability. There are no blue bars that exceed the WD Hill stability limit, as expected. The height of the green bars around the 3:2 commensurability demonstrates how it helps stabilize the simulated systems.

Also of interest is the evolutionary phase at which instability occurs. The collision of a planet with a star has been proposed to explain both the existence of extreme horizontal branch stars without stellar binary companions – as the envelope of the progenitor giant could be removed by the planet (Charpinet et al. 2011; Bear & Soker 2012) – and the enrichment of Lithium seen in a few per cent of stars at all parts of the RGB (Lebzelter et al. 2012). Of particular interest here is the planet candidate proposed orbiting the Lithium-rich giant BD+48 740, whose eccentric orbit suggests a past strong scattering interaction such as we are considering (Adamów et al. 2012).

Although the fraction of post-MS pre-WD instability can be deduced from Figs 9 and 10, we have created a separate figure, Fig. 12, which better visualizes the result. Fig. 12 displays the fraction of systems which are stable (green bars), and those which incur instability on the MS phase (purple bars), on the WD phase (orange bars) and in between the MS and WD phases (grey bars). Instability during a giant branch phase occurs relatively infrequently: [4.0, 5.4, 6.9, 6.7, 8.7, 4.6 per cent] of all unstable systems for each stellar progenitor mass. WD instability is not limited to systems with initial separations beyond the Hill stability limit.

### 5.5.4 Potential resonance behaviour

Mean-motion commensurabilities, shown on the upper x-axes of Figs 9–12, appear to play an important role in affecting stability. At these locations, stability may be either enhanced or additionally disrupted. A few suggestive instances of these commensurabilities making a contribution include the 5:4 location in the $M_*(0) = 8 M_\odot$ simulations and the 5:3 location in the $M_*(0) = 3 M_\odot$ simulations. The 3:2 commensurability demonstrably provides a protection mechanism for planetary systems in each plot. The 4:3 commensurability seems to yield all possible outcomes for $M_*(0) \ge 5 M_\odot$, but no stable systems for $M_*(0) < 5 M_\odot$.

Hence, exploring the resonant character of these systems is of potential interest. However, our 5 Gyr simulations are not well suited to determine if a given system is locked into mean-motion

<footnote n="6">⁶ The strong, first-order 2:1 mean-motion commensurability perhaps plays a role in establishing the Lagrange stability boundary (see, e.g., Barnes & Greenberg 2007) for our fiducial cases.</footnote>
<figure><img src="image_17.png" /></figure>

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

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

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

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

<figure><img src="image_22.png" /><figcaption>**Figure 11.** Types of instability. Blue indicates the fraction of systems (out of eight per bar) that went unstable because of a collision between the two planets. Orange represents instability due to a collision with the central star. Yellow indicates any other type of instability, which predominantly includes ejection. Green indicates no instability. Hence, Hill unstable systems are included in the blue bars only, and Lagrange unstable systems are included in the yellow and orange bars only. The black dashed and dotted lines are the MS and WD Hill stable boundaries, as in Fig. 9.</figcaption></figure>

resonance because our output frequency of 1 Myr is over four orders of magnitude greater than a typical orbital period. The sudden and drastic changes in eccentricity and inclination which can arise from purely three-body interactions (not including any type of dissipation nor external forces) may act well within 1 Myr (e.g. Naoz

et al. 2011), and hence disrupt and/or create resonances. Additionally, resonance behaviour may manifest itself only periodically due to repeated separatrix crossings, which yield different intervals of libration and circulation of one or more resonant angles (e.g. Farmer & Goldreich 2006). Recently, Ketchum, Adams & Bloch
<figure><img src="image_23.png" /></figure>

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

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

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

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

<figure><img src="image_28.png" /><figcaption>**Figure 12.** Phases of instability. Green indicates the fraction of systems (out of eight per bar) that were stable over 5 Gyr, as in Fig. 11. Purple indicates that instability occurred on the MS. Grey indicates that instability occurred between the MS and WD phases. Orange indicates that instability occurred during the WD phase. The plot illustrates that instability during a giant branch phase does occur, but infrequently. Also, WD instability can occur for Hill unstable planets which survive until the WD phase. The yellow dashed and dotted lines are the MS and WD Hill stable boundaries, as in Figs 9–11.</figcaption></figure>

(2013) have classified this behaviour as ‘nodding’ and analytically characterized it.

Despite these caveats, we have considered the evolution from selected resonant angles from our output. Identifying resonant systems requires defining a libration centre, maximum libration amplitude

about this centre and a duration. Veras & Ford (2009, 2010) demonstrated that fully characterizing potential resonant behaviour for two massive planets may require the sampling of several libration centres, as well as computing a mean absolute deviation or root-mean-squared deviation about each centre, for each resonant angle.
<figure><img src="image_29.png" /><figcaption>**Figure 13.** The geometric mean of instability times when at least one system out of eight per semimajor axis ratio goes unstable in each fiducial ensemble of simulations. The plot suggests that planets near strong mean-motion commensurabilities tend to survive for longer times before going unstable.</figcaption></figure>

Here, we do not pursue any such analysis, but instead simply point out that the majority of stable fiducial systems over 5 Gyr do not appear to exhibit resonant libration. These systems include planets with Hill unstable separations. Conversely, a smaller number of systems do appear to exhibit resonant libration, typically close to the strong first-order mean-motion resonant commensurabilities. This result is expected, as forming resonances of Jovian-mass planets from uniformly-sampled orbital angles should be infrequent (see, e.g., Veras 2007).

A perhaps better measure of the effect of mean-motion commensurabilities is by considering the geometric mean of instability times for each semimajor axis ratio set of eight simulations that produces at least one unstable system. We plot these mean times in Fig. 13, where the 4:3 commensurability is shown to have a clear effect.

### 5.5.5 Survivor orbit properties

Of potential interest to WD pollution investigators and WD planet hunters is the properties of planets undergoing Lagrange instability during the WD phase. In some cases, the inner planet simply collides directly with the WD, creating a direct pollutant. However, our simulations suggest that these occurrences are rare, occurring [1, 0, 2, 2, 0, 0] times for $M_\odot(0) = [8, 7, 6, 5, 4, 3 M_\odot]$.

In other cases, the outer planet is ejected and the inner planet survives on a bounded orbit. A bound, eccentric inner planet may induce pollution by scattering comets or asteroids close enough to the WD to be tidally disrupted and ingested by the WD. To gain insight into how an inner planet survives on an eccentric bounded orbit, we use conservation of energy and angular momentum. Although energy is not conserved in a system with mass-loss – and certainly not in our integrations – after the parent star has become a WD, mass-loss ceases and then energy is conserved for the future. Thus, we can compare the states at the beginning of the WD and at the moment of ejection (≡ $t_{ins}$). Furthermore, because the mass-loss is adiabatic, we can relate the semimajor axes of the planets on the MS and the WD phases through knowledge of how much mass is

lost. Therefore, we find that the semimajor axis of the bound planet should be at most

$$a_1(t_{ins}) \le \frac{M_*(0)a_1(0)a_2(0)M_1}{M_*(t_{WD})[a_2(0)M_1 + a_1(0)M_2]}, \quad (8)$$

where $M_*(t_{WD})$ is the mass of the WD.

Angular momentum is conserved throughout a planetary system's life, even under the effects of mass-loss. Thus, in principle, one can use conservation of angular momentum to determine the value of $e_1(t_{ins})$. However, doing so requires knowledge of the hyperbolic values of $a_2(t_{ins})$ and $e_2(t_{ins})$. These values are set by the ejection velocity, which is determined by the strength of the instability in each case. Our numerical simulations show that $e_1(t_{ins})$ varies considerably.

We plot the semimajor axes (blue squares) and periastra (orange dots) of the surviving planets in systems featuring ejections in the WD phase only, in Fig. 14. Superimposed on the plots through aqua lines are the analytically predicted maximum values of $a_1(t_{ins})$ through equation (8). Also plotted as a solid black horizontal line is the maximum stellar radius (see Figs 2 and 3) attained during the star's evolution. The presence of orange dots below the black line suggests the presence of a population of highly eccentric planets orbiting WDs whose present pericentres take them inside the maximum AGB radius. These planets cannot have been formed in situ or anywhere near their WD locations because otherwise they would have been destroyed or suffered radical orbital alterations on the AGB. The fractions of orange dots below the black line for each progenitor mass are 8.4 per cent (8 M$_\odot$), 10.9 per cent (7 M$_\odot$), 8.4 per cent (6 M$_\odot$), 5.1 per cent (5 M$_\odot$), 1.7 per cent (4 M$_\odot$) and 0 per cent (3 M$_\odot$).

## 5.6 Simulation results: additional cases

### 5.6.1 Terrestrial planet masses

Now we consider variations on the fiducial case. First, we set the planet masses such that $M_1 = M_2 = M_\oplus$. The results of those 632 simulations are summarized in Fig. 15, which include plots for instability time versus initial semimajor axis ratio, the types and phases of instability, and the geometric mean of instability time. The smaller planetary masses here shrink the chaos limit and the Hill stability limit – as well as the difference between the MS and WD Hill stability limits (see the bottom panel of Fig. 4) – and introduces a different set of potentially important commensurabilities, displayed on the top x-axis of all plots.

Earth-mass planets fail to go unstable beyond the MS Hill stability limit in all but a few cases; the blue stars continue in an unbroken chain out to $a_2/a_1 \approx 1.510$. Therefore, Hill stability and Lagrange stability appear to have almost identical boundaries in this case.

No system features a collision of a terrestrial planet with the WD. This result should not imply that this type of collision cannot occur, but rather that giant planet collisions are likely to be much more frequent.

Almost all instability involves collision between both planets. Interestingly, one planet–planet collision occurs just beyond the WD Hill stability boundary, meaning that the real boundary must differ from the line on the plot. This situation arises because planets arrive on the WD phase with slightly different osculating orbital parameters than they harboured at the start of the MS (primarily due to their mutual interactions, and slightly due to post-MS mass-loss). Hence, the real WD Hill stability limit for each individual
<figure><img src="image_30.png" /><figcaption>**Figure 14.** Orbits of inner planets that survive the ejection of the outer planet during the WD phase. The blue squares represent semimajor axes, and the orange dots represent periastra, all from the simulation outputs. The aqua curve represents the analytical estimate for where the maximum semimajor axis should be (equation 8). The black curve indicates the maximum stellar radius achieved during its evolution. Dots below this curve indicate the existence of WD planets orbiting inside the maximum AGB radius even though the planets were formed elsewhere.</figcaption></figure>

system is different, and differs from the line shown in the figure that was computed for $e_1 = e_2 = 0.1$ exactly.

Although only one system undergoes instability in the post-MS pre-WD phases, instability during the WD phase is common, and is

in fact greatest for separations close to the tightly packed 11:10 commensurability (corresponding to $a_2/a_1 \approx 1.066$). Other instances of WD instability occur around the 7:6 and 5:4 commensurabilities, and to either side of the 4:3 commensurability, on which all
<figure><img src="image_31.png" /><figcaption>**Figure 15.** Here, $M_1 = M_2 = 1 M_\odot$. The plots are of the same format as in Figs 9–13. Instability during the WD phase occurs almost exclusively at initial separations inside the MS and WD Hill stability limits. MS Hill unstable systems protected by mean-motion commensurabilities, many of which feature resonant behaviour, allow for survival during the MS before becoming unstable on the WD phase.</figcaption></figure>

eight systems are stable. The influence meted by these first-order commensurabilities therefore appears to be extensive. However, the bottom plot of the figure might suffer bias due to small number statistics. An investigation of individual systems indicates that resonant

behaviour appears to be more common for these terrestrial-mass planet systems than for the fiducial Jovian-mass systems, perhaps demonstrating the importance of the planet/star mass ratio in creating resonance behaviour from random initial conditions.
<figure><img src="image_32.png" /><figcaption>**Figure 16.** Here, $e_2(0) = 0.5$. The plots are of the same format as in Figs 9–14. This figure importantly demonstrates that Lagrange instability can occur readily in moderately eccentric systems.</figcaption></figure>

### 5.6.2 Different planetary eccentricities

Next, we consider systems with an initially moderately eccentric outer planet orbit ($e_2 = 0.5$) and the same fiducial inner planet eccentricity ($e_1 = 0.1$) in Fig. 16. In this case, the analytic MS Hill stability limit is much larger ($a_2/a_1 \approx 2.73$) than in the fiducial case, which is expected (see Fig. 4).

For these systems, instability beyond the MS and WD Hill limits is extensive. In fact, the Lagrange unstable region extends out nearly 7 au from the location of the MS Hill stability boundary for $a_1 = 10$ au. However, the Lagrange stability boundary lies at a distance that is 143 and 135 per cent of the MS and WD Hill stability limits; the former value is a much smaller factor than in the fiducial cases.
Collisions with the star are infrequent and ejections are common, perhaps because the inner planet has the much lower eccentricity. In fact, inner planets which survive outer planet ejections largely fail to attain high enough eccentricities to intrude within the maximum stellar radius achieved during the host star's evolution: only 0.85 per cent of the orange dots are below the black line in the bottom-right plot. In that same plot, the two outliers likely represent systems that experienced pre-WD instability that was missed by our low output frequency. Instability during the giant phases is uncommon, and is restricted to the region around the 4:1 commensurability.

# 6 DISCUSSION

## 6.1 Consequences for WD pollution

### 6.1.1 Background

One potential consequence of dynamical instabilities induced due to stellar mass-loss is WD pollution. WDs are surrounded by thin atmospheres of either hydrogen or helium. Heavier elements sink rapidly in such thin atmospheres and it is therefore puzzling that such a high fraction of WDs have evidence for metal pollutants in their spectra (25 per cent of single DA WDs; Zuckerman et al. 2003). Such metal pollution has been associated with excess emission in the infrared consistent with a close-in dust disc in more than a dozen cases (e.g. Zuckerman & Becklin 1987; Becklin et al. 2005; Kilic et al. 2005; Reach et al. 2005; Farihi, Jura & Zuckerman 2009) and gas discs in a handful of cases (Gänsicke et al. 2006, 2008; Gänsicke, Marsh & Southworth 2007). It has been suggested that these observations could be explained by planetary material, motivated by the similarity of the composition of the accreted material with planetary material (Klein et al. 2010; Girven et al. 2012).

In order for a WD to be polluted by material from an outer planetary system, comets, asteroids or planets must be scattered at least close enough to the star (at about $1R_\odot$) so that they are tidally disrupted. Changes to the dynamics of the planetary system following stellar mass-loss has been suggested as a potential cause of increased numbers of planetary bodies scattered on to star-grazing orbits (Debes & Sigurdsson 2002; Jura 2008; Bonsor et al. 2011; Debes et al. 2012). Even in planetary systems where the planets remain on stable orbits, Bonsor et al. (2011) and Debes et al. (2012) show that sufficiently many asteroids or comets can be scattered on to star-grazing orbits to explain some of the observations of polluted WDs. Such mechanisms, however, struggle to explain observations of high accretion rates in old polluted WDs (Koester et al. 2011; Girven et al. 2012).

Instabilities in planetary systems could provide a potential explanation for pollution in these, and other, WDs. Depending on the exact nature of the instability and structure of the individual planetary system, in many cases the number of planetary bodies scattered on to star-grazing orbits increases significantly following an instability. This means that the planetary systems presented in this work in which instabilities occurred during the WD phase have the potential to produce polluted WDs.

### 6.1.2 Our simulations

We qualify the following discussion by reminding the reader that we have not probed progenitor masses between 1 and $2M_\odot$, where the true WD population appears to peak. The reason for not considering this range is due to computational limitations, as tracking planetary orbits over such long MS lifetimes is prohibitive.

Fig. 9 showcases the different types of instability which can occur during the WD phase. First, many systems with giant planets that were Hill stable on the MS become Hill unstable due to mass-loss preceding the WD phase. Secondly, other systems that were technically Hill unstable on the MS, but were protected against instability by commensurabilities, become unstable during the WD phase. Thirdly, some of the systems that are Hill stable during the WD phase are actually Lagrange unstable, and experience instability at a late time.

We can now relate the WD instability to the cooling age of the WD, which is the time lapsed since the star became a WD. The cooling age of real (not simulated) WDs are readily determined from observations, and hence provide an opportunity for comparison with and motivation for numerical simulations. Fig. 17 displays the number of planetary systems that become unstable during the WD phase as a function of cooling age; each plot corresponds to a different progenitor mass. Any instance of instability could represent a trigger for a potentially polluting event. An ejection especially suggests that the inner planetary system will be significantly perturbed, and this perturbation can throw planetesimals or terrestrial planets that were originally in the inner system on to the star.

The figure contains several notable features: (1) planet–planet collisions tend to occur at short cooling ages ($\sim 10^4–10^7$ yr), (2) escape and stellar collisions becomes prevalent only after $\sim 10^7$ yr, but tails off after a few Gyr, (3) the cooling age of this type of instability increases as the progenitor mass increases and (4) direct collisions with the WD are rare. The striking delay of $\sim 10^7$ yr before escape instability becomes dominant reflects the (largely unexplored) dependences of Lagrange instability on number of orbits completed, the initial values of the semimajor axes and the mass ratios. In fact, perhaps the positive correlation of cooling age instability with progenitor mass is due to more orbits being required for smaller mass ratios for Lagrange instability to kick in, a feature also apparent on Fig. 9 on the MS. Alternatively, the correlation may result from the much wider orbits of the giant planets around higher mass WDs.

The few giant planets which hit the WD are particularly interesting as pollution sources, and imply that some polluted WDs may result from the disruption of giant planets rather than comets or asteroids. Planetary material that is accreted on to WDs might be composed of many small asteroids or comets (Jura 2003; Bonsor et al. 2011; Debes et al. 2012), but could also be the result of a single large object. Detailed compositional analysis of some objects concludes that the accreted material resembles more closely the bulk Earth in composition than chondritic material (e.g. Zuckerman et al. 2007; Klein et al. 2010); it is feasible that a disrupted planet would produce high levels of pollution in a WD. Such disrupted planets could provide the explanation for a handful of polluted WDs, in particular old, heavily polluted WDs. However, our simulations indicate that star-grazing giant planet collisions are too rare to alone account for the abundance of observed differentiated material. This result agrees with the implied low fraction of planetary collisions with WDs from the finding that at most 1–5 per cent of WDs have high accretion rates due to dust discs (Debes et al. 2011; Farihi et al. 2012a).

Observed polluted WDs can provide the inspiration for extensions to our work, with numerical simulations tailored to particular observational samples or campaigns. However, the purpose of this paper is not to make a detailed comparison with the observed polluted WD population, which is heavily biased and suffers from uncertainty of the distribution of their progenitor masses. Our computationally expensive simulations were set up to demonstrate that full-lifetime orbital integrations including realistic mass-loss are
<figure><img src="image_33.png" /><figcaption>**Figure 17.** Instability as a function of the time elapsed since the end of the of the AGB phase, when the WD was born. Planet–planet collisions occur quickly (~$10^4–10^7$ yr), whereas instances of escape or stellar collisions typically do not occur until after a few tens of Myr.</figcaption></figure>

now achievable, and to explore dynamical properties of the result-
ing instability. Nevertheless, we can make a crude comparison of
the potentially WD-polluting events in our simulations (Fig. 17) and
the observed distribution of polluted WDs (in Fig 18). Data for the
observed distribution of 78 WD ages were taken from Farihi et al.
(2012b) and Girven et al. (2012).

Both the observations and simulations show a broad consensus.
Observations confidently detect pollution at cooling ages from a

few tens of Myr (e.g. Gänsicke et al. 2012) to several Gyr (Koester
et al. 2011), but cannot detect pollution at earlier times when the
WD is too hot to rule out a primordial origin for metals<fnref n="7" /> nor at later
times when the WDs are too intrinsically faint. The simulations

<footnote n="7">⁷ These young WDs are able to radiatively levitate metals, masking what could be external pollution.</footnote>
<figure><img src="image_34.png" /><figcaption>**Figure 18.** Cooling ages for a sample of 78 observed WDs from Farihi et al. (2012b) and Girven et al. (2012). Although this figure may be compared with Fig. 17, the progenitor masses of these WDs are unknown, and the sample suffers from observational bias.</figcaption></figure>

begin to demonstrate ejection instability only at times greater than tens of Myr due to a true physical effect which may be explained by the time-scale for Lagrange instability to act on the WD phase. This instability appears to have largely run its course after a few Gyr have passed, so that relatively fewer systems are likely to become unstable beyond the 5 Gyr integration time. Observationally, comparing the pollution frequency at 5 Gyr with that at tens of Myr is difficult due to intrinsic biases, despite the suggestive nature of the distribution in Fig. 18. Nevertheless, Bonsor et al. (2011) found that asteroidal accretion on to WDs from particles scattered due to the overlap of mean-motion resonance exterior to the planet follows a similar trend of being present beyond a few tens of Myr after post-MS evolution before eventually tailing off after a few Gyr. Debes et al. (2012) later found that asteroidal accretion on to WDs from both exterior and interior resonances also follow this trend.

## 6.2 Phase space extrapolation

The results of this study raise several questions. One important question is how robust our conclusions are to variations of the initial conditions. Because of the computational expense of running simulation ensembles for 5 Gyr with the BS algorithm, we could not perform a wider phase space exploration. However, we can guess how the outcomes would vary in other circumstances.

Varying $a_1$ would change the number of orbits completed over the MS; the effect is similar to varying $M_*$, Hence, Fig. 9 demonstrates the likely consequence: if the planets complete enough orbits, they will be prone to long-term MS instability. If the initial planetary separation is beyond the Hill stability limit, then this instability must be Lagrange instability. Otherwise, the type of instability is unrestricted.

The consequence of varying planetary eccentricities and mutual inclinations is less clear. What is clear is that planets on eccentric orbits will feature in post-MS systems. Observations suggest that if giant planets are formed at several au on circular orbits, they are unlikely to retain their primordial eccentricities on the MS. According to the Extrasolar Data Explorer,<fnref n="8" /> as of 2012 October 31, only

20.4 per cent of exoplanets with $M_p \ge 0.1 M_J$ and $a \ge 1$ au have $e < 0.1$. This percentage shrinks to 10.4 per cent for $e < 0.05$.<fnref n="9" /> If one were to include only planets in multi-planet systems, these percentages become 18.0 and 14.0 per cent, respectively. Therefore, if the planets we currently observe predominately survive until the post-MS, the significant majority will enter those phases with non-zero eccentricities. Nevertheless, it is of interest to determine if MS Lagrange instability can occur for planets formed on perfectly circular orbits. Therefore, we have performed 48 additional simulations with $e_1(0) = e_2(0) = 0.0$ for $M_*(0) = 3 M_\odot$ at different locations beyond the Hill stability limit. We can confirm that a few of those systems become Lagrange unstable. A detailed statistical comparison is best left to a more comprehensive phase space study.

We can estimate the dependence on planetary mass by extrapolating from Figs 9 and 15. The trend suggests that as the test particle limit is approached, where there is no mutual interaction between secondaries, the difference between the Hill and Lagrange boundaries might tend to zero faster than the Hill boundary itself. In the opposite limit, for two brown dwarfs and a more massive evolving primary, we expect the stability boundaries to widen even further than in our fiducial case. In these instances, instabilities might be common, and could represent a significant catalyst for WD pollution.

Adding more planets to model real systems such as HR 8799 introduces a significant complication, but is a viable avenue of a future study with our numerical code (Mustill & Veras, in preparation). If additional planets are of similar masses, then we expect them to represent destabilizing influences, particularly in the post-MS stages.

## 6.3 Comparison with radial velocity exosystems

The majority of known exoplanets are likely to be engulfed during the post-MS phase, as the population of planets beyond 10 au is largely unknown. Nevertheless, scaled versions of many known systems could represent genuine test cases for our model, because of the observed close separations of pairs of exoplanets. Furthermore, our findings of MS instability warrant a closer look at these systems.

We have compiled a list of all pairs of planets in the same system with $a > 1$ au each that were detected by radial velocity measurements. We obtained the data from the Extrasolar Planets Explorer on 2012 November 6. Using equations (1)–(5), and assuming minimum masses and coplanarity, we determined if the systems are currently Hill stable.

We found that six pairs of planets are not (HD 181433 c,d; HD 204313 b,d; 24 Sex b,c; BD +20 2457 b,c; HD 128311 b,c; HD 200964 b,c). Hence, these planets are likely to be protected by their proximity to mean-motion commensurabilities inside small islands of stability, as in e.g. Wittenmyer, Horner & Tinney (2012). This behaviour is reflective of some of our Hill unstable systems in Figs 9, 15 and 16. BD+20 2457 is notable because its mass is about $2.8 M_\odot$, and may be related to our $3 M_\odot$ simulations. However, both planets in the system are likely instead to be classified as brown dwarfs ($M_1 \approx 22.7 M_J$ and $M_2 \approx 13.2 M_J$). The correspondingly high mass ratios in the system, coupled with their close separations ($a_2/a_1 \approx 1.39$), perhaps intuitively suggest that instability should be imminent. Hence, the existence of such a system is

<footnote n="9">⁹ Admittedly, there is a bias towards fitting higher eccentricity values than the true values. The disparity worsens with sparser data, which is often associated with the most distant radial velocity planets.</footnote>

<footnote n="8">⁸ http://exoplanets.org/</footnote>
stark confirmation that Hill unstable systems can remain stable for long periods. The system's future prospects for stability on the MS would require detailed modelling.

Three pairs of planets are Hill stable, but have a separation ratio which exceeds the Hill stability limit by no more than 30 per cent (HD 37124 c,d; 47 Uma b,c; HD 183263 b,c). Based on our simulations, the close proximity of these systems to the Hill stability limit suggests that they might not be Lagrange stable for the remainder of their MS lifetimes. Testing this suggestion would require a detailed suite of long-term simulations for each system. If Lagrange instability scales strongly with planet/star mass ratio, then the planets in HD 183263 are perhaps in the greatest danger: the mass ratios in those systems are about 9.5 times as great as any mass ratio that we considered in our simulations.

Three pairs of Hill stable planets have separations exceeding the Hill stability limit by between 45 and 65 per cent (HD 108874 b,c; HD 159868 c,b; HD 10180 g,h), and two planet pairs (HD 4732 b,c; mu Ara b,c) have separations that are over twice the Hill limit. Although the Lagrange stability boundary is likely dependent on several variables, we have performed additional simulations as proof of concept to show that in at least one case, this boundary can extend out to twice the Hill stability boundary. Hence, HD 108874, HD 159868 and HD 10180 are not guaranteed to be Lagrange stable without further detailed analyses.

If the planets in any of these observed systems are not coplanar, then they are more likely to be Hill unstable (see the middle panel of Fig. 4). A mutual inclination of just 12:3 would render the planets in HD 37124 Hill unstable. At the opposite extreme, for the widely separated planets in HD 4732, a mutual inclination of 61:0 would be required.

## 6.4 Non-adiabatic mass-loss

Adiabaticity in the two-body problem with mass-loss is well defined (e.g. Veras et al. 2011). If we were to assume the two-body definition of adiabaticity for each of the planets in each of our simulations, then we can claim that at no time did our stable planetary systems (with $a_1(0) = 10$ au) approach a regime that featured non-adiabatic mass-loss. However, if two planets were stably orbiting at separations of hundreds of au on the MS, then the mutual planet-planet interaction coupled with stellar mass-loss could yield unpredictable evolution during the post-MS.

Voyatzis et al. (2013) have recently explored non-adiabatic mass-loss in the three-body problem. By using Lypaunov characteristic numbers to create dynamical stability maps, they found that non-adiabatic mass-loss can cause a stably interacting pair of planets to shift to a chaotic region of phase space. The resulting instability (manifested by escape or collision) may be latent, sometimes not appearing for a time that exceeds the duration of the mass-loss by a factor of tens. In cases where the outer planet escapes, distinguishing whether the instability was triggered by Lagrange instability, non-adiabatic mass-loss or both might require detailed follow-up simulations.

## 6.5 Sharpness of the Hill stability limit

No violations of the Hill stability limit have occurred in our simulations, despite us using the two-body approximation for the energy in the analytical formulation (equations 1–5). Nevertheless, we can estimate the error in semimajor axis as a result. Fig. 19 plots the mean and median energy error incurred by using the two-body approximation for the energy of the system. In all cases, the energy

<figure><img src="image_35.png" /><figcaption>**Figure 19.** The error of using the two-body approximation to model the total energy of the system. The black curve and filled dots indicate the mean energy error, the blue dashed curve and filled squares the median energy error. These curves are drawn for fiducial initial orbital parameters described in Section 5.4.1. The orange filled diamond and upward-pointing triangle represent the mean and median energy error for the $e_2 = 0.5$ simulation. The purple filled downward-pointing triangle and open dotted circle represent the simulations with the pair of Earth-mass planets. These errors correspond to the semimajor axis differences of the order of $10^{-3}$ au for $a \approx 10$ au.</figcaption></figure>

error is between 0.02 and 0.06 per cent. This range corresponds to semimajor axis differences of several $10^{-3}$ au for planets with $a \approx 10$ au. Our gradient of semimajor axis values sampled for our fiducial simulations exceeded this value in each case, helping to confirm why the analytic limit was not violated and suggesting that the unknown integration error was comparably small, if not smaller.

## 7 CONCLUSION

Architectures of planetary systems during each stellar phase may represent historical tracers of formation and presage future evolutionary instability and death. We have performed 5 Gyr simulations that consistently treat the dynamics of two massive planets and every phase of stellar evolution for a wide range of progenitor stellar masses (3–8 $M_\odot$). These computationally demanding simulations suggest that stable MS systems are in danger of future instability. The zone of danger is wide, reaching out to 163–178 per cent of the MS Hill stability limit and 123–137 per cent of the WD Hill stability limit for our low eccentricity [$e_1(0) = e_2(0) = 0.1$] simulations. The consequences for WD pollution may be significant: For example, the inner planet can be perturbed on to a highly eccentric orbit which takes the planet close to the WD or hits the star directly.

## ACKNOWLEDGEMENTS

We thank the referee for carefully reading the manuscript and assistive and helpful suggestions. We also thank Jay Farihi for clarifying the observational biases, J. Richard Donnison for helping us confirm a typo in the previous literature about Hill stability, and
Boris T. Gänsicke and Mukremin Kilic for useful discussions. AJM is funded by the Spanish National Plan of R&D grant AYA2010-20630, 'Planets and stellar evolution'. AB acknowledges the support of the ANR-2010 BLAN-0505-01 (EXOZODI). This work made use of facilities funded by the European Union through ERC grant number 279973.

## REFERENCES

Adamów M., Niedzielski A., Villaver E., Nowak G., Wolszczan A., 2012, ApJ, 754, L15

Alcock C., Fristrom C. C., Siegelman R., 1986, ApJ, 302, 462

Barnes R., Greenberg R., 2006, ApJ, 647, L163

Barnes R., Greenberg R., 2007, ApJ, 665, L67

Bear E., Soker N., 2012, ApJ, 749, L14

Becklin E. E., Farihi J., Jura M., Song I., Weinberger A. J., Zuckerman B., 2005, ApJ, 632, L119

Binney J., Tremaine S. (eds), 2008, Galactic Dynamics, 2nd edn. Princeton Univ. Press, Princeton, NJ

Bonsor A., Wyatt M., 2010, MNRAS, 409, 1631

Bonsor A., Mustill A., Wyatt M., 2011, MNRAS, 414, 930

Bonsor A., Augereau J.-C., Thebault P., 2012, A&A, 548, A104

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

Charpinet S. et al., 2011, Nat, 480, 496

Darriba L. A., Maffione N. P., Cincotta P. M., Giordano C. M., 2012, preprint (arXiv:1205.0875)

Debes J. H., Sigurdsson S., 2002, ApJ, 572, 556

Debes J. H., Hoard D. W., Wachter S., Leisawitz D. T., Cohen M., 2011, ApJS, 197, 38

Debes J. H., Walsh K. J., Stark C., 2012, ApJ, 747, 148

Donnison J. R., 2006, MNRAS, 369, 1267

Donnison J. R., 2009, Planet. Space Sci., 57, 771

Donnison J. R., 2010a, Planet. Space Sci., 58, 1169

Donnison J. R., 2010b, MNRAS, 406, 1918

Donnison J. R., 2011, MNRAS, 415, 470

Farihi J., Jura M., Zuckerman B., 2009, ApJ, 694, 805

Farihi J. et al., 2012a, MNRAS, 421, 1635

Farihi J. et al., 2012b, MNRAS, 424, 464

Farihi J. et al., 2012c, MNRAS, 424, 519

Farmer A. J., Goldreich P., 2006, Icarus, 180, 403

Gänsicke B. T., Marsh T. R., Southworth J., Rebassa-Mansergas A., 2006, Sci, 314, 1908

Gänsicke B. T., Marsh T. R., Southworth J., 2007, MNRAS, 380, L35

Gänsicke B. T., Koester D., Marsh T. R., Rebassa-Mansergas A., Southworth J., 2008, MNRAS, 391, L103

Gänsicke B. T., Koester D., Farihi J., Girven J., Parsons S. G., Breedt E., 2012, MNRAS, 424, 333

Geier S., Edelmann H., Heber U., Morales-Rueda L., 2009, ApJ, 702, L96

Georgakarakos N., 2008, Celest. Mech. Dyn. Astron., 100, 151

Girven J., Brinkworth C. S., Farihi J., Gansicke B. T., Hoard D. W., Marsh T. R., Koester D., 2012, ApJ, 749, 154

Gladman B., 1993, Icarus, 106, 247

Gu P.-G., Lin D. N. C., Bodenheimer P. H., 2003, ApJ, 588, 509

Guillot T., Burrows A., Hubbard W. B., Lunine J. I., Saumon D., 1996, ApJ, 459, L35

Hadjidemetriou J. D., 1963, Icarus, 2, 440

Hadjidemetriou J. D., 1966, Z. Astrophys., 63, 116

Holberg J. B., Sion E. M., Oswalt T., McCook G. P., Foran S., Subasavage J. P., 2008, AJ, 135, 1225

Hurley J. R., Pols O. R., Tout C. A., 2000, MNRAS, 315, 543

Iben I., Jr, Renzini A., 1983, ARA&A, 21, 271

Jarosik N. et al., 2011, ApJS, 192, 14

Jura M., 2003, ApJ, 584, L91

Jura M., 2008, AJ, 135, 1785

Jurić M., Tremaine S., 2008, ApJ, 686, 603

Ketchum J. A., Adams F. C., Bloch A. M., 2013, ApJ, 762, 71

Kholshevnikov K. V., Kuznetsov E. D., 2007, Sol. Syst. Res., 41, 265

Kilic M., von Hippel T., Leggett S. K., Winget D. E., 2005, ApJ, 632, L115

Klein B., Jura M., Koester D., Zuckerman B., Melis C., 2010, ApJ, 709, 950

Koester D., Girven J., Gänsicke B. T., Dufour P., 2011, A&A, 530, A114

Kopparapu R. K., Barnes R., 2010, ApJ, 716, 1336

Kratter K. M., Perets H. B., 2012, ApJ, 753, 91

Kubala A., Black D., Szebehely V., 1993, Celest. Mech. Dyn. Astron., 56, 51

Kunitomo M., Ikoma M., Sato B., Katsuta Y., Ida S., 2011, ApJ, 737, 66

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

Lebzelter T., Uttenhalter S., Busso M., Schultheis M., Aringer B., 2012, A&A, 538, A36

Lee J. W., Kim S.-L., Kim C.-H., Koch R. H., Lee C.-U., Kim H.-I., Park J.-H., 2009, AJ, 137, 3181

Lee B.-C., Mkrtsichian D. E., Han I., Park M.-G., Kim K.-M., 2012, A&A, 548, A118

Lee B.-C., Han I., Park M.-G., 2013, A&A, 549, A2

Maercker M. et al., 2012, Nat, 490, 232

Mullally F., Winget D. E., De Gennaro S., Jeffery E., Thompson S. E., Chandler D., Kepler S. O., 2008, ApJ, 676, 573

Mullally F., Reach W. T., De Gennaro S., Burrows A., 2009, ApJ, 694, 327

Mustill A. J., Villaver E., 2012, ApJ, 761, 121

Mustill A. J., Wyatt M. C., 2012, MNRAS, 419, 3074

Namouni F., 2005, AJ, 130, 280

Namouni F., 2007, ApJ, 659, 1505

Namouni F., 2012, Ap&SS, 317

Naoz S., Farr W. M., Lithwick Y., Rasio F. A., Teyssandier J., 2011, Nat, 473, 187

Nordhaus J., Spiegel D. S., 2012, preprint (arXiv:1211.1013)

Omarov T. B., 1962, Izv. Astrofiz. Inst. Acad. Nauk. KazSSR, 14, 66

Parriott J., Alcock C., 1998, ApJ, 501, 357

Portegies Zwart S., 2013, MNRAS, 429, L45

Quillen A. C., Faber P., 2006, MNRAS, 373, 1245

Raymond S. N., Barnes R., Veras D., Armitage P. J., Gorelick N., Greenberg R., 2009, ApJ, 696, L98

Reach W. T., Kuchner M. J., von Hippel T., Burrows A., Mullally F., Kilic M., Winget D. E., 2005, ApJ, 635, L161

Sato B. et al., 2012, PASJ, 64, 135

Sato B. et al., 2013, ApJ, 762, 9

Schröder K.-P., Cuntz M., 2005, ApJ, 630, L73

Setiawan J., Klement R. J., Henning T., Rix H.-W., Rochau B., Rodmann J., Schulze-Hartung T., 2010, Sci, 330, 1642

Sigurdsson S., Richer H. B., Hansen B. M., Stairs I. H., Thorsett S. E., 2003, Sci, 301, 193

Silvotti R. et al., 2007, Nat, 449, 189

Smith A. W., Lissauer J. J., 2009, Icarus, 201, 381

Spiegel D. S., Madhusudhan N., 2012, ApJ, 756, 132

Vassiliadis E., Wood P. R., 1993, ApJ, 413, 641

Veras D., 2007, Celest. Mech. Dyn. Astron., 99, 197

Veras D., Armitage P. J., 2004, Icarus, 172, 349

Veras D., Evans N. W., 2013a, MNRAS, 430, 403

Veras D., Evans N. W., 2013b, Celest. Mech. Dyn. Astron., 115, 123

Veras D., Ford E. B., 2009, ApJ, 690, L1

Veras D., Ford E. B., 2010, ApJ, 715, 803

Veras D., Moeckel N., 2012, MNRAS, 425, 680

Veras D., Tout C. A., 2012, MNRAS, 422, 1648

Veras D., Wyatt M. C., 2012, MNRAS, 421, 2969

Veras D., Wyatt M. C., Mustill A. J., Bonsor A., Eldridge J. J., 2011, MNRAS, 417, 2104

Villaver E., 2011, AIP Conf. Ser. Vol. 1331, Planetary Systems Beyond The Main Sequence. Am. Inst. Phys., New York, p. 21

Villaver E., Livio M., 2007, ApJ, 661, 1192

Villaver E., Livio M., 2009, ApJ, 705, L81

Voyatzis G., Hadjidemetriou J. D., Veras D., Varvoglis H., 2013, MNRAS, preprint (arXiv:1301.5441)

Wickramasinghe D. T., Farihi J., Tout C. A., Ferrario L., Stancliffe R. J., 2010, MNRAS, 404, 1984
Wisdom J., 1980, AJ, 85, 1122

Wittenmyer R. A., Horner J., Tinney C. G., 2012, ApJ, 761, 165

Wolszczan A., 1994, Sci, 264, 538

Wolszczan A., Frail D. A., 1992, Nat, 355, 145

Zakamska N. L., Tremaine S., 2004, AJ, 128, 869

Zuckerman B., Becklin E. E., 1987, Nat, 330, 138

Zuckerman B., Koester D., Reid I. N., Hünsch M., 2003, ApJ, 596, 477

Zuckerman B., Koester D., Melis C., Hansen B. M., Jura M., 2007, ApJ, 671, 872

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