The post-main-sequence fate of the HR 8799 planetary system

The noteworthy four-planet HR 8799 system teeters on the brink of gravitational instability and contains an A-type host star which is characteristic of the progenitors of the majority of known white dwarf planetary system hosts. Gozdziewski and Migaszewski (2020) have demonstrated that the system can retain all four planets for at least 1 Gyr along the main sequence if the planets evolve within an externally unperturbed 8:4:2:1 mean motion resonance configuration. Here we propagate forward their most stable fit beyond the main sequence, and incorporate external effects from Galactic tides and stellar flybys. We find that (i) giant branch mass loss always breaks the resonance, and usually triggers the ejection of two of the planets, (ii) stellar flybys and Galactic tides rarely break the resonance during the main-sequence and giant branch phases, but play a crucial role in determining the final planetary configurations around the eventual white dwarf host star, and (iii) the meanderings of the surviving planets vary significantly, occupying regions from under 1 au to thousands of au. The ubiquitous survival of at least one planet and the presence of the debris discs in the system should allow for dynamical pathways for the white dwarf to be metal-polluted.


INTRODUCTION
Over a decade after its discovery (Marois et al. 2008(Marois et al. , 2010, the HR 8799 planetary system remains a benchmark of exoplanetary science. The multifaceted appeal of the systemthe first imaged multi-planet exosystem, and one that might be in the process of destabilizing -arises from both its observational and dynamical properties.

Planets
The HR 8799 system is comprised of a λ Boo A5V host star with close proximity to the Sun (41.29±0.15 pc, Gaia Collaboration 2018). Three co-moving companions were initially identified by Marois et al. (2008) at projected orbital separations of 68, 38 and 24 au over multiple epochs using adaptive optics and coronagraphy at the Gemini ⋆ E-mail: d.veras@warwick.ac.uk † STFC Ernest Rutherford Fellow and W.M. Keck Observatories. A host of age indicators for HR 8799, including its galactic space motion (indicating a high likelihood of being a member of the Columba moving group) and the placement of the star in a colourmagnitude diagram suggest a best-estimate age of ≈30 Myr (Zuckerman et al. 2011;Malo et al. 2013).
Based on this age estimate, as well as the measured near-infrared brightnesses of the co-moving companions, Marois et al. (2008) assigned mass estimates of 7 +4 −2 , 10 +3 −3 , and 10 +3 −3 MJup for the HR 8799bcd planets, respectively. Subsequent observations of this system at the more favourable observing wavelength of 3.8 µm (Marois et al. 2010) unveiled the direct detection of a fourth planet ("e") at a significantly smaller projected orbital separation of 14.5±0.4 au, and an estimated mass of 7 +3 −2 MJup, similar to that of the "c" and "d" companions.
The relatively wide angular separations on the sky (ranging from ≈ 0.4 − 1.7 ′′ ) and high intrinsic luminosity of the four planets has allowed for some of the most detailed characterization of the planets' atmospheres through direct spectroscopy using instruments that are dedicated to exoplanet imaging (e.g. GPI, SPHERE, Ingraham et al. 2014;Zurlo et al. 2016;Greenbaum et al. 2018). Also, multipurpose instruments equipped with much higher spectral resolution (λ/∆λ∼4000, Konopacky et al. 2013) have indicated the presence of vertical atmospheric mixing when combined with detailed atmospheric models (Baruteau et al. 2011). More recently, the GRAVITY instrument at the Very Large Telescope has revealed the ability of using interferometric measurements to measure a spectrum with resolution R ∼ 500, a value which is higher than that of typical dedicated exoplanet imaging instruments such as GPI and SPHERE. The GRAVITY measurements returned precise atmospheric composition measurements of HR 8799 e, allowing for the derivation of a remarkably precise value of the atmospheric Carbon-to-Oxygen ratio of 0.60±0.08 (Mollière et al. 2020), suggesting that the planet formed outside of the CO2 or CO ice line.
The 13-year time baseline since the discovery of the HR 8799 planets has allowed for careful astrometric monitoring of the planets' orbits with instruments that achieve high contrast (e.g. Maire et al. 2015;Pueyo et al. 2015;Konopacky et al. 2016;Wang et al. 2018). Nevertheless, some of the tightest astrometric constraints originate from 1998 Hubble Space Telescope observations of the system. Remarkably, when novel image processing techniques were applied to this 1998 dataset (Soummer et al. 2012), the locations of the three outermost planets HR 8799 bcd (Soummer et al. 2011) were revealed, providing a 20 year time baseline for astrometric monitoring. Most recently, GRAVITY has returned astrometry of the "e" planet relative to the star with precision of ≈100µas, placing strong constraints on the orbit, and disfavouring perfectly coplanar orbits (Gravity Collaboration et al. 2019).

Debris belts
In addition to containing four planets, the system also hosts belts of debris. Following the early identification with IRAS of HR 8799 being dust-rich (Sadakane & Nishida 1986), more recent observations of this system in the mid-infrared (Matthews et al. 2014) and submillimetre (Holland et al. 2017) have provided a more complete picture of the underlying structure of the population of dust and planetesimals.
Specifically, by using spectroscopy from the Spitzer Space Telescope, Su et al. (2009) were able to identify a debris structure comprised of two components. The first is a warm dust component (T ≈ 150 K) at 6-15 au, which is located interior to the orbit of the HR 8799 e planet. The second is a much more spatially extended zone (≈90 to ≈300 au) of cooler dust (T ≈ 45 K) extending outward from the orbit of the "b" planet. Hence, both debris components flank all four planets. Thus, with four giant planets with orbits residing in the cleared region between two dust belts, it has been suggested that the debris structure within the HR 8799 system resembles a "scaled-up" version of our our own solar system (e.g. Su et al. 2009;Hughes et al. 2018).
Importantly, careful characterization of the HR 8799 debris structure, combined with dynamical arguments, have placed constraints on the masses of the individual planets. Wilner et al. (2018) used 1.3 mm data from the Submillimeter Array along with archival ALMA data (Booth et al. 2016) to constrain the location of the inner edge of the outer debris belt to 104 +8 −12 au. This much higher precision on the location of the belt edge, along with the sensitive mass dependence of HR 8799 b on the spatial extent of the chaotic zone of the dust population, placed a strong constraint on the mass (5.8 +7.9 −3.1 MJup) of HR 8799 b. This mass is consistent with the measurements of Marois et al. (2008), who instead used evolutionary models of the luminosity of HR 8799 b.

Dynamical evolution of HR 8799
The architectures revealed by the observations have dynamical significance because of the very high masses of all four planets (≈ 5 − 10MJup), the wide planet-star separations of all four planets (≈ 15−70 au) and the relatively close planetplanet separations with respect to their stability boundaries. Initially, questions arose about how these planets were formed, particularly through fragmentation as opposed to core accretion (Dodson-Robinson et al. 2009;Kratter et al. 2010;Meru & Bate 2010). Further, even before the discovery of HR 8799 e, mean motion resonances between the planets were suspected as the driving forces behind the stability of the system (Goździewski & Migaszewski 2009;Reidemeister et al. 2009;Fabrycky & Murray-Clay 2010;Marshall et al. 2010).
The discovery of HR 8799 e (Marois et al. 2010) has prompted additional theoretical explorations of the formation and evolution of the system. Because disc fragmentation through gravitational instability is the favoured formation mechanism, this system has directly motivated several studies about the details of the process (e.g. Baruteau et al. 2011;Boss 2011;Vorobyov & Elbakyan 2018). The apparent fragility of the system has also been investigated in the context of stellar birth clusters; HR 8799 might have been lucky to survive the cluster phase with its four planets intact (Li et al. 2020). Post-cluster evolution investigations not only considered the interplay between resonant behaviour and stability amongst the four planets (Goździewski & Migaszewski 2014Götberg et al. 2016;Morrison & Kratter 2016), but also the interactions between the planets and the debris belts in the system (Moore & Quillen 2013;Goździewski & Migaszewski 2014Götberg et al. 2016;Morrison & Kratter 2016).
However, the fate of the system beyond the mainsequence phase has not yet been investigated in detail. One potential reason is because the system is a rare example of a known exosystem which is not guaranteed to remain stable until the end of the main sequence. In fact, identifying a long-term stable solution given the observational constraints has been challenging. In a series of papers, Goździewski & Migaszewski (2009 have been updating and refining such a solution. Only in the most recent iteration (Goździewski & Migaszewski 2020) have the authors obtained an exactly periodic configuration which guarantees stability in the absence of external perturbations (e.g. from the Galactic environment) or internal perturbations (e.g. from stellar evolution).
This solution allows us here to explore the post-mainsequence evolution of the system, albeit under a limiting "most stable" assumption. The primary motivation for this study is simply to consider this unexplored aspect of the fate of this benchmark system. We do, however, also have a secondary motivation, one which is relevant to the mounting discoveries of planetary systems orbiting white dwarfs.

An exemplar progenitor white dwarf planetary system
Over 1000 white dwarf planetary systems are known through planetary debris that is detected in the photospheres of the stars (Dufour et al. 2007;Kleinman et al. 2013;Kepler et al. 2015Kepler et al. , 2016Coutu et al. 2019). Given the frequency of this debris, accretion rates onto the white dwarfs, and planetary mass in the surrounding circumstellar environment (Koester et al. 2014;Vanderburg et al. 2015;Farihi 2016;Rappaport et al. 2016;Gurri et al. 2017;Manser et al. 2019;Guidry et al. 2020;Vanderbosch et al. 2020), we assume that the debris arises primarily from the destruction of minor planets like asteroids as opposed to major planets like exo-Jupiters or exo-Earths. Nevertheless, we now know that such large planets can reach separations of under 0.1 au Vanderburg et al. 2020).
Overall, major planets, minor planets, dust, gas and metallic debris have all been discovered around or in white dwarfs; for a recent review, see Veras (2021). The HR 8799 system is an exemplar progenitor white dwarf planetary system, for several reasons. First, the host star is an A-type main-sequence star, which represents the predominant progenitor stellar type for the known population of white dwarf planetary systems (Tremblay et al. 2016;Cummings et al. 2018;El-Badry et al. 2018;McCleery et al. 2020;Barrientos & Chanamé 2021). Further, the system contains multiple giant planets which will escape the reach of the host star's giant branch transformation (Mustill & Villaver 2012;Adams & Bloch 2013;Nordhaus & Spiegel 2013;Madappatt et al. 2016;Ronco et al. 2020). These planets are also on the verge of instability, and instability during the white dwarf phase is necessary to perturb objects close to and later onto the star (Debes & Sigurdsson 2002;Veras et al. 2013a). Finally, the HR 8799 system features multiple debris discs, and these could supply the metals that will be accreted by the white dwarf through interactions with the planets (Bonsor et

Plan of paper
In this paper, we model the full-lifetime post-cluster evolution of the four planets in the HR 8799 system by assuming that they currently reside in the stable, exact 8:4:2:1 mean motion resonance orbital solution given by Goździewski & Migaszewski (2020). We include both internal and external forces which could destabilize the system in the form of stellar evolution and Galactic forces. In Section 2 we describe the orbital solution of Goździewski & Migaszewski (2020). Basic stability considerations are discussed in Section 3. Our numerical simulation setup and results are then presented in, respectively, Sections 4 and 5. We then conclude with a discussion (Section 6) and summary (Section 7).

AN EXACT RESONANCE ALONG THE MAIN-SEQUENCE
We adopt the best-fitting orbital solution of Goździewski & Migaszewski (2020) -which ensures stability of the four planets in the absence of extra forcesand list the solution in Table 1. This configuration places the four planets in resonance immediately, such that at least one resonant angle between any two pairs of planets librates. In fact, the highlighted resonance achieved is a generalized 8:4:2:1 Laplace resonance, reflecting the approximate orbital period ratios of planets b:c:d:e. More technically, this resonance expresses that the angle λe − 2λ d − λc + 2λ b librates, where λ denotes mean longitude.
The resonance is "exact" in the sense of the configuration being strictly periodic with canonical elements rather than the amplitude of the Laplace resonance angle equalling zero. The libration amplitude is actually approximately 4 • . Goździewski & Migaszewski (2020) achieved the exact periodicity by appealing to the field of periodic orbit families, and in particular the four-planet case (Hadjidemetriou & Michalodimitrakis 1981). Goździewski & Migaszewski (2020) derived this solution by assuming that the mass of the star HR 8799 equals Main-sequence evolution 1.52M⊙ (Konopacky et al. 2016). For our study, this assumption may be crucial because the main sequence mass helps determine the duration and extent of giant branch mass loss, which in turn affects the stability of the system. Further, varying this mass, even by a few hundredths of a solar mass, may break the resonant configuration even on the main sequence. In this respect, the main-sequence lifetime then might affect the evolution.
Consequently, we have adopted three different stellar masses, straddling the uncertainties (1.52 ± 0.15M⊙ from Konopacky et al. 2016), although for most of our simulations we adopt the fiducial mass of 1.52M⊙. We use the SSE code to model the stellar evolution (Hurley et al. 2000). From this code, for main-sequence stellar masses of 1.37, 1.52, 1.67M⊙ respectively, we obtain main sequence lifetimes of (3.36, 2.62, 1.98) Gyr, giant branch lifetimes of (0.391, 0.284, 0.222) Gyr, and total fractions of mass lost during the giant branches of (0.591, 0.619, 0.643).
The current age of HR 8799 is uncertain, but as previously mentioned has been estimated to be on the order of tens of Myr. Our simulations are not highly sensitive to the current main-sequence age because the exactly periodic configuration is indefinitely stable in the absence of extra forces. Hence, we simply initiate HR 8799 at the zero-age main-sequence phase of evolution as a conservative limit.

BASIC STABILITY CONSIDERATIONS
Having established the physical and orbital parameters that we will adopt for our simulations, now we discuss potential triggers for instability. One trigger is stellar flybys, which represent a potential danger throughout the lifetime of the system. The frequency and impact parameters of each flyby is unknown and probabilistic, but can be estimated through either analytic arguments or numerical simulations.
For the HR 8799 system, the planet most susceptible to perturbations from flybys is HR 8799 b, the outermost planet. If this planet is sufficiently perturbed away from the exact periodic configuration, the result may be a (not necessarily immediate) scattering event. Another potential, but more unlikely, trigger for instability are Galactic tides. Nevertheless, we incorporate both stellar flybys and Galactic tides in our simulations during all phases of stellar evolution.
Those perturbations are external; the typically stronger perturbations are internal to the system, from the host star as it leaves the main sequence. The resulting increase in stellar radius and luminosity have just a negligible effect on the HR 8799 planets because of their wide orbits. However, the mass lost through stellar winds alters the gravitational potential of the system, expanding and stretching the planetary orbits (Omarov 1962;Hadjidemetriou 1963).
The resulting variations in orbital eccentricity and argument of pericentre for a single planet within about 10 3 au due to stellar mass loss can usually be considered negligible (Veras et al. 2011). However, the planet's semimajor axis would increase in proportion to the amount of mass lost: for our fiducial HR 8799 host star mass of 1.52M⊙, a mass loss fraction of 62 per cent would correspond to a semimajor axis increase factor of 2.6. If the mass loss is isotropic -an assumption that we adopt -then the inclination and longitude of ascending node of the orbit would remain fixed (Veras et al. 2013b;Dosopoulou & Kalogera 2016a,b).
In systems with more than one planet, stellar mass loss shifts stability boundaries such that a stable configuration on the main sequence becomes unstable on the giant branches or white dwarf phase (Debes & Sigurdsson 2002;Veras et al. 2013a;Voyatzis et al. 2013;Mustill et al. 2014 Figure 3. A representative evolution profile of semimajor axes and orbital pericentres, featuring the ejection of two planets (here planets b and e) just after the tip of the AGB, and the aftermath during the white dwarf phase. The abbreviations "MS", "GB" and "WD" stand for main sequence, giant branch, and white dwarf respectively.

Representative case
For systems of at least four planets (like HR 8799), instability in the entire system is easily triggered if any one pair of planets become unstable Veras et al. 2016;Maldonado et al. 2021). Only if each planet is sufficiently separated from one another -such as is the case with Jupiter, Saturn, Uranus and Neptunecan the system remain stable throughout the giant branch and white dwarf phases (Duncan & Lissauer 1998;Veras 2016a,b;Zink et al. 2020).
Although the four HR 8799 planets are spread out more than Jupiter, Saturn, Uranus and Neptune, the significant difference in masses between the two sets of planets suggests that the HR 8799 planets might actually be more tightly "packed" (e.g. Raymond et al. 2009). Hence, whether all four of these planets can survive giant branch evolution is not immediately clear, particularly if they are stabilized through a resonance. Now we perform numerical simulations to explore this evolution.

NUMERICAL SIMULATION SETUP
In order to perform our numerical simulations, we use the RADAU integrator within the N -body code developed in Veras et al. (2013a) and subsequently modified in Mustill et al. (2018). This code is originally based on the Mercury integration suite (Chambers 1999), but with the crucial difference that stellar evolution variations (with profiles from SSE, Hurley et al. 2000, as described earlier) are incorporated into the integrations. The code also includes effects from stellar flybys and Galactic tides.
The prescription we used for Galatctic tides is taken from Veras & Evans (2013), and assumes that the HR 8799 planetary system resides at an approximate distance of 8 kpc from the Galactic centre. For stellar flybys, we adopted the prescription from Veras et al. (2014a) but with an updated stellar mass function, stellar density and field encounter distribution from Martínez- Barbosa et al. (2017). The random introduction of stellar flybys provides a stochastic component to each simulation. Hence, although the initial masses and orbital parameters of each planet are equivalent across all simulations, the introduction of flybys creates a very slight perturbation which can change the fate of the system.
When considering full lifetime simulations of planetary systems, one must balance accuracy and speed. In order to keep the HR 8799 planets locked in resonance throughout the main sequence, the integrator accuracy must be sufficiently high. Hence, we sampled four different accuracy tolerance values of 10 −10 , 10 −11 , 10 −12 , 10 −13 both to tease out any potential dependence, and to create another initial differentiating factor amongst different simulations. For each of these values and an initial mass of HR 8799 of 1.52M⊙, we ran ten simulations for a total duration of 14 Gyr. We also ran ten simulations for each of the other two adopted initial stellar masses (1.37M⊙ and 1.67M⊙) assuming an accuracy tolerance of 10 −12 . Finally, in order to track resonant behaviour due to mass loss, we ran one set of ten simulations with a very high output frequency (2340 yr) between 2.739 Gyr and 2.973 Gyr for an initial stellar mass of 1.52M⊙ and an accuracy tolerance of 10 −10 .
In all simulations, we set the effective white dwarf physical radius to be 1R⊙. This value, while hundreds of times larger than the star's actual physical radius, is a representative value of the white dwarf's Roche sphere. Any planet entering this Roche sphere will be destroyed. For the gas giants of the HR 8799 system, this value is conservative by a factor of a few (Veras & Fuller 2020), but is large enough to prevent prohibitive timestep issues with the integration.

NUMERICAL SIMULATION RESULTS
We break down our results into both notable individual system cases and ensemble system statistics. First, we highlight the result that in no case modelled did the resonant configuration break before the giant branch phases. Although Galactic tides and stellar flybys always altered the planets' orbital parameters, these variations were at a small enough level to maintain the resonance. Even changing the initial stellar mass from 1.52M⊙ to 1.37M⊙ or 1.67M⊙ did not break the resonance along the main sequence. Although altering the stellar mass rescales all of the planet-to-star mass Three-survivor case Figure 4. A rare case where three planets survive. Although all planets initially reside on coplanar orbits, stellar flybys create slight non-coplanarities. These are then exacerbated due to angular momentum exchange during the scattering event. ratios by the same factor, maintenance of resonance lock is not necessarily guaranteed; for HR 8799, however, resonance lock appears to be a sufficiently weak function of the planet-to-star mass ratio.

Individual results
We first show some representative evolutions, before considering particularly interesting cases. Figure 1 illustrates the astrocentric geometry of the planetary orbits throughout the main sequence, at 10 Myr snapshots. The exactly periodic configuration assumes co-planarity, which is only slightly broken by perturbations from stellar flybys. This slight shift Role of stellar flybys is important because that can instigate a larger orbital plane divergence after a scattering event. Otherwise, because all of the orbital angular momentum for an exactly coplanar configuration would be normal to the orbital plane, in this case the planets' individual orbits would remain coplanar, and the incidence of planet-planet collisions would be artificially high.
Next, in Fig. 2, we show a representative case of the evolution of the resonant angle λe − 2λ d − λc + 2λ b as stellar mass is lost during the giant branch phases, for our fiducial initial stellar mass value (1.52M⊙). The time resolution of the plotted data points is 2340 yr, a small enough value to accurately assess resonant behaviour with respect to the giant branch timescale. The timespan of the plot covers the end of the red giant branch (RGB) phase, as well as the entire asymptotic giant branch (AGB) phase. The slight bump in the middle of the plot corresponds to the tip of the RGB, after which about 4 per cent of the stellar mass has been lost. The resonance breaks towards the right side of the plot, when instability sets in; this instability coincides with the significant and rapid mass loss that is characteristic of superwinds at the AGB tip (Vassiliadis & Wood 1993).
The most common outcome of the instability and resonant breaking is the ejection of two of the planets. The final planetary orbits of the survivors differ significantly from the final orbits which would be expected from mass loss-induced orbital expansion alone. Now we show some examples. In the following figures, we use osculating Jacobi coordinates to represent the orbital elements of the planets. We also highlight the orbital pericentre, which indicates the inner reach of a planet and hence its ability to and timescale for perturbing minor planets and debris into the eventual white dwarf.

Representative case
One representative orbital evolution is shown in Fig. 3. The semimajor axes of all four planets increase during the giant branch phases before an instability occurs at the tip of the AGB. The result of the instability is that planets b and e are ejected, and planets c and d swap orderings and Strongly nonuniform evolution are perturbed inward and outward. The semimajor axis of planet d becomes hundreds of au, making it more susceptible to Galactic tides and stellar flybys. The orbital pericentre of planet c becomes just 3-7 au, which is over 5 au smaller than the pericentre of planet e along the main sequence. The orbital evolution of both planets c and d remains relatively uniform throughout the white dwarf phase.

Three-survivor case
Next we show a rarer outcome, one where three planets survive giant branch evolution (Fig 4). In this case, a stellar flyby at 144 Myr (early on the main sequence) created a perturbation which was not strong enough to break the resonance, but strong enough to imprint a dynamical signature that manifests itself at the tip of the AGB. This signature allows planet c to be ejected, and all of the other planets to survive and maintain their ordering. Further, this flyby is strong enough to break the coplanarity of the orbits on the main sequence by a sufficient amount (although not noticeable on the plot), such that the instability triggers large (tens of degrees) oscillations in inclination and longitude of ascending node of the survivors.

Role of Galactic tides
If the surviving outermost planet in a system is sufficiently distant, then Galactic tides start to play a role. In Fig. 5, planet c is perturbed outward and subsequently harbours a semimajor axis of about 2,000 au. This value is high enough (Bonsor & Veras 2015) to create a gradual pericentre drift. The surviving inner planet, planet b, experiences secular nonuniform oscillations in eccentricity, which are manifested in the orbital pericentre plot. Because the oscillations are nonuniform, the planet may dynamically access different reservoirs of remnant planetary debris and minor planets throughout white dwarf evolution. The reason for the high amplitude secular oscillations is the von   Zeipel-Lidov-Kozai effect, which effectively exchanges angular momentum between the eccentricity and inclination (Hamers & Portegies Zwart 2016;Petrovich & Muñoz 2017;Stephan et al. 2017Stephan et al. , 2018Stephan et al. , 2020Muñoz & Petrovich 2020;O'Connor et al. 2021).

Role of stellar flybys
The strength of a stellar flyby cannot be characterized by the minimum close encounter distance with the star alone; in fact, the smallest closest encounter distance in any of our simulations (76 au) had no noticeable effect on the dynamics of that system. The more important measure is the direction-dependent closest proximity to any of the planets. Figure 6 illustrates a clear consequence of a stellar flyby, at a time of about 10 Gyr, on two surviving planets with semimajor axes at about 40 and 500 au. The flyby (with a closest encounter distance of about 785 au) increases the range of inclinations by one order of magnitude, albeit just from hundredths to a tenth of a degree. The figure also illustrates how much the initial inclination values are changed upon reaching the white dwarf phase.

Strongly nonuniform evolution
Example evolutions for the alternate host star mass cases are similar to the figures presented above. In one notable 1.67M⊙ case (Fig. 7), planets c and d survive the scattering event, and planet c is perturbed into an orbit with a semimajor axis of about 8,000 au and an eccentricity near unity. This planet is severely affected by Galactic tides, and the resulting perturbations on planet d, with a semimajor axis of about 30 au, produce a highly nonuniform orbital pericentre profile. The nonuniform radial incursions of both planets could potentially access and disperse reservoirs of minor planets towards the white dwarf at any time, including very late (≈ 10 Gyr) times. Accretion rates for such old white dwarfs have been measured (Hollands et al. 2018).

Ensemble results
The wide variety of outcomes for individual simulations can be represented through a few ensemble illustrations. One relatively common outcome is the ejection of two of the four planets (Fig. 8). The instabilities predominately produce ejections; rarely will a planet intersect the white dwarf's Roche radius, break up and be engulfed.
Due to angular momentum balance, the ejection event is likely to perturb a surviving planet much closer to the white dwarf than at any point during the main sequence or giant branch phases. Figure 9 illustrates the final orbital pericentre versus semimajor axis for all planets which survive for a 14 Gyr timespan. Although the semimajor axes of nearly all these planets exceed 10 au, their orbital pericentres extend inward to about 0.3 au. The most common survivor is planet d, which is also the most massive of the four planets. In two cases the initially outermost planet, planet b, is perturbed to within about 10 au of the host star and survives.
We also determine if integrator accuracy has a demonstrable effect on instability times in Fig. 10. The plot does not show any apparent correlation. The vast majority of instabilities occur within a cooling age (or time since the white dwarf was born) of 10 Myr. The ejection timescale is defined specifically as when the planet leaves the Hill ellipsoid of HR 8799 with respect to the Galactic centre (Veras & Evans 2013;Veras et al. 2014a) and is dependent on the direction of escape. Reaching this boundary typically takes on the order of Myrs, meaning that the gravitational scattering event for points on the plot could have occurred at or just after the tip of the AGB phase.
This plot also illustrates that the planets in HR 8799 do not typically experience late gravitational scattering events, i.e. for cooling ages longer than 1 Gyr. Nevertheless, as demonstrated in Fig. 8, the nonuniform orbital changes represent a continuous source of dynamical activity despite the lack of late ejections. Another tentative trend in the plot is the vertical scatter in the data points as a function of progenitor stellar mass: as this mass increases, so does the violence during the giant branch phases. More rapid and greater mass loss for the 1.67M⊙ case produces stronger scattering events, allowing planets to be more easily perturbed beyond 10 3 au. At these distances, Galactic tides and flybys are more likely to instigate instabilities.

DISCUSSION
Because our investigation considered only one set of initial conditions for the planets, our results do not necessarily represent the future evolution of the system. Further, even with this set of initial conditions, there is a small but finite probability that a very close stellar flyby could disrupt the system along the main sequence.
In fact, the striking variety of fates of HR 8799 from a single set of initial orbital parameters emphasizes the significant role external forces will play in the system's future evolution. Galactic tides could lead to qualitative secular behavioural changes for planets with post-scattered semimajor axes 10 3 au depending on the mutual inclination between the orbits of those planets and the Galactic plane (see e.g. Fig. 5 to see planet c's gradually increasing pericentre). Stellar flybys, however, affect systems differently, randomly generating short bursts of strong interactions, but only sometimes. Unfortunately, prospects for computing the future stellar dynamical properties of HR 8799's immediate neighbourhood are limited because of the uncertainties in the stellar dynamics of Gyr-scale Galactic evolution.
Nevertheless, the uncertainties about the minor planets in the system and their eventual locations are arguably greater. These bodies, currently unseen except for the dust they produce, are the potential metal polluters of the white dwarf. They can secularly or resonantly interact with the planets during giant branch evolution (Dong et al. 2010) and experience grind-down and radiation blowout during this violent phase (Bonsor & Wyatt 2010;Zotos & Veras 2020). In addition, these bodies will be subject to the radiative Yarkovsky effect, which can propel them either outward or inward by tens or hundreds of au (Veras et al. , 2019a. Some of asteroids closest to the star may also spin up to the point of rotational fission from the YORP effect (Veras et al. 2014b;Veras & Scheeres 2020), producing debris which then becomes more easily propelled by the Yarkovsky effect.
Wherever the minor planets and debris ultimately end up residing, we have shown the surviving planets typically and continually sweep out large areas encompassing the entire system from ∼ 0.1 au to ∼ 10 5 au. Based on this result, HR 8799 represents a viable future metal-polluted white dwarf. The pollution is likely to occur soon after the star becomes a white dwarf because the planets are so massive. Terrestrial-mass planets would allow minor planets to meander for multiple Gyr (Frewen & Hansen 2014;Mustill et al. 2018), whereas giant-mass planets would scatter the vast majority of minor planets within 1 Gyr (Veras et al. 2021). We then might expect HR 8799 to be polluted at cooling ages under about 1 Gyr.
Although HR 8799 appears to contain all of the ingredients for the eventual white dwarf to be polluted, planetary systems as massive as HR 8799 are rare. Only a few systems are known to contain multiple massive giant planets, such as PDS 70 (Keppler et al. 2018;Haffert et al. 2019) and TYC 8998-760-1 (Bohn et al. 2020a,b). Hence, although HR 8799 would likely represent a standard future polluted white dwarf, its surviving set of planets would not be typical.
In our simulations, only rarely did the planets themselves enter the white dwarf Roche radius, or reach close enough to it to be significantly influenced by starplanet tidal effects (Veras & Fuller 2019;Veras et al. 2019b;O'Connor & Lai 2020). Amongst the known population of four-planet systems (Maldonado et al. 2021), the fraction with close encounters between the white dwarf and planet is higher most likely because on average, those systems harbour smaller masses. Both  and Veras et al. (2016) demonstrated that in systems with at least four planets, the smaller the planet masses, the longer the planets will meander and the greater the chance that they will enter the white dwarf Roche sphere.
Our analysis of the HR 8799 system relies on a series of assumptions which we now discuss in more detail. The first is that the debris discs are not massive enough to alter the future stability of the four planets. Moore & Quillen (2013) investigated this situation in detail (but not with the Goździewski & Migaszewski 2020 orbital fit) and found that a sufficiently massive exterior disc may act as either a stabilizing or de-stabilizing influence on the planets. For discs to be effectual either way, they need to currently be at least as massive as Neptune. However, the disc masses in HR 8799 system remain unknown. The disc masses previously could have been larger, especially during and just after the protoplanetary disc phase.
Another one of our assumptions was that the parent star emits mass isotropically and does not experience a natal kick during the transition to a white dwarf. Anisotropic mass loss would change the equations of motion for the variable-mass two-body problem (Veras et al. 2013b;Dosopoulou & Kalogera 2016a,b), but not qualitatively change the primary trigger for the instability, which is a shifting of stability boundaries due to mass lost (Debes & Sigurdsson 2002;Veras et al. 2013a). Further, Fig.  2 indicates that the instability is triggered at the end of the AGB phase, which is also when the greatest anisotropy in stellar wind flow would occur. Overall then, we would not usually expect anisotropy to alter the instability timescale signficantly, but rather (if present) to change the fate of the system stochastically at a similar level to how weak random stellar flybys do so.
A third assumption is that no additional planets reside in the system. Hence, we now consider the prospects that as-yet unseen planets in HR 8799 exist. Following the initial discovery of the first three planets, Hinkley et al. (2011) used interferometric imaging to search for any massive inner companions residing within 10 au of the system. Although no new companions were found, these observations placed some of the strongest limits on companions in the innermost regions of the system, including the region within about 6 au that is relatively cleared of dust.
However, as more advanced coronagraphs have been deployed (e.g. Maire et al. 2015), and sophisticated image processing techniques have matured, more recent observations (e.g Currie et al. 2014;Wahhaj et al. 2021) have placed comparable limits on an inner fifth companion in the system. Similarly, the morphology of the outer circumstellar debris structure is perhaps the best indicator of where additional planets may reside at wide separations. Given that modern ground based observations easily achieve sensitivities in the HR 8799 system down to a few Jupiter masses within about 200 au, the masses of any additional companions in the system will likely be 1-2 MJup. Such low masses will be easily attainable given the sensitivity of the upcoming James Webb Space Telescope (Carter et al. 2021).

SUMMARY
We have modelled the fate of the HR 8799 planetary system under the assumption that all four planets will survive until the end of the main sequence in a resonant configuration. Because of the stochastic nature of the system, markedly different outcomes are generated by just tweaking the integrator accuracy or introducing external effects; stellar flybys and Galactic tides can play a significant role in the evolution during the white dwarf phase after one of the planets has been scattered out to a distance of hundreds or thousands of au. Because HR 8799 is an A-type star, it represents the typical progenitor of the white dwarfs currently observed in the Galaxy. With the star's four planets, debris belts, and dynamical activity during the white dwarf phase, HR 8799 contains all of the necessary ingredients to become a debrispolluted white dwarf.