Birth cluster simulations of planetary systems with multiple super-Earths: initial conditions for white dwarf pollution drivers

ABSTRACT

1 INTRODUCTION

Over 99 per cent of all known exoplanet host stars will eventually evolve into white dwarfs (WDs). This fact emphasizes the importance of being able to connect planetary architectures around WDs to their previous incarnations around giant branch and main-sequence stars and to the processes that occurred in their nebular birth clusters.

One way to pursue this connection is to consider the observations of known planetary systems around WDs. Veras (2021) partitions these observations into four classes: (i) major planets, (ii) minor planets (such as asteroids, comets, or moons, but also the remnants of larger planets), (iii) discs and rings, and (iv) chemical pollution by metals in the WD’s atmosphere from accreted planetary debris. The largest category is by far the last, which includes over 1000 systems (Dufour et al. 2007; Kleinman et al. 2013; Kepler et al. 2015, 2016, 2021; Coutu et al. 2019). The smallest category is the first, with just five examples of major planets known (Thorsett, Arzoumanian & Taylor 1993; Sigurdsson et al. 2003; Luhman, Burgasser & Bochanski 2011; Gänsicke et al. 2019; Vanderburg et al. 2020; Blackman et al. 2021). Nevertheless, both categories are connected, despite this observational gulf, because planets can dynamically drive this debris or their progenitor asteroids, moons, or comets into WDs.

The accreted debris is ubiquitous, appearing in between 25 and 50 per cent of Milky Way WDs (Zuckerman et al. 2003, 2010; Koester, Gänsicke & Farihi 2014). The debris has also been observed to occur at cooling ages (the time since becoming a WD) up to 8 Gyr (Hollands, Gänsicke & Koester 2018; Blouin & Xu 2022).

Identifying the planetary architectures that can allow planets to perturb smaller bodies towards the WD in a manner that mimics the distribution of accretion rate with cooling age is an ongoing challenge and is subject to a large number of degeneracies. The currently known exoplanets orbiting main-sequence stars will predominately be engulfed by their parent stars upon leaving the main sequence (Maldonado et al. 2020a,b, 2021), and sub-Saturn-sized planets remain largely hidden from view at separations where such planets would survive.

However, such distant super-Earths is a class of planets shown to be particularly efficient at polluting WDs at a wide variety of cooling ages (Frewen & Hansen 2014; Mustill et al. 2018), particularly because such planets meander with their orbits in continuous motion (Veras & Gänsicke 2015; Veras et al. 2016). In contrast, a Jupiter–mass planet acts more like a sledgehammer on a fixed orbit, dissipating the system and polluting the WD in short bursts (Veras et al. 2021). One key feature of these super-Earth polluters is that they are on eccentric orbits (as e.g. for β Pictoris; Beust & Morbidelli 1996), because a single planet on an exactly circular cannot perturb minor bodies sufficiently close towards a WD (Antoniadou & Veras 2016). Also helpful for pollution is when planets reside in multiple-planet systems, because otherwise perturbing minor bodies on to star-grazing orbits is challenging (Bonsor, Mustill & Wyatt 2011) and may require asteroid reservoirs several orders of magnitude more massive than the Solar systems (Debes, Walsh & Stark 2012). WD pollution due to secular chaos in multiplanet systems has been previously investigated in Smallwood et al. (2018), Smallwood et al. (2021), and O’Connor, Teyssandier & Lai (2021).

The vast majority of WD pollution investigations that contain perturbing planets (see fig. 6 of Veras 2021 for an extensive list) use initial conditions for their planetary systems, which are not outputs from birth cluster simulations. In a first attempt to bridge this gap, Veras et al. (2020) connected the outcomes of stellar cluster simulations involving outer Solar system analogues with their future evolution across different stellar phases. However, partly because their setup was limited to giant planets on nearly circular orbits that architecture is not necessarily representative of those found in chemically polluted WD systems.¹

Figure 1.

The simulated star cluster at the beginning of the simulation (⁠|$t=0\, \mathrm{Myr}$|⁠, |$r_\mathrm{vir}=1.0\, \mathrm{pc}$|⁠; left-hand panel) and at the time |$t=100\, \mathrm{Myr}$| (⁠|$r_\mathrm{vir}=3.9\, \mathrm{pc}$|⁠; right-hand panel). The star cluster has visibly expanded and is in the process of dissolving. Although the cluster’s centre of mass and its density centre move, which means that the plotted axes range must be adjusted with increasing simulation time, the physical scale (a total of |$12\, \mathrm{pc}$| on x- and y-axis) for both plots remains the same.

Open in new tabDownload slide

Figure 2.

Fractions of surviving planets for the 3P, 7PC, and 7PW model. The dotted black line represents the average survival fraction for each model.

Open in new tabDownload slide

Figure 3.

The a–e space for the 3P, 7PC, and 7PW model as well as for the different host star masses. The grey shaded area shows which planets would be engulfed by the star during the red giant phase.

Open in new tabDownload slide

Figure 4.

Cumulative, normalized histogram showing the distribution of eccentricities sorted by planetary model (bluish colours) and host star mass (reddish colours). The bin size is 0.02.

Open in new tabDownload slide

Figure 5.

The a–i space for the 3P, 7PC, and 7PW model. Planets above the dotted grey line are on a retrograde orbit.

Open in new tabDownload slide

Figure 6.

Cumulative, normalized histogram showing the distribution of inclinations sorted by planetary model (bluish colours) and host star mass (reddish colours). The bin size is 2.5°.

Open in new tabDownload slide

Here, we perform stellar cluster simulations with a wider variety of planetary architectures that are more likely to pollute the eventual WDs over long cooling times. We also use more representative progenitor WD masses (⁠|$1.5\, \mathrm{M}_{\odot }{-}2.5\, \mathrm{M}_{\odot }$|⁠; Tremblay et al. 2016; Cummings et al. 2018; McCleery et al. 2020) rather than |$1.0\, \mathrm{M}_{\odot }$| Sun-like stars. A key result of this study is publicly available sets of post-cluster initial conditions that modellers could use as starting points for their simulations of post-main-sequence planetary systems.

Given the computational expense and complexity of stellar cluster simulations that contain multiplanet systems, we devote Section 2 towards describing our methods. In Section 3, we report the results, and we state our conclusions in Section 4. Our output data tables are available as supplementary material in the online version of this paper.

2 METHODS AND INITIAL CONDITIONS

2.1 Computational approach

Stars form predominantly in groups, like stellar associations or star clusters (Lada & Lada 2003; Portegies Zwart, McMillan & Gieles 2010), and due to the close connection between star and planet formation, planets are accordingly born into these clustered environments. However, the simulation of multiplanetary systems in star clusters is challenging due to various reasons and requires different computational approaches depending on the underlying scientific question. One challenge for the numerical integration of the motions of the planets around the stars and the motion of the host star through the cluster is the completely different dynamical time-scales. While the dynamical evolution of planets takes place on time-scales of days and years, for star clusters it is typically in the range of several million years. Another aspect is the hierarchical nature of stars with (multi-)planetary systems. In principle, planetary systems can be treated and regularized similar to binary systems. Spurzem et al. (2009), who studied single-planetary systems in a star cluster in a fully coupled dynamical simulation, used this approach, as well as van Elteren et al. (2019), who studied multiplanetary systems in star clusters using the Nemesis module in AMUSE (Portegies Zwart 2011; Portegies Zwart & McMillan 2018). However, we want to be able to accurately trace resonant and secular effects in the dynamical evolution of the planetary systems. For this reason, we use a hybrid approach and simulate star cluster and planetary systems separately by using encounter information from the star cluster simulation for the integration of the planets. This approach is possible under the assumption that the motion of individual stars and the evolution of the whole star cluster can influence the dynamical evolution of the planets but not vice versa.

As a first step in this approach, the star cluster is simulated using NBODY6++GPU (Aarseth 2003; Wang et al. 2015, 2016, and references therein), where the motions of the stars are integrated using the Hermite scheme (e.g. Aarseth 2003; Aarseth, Tout & Mardling 2008). All necessary information is stored in high temporal resolution in the HDF5² format. Then, using the LonelyPlanets Scheme (LPS; Cai et al. 2015, 2017; Cai, Portegies Zwart & van Elteren 2018; Cai et al. 2019; Flammini Dotti et al. 2019; Stock et al. 2020), which is based on the AMUSE framework, all encounters of the selected host stars with each of the five nearest stars during the cluster simulation are calculated and stored, including the first and second time derivatives of the perturbers, in order to reconstruct the details of an encounter for the subsequent integration of the planetary systems. LPS uses REBOUND (Rein & Liu 2012) for the actual integration of the planetary systems, as well as additional features from REBOUNDx (Tamayo et al. 2020). For our simulations, we use REBOUND’s high-order, adaptive-step size integrator IAS15 (Rein & Spiegel 2015) to obtain accurate integration results of systems with close encounters between the planets.

2.2 Initial conditions for the star cluster simulation

We simulate an open star cluster consisting of 8000 stars whose masses follow a Kroupa (2001) initial mass function (IMF) in the mass range of |$0.08{-}20\, {\rm M}_\odot$|⁠. The total cluster mass is |$M_\mathrm{cl} = 4073.4\, \mathrm{M}_\odot$|⁠. The star’s initial positions and velocities in the cluster are drawn from a Plummer (1911) model. The star cluster is initially in virial equilibrium (|U| = 2T, where U is the total potential energy of the Plummer sphere and T is the total kinetic energy of the cluster stars). The virial radius, which is defined as |$r_\mathrm{vir}= \mathrm{G}M_\mathrm{cl}^2/(2|U|)$| (with G being the gravitational constant), is |$1\, \mathrm{pc}$| for our cluster, while the initial half-mass radius is |$r_\mathrm{hm} \approx 0.78\, \mathrm{pc}$|⁠. The star cluster is assumed to be on a Solar-like orbit around the Galactic Centre, which is why the tidal forces of the Galaxy acting on the cluster are assumed to be equal as for the solar neighbourhood (Heisler & Tremaine 1986). The cluster’s initial tidal radius r_tid = R_G(M_cl/M_G)^1/3 (with R_G being the distance to the Galactic Centre and M_G being the Galaxy’s mass contained inside R_G) is |$22.6\, \mathrm{pc}$|⁠. Stellar evolution is implemented (see e.g. Khalaj & Baumgardt 2015; Spera, Mapelli & Bressan 2015), but the mass loss is negligible for the host stars whose planetary systems are simulated only during a period when their host star is still on the main sequence. Primordial binary systems are not included; however, binaries can form during the course of the simulation. We also do not assume primordial mass segregation, but we observe the onset of mass segregation during the first few million years when the cluster experiences a short phase of core collapse.

According to equation (3) in Malmberg et al. (2007), encounters below |$r_\mathrm{min}=1000\, \mathrm{au}$| between our host stars and an average-massed star (⁠|$M_\star \sim 0.51\, \mathrm{M}_\odot$|⁠) in the cluster take place on time-scales of τ_enc ≈ 0.6–1.2|$\, \mathrm{Myr}$|⁠. This corresponds to encounter rates of 0.8–1.7 encounters per star and per Myr for the host star mass range used in our simulation.

We want to integrate the planetary systems until the cluster has sufficiently dissolved. Here, however, a compromise must be found between the computational costs of the planetary system simulations and the complete dissolution of the cluster. A good compromise for the 8000 star cluster we simulate is a period of |$100\, \mathrm{Myr}$|⁠. Although only about 20 per cent of the stars have completely escaped the cluster’s gravitational field by then, the cluster has already expanded significantly to |$r_\mathrm{vir} = 3.87\, \mathrm{pc}$|⁠, so that encounters between stars outside the dense core are very rare after more than |$100\, \mathrm{Myr}$|⁠. The median distance of the host stars to the cluster centre after |$100\, \mathrm{Myr}$| is around |$8\, \mathrm{pc}$| for the lowest and highest host star mass range, while the median distance to the nearest star is |$1.1\, \mathrm{pc}$| and |$1.4\, \mathrm{pc}$|⁠, respectively, and has increased by a factor of 11 and a factor of 15 compared to the beginning of the simulation. Thus, a maximum integration time of |$100\, \mathrm{Myr}$| for the planetary systems seems to be adequate. An optical comparison between the star cluster at the beginning of the simulation and after |$100\, \mathrm{Myr}$| is shown in Fig. 1.

2.3 Initial conditions for the planetary system simulations

In this work, we provide the results of 1224 planetary systems embedded in an open star cluster whose properties were described in Section 2.2. The host stars to be simulated are the typical progenitor stars of polluted WDs on the main sequence, which typically have masses of 1.5–|$2.5\, \mathrm{M}_\odot$| (see, for example, Tremblay et al. 2016; Cummings et al. 2018; El-Badry, Rix & Weisz 2018; McCleery et al. 2020; Barrientos & Chanamé 2021).

However, using a continuous mass spectrum for the host stars, as would be the case in a real star cluster, would reduce the comparability between the individual simulated planetary systems. For this reason, we divide the host stars into three different mass ranges, 1.25–|$1.75\, \mathrm{M}_\odot$|⁠, 1.75–|$2.25\, \mathrm{M}_\odot$|⁠, and 2.25–|$3.25\, \mathrm{M}_\odot$|⁠, and search for those stars whose masses lie in one of these mass ranges. Since the IMF drops off very steeply towards higher masses and the number of available stars is limited in the range around |$2.5\, \mathrm{M}_\odot$|⁠, the upper limit of the third mass range is deliberately chosen to be higher. We then calculate the encounters of all these stars with each of the nearest five stars in the cluster and store this information. For the subsequent integration of the planetary systems, the masses of the host stars are set to |$1.5\, \mathrm{M}_\odot$|⁠, |$2.0\, \mathrm{M}_\odot$|⁠, and |$2.5\, \mathrm{M}_\odot$| to ensure the comparability of the planetary systems within these three mass ranges and to be able to work out the pure effect of the cluster environment on the dynamical evolution of the individual systems. According to the number of stars present in the three mass ranges, we have a total of 408 host stars available (193 stars with |$M_\star = 1.5\, \mathrm{M}_\odot$|⁠, 114 stars with |$M_\star = 2.0\, \mathrm{M}_\odot$|⁠, and 101 stars with |$M_\star = 2.5\, \mathrm{M}_\odot$|⁠).

The planetary systems around these 408 host stars are then started in three different initial orbital configurations, while the host star and its trajectory through the cluster remain the same for all three different planetary system models. All planetary systems solely consist of super-Earths, each having a mass of 0.01 M_Jup (≈3.2 M_⊕). Due to the variety in multiplicity of actual observed planetary systems, we aim to simulate two bounding cases of systems consisting of three and seven planets. However, the compactness of a seven-planet system crucially determines its dynamical evolution, especially if it is externally perturbed (by stellar flybys).

Therefore, we simulate the following three scenarios in which all planets are equally separated in terms of mutual Hill radii (R_H,m = (a₁ + a₂)/2 · ((M₁ + M₂)/(3M_⋆))^1/3): (i) a system consisting of only three planets (called ‘3P model’) between 2.0 and 18.63 au, (ii) a rather tightly packed system of seven planets, called ‘7PC model’, within the same orbital boundaries as model 3P, and (iii) a wider system of seven planets, called ‘7PW model’, in which the five innermost planets are placed in the same orbital range as in the previous two cases, with two additional planets on wider orbits (the outermost planet has |$a=56.87\, \mathrm{au}$|⁠, resulting from the fixed number of mutual Hill radii).

The orbital configurations of the three models and the number of mutual Hill radii used for the orbital spacing are listed in Tables 1 and 2, respectively. In all three models, the planetary orbits are initially circular (e = 0) and coplanar (i = 0°). Furthermore, all systems are long-term stable if they are placed in isolation. The argument that the planetary systems should be stable in isolation over time was also decisive for outwardly increasing spacings between the planets, which we achieve by using mutual Hill radii instead of using similar orbital spacings according to the peas-in-a-pod theory (e.g. Millholland, Wang & Laughlin 2017; Weiss et al. 2018) based on findings from the Kepler mission.

The inner boundary of 2.0 au is chosen as a minimum semimajor axis to account for the potential engulfment of the innermost planet due to the expansion of the host star during the giant branch phase. As a basis for the outer boundary in the first two cases, we take into account that core accretion during planet formation becomes inefficient at larger semimajor axes, although it is not impossible for super-Earths to form at wider orbits (see, for example, fig. D.3 in Schlecker et al. 2021), which is why we additionally consider the possibility of a more extended planetary system in the third scenario.

3 RESULTS AND DISCUSSION

3.1 Fraction of surviving planets

If the eccentricity of a planet is excited to more than e > 0.99, we assume that the planet is close to ejection and remove it from the simulation. These ejected planets would, depending on their escape velocity, either continue to move through the cluster as free-floating planets, which may even allow re-capture by other cluster members, or they would directly escape not only the gravitational field of the host star but also that of the star cluster.

For technical reasons, we do not trace the motion of the planets through the cluster after their ejection from a planetary system. However, the fraction of ejected planets f_ej gives an estimate for the expected number of free-floating planets in open clusters with similar properties to the one we simulated. The fraction of surviving planets f_surv = 1 − f_ej is plotted in Fig. 2 for all three planetary models as a function of simulation time. The model with the highest |$\overline{f}_\mathrm{surv}$| (dotted black line in Fig. 2), averaged over all planets in the system, is the 3P model with a value of 0.76. The 7PC model shows little difference with |$\overline{f}_\mathrm{surv} = 0.74$|⁠, indicating that despite the higher planet density in this system compared to the 3P model, which in principle leads to more planet–planet interaction and to higher ejection rates, the orbit width of the outer planets is the more important factor for the averaged survival fraction. Consequently, the planets in the 7PW model have on average the lowest survival probability with a value of |$\overline{f}_\mathrm{surv.} = 0.71$|⁠. The fraction of escapers that arise in our simulations are slightly higher than, e.g. in van Elteren et al. (2019), who obtained |$\overline{f}_\mathrm{ej} \approx 0.14$|⁠. However, the difference can be explained by the significantly shorter simulation time and the considerably smaller number of stars in the host star cluster in van Elteren et al. (2019).

However, when considering the survival probability for the individual planets, the exact configuration of the planetary system, especially its multiplicity and consequently its compactness, does play a role. While in the 3P model the planet’s probability for being ejected correlates with its initial semimajor axis, this is not consistent with the 7PC and 7PW models, as can be seen in Fig. 2.

The planet with the highest survival fraction in the 7PC model is P3 (f_surv = 0.78), followed by P4 (f_surv = 0.77), P2 (f_surv = 0.76), and P1 (f_surv = 0.74), so the survival rate increases slightly for the middle planets and decreases only from the fifth planet towards the outer planets P7 (which has f_surv = 0.70). The spread in survivability for the 3P model is only slightly larger and ranges from f_surv = 0.68 to f_surv = 0.82. The reason for the changed order in the survival fraction of the planets in the 7PC model is the higher planetary density at constant orbital expansion of the system. Due to increased interactions among the planets after an external gravitational perturbation, the inner planets can experience delayed ejection from the system indirectly as a result of an earlier flyby of a neighbouring star. These delayed ejections have also been observed in van Elteren et al. (2019) and Stock et al. (2020).

As expected, for a system with wider orbits but the same number of planets, as in the 7PW model, the spread in the individual planet’s survival rate is larger than for the more compact case. Here, the values are between f_surv = 0.60 (P7) and f_surv = 0.79 (P2). As in the 7PC model, the second innermost planet in the 7PW model has a slightly larger survival fraction than the innermost planet P1 (f_surv = 0.76), but here the planets’ survival probability decreases beyond the second planet as expected with increasing initial semimajor axes.

3.2 Semimajor axis and eccentricity distribution and possible engulfment during red giant phase

The fraction of planetary systems whose dynamical evolution is considerably perturbed by passing stars (directly or indirectly by delayed planet–planet scattering) depends on the one hand on the planetary model used but also on the host star mass. For the 3P model and a host star mass of |$1.5\, \mathrm{M}_\odot$|⁠, we generally observe the lowest effect of the stellar environment on the dynamical evolution of the individual planets. In this scenario, 83 per cent of the planets remain largely unperturbed. As a criterion for a considerable perturbation, we look at whether the semimajor axis deviates by more than 5 per cent from the initial value by the end of the simulation, or whether the eccentricity increases to more than 0.1. The fraction of planets that are significantly perturbed in their dynamics increases for the models with more planets per star, the orbital separation of the outermost planet, and the host star mass. For the 7PW model and a host star mass of |$2.5\, \mathrm{M}_\odot$|⁠, the fraction of perturbed planets increases from 17 per cent to 29 per cent. The distribution in a–e space is wide for those planets that were considerably perturbed, as we demonstrate in Fig. 3, and comparable to results from previous studies (see e.g. fig. 10 in Malmberg, Davies & Heggie 2011). We observe inward migration for the perturbed planets in a few cases but outward migration in the vast majority of cases. Almost 1 per cent of all planets have wide orbits of more than |$100\, \mathrm{au}$| at the end of the simulation, and one planet out of a total of 6936 planets is even scattered to an orbit of more than |$1000\, \mathrm{au}$|⁠. More than 6 per cent of all planets are excited to high eccentricities (e > 0.5). The fractions vary depending on the planetary system model and host star mass and are listed in Table 3 for all the different scenarios.

For the question of whether the planetary system model used or the mass of the host star (and thus the stellar density in the vicinity of a planetary system) has a stronger influence on the formation of high eccentricities, we plot in Fig. 4 the cumulative distribution of eccentricities after |$100\, \mathrm{Myr}$|⁠. We distinguish between the three planetary system models (independent of the host star mass) and the host star mass (independent of the used planetary system model). For those systems that orbit a |$2.5\, \mathrm{M}_\odot$| host star, it can be clearly seen that the host star mass, and thus the position in the cluster, which in turn is related to the stellar density in the vicinity, plays a more important role in exciting planets to high eccentricities than the exact orbital configuration and multiplicity of the planetary system.

Many studies have investigated the critical engulfment distance during the giant branch phases at different levels of detail and with different underlying theories (Mustill & Villaver 2012; Adams & Bloch 2013; Nordhaus & Spiegel 2013; Villaver et al. 2014; Madappatt, De Marco & Villaver 2016; Privitera et al. 2016; Ronco et al. 2020). We use the critical engulfment distances along the asymptotic giant branch phases for Earth–mass planets from figs 2–4 in Mustill & Villaver (2012) and calculate for how many planets the periastron distance (r_p = (1 − e)a) would be below this limit. For the planets around a |$1.5\, \mathrm{M}_\odot$| star, the critical distance is about |$1.9\, \mathrm{au}$|⁠. About 5 per cent of our planets around such a star would be engulfed during the giant branch phases. The critical distance for 2.0 and |$2.5\, \mathrm{M}_\odot$| stars is 2.2 and |$2.3\, \mathrm{au}$|⁠, respectively. In these cases, 16 per cent and 15 per cent of the planets would be engulfed, respectively. Although the possibility of being engulfed by the host star mainly affects the innermost planet P1, it is not exclusive, since an external perturbation combined with internal planet–planet scattering can cause even the initially outermost planet to migrate to a small or highly eccentric orbit, causing its periastron distance to fall below the critical value. This can be seen in the bottom panels of Fig. 3, where we additionally plot the critical engulfment distance for each host star mass.

3.3 Inclination and retrograde orbits

The first exoplanets thought to have a polar or retrograde orbit (i ≥ 90°) were HAT-P-7 b (Winn et al. 2009) and WASP-17 b (Anderson et al. 2010; Bayliss et al. 2010). Since the number of confirmed retrograde planetary orbits is still small, the statistical abundance of these peculiar orbits is still uncertain. In addition to the expected clustering of prograde orbits, Albrecht et al. (2021) recently found a further clustering of polar orbits in a sample of 57 systems rather than a scattering over the entire range of possible obliquities. Since all planets in our simulations are initially co-planar, most planetary orbits are still only slightly inclined at the end of the simulations.

The distribution of planetary orbits in the a–i space for all three planetary system models at the end of the simulation is shown in Fig. 5.

The averaged fraction of planets with a retrograde motion at the end of the simulation is 1.4 per cent, with values ranging from 0.5 to 2.2 per cent for the different system models (the fraction of retrograde orbits for each scenario used in this work is listed in Table 3). This is somewhat higher but still in good agreement with the values in Stock et al. (2020), where we used variations of the Solar system around Sun-like host stars and found, depending on the initial planet configuration and the star cluster size, 0.1–1.6 per cent of all planets to be on a retrograde orbit after the same simulation time of |$100\, \mathrm{Myr}$|⁠. The subtle differences can be explained by the higher multiplicity in the 7PC and 7PW models as well as by the generally larger host star masses used in this study. This agreement, and the circumstance that we cover a wide range of possible planetary systems in this work and in Stock et al. (2020), leads us to the rough estimate that in open star clusters similar to the one simulated in this work and those simulated in Stock et al. (2020), about 1–2 per cent of all planets could be on stable retrograde orbits.

The number of planets that flip to a retrograde orbit for at least one integration step at some time during the simulation is significantly larger and gives an indication that unstable retrograde orbits are not uncommon in environments with frequent external gravitational perturbation. ‘Unstable retrograde orbit’ in this context means that the planet does not remain permanently on a retrograde orbit, either because it changes back to a prograde orbit or because it is ejected from the planetary system at a later time. In 33 per cent of all systems, we find at least one planet that flips to a retrograde orbit for at least one (stored) time-step of 1000 yr during the simulation. This fraction of systems is generally lowest for the 3P model and |$1.5\, \mathrm{M}_\odot$| stars and highest for the 7PW model and |$2.5\, \mathrm{M}_\odot$| stars.

Since inclined orbits, just like eccentric orbits, form through angular momentum exchange, and close stellar encounters are particularly good at introducing an angular momentum deficit into the planetary system, it is especially the systems around |$2.5\, \mathrm{M}_\odot$| stars that have inclined orbits, as can be seen in Fig. 6. Again, the mass of the host star and thus the frequency and strength of the encounters with other cluster members are more important than the exact planetary configuration for the formation of inclined planetary orbits.

3.4 Mean-motion resonances

Resonances in planetary systems can be a source for WD pollution, because the orbits of asteroids can increase in eccentricity due to planets near a secular or mean-motion resonance (MMR) (Debes et al. 2012; Smallwood et al. 2018, 2021; Antoniadou & Veras 2019; Veras et al. 2021). Finding MMRs in simulations is very challenging in view of the large number of simulated planetary systems, the long integration time, but especially because of the often chaotic dynamical evolutions of the planetary systems due to the steady external perturbation from the cluster. This type of dynamical evolution can lead to transitory resonances, which may endure over just a handful of output time-steps, or even within two consecutive outputs of the simulations.

In principle, we find MMR in our simulations at arbitrary times (except at the beginning), but they only rarely survive for longer times (i.e. several million years) due to the constant gravitational perturbation of the neighbouring stars. Only with increasing simulation duration and the expansion of the host star cluster, when the frequency and strength of the encounters with neighbouring stars decrease, can the planets actually remain in MMR for several million years. In particular, we therefore investigate how many and which MMRs are found in the last 1 million years of our simulations that are stable until the end of the simulation at |$t = 100\, \mathrm{Myr}$|⁠, since these can also persist once the cluster has completely dissolved. If we consider only planetary pairs that were direct neighbours at the beginning of the simulation, we find six planetary pairs that are in stable MMR at the end of the simulation and list them in Table 4.

For the 2:1 MMR in system 122 from the 7PC model with a |$1.5\, \mathrm{M}_\odot$|⁠, we show the FAIR plot for the time range t = 99–|$100\, \mathrm{Myr}$| in Fig. 7 and, in addition, also plot the orbital elements as well as the resonance angle of the 2:1 MMR as a function of the total simulation time in Fig. 8. This makes it possible to see whether a system was excited into this resonance by chance only at the end of the simulation or whether the planets were already close to resonance for a long time.

Figure 7.

2:1 MMR in planetary system 122 (7PC model, |$1.5\, \mathrm{M_\odot }$| host star) between planet 6 and 7 for t = 94.5–|$100\, \mathrm{Myr}$|⁠.

Open in new tabDownload slide

Figure 8.

Orbital elements a, e, and i for planetary system 122 (7PC model, |$1.5\, \mathrm{M_\odot }$| host star) and the resonance angle θ₂ with p = q = 1 (for planets 6 and 7) as a function of the simulation time. The thicker grey lines in the background represent the distance to the cluster centre (top panel), the distance to the closest stellar perturber (top middle panel), and the normalized angular momentum deficit (bottom middle panel), defined as |$\mathrm{NAMD} = \mathrm{AMD}/(\sum _k m_k \sqrt{\mathrm{G}m_\star a_k})$| (see Turrini, Zinzi & Belinchon 2020, and references therein).

Open in new tabDownload slide

As can be seen in Fig. 8, planets 6 and 7 both migrate outwards to eccentric crossing orbits with a ∼ 23–|$24\, \mathrm{au}$| as a result of a close encounter (⁠|$r_\mathrm{ p} \lt 231\, \mathrm{au}$|⁠) at |$t=8.4\, \mathrm{Myr}$| with a |$0.2\, \mathrm{M}_\odot$| star. After another encounter at time |$t=11.8\, \mathrm{Myr}$|⁠, both planets are thrown to orbits of |$a=66\, \mathrm{au}$| and |$a=91\, \mathrm{au}$|⁠, and thus happen to be near 2:1 MMR. Due to several further but weak perturbations, they first enter 2:1 MMR at |$t=54.2\, \mathrm{Myr}$|⁠, which does not completely resolve until |$t=66.2\, \mathrm{Myr}$| (however, the libration amplitude is very large in the interim). Due to subsequent weak perturbations, the planets re-enter 2:1 MMR at |$t=80\, \mathrm{Myr}$|⁠. The value around which the resonance angle librates changes during the remaining simulation time, but in principle, the planets remain in 2:1 MMR until the end of the simulation. That both planets are in stable resonance over the last 1 million years can be seen in Fig. 7.

Planet–planet scattering as a cause of MMRs is often disregarded. Yet, they can be particularly responsible for the higher-order resonances (Raymond et al. 2008), as these cannot arise so easily through migration (Papaloizou & Szuszkiewicz 2005; Tadeu dos Santos et al. 2015; Xu, Lai & Morbidelli 2018). In our simulations, however, star–planet interactions seem to be more important than planet–planet interactions. In addition to four first-order resonances, we find a third-order resonance and a seventh-order resonance, and plot the resonances (Figs A1, A3, A5, A7, and A9 in the appendix) and the overall dynamical evolution of the systems including perturber information (Figs A2, A4, A6, A8, and A10 in the appendix) for each of these systems. However, some librating resonance angles show a long-term trend that cannot be resolved in time. Whether the systems are in actual – and stable – resonance cannot be said with certainty in these cases. In all systems, one or several encounters with neighbouring stars lead to a migration of the planetary pairs near a certain MMR. Further weaker stellar perturbations, usually millions of years later, drive the planetary pair into actual resonance. Only in the case of system 192 from the 7PW model around a |$2.5\, \mathrm{M}\odot$| star (see Fig. A10 in the appendix), a resonance is created directly by the first strong stellar perturbation, but this resonance is repeatedly perturbed by stellar neighbours at later times without breaking it completely.

Our results are only partially comparable to Raymond et al. (2008), since the instability in our simulations has an external rather than an internal origin. In most cases, the external perturbation completely outweighs internal effects. In cases where planet–planet scattering may play an additional role in the origin of the resonance, the effects of external and internal perturbation cannot be separated clearly enough. However, our simulations confirm that most of the resonances that arise are low-order and that higher-order resonances can also arise in a few cases. Furthermore, an important difference from the simulations in Raymond et al. (2008) is the mass distribution within the planetary system. While we exclusively simulate planets of equal masses, Raymond et al. (2008) also use mixed systems, in which the interplay of large and small planet plays an important role in the formation of resonances due to planet–planet scattering.

3.5 Long-term stability

The focus of this investigation is to perform stellar cluster simulations including multiplanet systems and to present the results as initial conditions, which may be used for the simulation community. Nevertheless, we can make some preliminary rapid judgements about the stability of these systems along the main sequence by taking advantage of machine learning.

In order to estimate the long-term stability of each planetary system (with the host star mass remaining constant), we perform a SPOCK test (Tamayo et al. 2020) for each system. SPOCK is able to predict the long-term stability of compact multiplanet systems by using the statistics from machine learning training data sets in a fraction of the time compared to actual integration (see Tamayo et al. 2020, for a detailed description of the training and the model). Since SPOCK requires at least three planets in the system, for those systems where planets have been ejected, we replace the missing planets with massless particles at wide orbits (⁠|$\gt 100\, \mathrm{au}$|⁠) before performing the SPOCK test. We give the likelihood for the long-term stability of our planetary systems as additional feature in our result tables (see online or extracts in Tables B1 and B2 in the appendix) for the different planet system models.

4 CONCLUSIONS

We have presented the results from a total of 1224 simulations of planetary systems embedded in a cluster environment, differing in the planetary system model used and in the mass of the host star. We have aimed to publish a comprehensive data set that contains those planetary systems typically responsible for WD pollution, which, at the same time, have the dynamical imprint of a typical birth star cluster. This data set can now be used for further numerical integration beyond the main sequence to the WD phase.

The three different planetary system models should represent the extreme cases of possible planetary systems regarding their multiplicity and orbital spacing. For this reason, we simulate one model with only three planets (3P model), and two models with seven planets each in the system, which differ in their compactness (7PC and 7PW model). As host stars, we have chosen those stars from our 8000 star cluster that are most similar to the masses of |$1.5\, \mathrm{M}_\odot$|⁠, |$2.0\, \mathrm{M}_\odot$|⁠, and |$2.5\, \mathrm{M}_\odot$|⁠. All planetary systems are integrated to |$t=100\, \mathrm{Myr}$|⁠. By that time the star cluster has expanded considerably and perturbations from neighbouring stars are very rare.

In our simulations, it was not only the number of planets or the compactness of the system that played a role in the average survival rate but above all, the semimajor axes of the outermost planets. The 3P model has an average survival fraction of 76 per cent, the 7PC model has 74 per cent, and the 7PW model, with the widest orbits, has only 71 per cent. While the innermost planet has the highest survival probability in the 3P model, it is the third planet in the 7PC model and the second innermost planet in the 7PW model. The spread in survivability has been particularly small in the 7PC model, which is why we conclude that the compactness of the system and thus enhanced internal effects such as planet–planet scattering almost equalize the planets’ probability to survive the star cluster phase.

Given our initial conditions, we found that about 5 per cent of planets around |$1.5\, \mathrm{M}_\odot$| stars, roughly 16 per cent of the planets around |$2.0\, \mathrm{M}_\odot$|⁠, and approximately 15 per cent of the planets around |$2.5\, \mathrm{M}_\odot$| stars would be swallowed by the eventual asymptotic giant branch star’s envelope because the planets’ periastron distances would be below the critical engulfment distance.

The excitation in eccentricity and inclination correlates with the number of planets in the system and the initial semimajor axis of the outermost planet. In particular, it also correlates with the stellar density in the vicinity of the host star, which tends to be larger for higher mass stars due to the effect of mass segregation. On average, 1.4 per cent of all planets are on a retrograde orbit at the end of the simulation, which is, due to the higher host star masses and the higher multiplicity in the 7PC and 7PW model, somewhat higher but still in good agreement with the results from Stock et al. (2020).

Eccentric planets in the super-Earth mass regime are thought to be particularly efficient drivers for WD pollution over a wide range of cooling ages (Frewen & Hansen 2014; Mustill et al. 2018). 30 per cent of the planets in our simulations attained an eccentricity of e > 0.1 and 25 per cent have e > 0.17 after |$100\, \mathrm{Myr}$|⁠, showing that the birth environment of planetary systems can produce a sufficient distribution in eccentricity to help generate the architectures suitable for dynamical delivery of pollutants to WDs. Even if subsequent increases in eccentricity due to mutual perturbations and/or along the giant branch phases due to stellar mass loss alone (Veras et al. 2011) are negligible, the planets’ primordial eccentricities will persist into the WD phase.

Furthermore, we find planetary pairs in several planetary systems that are in resonance at the end of the simulation. These systems may also play a role in WD pollution, since asteroids near these resonances may be driven to eccentric orbits and subsequently be tidally disrupted by the WD (Smallwood et al. 2018, 2021).

SUPPORTING INFORMATION

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

We thank the reviewer for helpful comments that have improved the manuscript. We would also like to thank Dr. Martin Schlecker for the valuable discussions on the initial conditions of our simulations. KS and RS acknowledge the support of the DFG priority program SPP 1992 ‘Exploring the Diversity of Extrasolar Planets’ (SP 345/20-1 and SP 345/22-1). DV gratefully acknowledges the support of the STFC via an Ernest Rutherford Fellowship (grant no. ST/P003850/1). The simulations were carried out partly on the GPU accelerated cluster ‘kepler’, funded by Volkswagen Foundation grants 84678/84680. The authors acknowledge support by the High Performance and Cloud Computing Group at the Zentrum für Datenverarbeitung of the University of Tübingen, the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grant no. INST 37/935-1 FUGG.

DATA AVAILABILITY

The results of the planetary simulations are available as supplementary material in the online version of this paper. The data underlying this article are available on Zenodo at doi:10.5281/zenodo.5883613 in time-steps of 2000 yr. The data underlying this article will be shared in full temporal resolution (time-steps of 1000 yr) on reasonable request to the corresponding author.

Footnotes

REFERENCES

APPENDIX A: MEAN-MOTION RESONANCES AT THE END OF THE SIMULATIONS

Figure A1.

10:3 MMR between planets 2 and 3 in planetary system 177 (3P model, |$1.5\, \mathrm{M_\odot }$| host star) between 99 and |$100\, \mathrm{Myr}$|⁠.

Open in new tabDownload slide

Figure A2.

As Fig. 8 but for planetary system 177 (3P model, |$1.5\, \mathrm{M_\odot }$| host star). The resonance angle is shown for p = 3, q = 7, and planet pair 2/3.

Open in new tabDownload slide

Figure A3.

3:2 MMR between planets 5 and 6 in planetary system 22 (7PW model, |$2.0\, \mathrm{M_\odot }$| host star) between 99 and |$100\, \mathrm{Myr}$|⁠.

Open in new tabDownload slide

Figure A4.

As Fig. 8 but for planetary system 22 (7PW model, |$2.0\, \mathrm{M_\odot }$| host star). The resonance angle is shown for p = 2, q = 1, and planet pair 5/6.

Open in new tabDownload slide

Figure A5.

3:2 MMR between planet 5 and 6 in planetary system 45 (7PW model, |$2.0\, \mathrm{M_\odot }$| host star) between 98.5 and |$100\, \mathrm{Myr}$|⁠.

Open in new tabDownload slide

Figure A6.

As Fig. 8 but for planetary system 45 (7PW model, |$2.0\, \mathrm{M_\odot }$| host star). The resonance angle is shown for p = 2, q = 1, and planet pair 5/6.

Open in new tabDownload slide

Figure A7.

5:2 MMR between planets 4 and 5 in planetary system 38 (7PW model, |$2.5\, \mathrm{M_\odot }$| host star) between 99 and |$100\, \mathrm{Myr}$|⁠.

Open in new tabDownload slide

Figure A8.

As Fig. 8 but for planetary system 38 (7PW model, |$2.5\, \mathrm{M_\odot }$| host star). The resonance angle is shown for p = 2, q = 3, and planet pair 4/5.

Open in new tabDownload slide

Figure A9.

4:3 MMR between planets 3 and 4 in planetary system 192 (7PW model, |$2.5\, \mathrm{M_\odot }$| host star) between 99.5 and |$100\, \mathrm{Myr}$|⁠.

Open in new tabDownload slide

Figure A10.

As Fig. 8 but for planetary system 192 (7PW model, |$2.5\, \mathrm{M_\odot }$| host star). The resonance angle is shown for p = 3, q = 1, and planet pair 3/4.

Open in new tabDownload slide

APPENDIX B: ADDITIONAL MATERIAL