Abstract
Although 25–50 per cent of white dwarfs (WDs) display evidence for remnant planetary systems, their orbital architectures and overall sizes remain unknown. Vibrant close-in (≃1 R⊙) circumstellar activity is detected at WDs spanning many Gyr in age, suggestive of planets further away. Here we demonstrate how systems with 4 and 10 closely packed planets that remain stable and ordered on the main sequence can become unpacked when the star evolves into a WD and experience pervasive inward planetary incursions throughout WD cooling. Our full-lifetime simulations run for the age of the Universe and adopt main-sequence stellar masses of 1.5, 2.0 and 2.5 M⊙, which correspond to the mass range occupied by the progenitors of typical present-day WDs. These results provide (i) a natural way to generate an ever-changing dynamical architecture in post-main-sequence planetary systems, (ii) an avenue for planets to achieve temporary close-in orbits that are potentially detectable by transit photometry and (iii) a dynamical explanation for how residual asteroids might pollute particularly old WDs.
1 INTRODUCTION
Our Solar system is packed. In other words, it could quickly become dynamically unstable with the insertion of an additional planet.1 Mounting discoveries of packed extrasolar systems (e.g. Lissauer et al. 2011, 2013; Swift et al. 2013; Cabrera et al. 2014; Masuda 2014) suggest that the Solar system is not unusual, particularly as detection sensitivities improve and formerly inaccessible discovery space is becoming available for exploration.
The concept of packing has received increasing attention over the last decade. Planet formation often appears to be efficient, even in the extreme environment of the first confirmed exoplanetary system (PSR B1257+12), discovered over two decades ago (Wolszczan & Frail 1992; Wolszczan 1994) around a supernovae remnant. The ‘packed planetary systems’ hypothesis (Barnes & Raymond 2004; Barnes, Goździewski & Raymond 2008) advocates that efficient formation processes create planets with mutual separations that leave no room for additional companions. Dynamical scattering subsequent to formation represents another channel which can achieve tightly packed systems (e.g. Raymond et al. 2009a). The tally of packed systems continuously increases with observations that rely on transit photometry (Lissauer et al. 2011; Fang & Margot 2013), Doppler radial velocity measurements (e.g. Lovis et al. 2011; Anglada-Escudé et al. 2013) and even direct imaging (Marois et al. 2010). Now planetary packing has become the framework in which to analyse a variety of complex systems, such as those containing moons (Kane, Hinkel & Raymond 2013; Payne et al. 2013) or additional stars (Kratter & Shannon 2014). Theoretical interest in the formation and evolution of these packed systems has also been robust (Hansen & Murray 2012, 2013; Chatterjee & Tan 2014; Hansen & Murray 2014; Hands, Alexander & Dehnen 2014).
The frequent occurrence of packed systems in nature motivates analyses of their ultimate fates. Although we do not yet have a confirmed detection of a planet orbiting a WD within a few thousand au2 (Mullally et al. 2008; Hogan, Burleigh & Clarke 2009; Debes et al. 2011; Faedi et al. 2011; Steele et al. 2011; Fulton et al. 2014), we know that between 25–50 per cent of all white dwarfs (WDs), the end product of stellar evolution for over 90 per cent of all Milky Way stars, contain remnant planetary systems. This remarkable statistic follows from observations of abundant photospheric metal pollution in WD atmospheres (Zuckerman et al. 2003, 2010; Koester, Gänsicke & Farihi 2014). These metals must be accreted from circumstellar matter, which is detected in a number of cases in the form of dusty (Zuckerman & Becklin 1987; Becklin et al. 2005; Kilic et al. 2005; Reach et al. 2005; Farihi, Jura & Zuckerman 2009) and gaseous discs (Gänsicke et al. 2006; Gänsicke, Marsh & Southworth 2007; Gänsicke 2011; Farihi et al. 2012; Melis et al. 2012). High levels of variability in these discs (e.g. Gänsicke et al. 2008; Wilson et al. 2014; Xu & Jura 2014) on time-scales of years highlight the current dynamical activity occurring in these systems.
The debris discs, with a typical radial extent of about 1 solar radius, cannot have existed during the main-sequence (MS) or giant branch phases of stellar evolution. These discs are currently thought to originate from asteroids that are perturbed by planets on to highly eccentric orbits (Bonsor, Mustill & Wyatt 2011; Debes, Walsh & Stark 2012; Frewen & Hansen 2014), whereupon the asteroids disrupt close to the WD (Graham et al. 1990; Jura 2003; Debes et al. 2012; Bear & Soker 2013; Veras et al. 2014b). The disc then accretes on to the WD (Bochkarev & Rafikov 2011; Rafikov 2011a, b; Metzger, Rafikov & Bochkarev 2012; Rafikov & Garmilla 2012; Wyatt et al. 2014), creating signatures of metal pollution. Pollution in WD systems without observed orbiting discs still somehow must arise from rocky planetesimals (Bergfors et al. 2014).
The presence of planets orbiting WDs on tight orbits would help explain asteroid delivery, but also raise the possibility of harbouring life by residing in a WD habitable zone (Monteiro 2010; Agol 2011; Barnes & Heller 2013). Fossati et al. (2012) has demonstrated that photosynthetic processes associated with complex life can be sustained on these planets. Further, because habitable WD planets would reside within 0.1 au, these planets are in principle more easily detectable by transit (Agol 2011) and by polarized or reflected light (Fossati et al. 2012). Biomarkers in the atmospheres of habitable planets transiting WDs might be observable with the James Web Space Telescope (JWST; Loeb & Maoz 2013). Indeed, JWST is capable of detecting objects with masses three orders of magnitude less massive than the Moon (Lin & Loeb 2014). However, so far, transit searches have been unsuccessful (Faedi et al. 2011; Fulton et al. 2014), and one theoretical prediction based on a simulation suite of three Jovian-mass planets suggests that less than 1 per cent of WDs host close-in Jupiters (see section 4.3 of Mustill, Veras & Villaver 2014).
Motivated by both the robust signatures of remnant planetary systems and the possibility of detecting close-in planets, we here explore the link with packed planetary systems on the MS. First, however, we will review how long-term instability manifests itself in multiplanet systems, and comment on tidal effects between close-in planets and WDs.
1.1 Multiplanet instabilities
Before considering how multiple planets become unstable in post-MS systems, we should understand instability on the MS (Davies et al. 2013). During a star's MS lifetime, the star's mass and radius change negligibly. These changes translate into little dynamical excitation amongst orbiting bodies; Veras & Wyatt (2012) point out that the Solar system planets would increase their semimajor axes by at most about 0.055 per cent. Due to such small changes, the vast majority of all MS exoplanet studies treat the stellar mass and radius as fixed. However, after the star leaves the MS, then the mass-loss becomes extreme – causing orbital expansions by tens to several hundred per cent, or outright ejection (Veras et al. 2011; Adams, Anderson & Bloch 2013) – and the stellar radius increases by several orders of magnitude. Consequently, destabilization of planetary systems becomes more likely. In the Solar system, Mercury and Venus will be swallowed by the Sun, and Earth might too be destroyed (Schröder & Connon Smith 2008).
Here, we do not account for the gravitational effect of any residual smaller masses, such as planetesimals, asteroids or comets, on the planets: we treat our initial systems as already dynamically settled from formation. In reality, complex interactions between planetesimals and planets likely play important dynamical roles, primarily at the earliest MS ages; their incorporation enormously increases the available parameter space and computational cost of simulations (Tsiganis et al. 2005; Raymond, Armitage & Gorelick 2009b, 2010; Levison et al. 2011).
Stability boundaries are well defined for two-planet systems. If the distance between two planets exceeds the Hill stability boundary, then those planetary orbits will remain ordered (non-crossing) forever (for a review see Georgakarakos 2008), but not necessarily stable: the inner planet could collide with the star or the outer planet could escape. If that distance instead exceeds the Lagrange stability boundary, then the planets will forever remain bounded and ordered. The Hill stability boundary can be expressed in Jacobi coordinates for arbitrary eccentricities and inclinations (Donnison 2011) and is entirely analytical, except for a usually negligible truncation in the expression of the energy (fig. 19 of Veras et al. 2013a).3 Unfortunately, no such analytical formulation exists for Lagrange stability. However, empirical estimates (Barnes & Greenberg 2006; Veras & Mustill 2013) help establish this boundary. Also, two planets which are not Hill stable may still be bounded and ordered over long time-scales. This behaviour is especially connected to the overlap of mean motion resonances (Chirikov 1979; Wisdom 1980; Mardling 2008; Funk et al. 2010; Mustill & Wyatt 2012; Deck, Payne & Holman 2013; Giuppone, Morais & Correia 2013; Bodman & Quillen 2014).
Two examples of post-MS planetary unpacking with terrestrial planets. The top and bottom panels each correspond to one simulation from rows 1 and 6 of Table 1, respectively. The background colours refer to the main sequence (MS; yellow), giant branch (GB; pink) and white dwarf (WD; grey) phases of stellar evolution. In the top panels, the four planets remain packed for the entire MS and GB phases, and for about 3 Gyr of WD evolution. In the bottom panels, the 10 planets become unpacked during the GB phases (specifically, during the asymptotic giant branch). In both cases, planets make sweeping incursions well within the initial innermost semimajor axis and maximum stellar radius Rmax, and at WD cooling ages exceeding 5 Gyr. In the top-left panel, the initially outermost planet (red) spent Gyr of WD cooling time at close pericentres before being engulfed into the WD at 12.2 Gyr. The bottom-left panel features a particularly intrusive incursion: the red planet achieves a pericentre of 0.046 au – close to a typical WD habitable zone and beyond the WD Roche radius Rroche – at 10.68 Gyr, and survives. Depending on the internal properties of the planet, it might be subsequently tidally circularized – a process ignored in this work. The semimajor axis evolution showcases the complex energy exchanges triggered in these systems after the first instance of instability, and demonstrates that close pericentre passages are due to variations in orbital eccentricity.
For systems with more than two planets, stability estimates become decidedly more imprecise. No analytical criterion exists like for the two-planet Hill stability case. Consequently, several authors (Chambers, Wetherill & Boss 1996; Faber & Quillen 2007; Zhou, Lin & Sun 2007; Chatterjee et al. 2008; Smith & Lissauer 2009; Quillen 2011) have attempted to fit empirical instability time-scales to systems of three or more planets from numerical simulations. Because these estimates are primarily based on simulations with equal-mass planets, here we also adopt equal-mass planets in order to utilize these estimates. The rich dynamics in systems with a (more realistic) set of unequal-mass planets attests to its greater number of degrees of freedom (Veras & Armitage 2005, 2006; Ford & Rasio 2008; Jurić & Tremaine 2008; Raymond et al. 2011, 2012; Matsumura, Ida & Nagasawa 2013).
The dynamical complexities increase in post-MS systems, for which relatively few multiplanet studies have been performed. All such studies consider gravitational scattering which results from instability due to mass-loss. Duncan & Lissauer (1998) investigated the future evolution of collections of planets in artificial Solar systems for various planet–Sun mass scalings. Debes & Sigurdsson (2002) evolved systems with two and three equal-mass planets over 1000 orbits around a progenitor 1 M⊙ star that lost half of its mass in a uniform fashion. Veras et al. (2013a) evolved two-planet systems with equal planetary masses from the endpoint of formation to several Gyr into the WD phase for progenitor stellar masses of 3–8 M⊙; this work was recently extended with simulations of equal-mass planets in three-planet systems (Mustill et al. 2014). Voyatzis et al. (2013) instead considered a more analytical approach, tracking the regular and chaotic trajectories in phase space that results from the inclusion of mass-loss in multiplanet systems.
1.2 Tidal effects
Planets scattered close to their parent stars might change their orbits due to star–planet tides. MS and giant branch studies dominate investigations of star–planet tidal dynamics; few references to WD–planet interactions exist. Nordhaus & Spiegel (2013) pose in their section 4 that an Earth-mass planet driven to within half of a solar radius (about 0.002 au) of a WD might be circularized due to tidal friction, but do not elaborate further on this process. They do surmise that the resulting dissipation of heat would render the planet uninhabitable, even if detectable. Mustill et al. (2014) suggest that the dynamics of WD–planet tidal interactions would be nearly equivalent to that of MS stars, and hence feature circularization, because the dominant forces arise from tidal deformation of the planet and not the star.
If we make that assumption here, then all of the uncertainties which plague tidal theory on the MS carry over to WD systems. This problem is compounded by the newfound inoperability or extremely limited scope (Efroimsky & Makarov 2013) of classic tidal theories such as the constant geometric lag model (Goldreich 1966; Murray & Dermott 1999) and the constant time lag model (Mignard 1979; Hut 1981). Consequently, a process such as tidal circularization of a terrestrial planet orbiting a WD must be taken on a case-by-case basis, one that is highly dependent on the rheology of the planet (but independent of the WD). The circularization time-scale can easily vary by orders of magnitude depending on this rheology (Henning & Hurford 2014). Giant planets also suffer from this uncertainty, as they may contain solid cores, necessitating the combined modelling of fluid mechanics and solid mechanics. As observed by Ogilvie (2014), Earth's tidal dissipation is dominated by a fluid layer that represents just 0.023 per cent of the entire planet's mass.
Because our planets here represent point masses, we make no attempt to invent interior structures, neither compute maximum tidal reaches nor circularization time-scales. Nevertheless, we emphasize throughout the manuscript that we neglect these processes, which may strongly affect planets that achieve pericentres at orbital separations from the WD that are realistically within the reach of observational detection.
We now proceed to setup, execute and analyse our simulations in Section 2, discuss the results in Section 3 and conclude in Section 4.
2 SIMULATIONS
Our simulations incorporated self-consistently the evolution of both the host star and planetary orbits to model the progression of the planetary system throughout the MS, giant branch and WD phases. We performed integrations over the age of the Universe, which is approximately 14 Gyr.
2.1 Numerical code
We used an integrator with features from the simulators used in Veras et al. (2013a), Mustill et al. (2014) and Veras, Shannon & Gänsicke (2014d). The integrator contains a Bulirsch–Stoer algorithm, originally from Chambers (1999), that is modified with stellar mass and radius interpolation between substeps. The stellar mass and radius profiles are computed with the sse code (Hurley, Pols & Tout 2000), assuming a Reimers mass-loss prescription with coefficient 0.5 during early giant branch phases and a superwind prescription from Vassiliadis & Wood (1993) along the asymptotic giant branch phase. We assumed that the mass-loss is isotropic, an excellent assumption for the planet–star separations considered in this paper (Veras, Hadjidemetriou & Tout 2013b). When the star becomes a WD, the resulting (constant) radius is replaced with the WD's Roche radius assuming the planets all have a typical terrestrial planet density of 5.5 g cm−3 (see Section 2.2.2). When a planet collides with the star, the star is modelled to absorb the planet and increase mass accordingly, affecting the subsequent dynamics of the system. Given the long duration of our simulations, we outputted data every 1 Myr, but importantly continuously tracked the minimum pericentre that every surviving planet ever achieved. We, conservatively, adopted a Bulirsch–Stoer tolerance parameter of 10−13. Consequently, in every case, we have found that angular momentum was conserved to within 10−6; energy is not conserved in mass-losing systems.
2.2 Initial conditions
2.2.1 Stellar mass
Our choice of initial stellar masses in this work (1.5, 2.0 and 2.5 M⊙) is an improvement over previous studies (Veras et al. 2013a; Mustill et al. 2014). The average mass of a WD in the present-day WD population is 0.60–0.65 M⊙ (Liebert, Bergeron & Holberg 2005; Falcon et al. 2010; Tremblay et al. 2013), which corresponds to an A-type MS progenitor star with a mass around 2 M⊙ (Catalán et al. 2008; Kalirai et al. 2008; Casewell et al. 2009). Fig. 1 of Koester et al. (2014) illustrates that by considering the stellar mass range of 1.5–2.5 M⊙, we sample the progenitor mass range of nearly all known polluted WDs. This mass range corresponds to MS lifetimes between about 600 Myr and 3 Gyr. Consequently, the majority of the integration time was during the WD phase of stellar evolution. The final WD masses corresponding to these initial MS masses are 0.576, 0.637 and 0.690 M⊙; the vast majority of the mass-loss occurred on the asymptotic giant branch phase, as during the giant branch phase only 0.061, 0.002 and 0.001 M⊙ were shed.
2.2.2 Planetary mass and orbits
We adopted equal, Earth-mass planets (3.0035 × 10−6 M⊙) in our main suite of simulations; later, we perform a few additional full-lifetime simulations with Jupiter-mass planets (9.5460 × 10−4 M⊙) in order to showcase the differences. We focus on terrestrial planets because (i) at the separations probed so far, low-mass planets appear to be much more common than giant planets (Winn & Fabrycky 2015), (ii) observational evidence for their packing is greater and (iii) there already exist full-lifetime estimates for the critical survival separations throughout MS evolution for four or more terrestrial planets (Smith & Lissauer 2009).
The process of unpacking is largely scale invariant, i.e. independent of our choice of a1. However, there do exist physical constraints. Johansen et al. (2012) and Petrovich, Tremaine & Rafikov (2014) show an increased incidence of collisions with the central star as a1 is decreased, a phenomenon which can be explained through the Safronov number (Safronov & Zvjagina 1969). This number represents the square of the ratio of orbital speed to surface escape speed. As this ratio increases, the frequency of collisions increases.
Another constraint includes star–planet tides along the post-MS phases, as a star's radius will expand significantly and might swallow the planet. For the stellar mass range, we consider (1.5–2.5 M⊙), a terrestrial planet with a circular orbit and MS semimajor axis of at least about 3 au will survive engulfment from the expanding giant star (fig. 7 of Mustill & Villaver 2012). This bound remains robust even for different mass-loss prescriptions (Adams & Bloch 2013; Villaver et al. 2014). Another constraint is the age of the Universe, which limits the total number of orbits which could be completed as a function of a1, affecting Lagrange-type instability (e.g. Veras & Mustill 2013). However, as the number of orbits increase, computations take longer and may accumulate error in angular momentum conservation.
For these reasons, we choose a1 = 5 au for most of our simulations; β = 10 then gives (a2, a3, a4 = 5.67, 6.44, 7.30) au.
We do not include Galactic tides in our simulations, as planets in the solar neighbourhood are negligibly perturbed by these tides for separations up to thousands of au (Kaib, Raymond & Duncan 2013; Veras & Evans 2013a, b), and adiabatic mass-loss will increase the outermost planet's initial semimajor axis [a10 < 20 au] by a factor of just a few. The outcomes of our simulations justify this choice ex post facto. We also do not include flybys from passing stars, for both computational reasons and so that we can focus on the unpacking process. In reality, we should expect, over the age of the Universe, one Galactic field star to come within a few hundred au of another (Zakamska & Tremaine 2004; Veras & Moeckel 2012). However, post-MS mass-loss occurs for just a small fraction of this time. Therefore, during post-MS mass-loss, a stellar flyby would not be expected to penetrate to within 104 au in the solar neighbourhood (Veras et al. 2014a).
Our simulations also do not include the effect of general relativity (GR). During the course of each orbit, GR can displace a highly eccentric planet from its predicted pericentre location around a typical 0.6 M⊙ WD by up to 2.65 km (Veras 2014). This value is independent of the planet's semimajor axis or pericentre. Other changes to the orbit during one complete period are similarly minor. However, over many orbits, the pericentre may precess noticeably; the extent of this precession is in fact dependent on the semimajor axis and eccentricity. Equation (19) of Veras et al. (2014b) demonstrate that any highly eccentric object approaching the Roche radius of the WD can experience precession on the order of a few degrees within 105 yr. Consequently, any secular behaviour which is strongly dependent on the argument of pericentre, such as the Lidov–Kozai mechanism, could be affected by GR.
Unlike Smith & Lissauer (2009), we initialized our planets on slightly non-coplanar orbits. We drew inclination values from random uniform distributions with ranges of [−0.5°, 0.5°] in order to avoid an artificially high number of collisions, given that real systems are not exactly coplanar but tend to harbour planetary orbits which are mutually inclined on the order of degrees. 4 We performed a similar random selection for all orbital angles (mean anomaly, argument of pericentre, longitude of ascending node) over their entire ranges. We set the planets on initially circular orbits for simplicity and because Jurić & Tremaine (2008), among others, suggest that qualitative scattering outcomes are independent of initial eccentricity.
2.3 Simulation results
2.3.1 Terrestrial planet suite
Here we discuss the results of our main set of simulations, summarized in Table 1, and show two illustrative evolution examples in Fig. 1. Fig. 2 summarizes the results from the majority of the simulations in the table. We performed four simulations with slightly different initial conditions (Section 2.2.2) for each set of parameters, i.e. for each row in Table 1. Within each row, the number of planets as well as the values of M⊙, β and a1 were fixed.
Minimum pericentres (in au) which are less than 5 au for systems in Table 1. Each numberline corresponds to a different (labelled) row of the table (rows 10 and 14 did not unpack), and each numberline displays all planets in all four runs of a given row; each run is colour-coded. Note that even for rows 7 and 8, where a1 = 10 au, at least one planet achieves a minimum pericentre within 5 au. This figure demonstrates that an endemic feature of post-MS unpacking is inward incursions of one or more planets to values which are just a fraction of their original MS semimajor axes.
Summary of simulations, four per row. Each simulation contains either Earth-mass planets or Jupiter-mass planets, separated by a mutual Hill separation indicated by β through equation (1). If β is too large, instability will never occur. If β is too small, the first instability will occur quickly (i.e. on the MS).
| Terrestrial planets | |||||
|---|---|---|---|---|---|
| Row | No. of | Progenitor | β | a1 | Per cent unpacked |
| No. | planets | mass ( M⊙) | (au) | on post-MS | |
| 1 | 4 | 1.5 | 10 | 5 | 75 |
| 2 | 10 | 1.5 | 10 | 5 | 100 |
| 3 | 4 | 2.0 | 10 | 5 | 100 |
| 4 | 10 | 2.0 | 10 | 5 | 100 |
| 5 | 4 | 2.5 | 10 | 5 | 100 |
| 6 | 10 | 2.5 | 10 | 5 | 100 |
| 7 | 4 | 1.5 | 10 | 10 | 25 |
| 8 | 10 | 1.5 | 10 | 10 | 100 |
| 9 | 4 | 1.5 | 12 | 5 | 25 |
| 10 | 4 | 1.5 | 15 | 5 | 0 |
| 11 | 4 | 2.0 | 12 | 5 | 100 |
| 12 | 4 | 2.0 | 15 | 5 | 25 |
| 13 | 4 | 2.5 | 12 | 5 | 100 |
| 14 | 4 | 2.5 | 15 | 5 | 0 |
| Giant planets | |||||
| 15 | 4 | 1.5 | 8 | 5 | 0 |
| 16 | 4 | 2.0 | 8 | 5 | 0 |
| 17 | 4 | 2.5 | 8 | 5 | 0 |
| 18 | 4 | 1.5 | 6 | 5 | 100 |
| 19 | 4 | 2.0 | 6 | 5 | 100 |
| 20 | 4 | 2.5 | 6 | 5 | 100 |
| 21 | 4 | 1.5 | 4 | 5 | 0 |
| 22 | 4 | 2.0 | 4 | 5 | 0 |
| 23 | 4 | 2.5 | 4 | 5 | 0 |
| Terrestrial planets | |||||
|---|---|---|---|---|---|
| Row | No. of | Progenitor | β | a1 | Per cent unpacked |
| No. | planets | mass ( M⊙) | (au) | on post-MS | |
| 1 | 4 | 1.5 | 10 | 5 | 75 |
| 2 | 10 | 1.5 | 10 | 5 | 100 |
| 3 | 4 | 2.0 | 10 | 5 | 100 |
| 4 | 10 | 2.0 | 10 | 5 | 100 |
| 5 | 4 | 2.5 | 10 | 5 | 100 |
| 6 | 10 | 2.5 | 10 | 5 | 100 |
| 7 | 4 | 1.5 | 10 | 10 | 25 |
| 8 | 10 | 1.5 | 10 | 10 | 100 |
| 9 | 4 | 1.5 | 12 | 5 | 25 |
| 10 | 4 | 1.5 | 15 | 5 | 0 |
| 11 | 4 | 2.0 | 12 | 5 | 100 |
| 12 | 4 | 2.0 | 15 | 5 | 25 |
| 13 | 4 | 2.5 | 12 | 5 | 100 |
| 14 | 4 | 2.5 | 15 | 5 | 0 |
| Giant planets | |||||
| 15 | 4 | 1.5 | 8 | 5 | 0 |
| 16 | 4 | 2.0 | 8 | 5 | 0 |
| 17 | 4 | 2.5 | 8 | 5 | 0 |
| 18 | 4 | 1.5 | 6 | 5 | 100 |
| 19 | 4 | 2.0 | 6 | 5 | 100 |
| 20 | 4 | 2.5 | 6 | 5 | 100 |
| 21 | 4 | 1.5 | 4 | 5 | 0 |
| 22 | 4 | 2.0 | 4 | 5 | 0 |
| 23 | 4 | 2.5 | 4 | 5 | 0 |
Summary of simulations, four per row. Each simulation contains either Earth-mass planets or Jupiter-mass planets, separated by a mutual Hill separation indicated by β through equation (1). If β is too large, instability will never occur. If β is too small, the first instability will occur quickly (i.e. on the MS).
| Terrestrial planets | |||||
|---|---|---|---|---|---|
| Row | No. of | Progenitor | β | a1 | Per cent unpacked |
| No. | planets | mass ( M⊙) | (au) | on post-MS | |
| 1 | 4 | 1.5 | 10 | 5 | 75 |
| 2 | 10 | 1.5 | 10 | 5 | 100 |
| 3 | 4 | 2.0 | 10 | 5 | 100 |
| 4 | 10 | 2.0 | 10 | 5 | 100 |
| 5 | 4 | 2.5 | 10 | 5 | 100 |
| 6 | 10 | 2.5 | 10 | 5 | 100 |
| 7 | 4 | 1.5 | 10 | 10 | 25 |
| 8 | 10 | 1.5 | 10 | 10 | 100 |
| 9 | 4 | 1.5 | 12 | 5 | 25 |
| 10 | 4 | 1.5 | 15 | 5 | 0 |
| 11 | 4 | 2.0 | 12 | 5 | 100 |
| 12 | 4 | 2.0 | 15 | 5 | 25 |
| 13 | 4 | 2.5 | 12 | 5 | 100 |
| 14 | 4 | 2.5 | 15 | 5 | 0 |
| Giant planets | |||||
| 15 | 4 | 1.5 | 8 | 5 | 0 |
| 16 | 4 | 2.0 | 8 | 5 | 0 |
| 17 | 4 | 2.5 | 8 | 5 | 0 |
| 18 | 4 | 1.5 | 6 | 5 | 100 |
| 19 | 4 | 2.0 | 6 | 5 | 100 |
| 20 | 4 | 2.5 | 6 | 5 | 100 |
| 21 | 4 | 1.5 | 4 | 5 | 0 |
| 22 | 4 | 2.0 | 4 | 5 | 0 |
| 23 | 4 | 2.5 | 4 | 5 | 0 |
| Terrestrial planets | |||||
|---|---|---|---|---|---|
| Row | No. of | Progenitor | β | a1 | Per cent unpacked |
| No. | planets | mass ( M⊙) | (au) | on post-MS | |
| 1 | 4 | 1.5 | 10 | 5 | 75 |
| 2 | 10 | 1.5 | 10 | 5 | 100 |
| 3 | 4 | 2.0 | 10 | 5 | 100 |
| 4 | 10 | 2.0 | 10 | 5 | 100 |
| 5 | 4 | 2.5 | 10 | 5 | 100 |
| 6 | 10 | 2.5 | 10 | 5 | 100 |
| 7 | 4 | 1.5 | 10 | 10 | 25 |
| 8 | 10 | 1.5 | 10 | 10 | 100 |
| 9 | 4 | 1.5 | 12 | 5 | 25 |
| 10 | 4 | 1.5 | 15 | 5 | 0 |
| 11 | 4 | 2.0 | 12 | 5 | 100 |
| 12 | 4 | 2.0 | 15 | 5 | 25 |
| 13 | 4 | 2.5 | 12 | 5 | 100 |
| 14 | 4 | 2.5 | 15 | 5 | 0 |
| Giant planets | |||||
| 15 | 4 | 1.5 | 8 | 5 | 0 |
| 16 | 4 | 2.0 | 8 | 5 | 0 |
| 17 | 4 | 2.5 | 8 | 5 | 0 |
| 18 | 4 | 1.5 | 6 | 5 | 100 |
| 19 | 4 | 2.0 | 6 | 5 | 100 |
| 20 | 4 | 2.5 | 6 | 5 | 100 |
| 21 | 4 | 1.5 | 4 | 5 | 0 |
| 22 | 4 | 2.0 | 4 | 5 | 0 |
| 23 | 4 | 2.5 | 4 | 5 | 0 |
The table reveals that our choice of 10 mutual Hill radii successfully keeps the planets ordered and packed along the MS, before unpacking occurs along the post-MS phases. Rows 9, 11 and 13 indicate that increasing β to 12 reproduces similar behaviour. When β = 15, the planets are too widely separated to ever unpack, except in a single case (see rows 10, 12 and 14). The finding that unpacking ceases when 12 < β < 15 is in line with the prediction from fig. 11 of Mustill et al. (2014) after their single-planet Hill radii is converted into the mutual-planet Hill radii used here (originally from Smith & Lissauer 2009). Doubling the value of a1, as in rows 7 and 8, still produces unpacking in all four 10-planet cases, but in just one of the four-planet cases. The likely reason is because these planets have not completed as many orbits (see fig. 6 of Veras et al. 2013a and Veras & Mustill 2013). Altering the stellar mass does not have a significant effect on the final unpacking statistics, but dramatically changes the MS lifetimes and therefore the possible extent of the WD cooling ages. Nevertheless, comparing the last column in the first six rows suggests that the amount of mass lost from 1.5 M⊙ progenitor mass stars was just not high enough to produce unpacking in every case (row 1).5
Fig. 1 provides two examples of the time evolution of individual simulations from rows 1 and 6. The plots visually demonstrate that the MS lifetimes of the 1.5 and 2.5 M⊙ stars are roughly 3 and 0.6 Gyr. Along this phase, the planets remain ordered and packed.
Now consider the upper panels in more detail. The four planets here remain in this effectively stable state throughout both the giant branch phases, when their semimajor axes are more than doubled, and through nearly 3 Gyr of WD cooling. Only after this time do the planets unpack, crossing one another's orbits. Subsequently, the planets’ orbital eccentricities are excited, and some planets sweep inwards towards the WD at distances much less than the original innermost orbit (a1 = 5 au), or even the maximum stellar radius (Rmax = 1.42 au). Eventually, the initially outermost planet enters the Roche radius of the WD (Rroche = 0.0035 au) at 12.2 Gyr, corresponding to a cooling age of about 9 Gyr.
Contrastingly, in the bottom panels, the initially packed system of 10 planets becomes unpacked during the giant branch phases, and specifically during the asymptotic giant branch phase.6 With 10 planets and a high maximum amount of total mass lost, the stability boundary is more easily critically shifted than in the four-planet case with a smaller maximum amount of mass lost. During the entire WD phase, lasting over 13 Gyr, the inner system is swept through by planets on eccentric orbits. Three of these planets hit the Roche radius of the WD, at 5.61, 6.19 and 9.72 Gyr, at which point the planets are absorbed into the WD and increase its mass. One of the planets (red) which survives for the age of the Universe drifts within 0.05 au of the WD at 10.68 Gyr. A planet with such a close pericentre has a non-negligible probability of being detected while transiting the WD. If the pericentre is close enough to the WD for tidal circularization to occur, then prospects for detection increase dramatically. As the orbit circularizes, the geometric transit probability remains the same, while the transit duration increases and the orbital period decreases (Barnes 2007; Socrates et al. 2012); the latter effects provide the enhancement in detectability.
All instabilities manifested themselves as planet–planet collisions or planet–star collisions; no planets were ejected. Ejections require planets to be scattered out to large distances from the star, and these distances were not attained with the terrestrial planets simulated here. Our integrator faithfully tested for ejections by incorporating the true Hill escape ellipsoid in the solar neighbourhood, with semimajor axes on the order of 105 au (Veras & Evans 2013b). Therefore, planets flung out on orbits of several hundred au are retained (see e.g. the bottom-right panel of Fig. 1). Further, Ford & Rasio (2008) suggest that the smaller the planets, the less likely ejection is to occur. Finally, stellar mass-loss by itself is not fast and great enough to cause escape due to both the physical properties of the stars simulated and the planetary orbits considered (Veras et al. 2011). Of the instabilities across all terrestrial planet simulations, three quarters were engulfments into the star, and one quarter were planet–planet collisions. Planet–planet collisions are a plausible avenue to generate large amounts of fresh small body debris, which could explain strong metal pollution detected in a number of old, cool WDs (Koester et al. 2011).
We can improve our understanding of the dynamical evolution in both panels by focusing on the planetary semimajor axis changes throughout the star's lifetime (right-hand panels of Fig. 1). Although a system's total orbital energy includes eccentricity-dependent and inclination-dependent planet–planet interaction terms (see e.g. equation 2.27 of Veras 2007), we can use the other terms, whose only orbital parameters are semimajor axes, as excellent proxies (see fig. 19 of Veras et al. 2013a) for energy transfers in the system. Further, after the star has become a WD, and mass-loss ceases, energy is conserved for the remainder of the system's existence. In Fig. 1, the semimajor axes ‘jump’ around in a vaguely sawtooth-like manner, indicating shifts in orbital energy due to instances of strong interactions with the other planets.
In the upper-right panel, after an initial period of orbit crossing lasting a few Gyr, the planets’ orbits remain largely ordered, except for one notable orbit crossing at just over 11 Gyr. In no instance is the semimajor axis of a planet perturbed down to values within a few au. This result is expected because when the planets’ orbits are ordered and not undergoing strong interactions, they are evolving secularly, and secular interactions do not alter semimajor axes. Also, during scattering the apocentre must remain high in order to interact with the other planets. Hence, the semimajor axis cannot be reduced by much more than 50 per cent. Therefore, the small pericentres achieved are due primarily to changes in orbital eccentricity and not semimajor axis. Nevertheless, this analysis neglects the effects from tidal circularization, which may effectively ‘grab’ a planet at its pericentre and keep that pericentre to within a factor of 2 while shrinking the orbital semimajor axis.
Regardless of semimajor axis, the minimum planetary pericentre is the value of greatest interest here because of (i) how close-orbit planets might potentially perturb asteroids into the Roche radius of the WD and (ii) the possibility of detecting these planets (Agol 2011; Fossati et al. 2012). Therefore, we present a collection of these minimum pericentres in Fig. 2 for all simulations that became unpacked and for all planets that have survived for the age of the Universe. For each of seven different numberlines (sets of initial conditions), at least one simulation featured an incursion within 0.1 au of the star. In one 10-planet simulation (row 2), five of the surviving planets each individually drifted to within 5 au of the star. Each numberline demonstrates a wide variation of behaviour across the simulations, perhaps emphasizing the highly chaotic nature of unpacked planetary systems.
2.3.2 Selected giant planet simulations
We performed some additional simulations, this time with packed giant planets, in order to determine how the process of unpacking changes for a higher planet/star mass ratio. In these simulations, we assumed a WD Roche radius based on a typical giant planet density of 1.0 g cm−3. There are two particularly strong cases of packed giant planets which act as observational motivation: our Solar system and the HR 8799 system (Marois et al. 2008, 2010).
We found that closely packed giant planets which remain stable throughout the MS before undergoing instability require β < 8. These lower values differ from the terrestrial case probably because the Hill sphere is no longer the relevant parameter governing the dynamics; see equation 2 of Faber & Quillen (2007) for the explicit dependence on planet mass. The predictions for stability of post-MS systems based on this equation seem to agree with the results presented in this paper for giant planets too (see fig. 11 of Mustill, Veras & Villaver 2014). We also find that giant planet unpacking is more violent than terrestrial planet unpacking, resulting in less orbit meandering, speedy dynamical settling following instability and escaping planets.7
Fig. 3 illustrates two specific examples. For both cases, β = 6, and the instabilities in the system occur during the WD phase, after the stellar radius has been reduced to approximately Earth-size. In both cases, two planets are ejected, and one remaining planet is perturbed on to a wide orbit. The other planet experiences eccentric excursions close to the WD, and may be tidally circularized depending on the internal properties of the planet and its orbital pericentre.
Two examples of post-MS planetary unpacking with Jovian planets, from rows 19–20 of Table 1. Both panels feature four planets with initial semimajor axes between 5 and 17.38 au (left-hand panel) or 5 and 15.86 au (right-hand panel). In both cases, instability occurs only after the star has become a WD and already shrunk its radius down from its maximum value Rmax to approximately that of the Earth. In the left-hand panel, the inward incursions of the initially innermost planet are periodic with increasing penetration. Alternatively, in the right-hand panel, the incursions are more frequent and irregular. Neither GR nor tidal forcing is included in these simulations, which may affect the evolution of the simulation in particular in the left-hand panel.
The manner in which the incursions take place can be periodic and well defined (left-hand panel) or irregular and frequent (right-hand panel). The oscillations in the left-hand panel have an increasing amplitude, suggesting a way for giant planets to potentially collide with a WD after 10 Gyr of evolution. An event of this kind would give rise to a noticeable optical (and probably X-ray) transient, and result in a dramatic amount of metal pollution of the WD photosphere. Such large and periodic eccentricity variations are classic signatures of Lidov–Kozai oscillations (see section 7 of Davies et al. 2013 for a review); Fig. 4 confirms this suspicion with accompanying inclination oscillations. The changing amplitude may be indicative of the role of the octupole term in Lidov–Kozai evolution (see fig. 10 of Naoz et al. 2013). Of note is that when the orbit pericentre is at a minimum, its inclination is also at a minimum. Then, we might expect a detectable planet in this case to harbour an inclination at the critical Lidov–Kozai value, around 40°. Because this mechanism is strongly dependent on the evolution of the argument of pericentre, the evolution of this particular system may have been different had we included GR in our simulations.
Lidov–Kozai oscillations of inclination from the left-hand panel of Fig. 3. The interplay between the inclination and eccentricity variations results in close pericentre passages around the WD in this one case when two planets survive instability.
If tidal circularization is an active process in WD systems, then the inclinations of highly eccentric planets will be preserved as they are circularized. The reason is that inclination damping is primarily due to tides raised on the star, which are weak for WDs (as they are for MS stars). Consequently, close-in planets orbiting WDs might harbour a broad range of inclinations.
3 DISCUSSION
Our results demonstrate that planets can be perturbed on to orbits taking them close to the WD at late ages, after an entirely quiescent existence. This finding is particularly important in the context of observations of old WDs with extreme metal pollution (Koester et al. 2011 and fig. 8 of Koester et al. 2014) or orbiting debris discs (e.g. Farihi et al. 2011). If asteroids are the progenitors of the pollution and the discs, then the asteroids must be flung towards the WD, and probably by a planet.
One effective way to deliver asteroids is with a planetary mean motion resonance (Debes et al. 2012). Inside of a resonance, asteroid eccentricities can be pumped high enough to eventually achieve a WD-skimming orbit. This physical process, when applied to a single planet and an asteroid or Kuiper belt, produces a gradually decreasing stream of asteroids due to a depleting reservoir. Only a massive-enough belt would be able to sustain a high-enough delivery rate to explain accretion in WDs which are many Gyr old.
An alternative to resonant diffusion is direct scattering. However, a single planet will never change its orbit, except during giant branch mass-loss, due to interactions with a relatively massive planetesimal belt, or due to a flyby star. Therefore, scattering between a single planet and an asteroid belt will occur at young WD cooling ages (Bonsor et al. 2011; Frewen & Hansen 2014). Penetrative stellar flybys can scatter the planet, thereby refreshing planet–asteroid dynamics in the system, but only for favourable geometries, and infrequently (Veras & Moeckel 2012).
Neither single-planet scattering nor resonant diffusion commencing when the WD is born can easily explain pollution or discs which exist around WDs that are many Gyr old. Multiplanet systems, in concert with smaller bodies, may represent a solution. Two-planet systems are likely to be inadequate, because after an instance of instability, the surviving planet will have a fixed orbit, akin to the one-planet case. Although this orbit can be highly eccentric and penetrate to within the maximum radius of the star (Veras et al. 2013a), the orbit is periodic, and will too quickly exhaust its supply of fortuitously located asteroids. Three-planet systems are more promising because after instability, two surviving planets can exist in ever-changing orbits for the entire WD cooling age (Mustill et al. 2014), similar to the behaviour here in the right-hand panel of Fig. 3. Nevertheless, if these two surviving planets are secularly evolving, they will experience predictably periodic changes in their orbits (left-hand panel of Fig. 3).
Systems with at least four planets appear to harbour dynamics which are rich enough to potentially generate dynamical interactions with asteroids over long timespans and over an extensive range of physical space. Given the four terrestrial planets in the Solar system, and the mounting discoveries of packed exosystems, such configurations are perhaps common. Indeed, over half of all stars in the Milky Way are thought to contain at least one terrestrial planet (Cassan et al. 2012), and perhaps an average of at least five (Ballard & Johnson 2014). This average may be even higher considering that Ballard & Johnson (2014) focused on M dwarfs, which have much lower masses than the stars considered here. Further, the asteroidal material in the remnant WD system may be widely dispersed and have a size distribution skewed towards small bodies due to rotational fission from giant branch radiation (Veras, Jacobson & Gänsicke 2014c). Consequently, rocky material may be available in enough locations to be accessed at different times during WD cooling.
Our simulations also provide a potential evolutionary pathway – a dynamical history – for eccentric planets which could be detected transiting WDs within 0.1 au. If these planets fail to be tidally circularized, then they would not remain in a secularly stable compact configuration for Gyr. Instead, their pericentres would drift by orders of magnitude, greatly reducing their prospects for habitability. Consequently, observations would have to ‘catch’ a system at just the right point in its secular evolution (e.g. Veras & Ford 2009). We suggest that detectable eccentric planets would likely be accompanied by an unseen outer companion. The prospects for detection will increase with current and upcoming missions such as Gaia, LSST and PLATO, which will survey 105–106 WDs tens to thousands of times.
4 CONCLUSIONS
We find that the post-MS peregrinations of terrestrial planets which were packed and quiescent throughout the MS cover a large volume of space at both early and late WD cooling ages. These planets often but irregularly invade regions which are well interior to their nascent orbits and potentially contain asteroidal material that remained largely unperturbed for the entire MS and giant branch phases of stellar evolution. Consequently, these disturbances could represent catalysts for dynamical interactions which lead to the formation of debris discs and/or atmospheric pollution in old WDs. The planets themselves, when travelling through the pericentres of their orbits during epochs of high eccentricity, should be detectable and might be tidally circularized, but are unlikely to harbour life.
We thank the referee for a thorough and thoughtful review, explicit and constructive comments, and an expeditious response to our submission. We also thank Alex Bowler and Roberto Silvotti for pointing out typographic errors. 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 n. 320964 (WDTracer).
Theoretical models constrain the mass of object WD 0806−661 b, which orbits a WD at a distance of approximately 2500 au, to the planetary mass regime (Luhman, Burgasser & Bochanski 2011).
Perhaps related is Marzari (2014)'s frequency map analysis finding that the dependence of the Hill boundary on planet mass has not yet fully been accounted for.
See the Exoplanet Data Explorer at http://exoplanets.org
More precisely, in one case the mass-loss failed to shift the stability boundary sufficiently enough to cause orbit crossing either during a giant branch or WD phase.
None of our simulations were destabilized along the red giant branch phases.
The stronger close encounters likely decrease the accuracy of the simulations; here, typically total angular momentum is conserved to only 10−4–10−3. These values indicate that future studies modelling full-lifetime simulations containing strong interactions amongst many giant planets would benefit from either tiny Bulirsch–Stoer tolerances (<10−13) – prolonging the simulations even further – or by utilizing other codes. One may combine planetary and stellar evolution codes in the amuse package (Pelupessy et al. 2013), or modify an existing planetary integrator that emphasizes greater speed and accuracy (Grimm & Stadel 2014; Rein & Spiegel 2015).
