Disentangling the parameter space: The role of planet multiplicity in triggering dynamical instabilities on planetary systems around white dwarfs

In this work, we test four different planet masses: 1, 10, 100, and 1000 $\mathrm{M_{\oplus}}$. All planets have the same mass in each multiple-planet template. When the radius is needed we calculate it using the mass–radius relation proposed by Chen & Kipping (2017), as in \al@maldonado2020,maldonado2020b,maldonado2021; \al@maldonado2020,maldonado2020b,maldonado2021; \al@maldonado2020,maldonado2020b,maldonado2021. To give an overview on how the parameters used in this work compare with published studies, in the upper panel of Fig. 1 we show the planet mass as a function of the number of planets, where the green squares represent the planet masses tested in this work. The values adopted in previous works in the literature are shown with different symbols (see the legend and the references given in the figure’s caption). The box plots correspond to the mass distributions used in \al@maldonado2020,maldonado2020b,maldonado2021; \al@maldonado2020,maldonado2020b,maldonado2021; \al@maldonado2020,maldonado2020b,maldonado2021. The medians of each box are displayed as gray lines inside the boxes, the upper whisker extents from the third quartile Q3 to the lower datum below Q3+1.5(Q3-Q1), with Q1 the first quartile and the lower whisker goes from Q1 to the lowest datum above Q1-1.5(Q3-Q1) in each box. The outliers are shown with the small gray dots. The gray triangles in the upper part of each box distribution depict that there are mass values larger than 5000 $\mathrm{M_{\oplus}}$ in the samples.

In our simulations, planets are separated keeping equal distances (measured in mutual Hill radii) between adjacent pairs. The planet separation in terms of mutual Hill radius units is defined as $\Delta=(a_{j}-a_{i})/\mathrm{R_{m,Hill}}$, where $a_{i}$, $a_{j}$ are the semimajor axes in planets $i$, $j$, with $a_{j}>a_{i}$, and the mutual Hill radius is defined as:

where $m_{i}$, $m_{j}$ are the planet masses and $M_{*}$ is the mass of the central star. In this work we test 29 $\Delta$ values, ranging from 2 to 30 in steps of 1 mutual Hill radius. In the middle panel of Fig. 1, we display the $\Delta$ values explored (the green squares) as a function of the number of planets for the smallest planet mass. $\Delta$ values used in previous works are shown as grey symbols. References are provided in the figure caption. As in the upper panel, the box plot show the $\Delta$ distributions used by us in \al@maldonado2020,maldonado2020b,maldonado2021; \al@maldonado2020,maldonado2020b,maldonado2021; \al@maldonado2020,maldonado2020b,maldonado2021.

In each planetary system, we place the innermost planet at $a_{0}$= 10 au. This distance is chosen in order to avoid the influence of tidal forces on the planets and to ensure the survival of the planet to engulfment on the AGB phase of evolution of the host star (Villaver & Livio, 2009; Mustill & Villaver, 2012; Villaver et al., 2014). The consecutive planets are placed at distances where $\Delta$ is kept the same between the adjacent planet pairs in the range between 2 and 30 mutual Hill radii. Nevertheless, in systems with four, five and six planets with planets of 1000 $\mathrm{M_{\oplus}}$, the outermost planet for some of the largest $\Delta$ values can reach a larger distance than the semimajor axis at which we have set up MERCURY to consider a planet ejection. For this reason, in the following and for the statistics the maximum $\Delta$ that we consider in the four-, five- and six-planet systems with 1000 $\mathrm{M_{\oplus}}$ planets is 21, 19 and 17, respectively.

Furthermore, we would like to mention that there is a maximum planet separation $\Delta_{\mathrm{max}}$, since the semimajor axis $a_{j}$ $\rightarrow$ $\infty$ as $\Delta_{\mathrm{max}}\rightarrow 2(\frac{m_{i}+m_{j}}{3M_{*}})^{-1/3}$ (cf. Mustill et al., 2014). For planetary systems involving 1, 10, 100 $\mathrm{M_{\oplus}}$, testing $\Delta\leq 30$ is completely feasible. However, for 1000 $\mathrm{M_{\oplus}}$ planets, $\Delta\to\Delta_{\mathrm{max}}=22.88$, which means that for these simulations, the maximum $\Delta$ we can test is 22.

Besides the mass and planet separation $\Delta$, eccentricities, and orbital inclinations are also required for the initial set-up of the simulations. Moorhead et al. (2011); Van Eylen & Albrecht (2015); Pu & Wu (2015) have demonstrated that the eccentricity and inclination in observed multiple-planet systems are well described by Rayleigh distributions. Thus, in this work, we select randomly for each planet the eccentricity and inclination from a Rayleigh distribution with parameters $\sigma=0.03$ and $\sigma=1.12^{\circ}$ (Xie et al., 2016). In the lower panel of Fig. 1, we show the eccentricity distributions used in this work as green box plots. The outliers correspond to eccentricities in the tail of the Rayleigh distribution (values larger than 0.15 were not generated in the random selection from the Rayleigh distribution). For comparison, eccentricity values from previous studies are also displayed with gray symbols (references are given in the figure’s caption). We can see that previous two-planet studies have performed simulations with several eccentricity values, but works that simulated planetary systems with three and more planets have mainly tested the circular case. The only exception is our \al@maldonado2020,maldonado2020b,maldonado2021; \al@maldonado2020,maldonado2020b,maldonado2021; \al@maldonado2020,maldonado2020b,maldonado2021 where we simulated a distribution of eccentricities based on parameters obtained from observations of planetary systems (the white box symbols). We see a slight difference between the white and green box plots: it is due to the fact that the eccentricities in the white box plots were done using a $\sigma=0.02$ (Pu & Wu, 2015) for transiting planets and in this work both eccentricity and inclination $\sigma$ parameters are taken from the mean values derived from the $\it{Kepler}$ multiple-planet sample analyzed by Xie et al. (2016).

We have randomly selected the angular phases in the planet’s orbit, i.e. argument of the pericentre, mean anomaly and ascending node, from uniform distributions from 0 to 360^∘.

There is a fraction of simulations where orbit crossing and/or planet losses happen during the MS phase of the host star. In numbers, we find that 172/1080 (15.9$\%$) in the two-planet, 344/1080 (31.9$\%$) in the three-planet, 367/1070 (34.3$\%$) in the four-planet, 394/1050 (37.5$\%$) in the five-planet and 416/1030 (40.4$\%$) in the six-planet simulations are unstable on the MS. The early instabilities occur because a fraction of planet pairs in these systems are separated by distances where mean motion resonances (MMR) overlap. Two planets are in MMR when the ratio of orbital periods between the planets is expressed as the ratio of two integers $(p+q)/p$, with $p$ the integer and $q$ the order of the resonance.

We estimate the planet separation of the chaotic zone by using $(a_{j}-a_{i})/a_{i}=1.8e^{1/5}(m_{p}/M_{*})^{1/5}$, with $m_{p}=m_{i}=m_{j}$ the planet mass and $e$ the planet eccentricity (c.f. Mustill & Wyatt 2012). Despite the fact that the latter equation predicts better the stability of systems with two equal-mass planets, it give us clues where systems with a higher number of planets are prone to have early instabilities when the adjacent planet pairs are separated less than such limiting distance. By using the eccentricity values randomly selected from the Rayleigh distribution in this work, the average of the limiting distance $\Delta$ in terms of mutual Hill radius for 1, 10, 100 and 1000 $\mathrm{M_{\oplus}}$ is respectively 6.5, 4.7, 3.4, 2.4, implying that simulations with $\Delta$ lower than the limiting distances previously calculated imminently lose planets or have orbit crossing on the MS. The rest of unstable systems in the MS have adjacent pairs of planets whose semimajor axis ratio is found within the 5$\%$ from the corresponding semimajor axis ratio in some first and second order MMR, mainly the 3:1, 2:1, 5:3, 3:2, 4:3, 5:4, 6:5, influencing the dynamics of the system (with the exception of 12 simulations in systems of four, five and six planets with 1000 $\mathrm{M_{\oplus}}$ having adjacent planet pairs within the 8 $\%$ of the 2:1 resonance).

Due to the fact that the resonant systems need a specific orbital configuration in order to be dynamically stable and that the finding of such configuration is beyond the scope of this work, we removed the systems in MMR with unstable outcomes in the MS for the final statistics. However, we do consider for the statistics the resonant systems in which the random selection of the orbital angles prevents the early instability on the MS.

The adiabatic approximation is defined as the condition fulfilled when the stellar mass-loss time-scale is much larger than the orbital time-scale of a planet. In this case, the semimajor axis evolution is given by

with $\mu$ the planet–star mass ratio. In this regime, the eccentricity remains almost unchanged.

The effects of adiabatic and non-adiabatic mass-loss on the dynamical evolution of single planets have been explored by e.g. Veras et al. 2011. In general, they find that planets located at orbital distances $\geq$ 1000 au with different eccentricities are likely to suffer an ejection from the system during the mass-loss of the star: the planet’s eccentricity will increase drastically, and as a consequence the orbit of the planet becomes hyperbolic and the planet is imminently ejected during or just after the major mass-loss event.

The non-adiabatic regime produces planet ejections due to the excitement of the planet’s eccentricity and we consider that it is important to know when the non-adiabatic effects become important in our simulations. Therefore, we perform additional simulations of single planets orbiting a 3$\mathrm{M_{\odot}}$ and evolving the systems for 500 Myr (from the MS to a few Myr after the beginning of the WD phase). We test one-planet systems with the four planet masses studied in this work (1, 10, 100 and 1000 $\mathrm{M_{\oplus}}$). The eccentricity, inclination and the phase angles in the planet’s orbit are randomly selected from the Rayleigh and uniform distributions described in Section 2 for each template.

In the first trial, we test different initial semimajor axis ($a_{\mathrm{initial}}$) values, from 50 to 1000 au in steps of 50 au and we find that the eccentricity increases as the planet is located at further distances from the host star; we do not find any difference in the eccentricity change among the simulations with different planet masses. We find that around 250 au the final eccentricity ($e_{\mathrm{final}}$) increases to $\sim$ 0.2.

In the second trial, we make additional simulations in order to refine $a_{\mathrm{initial}}$ from 50 to 350 au, with steps of 10 au. This time, we test the four planet masses with the orbital parameters randomly selected in each run and we make 10 runs per planet mass and per $a_{\mathrm{initial}}$ value. In Fig.3, we show eccentricity as a function of the initial semimajor axis in each one-planet simulation. Due to the fact that we do not find any difference among the eccentricities for different planet masses, we display in the blue dots the final eccentricities in each system taken after the formation of the WD (478 Myr). For comparison, we also show with the red dots the initial eccentricities randomly selected from the Rayleigh distribution in each simulation. The horizontal black dashed line depicts $e_{\mathrm{final}}=0.2$ for better visualization along the x-axis. We see that the eccentricity is increasing as the planet is located further from the central star and for $a_{\mathrm{initial}}<200$ au, $e_{\mathrm{final}}<0.2$. Moreover, we see that the initial and final eccentricities converge as $a_{\mathrm{initial}}$ goes to 50 au, which it is expected in the adiabatic regime.

In this work, the planetary templates with large $\Delta$ values have planets at orbital distances where non-adiabatic ejections may affect the dynamical evolution. Indeed, we find for planetary systems with 100 $\mathrm{M_{\oplus}}$ planets, 27 simulations in the five-planet case ($\Delta\geq$ 28) and 70 simulations in the six-planet case ($\Delta\geq$ 24) where planets are ejected or have hyperbolic orbits when the star becomes a WD (at 478 Myr). Moreover, in systems with 1000 $\mathrm{M_{\oplus}}$ planets, we have 30 simulations in the three-planet case ($\Delta\geq$ 20), 58 simulations in the four-planet case ($\Delta\geq$ 16), 67 simulations in the five-planet case ($\Delta\geq$ 13) and 68 simulations in the six-planet case ($\Delta\geq$ 11) where non-adiabatic ejections happen. Since the ejections in all these simulations are produced by non-adiabatic effects and not by dynamical interactions among the planets, we decide to remove these simulations for the final statistics.

We take into account simulations dynamically active on the WD phase (having planet losses, orbit crossing and orbital scattering in at least one planet, as the examples depicted in Fig. 2). As described in sections 3.1 and 3.2, for the final statistics of unstable systems we remove simulations where planet losses and orbit crossing happens on the MS. We do not consider simulations where at least one planet is ejected by non-adiabatic mass-loss.

Due to the fact that many observed protoplanetary disks have extensions $\lesssim$ 200 au (e.g. Tazzari et al. 2017; Garufi et al. 2020; Andrews 2020) and based on the previous results from the one-planet tests performed in section 3.2, we check in our multiple-planet templates with the highest number of planets which planet separation $\Delta$ places all planets within 200 au in each system. For six-planet systems with 1 and 10 $\mathrm{M_{\oplus}}$, the outermost planet is always located at $a<200$ au with $\Delta\leq 30$. However, for 100 and 1000 $\mathrm{M_{\oplus}}$, the sixth planet fulfils the condition only with $\Delta\leq 14$ and $6$, respectively. We choose $\Delta=14$ as one of the threshold limits in the planet separation space in order to analyze our results. We make this choice because for $\Delta=6$, we would have very few or no simulations with dynamical instabilities on the WD phase to count after removing the MS unstable simulations for the statistics.

We consider $\Delta=$ 17 and 30 as additional $\Delta$ thresholds in order to analyze the fraction of unstable simulations. The former considers all the planet masses tested in this work with the largest planet separation simulated ($\Delta$=17 in the six-planet systems with 1000 $\mathrm{M_{\oplus}}$ planets). The latter considers only planetary systems with 1, 10 and 100 $\mathrm{M_{\oplus}}$ planets since $\Delta=30$ is the largest planet separation considered in this study and is only tested for those planet masses ($\Delta_{\mathrm{max}}=22.88$ for 1000 $\mathrm{M_{\oplus}}$ planets; see section 2.2).

We show in Fig. 4, 5 and 6 several pie charts with the fraction of stable and unstable simulations considered for the statistics ordered in rows (each one per planet mass) and columns (each one referring to the number of planets in the system). The different colours depict different dynamical outcomes. In green, we mark simulations that do not lose any planet in the 10 Gyr but show one or more orbit crossings of the planets (e.g. see as an example of such behaviour the middle panel in Fig. 2). In brown, we show simulations where at least one planet have orbital scattering in the semimajor axis without orbit crossing and no planets are lost in the entire simulated time (e.g. the upper panel in Fig. 2). Regarding planet losses, either by Hill or Lagrange instabilities, there are cases where multiple planets are lost in an individual simulation (e.g. the lower panel in Fig. 2), specially in systems with five and six planets. Therefore, the pink, yellow, blue, orange and light-blue colours refer to simulations that loss 1, 2, 3, 4, 5, planets respectively. Lastly, the gray colour marks completely stable systems that only show the adiabatic expansion of the planets’ semimajor axes due to the stellar mass loss when the star becomes a WD.

We see different dynamical behaviours regarding the planet mass, the planet separation, and the number of planets considered in the systems. In Fig. 4, for $\Delta\leq$ 14, we first see that there are no pie charts for the five- and six-planet systems on the first row. This is due to the fact that all the simulations in these cases experience dynamical instabilities on the MS. Regarding orbit crossing without planet losses, we observe that this instability is happening in all but restricted to simulations with 1 and 10 $\mathrm{M_{\oplus}}$ planets (low-mass planets). Orbital scattering without planet losses is present in a fraction $\sim$ 20 $\%$ of the two-planet systems with low-mass planets. A much lesser fraction is obtained in the five-planet systems with 100 and 1000 $\mathrm{M_{\oplus}}$ (high-mass planets), in which the outermost planet shows the scattering. Moreover, simulations losing one planet occur in all but 4 cases, where the highest fractions are obtained in five- and six-planet systems with 10 $\mathrm{M_{\oplus}}$. Planetary systems losing two and more planets are restricted to templates with four, five and six planets with 10, 100 and 1000 $\mathrm{M_{\oplus}}$. We do not have any simulation that loses all the planets during the 10 Gyr. The number of stable simulations decreases in all the cases as the number of planets increases.

Regarding simulations with $\Delta\leq$ 17, in Fig. 5 we have all the rows and columns with pie charts. In general, we observe a similar dynamical behaviour with respect to Fig. 4. The fraction of stable simulations decreases as the number of planets increases. Simulations with orbit crossing without planet losses are only present in systems with low-mass planets and we see an increasing trend in systems with 1 $\mathrm{M_{\oplus}}$ planets, as the number of planets increases. Besides having the same cases with simulations showing only orbital scattering without planet losses as $\Delta\leq$ 14, we have a new small percentage of simulations with orbital scattering and no planet losses appearing in the four-planet systems with planets of 1000 $\mathrm{M_{\oplus}}$, where the fourth planet has the scattering. We note that the fraction of stable and unstable simulations in systems with five and six planets having 1000 $\mathrm{M_{\oplus}}$ remains constant. It happens because non-adiabatic mass-loss ejections happen in these simulations with $\Delta$ $\geq$ 13 and 11, respectively, they are removed from the statistics and no extra simulations are counted for $\Delta\leq 17$ (see the analysis of section 3.2 and Fig. 4).

Finally, in Fig. 6, we display only planetary systems with planets having 1, 10 and 100 $\mathrm{M_{\oplus}}$ in which we perform simulations with $\Delta\leq$ 30. We see that in all the cases the fraction of stable simulations is in the range between 50 $\%$ and 93 $\%$. We confirm that the fraction of simulations with orbit crossing and no planet losses increases as the number of planets increases in simulations with planets of 1 $\mathrm{M_{\oplus}}$. It does not happen for planetary systems with 10 $\mathrm{M_{\oplus}}$ planets, instead, the fraction of planet losses increases, while the overall fraction of unstable systems increases only slightly. In addition to these low-mass cases, we have another fraction of simulations with orbit crossing and no planet losses in the six-planet systems with 100 $\mathrm{M_{\oplus}}$ planets. Regarding simulations with only orbital scattering without planet losses, we obtain the same cases as described in $\Delta\leq$ 14, and 17; however, for the simulations with high-mass planets, there is an increase in the fraction of five-planet systems experiencing weak scattering, and the six-planet systems have only a small fraction. All the planetary systems lose at least 1 planet and for systems with four and more planets, the higher the planet masses are, the more planets in a single planetary configuration get lost. The six-planet systems with planet mass of 100 $\mathrm{M_{\oplus}}$ have simulations that lose from 1 to 5 planets in a single template.

As mentioned in section 3.2, the planet’s eccentricities are excited due to non-adiabatic effects and the further the planet is located from the central star, the greater the eccentricity excitation induced during the mass-loss of the stellar host, even if non-adiabatic ejections do not happen. For the case of six-planet simulations with 100 $\mathrm{M_{\oplus}}$ planets, we obtain that non-adiabatic ejections happen for $\Delta\geq$ 24. However, for $\Delta=$ 22 and 23, the sixth planet’s eccentricity ranges from 0.72 to 0.93 and the fifth planet’s eccentricity is in the range $0.29<e<0.44$. As a consequence, orbit crossing between the sixth and fifth planets is observed in the 7.9 $\%$ of simulations without planet losses observed in Fig. 6 and starts from the beginning of the WD phase. A similar eccentricity evolution is observed in the five-planet systems with 100 $\mathrm{M_{\oplus}}$ for $26\leq\Delta\leq 28$, where the fifth planet has $e\geq$ 0.75, not enough to produce orbit crossing with the fourth planet ($0.17<e<0.35$); however, the perturbations induced by the other planets to the fifth planet produce the 8.8 $\%$ of simulations having orbital scattering without losing the planet observed in the five-planet column in Fig. 6. Having dynamical instabilities in the outer planets of a planetary system could perturb a putative outer asteroid belt and some of its components may be launched towards the inner planetary system, where the inner planets might continue sending the material toward the WD, as suggested in Bonsor et al. (2011); Mustill et al. (2018); Smallwood et al. (2018).

By considering the total fraction of unstable simulations (summing the simulations that lose planets, have orbit crossing and orbital scattering), we present in Fig. 7 three panels with the total fraction of simulations dynamically active on the WD phase with respect to the total number of simulations considered for the statistics as a function of the number of planets in the systems. We refer to systems involving a specific planet mass by different colours (brown, green, blue and yellow for 1, 10, 100, 1000 $\mathrm{M_{\oplus}}$ planets, respectively). Error bars are calculated assuming a binomial sampling distribution (Jaynes & Bretthorst, 2003) and are equivalent to 1$\sigma$ in Gaussian distributions. Each panel from left-hand to right-hand displays the fraction of simulations calculated within the planet separation $\Delta\leq$ 14, 17 and 30. Globally, we can see that planetary systems with 1 $\mathrm{M_{\oplus}}$ are more dynamically active on the WD phase than the systems with higher mass planets, with the exception of systems with 1000 $\mathrm{M_{\oplus}}$ planets, where the five- and six-planet systems with $\Delta\leq$ 14 and 17 have an increase of unstable simulations reaching fractions higher than 70 $\%$. This increment is due to the fact that simulations with 1000 $\mathrm{M_{\oplus}}$ planets are dynamically active on the WD phase with $\Delta$ values smaller than lower-mass planet simulations and non-adiabatic ejections start to appear at $\Delta$ values smaller than for the other planet masses (the $\Delta$ ranges considered for the statistics without the MS unstable simulations and after removing the non-adiabatic ejections are 6–13 and 6–11 for the five- and six-planet systems, respectively). The same happens to a lesser extent in five- and six-planet simulations having 100 $\mathrm{M_{\oplus}}$ planets for $\Delta\leq$ 30 ($\Delta$ values considered for the statistics: 8–28 and 8–23 for systems with five and six planets, respectively). Finally, we see a general tendency to have a larger fraction of unstable simulations as the number of planets increases in the systems, independently of the planet mass and the limits in $\Delta$.

After adding all the unstable simulations in the four planet masses but keeping separated the $\Delta$ thresholds, in Fig. 8 we show the total fraction of simulations dynamically active as a function of the number of planets in the system. Each colour refers to the $\Delta$ limits according to the figure’s legend. The error bars are calculated with the binomial sampling distribution (Jaynes & Bretthorst, 2003). We observe that higher fractions of unstable systems are obtained for simulations with $\Delta\leq$ 14 than systems with $\Delta\leq$ 17 and 30, which is an expected result since planets separated by large mutual Hill radius are less likely to interact each other and remain completely ordered the 10 Gyr of evolution. On the other hand, we can see the clear tendency of having more dynamically active simulations as the multiplicity of planets increases for the three $\Delta$ thresholds. The fraction goes from $10.2^{+1.2}_{-1.0}$ $\%$, $19^{+1.9}_{-1.7}$ $\%$ and $25.2^{+2.5}_{-2.2}$ $\%$ in the two-planet systems to $33.6^{+2.3}_{-2.2}$ $\%$, $60.8^{+3.5}_{-3.8}$ $\%$ and $74.1^{+3.7}_{-4.6}$ $\%$ in the six-planet systems, for $\Delta\leq$ 30, 17 and 14, respectively.

We find that the fraction of dynamically active simulations on the WD phase becomes larger as the number of planets within the planetary system increases. This result has been probed in a sample of multiple-planet systems built ad-hoc so the $\Delta$ thresholds and planetary masses and eccentricities are the same among the multiple systems tested. This result confirms the main conclusion from Paper III, in which scaled versions of observed multiple-planet systems were simulated but in which multiplicity was just one of the parameters that could be hidden among a mix of planet masses, eccentricities and planet separations of the simulated systems. We find that the fraction of instabilities increases with the multiplicity of the system independently of the considered planet separation $\Delta$ thresholds (see Fig. 8) and it is consistent for different planetary masses. The driving force of the increased instability at higher multiplicity is probably mainly a consequence of the intrinsically lower stability of high-multiplicity systems (Chambers et al., 1996; Funk et al., 2010), which is found both in previous works without mass loss and in our own simulations when we look at stability along the MS before mass loss. We speculate that this is a result of a larger number of multiple-planet resonances. However, non-adiabatic expansion may play a role for the systems of more massive planets where orbital radii can be large.

We also find that a different dynamical behaviour is observed on the WD phase depending on the planet mass. Here, we refer to planet masses of 1 and 10 $\mathrm{M_{\oplus}}$ as low-mass planets and planet masses of 100 and 1000 $\mathrm{M_{\oplus}}$ as high-mass planets. Simulations having orbit crossing without planet losses in the 10 Gyr tested are mainly restricted to systems involving low-mass planets. In addition, low-mass planets experience close encounters which are not strong enough to trigger a planet loss. The same dynamical outcome is found by Veras et al. (2016); Mustill et al. (2018) in their simulations of terrestrial planet mass. Indeed, the planet losses in our 1 $\mathrm{M_{\oplus}}$ simulations are restricted only to planet–planet collisions in all the multiple-planet configurations (see as well e.g. Veras et al. (2013); Veras & Gänsicke (2015)). In addition, planet–star collisions happen more frequently than planet–planet collisions but much less frequently than ejections, and the higher the mass of the planets, the smaller the number of collisions with the star obtained.

On the other hand, simulations involving high-mass planets preferably lose planets and typically more than one planet is lost in a single simulation in the five- and six-planet systems (see Figs. 4, 5, 6). We point out that high-mass systems involve close encounters among the planets in which the orbit crossing strongly induces a planet loss afterwards, mainly by ejection. In our simulations, the largest fractions of ejections occur in the systems with high-mass planets being the six-planet systems with 1000 $\mathrm{M_{\oplus}}$ planets where $\sim$ 60 $\%$ of the planets are ejected for $\Delta\leq$ 14, 17. This tendency of a higher rate of ejections for higher mass planets has been previously reported by simulations in the literature (see e.g. Veras & Mustill 2013; Veras et al. 2013; Mustill et al. 2014; Mustill et al. 2018; Veras & Gänsicke 2015; Veras et al. 2016; Paper I; Paper II; Paper III) and it is well documented analytically by the Safronov number $\Theta$ c.f. Safronov & Zvjagina 1969; Mustill et al. 2014) that states that it is more likely that a planet scatter out other bodies as the rate of the escape velocity with the orbital velocity increases. The Safronov number using the mass of the WD and the innermost planet’s orbit after the adiabatic expansion due to the mass-loss, is 3.8, 12.1, 30.9, 290.0 for the four planet masses (1, 10, 100 and 1000 $\mathrm{M_{\oplus}}$) respectively. Thus, our numerical results confirm the expected scatter out of more massive bodies as the planet/star mass ratio increases.

Furthermore, we find that simulations with high-mass planets tend to have the first orbit crossing at a WD cooling time $<$ 100 Myr, with the planet losses starting a few Myr later. The behaviour is slightly different for low-mass planets that have the first orbit crossing at cooling times $>$ 100 Myr and the first planet losses happen a few Gyr afterwards¹¹1We have calculated the geometric mean of four samples having the instability times converted to the WD’s cooling time, two of them (1 and 10 $\mathrm{M_{\oplus}}$ simulations together and the other with 100 and 1000 $\mathrm{M_{\oplus}}$ together) with the information of the first orbit crossing in all the multiple-planet simulations dynamically active on the WD phase. The other two samples have the time when a planet is lost or when the first planet loss happens in case that multiple planets are lost in the same simulation. The geometric means for orbit crossing: high-mass – 68 Myr, low-mass – 351 Myr; planet losses: high-mass – 331 Myr, low-mass – 2869 Myr.. Previous studies have described a similar tendency on the instability times between low-mass and high-mass planets (e.g. Veras & Mustill 2013; Mustill et al. 2014; Paper II; Paper III) and Veras et al. (2016); Mustill et al. (2018) also concluded that planets from 1 to 30 $\mathrm{M_{\oplus}}$ kept being dynamically active for Gyr time-scales. We should note that although we have stated that at high planet masses, one or more planets are rapidly ejected, this does not necessarily lead to long term WD pollution since simulations that include asteroids show that the resulting accretion rate decays much more quickly (e.g. Mustill et al. 2018), and they will produce a burst of accretion but little long-term delivery.

It can be clearly seen in Figs. 7 and 8 that the fraction of unstable simulations on the WD phase decreases as we consider systems with a higher planet separation $\Delta$, i.e. the fraction of dynamically active simulations is higher in systems with $\Delta\leq$ 14 than in systems with $\Delta\leq$ 17 and 30, respectively. This is an expected result given that it is the planet separation that provides the Hill and Lagrange stability limits at which planets are separated enough to have a dynamical evolution completely ordered and stable. This behaviour is also seen by comparing the fraction of stable systems obtained in our simulations in Figs. 4, 5, 6, which increases as $\Delta$ also increases. In general, we find that the required planet separation $\Delta$ to have fully stable systems during the entire simulated time increases when the system has more planets, in agreement with (Funk et al., 2010). The stability limits in our multiple-planet simulations are in agreement with the limits found in previous studies with two planets (Veras & Mustill, 2013; Veras et al., 2013), three planets (Mustill et al., 2014; Mustill et al., 2018) and four planets (Veras et al., 2016) by comparing the equivalent planetary masses and architectures. Although we can certainly know when a pair of planets are prone to have close encounters by the Hill stability formalism (Hill, 1878; Donnison, 2011), there is not any analytical formulation to predict when a planetary system might be Lagrange unstable and it can only be known by performing numerical simulations.

We highlight that orbital distances scaled close to planet mass to the 1/4th power should better suit the stability criterion in multiple-planet systems than planet distances scaled to planet masses to the 1/3rd power, as the Hill radius unit (see e.g. Petit et al. 2020; Tamayo et al. 2021; Rath et al. 2021). Therefore, in order to see if there is a significant difference in the results adopting one or the other scaling, we compare the fractions of unstable simulations on the WD phase that we have found within the threshold distances in mutual Hill radius with respect to their corresponding fractions found within their equivalent thresholds that scale to 1/4th power of the planet mass ²²2We have transformed our planet distances $\Delta$ in mutual Hill radius units to dynamical distances defined by equation 6 in Tamayo et al. (2021).. We find small, but not significant, differences in the number of unstable simulations when comparing both scalings in their respective threshold distances; the general tendency of increasing fraction of unstable systems with the number of planets is kept. Perhaps with more simulations a significant difference might be found when comparing both scalings; however, it is hidden within the error bars in our sample.

We note that we have chosen the eccentricity from the same Rayleigh distribution with $\sigma$ = 0.03 for all the planet masses. As a consequence, as the planet mass increases, so does the planet separation, thus, the eccentricity will become a fraction of the orbit crossing eccentricity. A caveat of this choice is that systems involving lower mass planets are, a priori, more likely to be dynamically unstable than systems having higher mass planets with the same planet separation. We defer to a future work to build a more dynamically comparable sample in terms of eccentricity, that is to re-scale the eccentricity with respect to the 10 $\mathrm{M_{\oplus}}$ simulations to be $\sigma_{e}$ $\sim$ 0.03.

It is important to identify the link between the results obtained in this study and what is expected to be observed in real multiple-planet systems. However, we can not make a direct comparison because the real number of systems with multiple planets is hard to know (Johansen et al., 2012; Tremaine & Dong, 2012; Zhu et al., 2018; Zink et al., 2019), and it is likely that the observed planetary systems with one or two planets have additional planets but they remain hidden due to the planetary architectures (for example, high mutual inclinations: Lissauer et al., 2011; Fang & Margot, 2012; Fabrycky et al., 2014; Ballard & Johnson, 2016). It is important to note however that it has been estimated that the number of planets per G and K star is 5.86 $\pm$ 0.18 as a correction to the number of planets in the Kepler parameter space (Zink et al., 2019). An additional reason why we can not compare our results with observed multiple-planet systems is due to the fact that the majority of these systems have all their planets located at semimajor axes $<$ 0.5 au and have sizes between Super-Earth to Sub-Neptune (Pu & Wu, 2015; Weiss et al., 2018). Only very few observed systems have planets with orbital distances comparable to the ones studied in this work, such as the planetary system HR 8799 with four giant planets (Marois et al., 2008; Marois et al., 2010). Moreover, the observed multiple-planet systems have a mix of planet masses, eccentricities, planet separations (Schneider et al., 2011; Akeson et al., 2013) and the planetary multiplicity is entangled with this mix in the physical and orbital parameter space of real systems, as it is shown in \al@maldonado2020,maldonado2020b,maldonado2021; \al@maldonado2020,maldonado2020b,maldonado2021; \al@maldonado2020,maldonado2020b,maldonado2021.

Systems may also provide pollution to the WD even if they do not undergo instability, for example through secular chaos (O’Connor et al., 2021). The increase in orbital eccentricities we find going from the MS to the WD will enhance such secular chaos. The mechanism will also be more significant in high-multiplicity systems where secular resonances are more numerous. Thus, our approach of looking at direct instability among the planets gives a conservative estimate of the fraction of multiple-planet systems that may be responsible for WD pollution.

Finally, we point out that a hint about the trend of finding more dynamical instabilities with more planets in a planetary system was found in a more limited analysis by Veras & Gänsicke (2015); Veras et al. (2016) when comparing their four-planet unstable simulations with the few tests involving six-, eight- and ten-planet systems. Furthermore, the fractions of dynamically active simulations in this work, regarding the two-, three- and four-planet systems with $\Delta\leq$ 14 on the WD phase, are comparable to the fractions found by Veras & Mustill (2013); Veras et al. (2013); Mustill et al. (2014); Mustill et al. (2018); Veras & Gänsicke (2015), respectively.