Long-term evolution of three-planet systems to the post-main sequence and beyond Free

Abstract

1 INTRODUCTION

The study of planetary systems as their host stars evolve beyond the main sequence (MS) is now becoming an important area of research. The detection of the pollution of the atmospheres of white dwarfs (WDs) by comets and asteroids has allowed the composition of extrasolar planetesimals to be determined (e.g. Zuckerman et al. 2007; Gänsicke et al. 2012; Xu et al. 2013). Binaries comprising a WD and a low-mass substellar companion provide a calibration for atmospheric models of brown dwarfs (BDs) and giant planets (Pinfield et al. 2006; Steele et al. 2009; Day-Jones et al. 2011). The population of planetary, sub-stellar and low-mass stellar companions that survive the star's giant phases can help to constrain poorly understood aspects of tidal interactions (Rasio et al. 1996; Villaver & Livio 2009; Mustill & Villaver 2012) and common envelope evolution (Schreiber & Gänsicke 2003; Zorotovic et al. 2010; Casewell et al. 2012; Mustill et al. 2013; Nordhaus & Spiegel 2013; Portegies Zwart 2013).

Many physical processes are involved in leading planetary systems from the MS through to the star's WD stage. Planets and planetesimals see their orbital radii increased by stellar mass loss, which in extreme cases can also cause eccentricity excitation and even ejection of distant planets and comets (Omarov 1962; Hadjidemetriou 1963; Parriott & Alcock 1998; Veras et al. 2011). This mass-loss can also affect the system's many-body dynamics, as planet:star mass ratios increase (Debes & Sigurdsson 2002; Bonsor, Mustill & Wyatt 2011; Debes, Walsh & Stark 2012; Adams, Anderson & Bloch 2013; Veras et al. 2013a; Voyatzis et al. 2013). Planets closer to the star will feel the effects of photoevaporation (Villaver & Livio 2007) and tidal forces (Villaver & Livio 2009; Kunitomo et al. 2011; Mustill & Villaver 2012; Nordhaus & Spiegel 2013), while smaller bodies may be sublimated or destroyed by stellar wind drag (Jura 2008; Bonsor & Wyatt 2010; Jura & Xu 2010). Achieving a full understanding of the behaviour of planetary systems around evolving stars is therefore a complex problem.

Around 1 in 5 of all known planetary systems has more than one known planet,¹ a fraction that may increase as more difficult-to-detect companions to existing planets are found. These systems range from compact systems such as Kepler-11 (Lissauer et al. 2011), with 6 known planets within 0.5 au of the host star, to wider ones such as HR 8799, with four known planets beyond 10 au (Marois et al. 2010). Many-body dynamics is therefore important for a significant fraction and wide diversity of planetary systems. Moreover, many multiplanet systems are dynamically ‘packed’, i.e. planets cannot be much closer together, nor can more planets be introduced into gaps between them, without the systems being violently unstable (Barnes & Quinn 2004; Barnes & Raymond 2004; Raymond & Barnes 2005; Fang & Margot 2013). The seminal studies of Duncan & Lissauer (1998) and Debes & Sigurdsson (2002) showed that mass-loss can destabilize systems that are stable in the absence of mass-loss, or induce instability more rapidly in unstable systems, as planet–planet interactions strengthen when the star loses mass.

Debes & Sigurdsson (2002) studied both two- and three-planet systems; however, they did not self-consistently study the dynamical evolution of systems before, during and after mass-loss. In Veras et al. (2013a, hereafter Paper I), we set out to rectify this, incorporating the stars’ mass-loss history into N-body integrations and following the evolution of two-planet systems from the zero-age MS through to WD cooling ages of up to 5 Gyr.

In this paper, we extend our previous work to systems of three planets, where the stability properties are more complex and less well understood. In common with our earlier work, we focus on stellar masses of 3 M_⊙ and above, for computational reasons, to giant planets, as these are the easiest to detect, and to planets beyond 10 au, to avoid complicating the evolution with tides, photoevaporation, or common envelope evolution.

The paper is laid out as follows. In Section 2, we discuss semi-analytical insights into the expected stability of planetary systems at different stellar evolutionary phases. In Section 3, we describe our N-body simulations of three-planet systems. In Section 4, we discuss our results. We conclude in Section 5.

2 STABILITY ESTIMATES

Because the parameter space for the initial conditions of arbitrary three-planet systems is so large, we seek some analytical or semi-analytical guidance to help to inform our choice of initial conditions for our numerical integrations. The problem we face is one of determining the stability of a given planetary system orbiting a given star, over a given time-scale, in particular those critical time-scales determined by the stellar mass such as the MS lifetime.

The stability of two-planet systems throughout the lifetimes of their 3–8 M_⊙ host stars was studied in detail in Paper I. In this case, there is an analytical result on stability, the Hill stability criterion, stating under what conditions a system may see a collision between the two planets (Gladman 1993). However, in Paper I and Veras & Mustill (2013), we showed that the Hill criterion can underestimate instability which can occur on long time-scales comparable to a star's MS lifetime, since it does not account for Lagrange instability, which typically results in the ejection of a planet. Moreover, the four-body problem is complex enough that often one cannot apply two-planet stability criteria pairwise amongst the three planets. In fact, three planets which are well-separated enough so that adjacent pairs of planets satisfy the two-body stability criteria may eventually suffer collisions and close encounters. The reasons for this behaviour are due to the properties of Kolmogorov–Arnol'd–Moser tori embedded in a chaotic sea, and are explored in detail in Shikita, Koyama & Yamada (2010).

2.1 Instability in MS planetary systems

Figure 1.

Approximate stability behaviour of three-planet systems as a function of stellar mass and separation on the MS. The planet mass and the innermost semimajor axis are fixed at 1 M_J and 1 au, respectively. Systems below the black solid line will be unstable on the MS. Systems in the tiny strip between the black solid and the black dotted lines will become unstable at some point between the end of the MS and the end of the AGB. Systems below the red solid line will be unstable during the WD phase, assuming that the star formed at the big bang; systems above this line must be stable for the age of the Universe. The red dotted and dashed lines show the stability boundaries for WD cooling ages of 1 and 5 Gyr. This plot demonstrates that after the star has become a white dwarf, many more systems can become unstable.

Open in new tabDownload slide

2.2 Instability in post-MS planetary systems

We can now make estimates of planetary systems’ stability as a function of stellar mass and evolutionary stage. We use the sse code (Hurley, Pols & Tout 2000) to find the stellar lifetimes and mass evolution of stars with metallicity Z = 0.02 and Reimers’ η = 0.5. The sse code applies Reimers mass-loss at early evolutionary stages, but during the AGB it applies the semi-empirical mass-loss rate of Vassiliadis & Wood (1993), reaching a maximum during the superwind of 1.36 × 10⁻⁹(L_*/L_⊙) M_⊙ yr⁻¹, where L_* is the stellar luminosity. An example of the mass-loss rate history on the AGB is shown in Fig. 2, for an 8 M_⊙ star. Mass-loss from the lower mass stars is less rapid, but still attains a maximum of 3 × 10⁻⁵ M_⊙ yr⁻¹ for a 3 M_⊙ star. The fraction of mass lost by the WD stage is 75 per cent for a 3 M_⊙ star and 82 per cent for an 8 M_⊙ star.

Figure 2.

Mass-loss rate during the AGB for our 8 M_⊙ stellar model.

Open in new tabDownload slide

Up to the end of the AGB, we take t_PMS to be the elapsed lifetime of the star. For the WD phase, we assume that the system clock is reset following mass-loss on the AGB, and take the relevant age to be the cooling age of the WD; this may underestimate the extent of instability if systems become mildly excited on the MS. In Fig. 1, we show the stability limits estimated from equations (5) and (6), for the MS lifetime, the end of the AGB, WD cooling ages of 1 and 5 Gyr, and a total age of 13.7 Gyr. From the figure, we make the following predictions.

• On the MS, planets with Δ ≲ 9r_H will be susceptible to instability (black solid line). This separation is remarkably insensitive to stellar mass, and hence the very different MS lifetimes, since the stellar lifetime only enters equation (5) logarithmically.

• The relatively short time between the end of the MS and the end of the AGB will not allow many more systems to become unstable (black dotted line).

• After the star has become a WD, many more systems can become unstable (red solid line), which is largely attributable to the power-law dependence of the critical separation on the final:initial stellar mass ratio (the increase in Hill sphere size). More systems become unstable around higher mass progenitor stars, because the mass ratio becomes more unequal. The stability boundary increases to around 16r_H for 8 M_⊙ stars. Almost all of the increase is accounted for by the change in the Hill radius, with the power-law dependence on |$M_\ast ^\mathrm{f}/M_\ast$| in equation (6) accounting for an increase in critical separation to ∼14r_H at 8 M_⊙, the logarithmic terms making up the difference.

• The 1 and 5 Gyr lines (red dotted and dashed) illustrate that most WD instabilities will occur early in the WD lifetime.

These stability limits exhibit a dependence on both planet mass and semimajor axis. As either mass or semimajor axis is reduced, the critical separation measured in Hill radii increases. This dependence is however weak: For a 5 M_⊙ MS star, reducing the innermost planet's semimajor axis from 1000 to 0.1 au, or planet mass from 10 M_J to 10⁻⁴ M_J, each results in an increase in the critical separation of about 4 Hill radii. We discuss this further in Section 4.1.

3 N-BODY SIMULATIONS

Our choice for the inner planet's orbital radius should ensure that its orbit will be unaffected by tidal forces on either the red giant branch (RGB) or AGB. Around lower mass stars, Villaver & Livio (2009), Mustill & Villaver (2012) and Nordhaus & Spiegel (2013) showed that giant planets initially at 10 au will not feel significant tidal forces. Note however that the fate of planets around the more massive AGB stars (>5 M_⊙) has not been studied; extrapolating the results of Mustill & Villaver (2012) and Nordhaus & Spiegel (2013) to higher masses suggests that Jovian planets at 10 au may feel the star's tidal forces on the AGB. Also, if scattering brings planets on to orbits with smaller pericentres, they will feel tidal forces. In this paper, we ignore any tidal effects, leaving the coupling of tidal and N-body dynamics for future work.

We run two sets of simulations, one with 1 M_J planets and one with 10 M_J planets. The planets’ densities used to calculate the Roche limit for tidal disruption are 1.3 and 13 g cm⁻³, respectively. In all cases, the inner planet is initially placed at 10 au. The second planet is placed at a separation Δ = 3, 3.5, 4, …, 17.5, 18 in units of the inner planet's Hill radius, a spacing that corresponds to semimajor axis ratios of 1.10-1.85 for the 1 M_J planets, depending on stellar mass, and 1.31-2.84 for the 10 M_J planets. The third planet is then placed at the same semimajor axis ratio from the second planet. As discussed above, at a separation of 18 Hill radii, systems should be stable throughout the WD lifetime. Our closest separation of 3 Hill's radii is well within the MS stability limit, and very closely spaced systems may even be destabilized while the gas disc is still present (Lega, Morbidelli & Nesvorný 2013); however, we choose to go so close to see whether mass-loss on the AGB can trigger a second round of instabilities in systems than have already relaxed during their MS lifetime. We chose the closest separation post hoc by finding the first separation for which all systems integrated were unstable on the MS. Planets are initially on circular orbits, measured in Jacobi coordinates, and are assigned small inclinations of up to 1° in order to avoid unrealistic numbers of planet–planet collisions. We integrate each system up to total age of 5 Gyr, encompassing the entire MS evolution and several Gyr of WD cooling. Planets’ orbital elements are recorded every 1 Myr, and each instance of a planet being lost from a system is recorded, where the loss can be due to ejection beyond 10⁴ au, or collision with another planet or the star. We use the loss of one or more planets as a proxy for saying a system is unstable.

In total, with six values of M_*, 31 values of Δ, and 8 systems with randomly chosen orbital phases at each separation, there are 1488 systems for our 1 M_J simulations. For the 10 M_J simulations, we only considered the lowest stellar mass, with 248 systems integrated. A sample size of the order of 1000 allows the distribution of final orbital parameters to be characterized with good accuracy (Chatterjee et al. 2008).

We now describe in detail the outcomes of our 1 M_J integrations (Section 3.1), then the outcomes of the 10 M_J integrations (Section 3.2).

3.1 Fiducial case: M_pl=1 M_J

First, we discuss the instabilities in our 1 M_J systems. The times at which a planet is lost, as well as the nature of the loss – be it through collision with the star or another planet, or ejection from the system – are plotted in Fig. 3 as a function of initial separation for each stellar mass. We also mark on this plot the critical separations from the two-planet Hill stability criterion on the MS (vertical solid line), as well as the estimates for three-planet stability limits on the MS and WD phases (vertical dashed and dotted lines). Important time-scales are marked on the figure with horizontal lines: the end of the MS lifetime of the star, the formation of the WD and the 5 Gyr duration of our integrations. We treat stellar evolution in three steps: first the MS; then the post-MS up to the formation of the WD, during which interval the stellar radius becomes very large and the star loses most of its mass; and finally the WD stage. A breakdown of the numbers of planets lost, by type of loss, stellar evolutionary stage and stellar mass, is given in Table 1 and displayed graphically in Fig. 4.

Figure 3.

Instability times in systems of three 1 M_J planets. Symbols are plotted every time a planet is lost; symbol types show whether planets were lost due to ejection or due to collision with another planet or the star. Horizontal lines mark the end of the MS, the beginning of the WD cooling track, and the 5 Gyr duration of the integrations. Solid vertical lines mark the two-planet Hill stability limit on the MS, and the dashed and dotted lines mark the three-planet stability limits for MS and WD stars determined in Fig. 1. Important MMRs are marked with triangles.

Open in new tabDownload slide

Crudely, the outcomes of instability as a function of stellar age can be summarized as follows. The systems around MS stars see planets lost due to collisions with the star and between planets, as well as ejections; all these loss mechanisms contribute in roughly equal measure. Around post-MS stars, from the subgiant (SG) to the AGB, planets are lost overwhelmingly through collision with the expanded star, with very few ejections or planet–planet collisions. Around WDs, planets are lost primarily through ejections, with few being lost due to collisions. We now describe our simulations in more detail, going through each stage of stellar evolution.

3.1.1 The main sequence

Instability on the MS in our integrations occurs for systems separated by up to 9r_H, as predicted by equation (5). At all these separations, planets are lost through ejection and collision with the star or another planet. The two-planet Hill stability criterion predicts that planetary collisions cannot occur at separations beyond ∼4.6r_H, and hence fails to accurately describe the stability properties of these three-planet systems.

We also see from the integrations that, at early times on the MS, instability is dominated by planet–planet collisions, while at later times collisions with the star and ejections become more important. The number of planets lost on the MS decreases with stellar mass, from 105 for the 3 M_⊙ case to 66 for the 8 M_⊙ case, as a consequence of the shorter MS lifetime. Meanwhile, the fraction of ejections³ falls from |$44.8^{+4.9}_{-4.7}$| per cent of planets lost to |$21.2^{+5.9}_{-4.2}$| per cent, for the 3 M_⊙ and 8 M_⊙ stars respectively, due to the star's deeper potential well.

Systems may lose more than one planet on the MS. Fig. 5 shows whether each planet lost was the first, second or third from that system, for the 3 M_⊙ star. Often, the loss of the second planet occurs considerably later than the loss of the first. Considering all stellar masses, of the 437 systems losing a planet on the MS, a further 60 (⁠|$13.7^{+1.8}_{-1.5}$| per cent) went on to lose a second planet on the MS as well. We note that in no cases did a star lose a third planet on the MS (<0.2 per cent); i.e. scattering events never resulted in the destruction or ejection of all of a system's planets.

Figure 4.

Number of planets lost in the three-1 M_J runs, broken down by initial stellar mass, type of loss and stellar evolutionary phase.

Open in new tabDownload slide

Figure 5.

Time of instability as a function of planet separation, showing whether each planet lost is the first, second or third to be lost from the system. Ancillary markings are as in Fig. 3.

Open in new tabDownload slide

3.1.2 Early post-MS evolution

After the end of the MS, the stellar radius increases to much larger values, with the 8 M_⊙ star attaining 1 au on the RGB and 8 au at the tip of the AGB. The consequence of this is a large increase in the number of planets being engulfed by the star: between the end of the MS and the end of the AGB, an additional 150 planets were lost due to stellar collision (⁠|$90.9^{+1.8}_{-2.7}$| per cent of all losses at this stage), compared to 13 more ejections (⁠|$7.9^{+2.6}_{-1.6}$| per cent) and two planet–planet collisions (⁠|$1.2^{+1.6}_{-0.4}$| per cent). In contrast to planet loss on the MS, planets lost during these stages were primarily from systems that had already lost a planet: 129 of the planets lost (⁠|$78.2^{+2.8}_{-3.5}$| per cent) were the system's second, and 9 (⁠|$5.4^{+2.4}_{-1.2}$| per cent) were the system's third. Destabilization of new systems is rare because the short duration of the post-MS phases means that there is not much time for hitherto stable systems to experience an instability, as predicted from equation (6) and Fig. 1. A system may lose all its planets either when the sole survivor has a small pericentre that results in its engulfment by the expanding stellar envelope, or when it has a large apocentre which makes the stellar mass loss non-adiabatic and results in a ‘Great Escape’ scenario (Veras et al. 2011). The latter occurs infrequently in our integrations however, with only four of the instances of third-planet loss being through this mechanism.

Mass-loss on the AGB does not typically result in immediate instability and loss of a planet: the results of AGB mass-loss are played out during the host star's WD lifetime, as described next.

3.1.3 The white dwarf stage

The range of separations vulnerable to instability on the WD stage obtained from the numerical integrations is broadly in line with that predicted from equation (6) and Fig. 1, which is shown as the dotted vertical line in Figs 3 and 5. However, there is variation about this value. In particular, particles just interior to the 2:1 resonance have noticeably shorter lifetimes than their neighbours.

The stellar mass loss just before the formation of the WD causes both instability in previously stable systems, and renewed instability in systems that had previously lost a planet. Instabilities occur at about the same rate in systems that have and have not already seen an instability before the end of their AGB evolution: of 1024 three-planet systems surviving to the formation of the WD, 496 (48.4 ± 1.6 per cent) experienced a subsequent instability, while of 275 two-planet systems at the formation of the WD, 149 (54.2 ± 3.0 per cent) experienced a subsequent instability.

However, the outcome of early instabilities can play a role in the systems’ future evolution. The effects of MS instabilities on the subsequent fates of planetary systems is shown in Fig. 6. For systems which lose one planet during the MS, the average a and e of their surviving planets for the remainder of the MS are shown. The subsequent fate of the systems is also shown: whether they lose a second planet during the SG through to early AGB phases, during the thermally pulsating AGB (TPAGB) phase, or as a WD, or whether they subsequently remain stable. The orbital elements of surviving planets following the first instability clearly affects the subsequent evolution: during the post-MS stages, systems with one planet on an orbit with a small semimajor axis lose it due to engulfment, as the stellar radius expands, the higher eccentricity systems being swallowed first when the stellar radius is smaller. During the WD phase, the planets lost come from more closely packed systems with a higher e for a given a on the MS.

Figure 6.

For systems which lose precisely one planet on the MS, the average a and e of the two surviving planets for the remainder of their MS lifetime is shown. Symbols show their subsequent fate: whether they lose a second planet during early post-MS evolution before the TPAGB, during the TPAGB or the WD stage, or remain stable. The vertical solid line shows a semimajor axis of 5 au. Dashed lines show an apocentre of 10 au and a pericentre of 15 au. All stellar masses are considered, with planets of 1 M_J.

Open in new tabDownload slide

The distribution of a and e in Fig. 6 is similar to that seen in other scattering studies (e.g. Chatterjee et al. 2008). A large population of planets ends up at ∼4 au. This implies a change in orbital energy, relative to an initial orbit at 10 au, just greater than that required to unbind a co-orbital planet, which would move a surviving planet in to 5 au. The rest of the population falls in two tails: one with an apocentre at 10 au, and another, less well-defined, with a pericentre at ∼15 au. These are respectively the initial semimajor axis of the inner planet, and the average initial semimajor axis of all planets. The region in between is underpopulated, since planets here in general are still experiencing encounters and hence systems are unstable (Chatterjee et al. 2008).

Instabilities during the WD phase tend to be delayed. Fig. 7 shows a histogram of the cooling ages of the WD at which planets are lost, summed across all stellar masses but broken down by the result of the instability. The ejections, forming the bulk of the instabilities, mostly occur after 1 Myr and continue until the end of the integrations. It appears that ending the integrations at 5 Gyr has truncated the tail of the distribution and instabilities would still occur, albeit less frequently, at still later ages. The planet–planet collisions begin much earlier, and their rate per logarithmic time unit is roughly constant. Collisions with the WD begin much later, at 10 Myr, and there is a hint of bimodality in the distribution of the logarithm of their ages. There is also a difference in the times at which the systems that did or did not experience pre-WD planet losses experience instability, with a median of 10^7.9 yr for the first WD instability in hitherto stable systems, versus 10^8.4 yr for the first WD instability in previously unstable systems. A KS test shows that the two distributions are significantly different (p = 5 × 10⁻⁸). Hence, systems which had previously experienced instability see their WD instabilities occurring a little later than the stable systems.

Figure 7.

Cooling ages at which WD instabilities occur in systems of three 1 M_J planets, and their instability outcome.

Open in new tabDownload slide

In Fig. 8, we show the final orbital elements of all systems at 5 Gyr, as plots of eccentricity e and pericentre q against semimajor axis a. Planets in stable systems, marked in red, have semimajor axes in the range 40–150 au, as expected from the adiabatic expansion of their initial orbits, and low eccentricities which have been somewhat excited by planet–planet interactions. Unsurprisingly, planets in unstable systems have wildly different orbital elements, with semimajor axes up to several thousands of au and eccentricities up to 1. The distribution of elements in unstable systems is similar to that in Fig. 6, but the structure is somewhat washed out by the different amounts of orbit expansion on the AGB. A particularly dense concentration can be found with semimajor axes of 20–30 au, smaller than the minimum adiabatically expanded initial value. Few semimajor axes are lower than these values. This is in accordance with the finding of Chatterjee et al. (2008) that scattering does not reduce semimajor axes below about 40 per cent of their initial value. Nonetheless, some of these planets can approach closer to the star, since they have very high eccentricity and hence very small pericentres. The survivors’ pericentres at 5 Gyr are shown in the lower panel of Fig. 8, where we see that some planets come within 10 au of their host star, a region expected to be cleared during the AGB phase (Mustill & Villaver 2012). The structure of the a-e distribution will be somewhat sensitive to initial conditions, since the semimajor axis of the inner planet sets the scale for the major concentrations of points, as discussed for Fig. 6 above. However, the existence of a hard inner cut-off in a will be robust, since if planets are started initially too close to the star, they will be swallowed by the stellar envelope. We note that three planets classed as being in ‘stable’ systems clearly belong to the unstable population from their orbital elements, with a ≈ 30, 350 and 400 au. These planets all belong to one system, which underwent a scattering event at ∼4.975 Gyr, just before the end of the integration. As we have no way of telling what the outcome of this instability will be, we choose to keep the system classified as ‘stable’.

Figure 8.

Final (at 5 Gyr of total evolution, i.e. more than 4 Gyr into the WD stage) eccentricities (top) and pericentres (bottom) as a function of final a. Planets in systems that have remained stable for 5 Gyr are marked in red; those in systems that have experienced an instability, in black.

Open in new tabDownload slide

We have seen that systems may lose several or even all three planets during the course of stellar evolution. Table 3 shows the number of surviving planets at 5 Gyr in all our simulation sets. Nearly all WDs that initially had systems of giant planets on wide orbits will keep at least one. In general, the more massive a star, the more systems will be destabilized and the fewer planets will be found around WDs. We discuss the observational consequences of this in Section 4.

The planetary systems of WDs are more unstable than those around MS stars for one simple reason: the increase in the planet:star mass ratio, causing a decrease in the planetary separation when measured in Hill radii. Several terms in equation (6) have a smaller effect. The decrease in stellar mass directly causes an increase in the orbital time-scale, and hence any process that takes a specified number of orbits would take a larger amount of real time. The increase in the planets’ semimajor axes causes a further increase in the orbital period, compounding this effect. However, since time and semimajor axis only affect the critical separation logarithmically, the power-law dependence of the Hill radius on the stellar mass is the dominant effect.

As we showed in Paper I, systems of two planets are also susceptible to instability following mass-loss. A notable difference when considering three-planet systems is the potential to undergo two rounds of instability: a first round on the MS which causes the loss of one planet while the survivors are left in a stable state, and then a second round as a two-planet system as a WD. Indeed, systems that had been already destabilized on the MS were just as likely in our simulations to then lose a second planet as a WD than those systems that survived from the MS intact. While in both cases ejections dominate as an outcome of WD instability, in our three-planet simulations we see far more collisions with the star, and far fewer planet–planet collisions, than in Paper I. Rather than being a fundamental change in the underlying dynamics, this likely reflects our attempts to improve the realism of our simulations: by assigning the planets a small initial inclination we greatly reduced the chances of a planet–planet collision compared to the coplanar simulations of Paper I, while increasing the radius for stellar collisions from the physical WD radius to the Roche limit has greatly increased the number of planets classed as ‘colliding with the star’. The precise fate of such bodies will be discussed below.

3.2 High planet mass case, M_pl = 10 M_J

We now discuss the case of three 10 M_J planets. The times and types of instability are shown as a function of initial separation in Fig. 9, while the planet losses by evolutionary stage are listed in Table 4. Ejections dominate the outcomes during the MS (62 of 81 planets lost, |$76.5^{+4.0}_{-5.3}$| per cent). On the MS, systems are unstable out to ∼6.5r_H, while the semi-analytical estimate from equation (5) predicts that the stability boundary should be at 8r_H. There is also a small region of instability associated with the 3:1 MMR at ∼10.5r_H. Hence, the predictions of the semi-analytical estimates are not so good as in the lower planet mass case, possibly because the formulae have not been calibrated at such high mass ratios. In contrast, the Hill criterion applied pairwise predicts that planet–planet collisions can occur out to 5r_H, which is indeed borne out by the integrations.

Figure 9.

As Fig. 3 for the case of three 10 M_J planets.

Open in new tabDownload slide

As in the case of the 1 M_J systems, planet losses between the end of the MS and the end of the AGB are primarily collisions with the expanded stellar envelope (19 out of 26, |$73.1^{+6.8}_{-10.2}$| per cent). Again, losses primarily occur in systems that had previously lost a planet, with three planets lost being the systems’ first (⁠|$11.5^{+9.2}_{-3.6}$| per cent), 17 being their second (⁠|$65.4^{+7.8}_{-10.1}$| per cent) and 6 their third (⁠|$23.1^{+10.0}_{-6.2}$| per cent).

When compared with the 1 M_J case, destabilization of systems round WDs remains common. Numerically, we find that systems can be unstable out to a separation of 11r_H, close to the 12r_H predicted from equation (6). However, this may be due to the enhanced destabilizing effects of the 3:1 MMR.

The nature of planet loss around WD stars in the high planet mass case is not the same as for lower planet masses. Of 182 three-planet systems entering the WD stage, 61 experienced instability (⁠|$33.5^{+3.7}_{-3.3}$| per cent). However, only 4 of the 31 systems of two planets at the start of the WD stage lost a further planet (⁠|$12.9^{+8.4}_{-3.9}$| per cent). This contrasts with the overall results from the lower planet mass case, where already destabilized systems were as likely to lose an additional planet as intact ones; although, when considering only 3 M_⊙ primaries with 1 M_J planets, the difference is less noticeable (here, |$31.5^{+3.8}_{-3.4}$| per cent of intact systems lost a planet, compared to |$20.0^{+7.7}_{-4.9}$| per cent of previously destabilized ones). The times of instability are also shifted to much earlier ages (Fig. 10), with the median age of loss at 10^6.5 yr, compared to 10^8.3 yr, the distributions differing significantly (KS test p-value of 3 × 10⁻¹⁴). As in the lower mass case, planet–planet collisions tend to occur earlier and planet–WD collisions later, although the numbers of such collisions are small. Indeed, ejection is the most common outcome by far, with 77 of the 82 planets lost (⁠|$93.9^{+1.7}_{-3.8}$| per cent) being ejected.

Figure 10.

Histogram of cooling ages at which 10 M_J planets are lost around WDs.

Open in new tabDownload slide

4 DISCUSSION

4.1 Generality and reality of our simulations

As the expensive nature of N-body integrations only permits a partial exploration of parameter space, we need to address two questions. First, how general are our results? And secondly, how representative are our systems of those occurring in reality?

Can we generalize the results of our simulations to other values of the planet and stellar mass, and initial semimajor axes, using the semi-analytical scalings presented in Section 2? We restricted our integrations to the study of 1 M_J planets around 3–8 M_⊙ progenitor stars and 10 M_J planets around 3 M_⊙ progenitor stars, with the inner planet initially at 10 au. For the lower mass planets, for all the stellar masses which we explored in our N-body integrations, the estimated stability limits from equations (5) and (6), which were based on the work of Faber & Quillen (2007) who studied planet:star mass ratios in the range 10⁻⁷–10⁻³, gave a good description of where systems can be unstable for MS and WD primaries, as can be seen in Fig. 3. When counting the star's mass-loss, this covers more than a factor of 10 in planet:star mass ratio and 5 in the inner planet's semimajor axis. While these ranges are small compared to the full possible range of planetary masses and semimajor axes, they do inspire confidence that the estimates in Section 2 can be extrapolated to lower planet and stellar masses and to different semimajor axes. We do note that the 10 M_J simulations are rather more stable than the analytical estimates predict. However, below we only consider extrapolation to lower planet masses.

Assuming that equations (5) and (6) accurately reflect actual stability properties, we can predict the range of stable and unstable planetary separations as a function of the inner planet's semimajor axis and the planetary and stellar mass. These ranges are shown in Fig. 11, for planets around 1, 2 and 5 M_⊙ stars. Less massive planets, and those at smaller radii, require a larger separation in Hill radii in order to be stable over any given time period. However, considering the ranges of parameters considered, the critical separations are remarkably insensitive to the planets’ masses and semimajor axes. In contrast, reducing the stellar mass to smaller values has a large effect: the increase in the number of systems that are unstable during the WD phase will be significantly smaller at lower stellar masses, as a consequence of the smaller fraction of mass lost from the star. There will also be an increase in the number of unstable systems on the MS at lower stellar masses, a trend which is seen in our integrations (see Table 1). We note that for solar mass stars significant mass-loss can occur on the RGB as well as the AGB, and this may induce orbital instability during the core Helium burning phase.

Figure 11.

Dependence of the estimated stability limits on inner planets’ semimajor axes (upper panel) and masses (lower panel) for different progenitor stellar masses, shown in different colours. Systems below the solid lines are unstable on the MS; between the solid and the dotted line they are unstable between the end of the MS and the onset of AGB mass-loss; between the dotted and dashed lines they are unstable around WDs; and above the dashed lines they are stable for the age of the Universe. In the upper panel, the planetary masses are each 1M_J. In the lower panel, the innermost semimajor axis is 10 au. Upper and lower solid lines show the stability limits at the end of the MS and at a total age of 13.7 Gyr. Dotted lines show the stability limits just before mass-loss on the AGB. In the upper panel, the WD stability limits are truncated at the moriturus ultimus radius of Mustill & Villaver (2012), within which the inner planet will be engulfed by the AGB star.

Open in new tabDownload slide

We have considered an idealized case where three planets of the same mass form on orbits separated by a constant semimajor axis ratio. In general, systems will not form in such a tidy manner. However, Chambers et al. (1996) showed that including a modest spread of planet masses or semimajor axis ratios does not significantly change the time-scale for instability in a system, although some features such as the effects of MMRs become less pronounced. We do note, however, that the outcomes of instability, including the effects of a second round of instability following mass-loss, will be affected. For example, reducing the mass of one planet in a system will increase the likelihood of its being ejected, as its binding energy is lower.

More generally, we acknowledge that the statistics on instability occurrence and outcomes given in Section 3 may not exactly relate to what will happen in real systems. Even if all systems were indeed of equal-mass, equally spaced planets as we have assumed, the relative frequencies of instabilities around stars at different ages would still depend on the initial separation of the planets. This must be borne in mind when interpreting our statements such as ‘The number of planets lost on the WD exceeds the number lost on the MS in all but the 3 M_⊙ case.’ If planets formed preferentially on very packed orbits, the frequency of MS instability would rise relative to that seen in our simulations.

We must therefore consider more carefully the relationship of our simulations to real planetary systems. Unfortunately, the statistics on systems similar to the ones we study are poor. First, most host stars are smaller than those we consider, with few having masses above 3 M_⊙. Secondly, we are considering relatively large semimajor axes and most detection techniques are biased towards small semimajor axes, with only direct imaging being able to probe the separations we consider . Several planetary systems have now been imaged around stars more massive than the Sun, at wide separations. Vigan et al. (2012), in a survey of 42 MS AF stars, find two with super-Jupiter planets (above 3 M_J), implying a fraction of 5.9–18.8 per cent having such planets between 5 and 320 au after correcting for sensitivity. This limit is consistent with the determination by Nielsen et al. (2013) that ≲20 per cent of 1.5–2.5 M_⊙ stars host 4 M_J planets, and ≲10 per cent host 10 M_J planets, at separations of tens to hundreds of au. One of the Vigan et al. (2012) systems is the famous multiple-planet system HR 8799 (Marois et al. 2008, 2010), which is the closest match to our simulated systems amongst those known. This suggests that multiplanet systems such as we have considered do indeed exist, although they might be uncommon.

Further estimates of the occurrence of the planetary systems herein considered must rely on extrapolations from closer-in planets detectable by radial velocity (RV) measurements or on models of planet formation. From an RV survey of 31 SGs, Bowler et al. (2010) find the occurrence rate of giant planets orbiting within 3 au of evolved intermediate-mass (1.5–2 M_⊙) stars to be |$26^{+9}_{-8}$| per cent. Maldonado, Villaver & Eiroa (2013) find that planet-hosting giant stars of ≳1.5 M_⊙ show metal enrichment similar to MS planet hosts, suggesting similar formation mechanisms around stars of different masses. From a theoretical perspective, Kennedy & Kenyon (2008) argue that giant planet formation through core accretion should be at its most efficient at stellar masses around 3 M_⊙: protoplanetary disc mass may increase with stellar mass, encouraging core formation, but the discs of more massive stars are short lived, decreasing the time available to form cores. However, planet formation through gravitational instability will not be subjected to the same time-scale restrictions and hence may not be disfavoured around higher mass stars.

4.2 WD planets on wide orbits

We now turn to discuss the implications of our study for searches for planets around WDs. Unless planets can form from material ejected from the AGB star, planets around WDs must have survived the full evolution of the star since they formed in the protoplanetary disc. Planets initially within a few au of the star will be engulfed and probably destroyed during the star's AGB evolution (Villaver & Livio 2007, 2009; Mustill & Villaver 2012), although bodies in the BD mass range may survive common envelope evolution to end up on very tight orbits around the WD (Maxted et al. 2006; Parsons et al. 2012; Nordhaus & Spiegel 2013). Mustill & Villaver (2012) showed that giant planets that survive engulfment by the star must be found at distances beyond 2 au, for a 1 M_⊙ progenitor, rising to 10 au for a 5 M_⊙ progenitor. How does the incorporation of multiplanet interactions change the expected distribution of WD planets?

First, a very few systems lost all planets during the course of their evolution. While all single-planet systems where the planet is in an orbit sufficiently wide to avoid engulfment in the envelope will survive to the WD, 15 of our 1736 multiplanet systems simulated ended up devoid of planets. However, neglecting this small rate of total planet loss, it will be the case that every intermediate-mass MS star that hosts at least one wide-orbit (≳10 au) giant planet should give rise to a WD which also hosts at least one wide-orbit giant planet. Thus, the fraction of wide-orbit planets around MS stars and WDs should be identical.

While there have been several surveys for planets around single WDs, there are currently no firm detections, and even the upper limits on the occurrence rate are as yet fairly weak. Direct detection of planets’ thermal emission is aided by the low luminosity of WDs compared to MS stars, and by the shift of the stellar spectrum towards shorter wavelengths (Burleigh, Clarke & Hodgkin 2002). However, the age of systems, typically >1 Gyr, means that planets and BDs have cooled significantly, making only the most massive planets detectable (although planets can experience significant heating during the planetary nebula phase which enhances the luminosity of planets orbiting young WDs; Villaver & Livio 2007). In searches for unresolved companions, Farihi, Becklin & Zuckerman (2008) and Kilic, Gould & Koester (2009) found no evidence for IR excesses consistent with unresolved 10 M_J companions in a sample of 40 WDs, implying a 2σ upper limit of <7 per cent. Unfortunately, the limits for Jovian-mass objects are very weak. Both of these surveys focused on WDs with progenitor masses ≳3 M_⊙. Debes et al. (2011) reported a low fraction (1–5 per cent) of candidate BD companions to WDs, but background confusion prevented secure identification of these sources. While these studies focused on unresolved companions, Hogan, Burleigh & Clarke (2009) sought common proper motion companions to 23 WDs, putting a limit of ≲9 per cent on the fraction of WDs hosting 10 M_J planets at separations of 60–200 au.

Kilic et al. (2009) compared their upper limit on planets orbiting WDs to the detection rate of closer-in RV-detected planets round intermediate-mass (1.3–1.9 M_⊙) SG stars (8.9 per cent, reported in Johnson et al. 2007), finding that the non-detections around WDs are consistent with the occurrence rate around MS stars, but that they are suggestive of a lower frequency. With the more recent results of Bowler et al. (2010) described above, this discrepancy becomes larger. Kilic et al. (2009) suggest three reasons for this difference: first, stars of more than 3 M_⊙ may be inefficient at forming planets (Kennedy & Kenyon 2008); secondly, planets may be lost due to the direct effects of stellar evolution, particularly on the AGB (e.g. Villaver & Livio 2007); thirdly, planets may be lost through instabilities induced through stellar mass loss (Debes & Sigurdsson 2002). They discount the second possibility being important for their sample as planets on wide orbits are safe from the effects of photoevaporation and tides. As a result of our simulations, we have shown that the third reason should not cause a large difference between occurrence rates of planets around MS and WD stars: In only 15 of our integrations, just 1 per cent, did the system lose all of its planets during stellar evolution, despite instability being so common, and in none of these cases was the final planet lost as a result of instability during the WD stage. Hence, we should find wide-orbit planets at similar rates around both MS and WD stars.

We caution that the IR surveys relying on IR excess are biased against wider companions, and comparisons of planet occurrence rates will have to take into account the often large expansion of planetary orbits. In Table 5, we show the number of systems in our simulations whose only surviving planets’ semimajor axes were outside a specified radius. In the most extreme case, 47.8 ± 1.3 per cent of the 1 M_J systems had no planets surviving within 50 au, which may not show up as unresolved IR excesses and would require multi-epoch imaging. The exact numbers given in this table will be rather sensitive to our initial planet masses and semimajor axes, but it will remain the case that in some systems the surviving planet(s) end up on fairly wide orbits around the WD. While planets from single-progenitor systems are not expected to be found within a few tens of au of massive WDs (Villaver & Livio 2007; Mustill & Villaver 2012), scattering in multiple systems can clear planets from much wider regions.

4.3 WD planets on close orbits

There have also been searches for planets close to the WD. Faedi et al. (2011) searched for transits of 194 WDs in the WASP survey, finding an upper limit of 10 per cent for the fraction of WDs hosting very short-period (<0.2 d) giant planets, and progressively weaker limits at larger orbital distances. It is unlikely that planets below ∼10 M_J can survive engulfment in the stellar envelope during the AGB (Villaver & Livio 2007; Nordhaus & Spiegel 2013), although the existence of planets around a horizontal branch star (Charpinet et al. 2011) suggests that in some cases planets can survive, or reform following tidal disruption (Bear & Soker 2012). Hence, close-in planets around WDs must have been brought to the WD vicinity after the WD formed, likely through dynamical processes followed by tidal circularization of their orbit. Our simulations see a small fraction of systems having planets which pass within the WD's Roche limit: 85 planets, out of 755 lost during the WD evolution, in our fiducial 1 M_J case, in 82 separate systems (⁠|$5.5^{+0.6}_{-0.5}$| per cent of all systems). We can make a simple physical argument as to why the delivery of planets close to the WD from these large distances is unlikely. Consider a simple two-planet system where both planets are at approximately the same semimajor axis a⁽ⁱ⁾ and have the same mass. If one planet acquires just enough angular momentum to become unbound, its angular momentum increases by a factor of |$\sqrt{2}$|⁠. The survivor's angular momentum after ejection is therefore |$L^\mathrm{(f)}=(2-\sqrt{2})L^\mathrm{(i)}$|⁠, so it has lost angular momentum during the encounters. However, its apocentre Q^(f) will be at approximately the original semimajor axis a⁽ⁱ⁾, and so the final eccentricity would be around 0.66, meaning a pericentre of 8 au for a planet with an apocentre at 40 au. In numerical simulations, Ford & Rasio (2008) found a maximum eccentricity of ≈0.8 for the survivor of a two-planet system. More complicated exchanges of angular momentum in a three-planet system can of course take place, but it remains challenging to drive a planet's eccentricity to the very high values (≳0.99) required to be able to tidally circularize it.

The actual fate of these planets may vary. In the next subsection, we discuss the effects of their tidal disruption. Here, however, we ask what happens if, before the planets are driven inside the Roche limit, they first pass the star at a slightly greater distance and their orbits become tidally circularized. This would result in their being dynamically decoupled from the source of perturbations that would have eventually driven them into the star, while eventually shrinking the orbit into a circular one at just beyond 2R_Roche. This mechanism is proposed to explain the distribution of hot Jupiters around MS stars (Fabrycky & Tremaine 2007; Nagasawa, Ida & Bessho 2008), and the dynamics around WDs will be the same, since the dominant tidal forces occur in the planet and do not depend on the evolutionary state or structure of the star.

As shown in Fig. 8, surviving planets at the end of our integrations in general do not have pericentres close to the Roche limit. Because the pericentres evolve with time, some may come closer-than recorded at the end of the simulation. Over the whole WD stage, we found 76 planets with a pericentre recorded within 0.1 au of the WD, of which 48 later collided with the star. Although our coarse time sampling (every 1 Myr) means that we may have missed some planets that have close approaches and may be tidally circularized, these figures suggest that the bulk of potentially close-in planets come from those that would suffer a stellar collision in the absence of tides. Hence, we assume the maximum fraction of systems forming ‘hot Jupiters’ around WDs would be given by the fraction that eventually enter the Roche limit, ∼5.5 per cent. Since ≲20 per cent of MS systems host the systems of massive planets that we consider in this paper, this would set an upper limit of ≲1 per cent of WDs hosting close-in Jupiters, which then would require a favourable geometric alignment to be observed as transiting systems. This is consistent with the upper limit of 10 per cent of WDs hosting close-in Jupiters from Faedi et al. (2011), and we can predict that we can only expect detections of such planets with greatly increased sample sizes of several thousands of WDs.

The number of planets scattered-in does of course depend on the system architecture. In our 10 M_J simulations, only |$0.8^{+1.0}_{-0.3}$| per cent of systems saw a planet driven to the Roche limit, meaning the number of transiting super-Jupiters will be even smaller. However, reducing the planet mass considerably will not result in many bodies being scattered in: in our two-planet simulations in Paper I, for the case of two Earth-mass planets with the inner at 10 au, we saw no collisions with the star over the whole stellar lifetime, including during the giant stages (<0.6 per cent). Although there was no strong trend in the number of collisions with the star as a function of stellar mass, it is the case that smaller stellar masses allow survival closer to the star during the AGB, and scattering from tighter final orbital radii may increase the number of close-in planets slightly. On the other hand, as discussed above, the extra parameter space destabilized is likely to be smaller. Hence, lower mass WDs with progenitors in the range 1–3 M_⊙ may offer a higher probability of hosting transiting planets, but this will still likely be small.

4.4 WD pollution

If a planet is not tidally circularized before reaching the Roche limit, it will suffer partial or complete disruption (Zuckerman et al. 2007; Klein et al. 2010; Debes et al. 2012; Gänsicke et al. 2012). The dynamics of planets passing within or close to the Roche limit are complex and depend on the planet's internal structure and the distance of approach to the star. If total disruption does not occur, mass is stripped from the outer layers of the planet, the remainder changing its orbit to become either more or less bound to the star, or even in some cases unbound (Faber, Rasio & Willems 2005; Guillochon, Ramirez-Ruiz & Lin 2011; Liu et al. 2013). A delivery of all of a planet's mass to pollute a WD would represent only an extreme case.

Direct pollution of a WD in this way is however rare: up to 5.5 per cent of massive WDs with multiplanet systems may be polluted directly by a gas giant. The accretion of a 1 M_J planet on to a WD would lead to a peak luminosity of order 10⁴ L_⊙, using the scalings in Bear & Soker (2013). However, such events would occur infrequently. Assuming a Galactic population of ∼10¹⁰ WDs (Napiwotzki 2009), 10 per cent of which have multigiant planet systems, and considering that 5 per cent of these will send a planet into the WD over a 5 Gyr time period, suggests that the rate of WDs accreting Jovian planets should be about one per century in the Galaxy.

Observed signatures of metal pollution in WDs, which occurs in around 25 per cent of DA WDs (Zuckerman et al. 2003), are generally consistent with the accretion of rocky asteroids. Debes & Sigurdsson (2002) proposed that the destabilization of giant planet systems could trigger the destabilization of planetesimals in the system and their subsequent accretion on to the WD. Our integrations confirm that instability is a common outcome for closely packed multiplanet systems around WDs. The time at which the instabilities occur in these integrations is at cooling ages of several 100 Myr, which was also the case in the two-planet simulations we studied in Paper I. As discussed in that paper, this is similar to the ages at which metal pollution is observed. However, the low fraction of MS stars and WDs that host these multigiant planet systems shows that such instabilities cannot occur around a large fraction of WDs. Moreover, around a quarter of planetary systems that experience instability around WDs have previously experienced instability, which would have already depleted their planetesimal reserves, although such systems may still see low levels of pollution from residual planetesimals. If we take an observational upper limit of 20 per cent of stars hosting systems of multiple giant planets on wide orbits, as discussed in Section 4.1, and consider that in our integrations 43 per cent of WD systems are unstable, we find that giant planet instability alone would predict pollution of ≲9 per cent of WDs. This falls somewhat short of the observed pollution occurrence fraction; furthermore, studies of planetesimal scattering in the wake of planetary instability around MS stars show that a high planetesimal flux to the regions near the star lasts only a short time (Booth et al. 2009; Bonsor, Raymond & Augereau 2013). Pollution can be achieved with a single planet, although the required planetesimal reservoirs are often orders of magnitude larger than the Solar system's asteroid belt (Debes et al. 2012). We note that a collision between terrestrial planets may create a fresh reservoir of debris during the WD phase that could then be perturbed by other planets (Gänsicke, private communication). While planet–planet collisions were not a common event during the WD phase in the simulations described in this paper, we did not study terrestrial-mass planets. In Paper I, we found that planet–planet collisions were more common in systems of two Earth-mass planets than in systems of two giant planets. Further studies of systems including low-mass planets would be needed to quantify this process and its sensitivity to the system set-up. It is likely that a variety of mechanisms plays a role, commensurate with the diversity of planetary system architectures.

4.5 Free-floating planets

A purportedly vast population of free-floating giant planets is travelling between and outnumbers the stars of the Milky Way, as revealed by microlensing observations (Sumi et al. 2011). If sufficiently young and nearby, free-floating planets may be imaged. Although the mass of these free-floaters is uncertain due to a degeneracy in their age, luminosity and mass, Delorme et al. (2012) have recently provided a near-confirmation of the mass for one of these objects, helping us to confirm the existence of free-floating giant planets.

The source of such a large free-floating planet population is unlikely to be dynamical instabilities arising from gravitational scattering alone because otherwise the average number of giant planets formed per planetary system would have to be unrealistically high (Veras & Raymond 2012). However, that study considers scattering only on the MS, and claims that scattering during the WD phase contributes to the total free-floating giant planet population of the order of 1 per cent based solely on the currently observed space density of WDs.

Our work can finally help assess the relative contributions to the free-floating population from MS scattering versus WD scattering. Although the fraction is dependent on our initial conditions, our initial condition choices encompass the entire phase space of separations from near-guaranteed instability to near-guaranteed stability. Furthermore, our assumption of initially circular orbits essentially controls the time-scale for instability to occur, and hence might have a weak dependence on the resulting dynamics for ensembles of systems (e.g. Jurić & Tremaine 2008).

We find that, summed over all simulations reported in Table 1, an ejection is 3.9 times more likely to occur while the parent star is a WD than when the star is an MS star. This difference is most pronounced for the highest mass stars (with a factor of 8–10 for 6-8 M_⊙ stars) and least pronounced for the least massive stars (with a factor of 1-2 for 3–4 M_⊙ stars). Because of the rarity of Milky Way stars with masses greater than or equal to 3 M_⊙, we do not expect the increased rate of ejections during the WD phase to provide a major source for the current free-floating planet population. However, if future work can show that the potentially complex orbital dynamics accompanying M ≲ 1 M_⊙ stellar evolution yields a planetary ejection rate around WDs that is several factors higher than the MS rate, then the Galaxy's free-floating planet population may be a strong function of time.

5 CONCLUSIONS

We have carried out numerical integrations of systems of three giant planets orbiting 3–8 M_⊙ stars up to a total system age of 5 Gyr, covering all the star's evolutionary stages from the beginning of the MS and for several Gyr of WD cooling. We have studied systems of planets on wide orbits, initially beyond 10 au, which are unaffected by stellar evolutionary effects such as tides, feeling only the changing stellar mass. By following systems through the entirety of the star's lifetime, we ensure that the systems we simulate on the post-MS have survived MS evolution.

In common with previous scattering studies, we find that on the MS stability is dependent on the planets’ initial spacing, with planets separated by less than 9 (single-planet) Hill radii prone to disruption. The outcome of these instabilities is the ejection of a planet or the collision of one planet with another or with the star. In some cases, two planets are lost from the systems, leaving only one survivor. Loss of all planets on the MS does not occur in our simulations.

When the host star evolves off the MS, the increase in stellar radius can engulf planets in systems that experienced instability on the MS, where these planets are left on eccentric orbits with pericentres close to the star. Instability can also leave a surviving planet with a large apocentre, in which case the planet can become unbound as a result of the star's mass-loss. In a very few cases, one of these effects succeeds in removing the system's third and final planet, leaving the star totally denuded of its planetary entourage. This happens but rarely however, in about 1 per cent of all our integrations. This implies that the fraction of MS stars hosting planets on orbits of ∼10 to a few hundreds of au should be almost identical to the fraction of their WD descendants hosting planets at similar radii.

Few hitherto-stable systems manifest instability between the end of the MS and the start of the WD stage, due to the relatively short duration of this evolutionary phase. Nor does the mass-loss at the end of the AGB immediately trigger destabilization of previously stable systems. Instead, instability around WDs becomes manifest typically at cooling ages of several hundred Myr. Instability in WD systems strikes both systems that survived previous evolution unscathed, and systems that previously lost a planet and have two remaining on eccentric orbits. In the latter case, instability occurs at slightly later cooling ages, with a median of 10^8.4 yr as opposed to 10^7.9 yr in our 1 M_J integrations.

The overwhelming outcome of WD instability is the ejection of at least one planet. Collisions between planets or between a planet and a star are relatively rare. We find that very few planets are scattered on to orbits with small pericentres, where they may then become tidally circularized and discovered as WD-transiting planets. Combining our simulation results with planet detection rates around MS stars, we expect ≲1 per cent of WDs to host close-in planets.

Similarly, we expect that the pollution of WDs following instability in systems of giant planets, whether directly through the accretion of a planet itself or indirectly through the destabilization of small bodies following planetary instability, would occur in ≲7 per cent of WDs, insufficient to explain observed pollution rates.

Although we restricted our integrations to a particular range of parameters, we have made predictions for behaviour for other parameter values. In particular, we expect the number of systems destabilized around low-mass WDs with progenitors <3 M_⊙ to be much lower, as a consequence of the longer MS lifetimes and particularly the smaller amount of mass-loss.

AJM and EV are supported by Spanish grant AYA 2010/20630. EV also acknowledges the support of the Marie Curie grant FP7-People-RG268111. Additionally, the research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement no. 320964 (WDTracer). We thank the anonymous referee for inspiring a number of points of clarification in the paper. We thank Boris Gänsicke and Amy Bonsor for comments on the manuscript.

REFERENCES