planet-disc interactions – planets and satellites: dynamical evolution and stability – circumstellar matter

1 Introduction↩︎

Debris disks are belts of particles, from planetesimals to small dust grains, surrounding many main-sequence stars (e.g. [1]–[4]). Our own Solar System has two debris disks, the Asteroid Belt and Kuiper Belt, with the latter located beyond Neptune’s orbit.

Extrasolar debris disks are detected through observations of micrometre- to millimetre-sized dust grains, in either thermal emission (e.g. [5]) or scattered light (e.g. [6]). However, where does this dust come from? A dust grain in a debris disk has a relatively short lifetime compared to the host star, before being removed or destroyed by various processes. Dust therefore needs to be continuously replenished from some source. It is thought that most dust originates from larger planetesimals with longer lifetimes, which collide with each other (e.g. [7]–[9]). This manifests itself as a collisional cascade, where larger bodies in a debris disk collide and break apart into smaller pieces, which in turn collide and break into even smaller pieces, and so on (e.g. [10], [11]). For debris particles to reach sufficiently high velocities to break apart in collisions, their eccentricities and/or inclinations must somehow be excited from the almost circular, unstirred orbits they are assumed to have formed with; this excitation process is known as “stirring” ([1]).

There are several commonly considered mechanisms for how debris disks get stirred. The first, secular-planet-stirring, is where a nearby, eccentric planet excites debris through secular interactions [12]. The efficiency of this process depends on the planet mass, eccentricity, and semimajor axis, but is only possible if the planet’s orbit is eccentric.

The most commonly considered alternative mechanism is self-stirring, where larger planetesimals within the disk excite smaller debris (e.g. [13]–[17]). This process does not require planets to be present. [13], [15] proposed that large debris bodies of size \(\sim\)​500-1000 km would be sufficient to stir debris disks. Later, [17] concluded that such large objects are not necessary for stirring (they take a long time to form, even though they would stir the disk quickly), and that bodies of \(\sim\)​200 km in size, which form much quicker, are sufficient to self-stir a debris disk.

However, these stirring models do not work in all cases. For secular-planet-stirring, not only must the planet’s orbit be sufficiently eccentric, but the self-gravity of the disk can significantly resist or even cancel the stirring by the planet (Löhne et al., in prep.). For self-stirring, the required disk mass is in some cases unphysically high; [18] concluded that if all debris disks contained objects up to 200 km in size, the masses of some disks would be in the range of \({10^3 - 10^4 \; {\rm M_\oplus}}\), which exceeds the plausible maximum of \({10^2 - 10^3 \; {\rm M_\oplus}}\)thought to be inheritable from protoplanetary disks. Similarly, [19] analyzed 67 resolved debris disks and concluded that 26 of these would need masses greater than \({10^2\; {\rm M_\oplus}}\) to self-stir if considering a conventional debris-size distribution, with 6 of these requiring masses greater than \({10^3 \; {\rm M_\oplus}}\). Therefore, in some cases self-stirring appears to require unphysically high debris-disk masses to operate (though also see [20]).

Several other stirring mechanisms have also been proposed, including “pre-stirring” (where debris is already stirred during the protoplanetary disk phase; e.g. [21]–[23]) and “flyby stirring” (where another star passes close to the system and stirs the debris disk; e.g. [24]–[26]). However, the prevalence of these processes is unclear.

The aim of our study is to look at other possibilities of stirring, which could help explain some cases where secular-planet or self-stirring falter. We identify two new mechanisms that can significantly increase the efficiency of planet-stirring: “projectile stirring” and “resonant stirring”. We define projectile stirring to be where a planet scatters one or more large bodies (“projectiles”) from an initially unstirred debris disk onto eccentric orbits. These newly eccentric projectiles can in turn sufficiently increase the eccentricities of debris within the disk through secular effects, before the planet eventually ejects the projectiles from the system. We will demonstrate that even low-mass projectiles can significantly increase planet-stirring efficiency.

We define resonant stirring to be where a planet stirs a wide region of a debris disk via its mean-motion resonances (MMRs). This effect has historically been overlooked as a stirring mechanism; whilst it is well known that MMRs can excite debris in narrow regions around nominal MMR locations (e.g. [27]–[30]), it is often overlooked that individual MMRs can be very broad, particularly for high-mass planets (e.g. [31]). This means that MMRs can significantly excite debris not only at the nominal resonance locations, but between those locations as well. We will demonstrate that broad MMRs can stir debris across a very wide region of a debris disk, far wider than simply around the nominal MMR locations as might naı̈vely be assumed. We will also show that this effect is efficient even for planets with very low eccentricities, in contrast to the traditional, secular model of planet-stirring.

The concept of projectile stirring is motivated by the increasing prevalence of suspected scattered populations in debris disks. The Kuiper Belt hosts a scattered population, comprising bodies that previously had close encounters with Neptune [32]. This population is highly eccentric, and these scattered bodies cross the orbits of other, dynamically colder Kuiper-Belt Objects. The scattered population includes dwarf planets such as Eris, demonstrating that high-mass bodies can be scattered onto disk-crossing orbits by planets located at the inner edges of debris disks. There are also tentative indications that even larger bodies may once have existed in this region, which got scattered and subsequently ejected from the Solar System (e.g. [33]). These arguments motivate questions about the dynamical impact of scattered bodies on the excitation level of debris disks. In addition to the Solar System, there is also evidence of scattered material in extrasolar debris disks; [34] found evidence for a scattered disk in \({{\rm HR} \;8799}\), and [35] identified both “hot” and “cold” dynamical populations in the \(\beta\) Pictoris debris disk. In addition, the sharp inner edges of some debris disks could indicate the presence of sculpting planets at those locations [36], [37], and such planets would generate scattered populations. Based on both the Solar System and extrasolar systems, scattered populations could be a common feature of debris disks; given that observed disks must be dynamically excited, and that there are problems associated with known stirring mechanisms, it is pertinent to ask whether scattered material could be helping to stir observed debris disks.

Both projectile and resonant stirring are new ways to increase planet-stirring efficiency beyond the original secular effect of [12]. [38] have also recently investigated a similar topic that they refer to as “mixed stirring”. In their work, they simulated systems with both a giant planet and a debris disk comprising massive bodies, in order to investigate stirring through a combination of effects (see also [39]). In contrast, we focus on two specific stirring processes (projectile and resonant stirring) and isolate them from other interactions to assess their individual effectiveness. These effects likely also occur in the simulations of [38].

In this paper we run \(N\)-body simulations across a wide parameter space, to quantify the effectiveness of projectile and resonant stirring. We first describe how we quantify stirring, including providing a public code to quantify the stirring level in a rebound simulation (Section 2). We then present the setup and results of our N-body simulations (Section 3). In Section 4 we investigate projectile stirring in more detail, including deriving an analytic criterion for when it will occur (Equation 10 ), and in Section 5 we further investigate resonant stirring. In Section 6 we discus our results, and we conclude in Section 7. The Appendix contains more details about our stirring analyses.

2 Stirring definitions↩︎

In Section 3, we will run N-body simulations to investigate the effectiveness of projectile and resonant stirring. However, before proceeding we must first define quantitatively what we mean by stirring. For this paper we consider two different stirring definitions; first, a simple analytic form that quantifies the minimum eccentricity that two adjacent bodies would need to collide destructively, and second, a detailed numerical approach that properly accounts for the specific orbits and orientations. The former provides a rough idea of whether a disk is stirred or not, whilst the latter provides a much more accurate picture and is the method that we use to assess our simulations. We outline both analyses in this section, and give more details in Appendix 8.

2.1 Rough criterion for debris to be considered stirred↩︎

We first perform a simple, rough calculation of the minimum eccentricity that a debris particle must attain to be considered stirred. We consider a system containing a star and two massless debris bodies. We assume that these particles have co-planar orbits, and that the two orbits are rotationally offset but otherwise identical (see Figure 12). In Appendix 8.2 we show that such bodies are unstirred if their eccentricities are below some value:

\[e_{\rm unstirred} < \frac{v_{\rm frag}}{2}\sqrt{\frac{a_{\rm deb}}{G m_*}}, \label{eq:32unstirred95ecc}\tag{1}\]

where \(v_{\rm frag}\) is the fragmentation speed of a debris body, which is dependent on its size and material, \(G\) is the gravitational constant, \(a_{\rm deb}\) is the semimajor axis of the debris particles, and \(m_*\) is the star mass. We can also relate this to \(v_{\rm circ}\), the speed of a circular orbit at the debris semimajor axis; since \({v_{\rm circ}^2 = Gm_*/a_{\rm deb}}\), Equation 1 can alternatively be written as

\[e_{\rm unstirred} < \frac{v_{\rm frag}}{2v_{\rm circ}}. \label{eq:32unstirred95ecc95v95circ}\tag{2}\]

For this paper we want to evaluate Equation 1 , to roughly constrain whether debris could be stirred as a function of its semimajor axis. This requires us to evaluate the fragmentation speed, \(v_{\rm frag}\), which means we must make assumptions about debris size and composition. First, in line with previous studies, we approximate the material to be basalt using the prescription of [11]. This yields a size-dependent form for \(v_{\rm frag}\), given by Equation 14 . Then, we assume the debris particles to have radii of \({1\; {\rm cm}}\); this is because debris-disk observations are sensitive to grains of millimetre size or smaller, so larger particles must be undergoing destructive collisions to produce the observed dust. Hence a pair of centimetre-sized, basalt grains on identical, azimuthally rotated orbits are unstirred if their eccentricity is less than

\[e_{\rm unstirred}(1~{\rm cm}) < 8.9 \times 10^{-4} \left( \frac{a_{\rm deb}}{\rm au}\right)^{1/2} \left( \frac{m_*}{\; {\rm M_\odot}}\right)^{-1/2}. \label{eq:32final95unstirred95ecc}\tag{3}\]

Note that having an eccentricity above this value is not a sufficient condition to be stirred, since collisional speed depends on relative orientation to adjacent-particle orbits; nonetheless, debris whose eccentricities never exceed this critical value are almost certainly unstirred. The Equation 3 criterion will be displayed on figures throughout our paper, to provide a rough idea of the stirring level in N-body simulations.

2.2 Full, numerical stirring analysis↩︎

Equation 3 is only a rough estimate of whether debris is stirred, because it assumes that nearby debris are on identical, optimally misaligned orbits. A real disk may contain debris with a wide range of orbits and orientations. To properly assess stirring in our N-body simulations, we will therefore use a numerical method instead.

The principle of our numerical stirring analysis is to take the results of an N-body simulation, and post-process it to determine the stirring level. We do this by numerically analysing pairs of debris particles with instantaneously overlapping orbits, to calculate the collision speed at the orbit-intersection points. If this speed is greater than the fragmentation speed, again assuming \({1\; {\rm cm}}\) basalt grains, then both particles are considered stirred and are not included in future stirring calculations. Conversely, if the collision speed is below the fragmentation speed, then a different pair of particles is considered instead. Only pairs of particles that were initially near each other are considered, so the stirring level at the outer edge of a broad disk would not be contaminated by high-eccentricity debris scattered by a planet at the inner edge. The process continues until all suitable particle pairs have been checked, and is then repeated across multiple simulation snapshots at different times. This yields all debris particles that have ever been stirred in a simulation, and is much more accurate at assessing the disk’s stirring level than Equation 3 , because it properly accounts for different orbit orientations and semimajor axes. The full details of this numerical analysis are given in Appendix 8.3, and we publicly release the code as a python programme to quantify stirring in rebound N-body simulations³. We will use this numerical analysis to assess stirring in all simulations in this paper.

3 N-body simulations↩︎

To assess the effectiveness of projectile and resonant stirring, we run a large suite of N-body simulations, covering a wide range of system setups. Individual setups are randomly generated, but have similar qualitative characteristics as described in the following section; each comprises at least a star, a planet, and a massless, external debris disk. For each setup, we run two simulations; the first comprises only the star, planet, and disk, whilst the second also contains several massive debris bodies (projectiles). This approach lets us distinguish the stirring effects of the planet and projectiles. Once each simulation has finished, we apply our numerical analysis code (Section 2.2) to quantify the stirring level. In Section 3.1 we describe the setups of our N-body simulations, in Section 3.2 our simulation strategy, and in Section 3.3 we present the results of an example pair of simulations.

3.1 Simulation setup↩︎

Our aim is to investigate a scenario where a planet scatters the inner region of an external debris disk, whilst leaving the outer region unscattered. By naturally generating a scattered population, we can investigate the dynamical effect of scattered bodies (projectiles) on the rest of the disk. Having a sufficiently wide disk also lets us investigate the effect of the planet on distant debris, well beyond the scattering region. To this end, we set up our simulations with the debris disk initially extending from the planet’s orbit out to at least 5 Hill radii beyond this. We provide a setup diagram on Figure 1, and describe the setup and parameters in this section.

3.1.1 Integrator↩︎

Our N-body simulations are performed using the python version of rebound [40], with the mercurius hybrid-symplectic integrator [41]. This uses a symplectic Wisdom-Holman integrator (whfast) when particles are far apart, and switches to a high-order integrator (ias15) during close encounters. This lets us efficiently and accurately model both individual scattering events and long-term evolution.

3.1.2 General setup↩︎

Our randomised system setups are designed to cover a wide parameter space, but also be roughly representative of real debris-disk systems (e.g. those in Table 1 of [42]). Our parameter space is listed in Table 1. When initialising a system, we first draw the star mass \(m_*\) from the range \({0.08 \leq m_*\leq 2\; {\rm M_\odot}}\), using a uniform-logarithmic distribution. We then draw the planet mass \(m_{\rm plt}\) from the range \({1\; {\rm M_{Neptune}}\leq m_{\rm plt}< 10\; {\rm M_{Jup}}}\) (equivalent to \({17 \leq m_{\rm plt}/{\rm M_\oplus} < 3200}\)), again using a uniform-logarithmic distribution. The planet’s semimajor axis is uniformly drawn between 12.5 and \({100\; {\rm au}}\), and its eccentricity⁴ set to \({10^{-4}}\).

The massless debris disk is initialised with its inner edge, \(r_{\rm inner}\), equal to the planet’s semimajor axis \(a_{\rm plt}\). The disk’s outer edge is drawn to be at least 5 Hill radii (\(5r_{\rm Hill, plt}\)) beyond this, to ensure that some of the disk remains unscattered by the planet⁵. To avoid a bias towards generating wide disks, the outer-edge assignment is performed by first drawing the disk’s width-to-radius ratio, \({\Delta R / R}\). This is uniformly drawn in the range \({(\Delta R / R)_{\rm min}}\) to 1, where \({(\Delta R / R)_{\rm min}}\) is

\[\left(\frac{\Delta R}{R}\right)_{\rm min} = 2 \times \frac{a_{\rm plt}+ 5r_{\rm Hill, plt}- r_{\rm inner}}{a_{\rm plt}+ 5r_{\rm Hill, plt}+ r_{\rm inner}}. \label{eq:32width95ratio}\tag{4}\]

Once \(\Delta R/R\) is drawn, the disk’s outer edge, \(r_{\rm outer}\), can be calculated as

\[r_{\rm outer}= r_{\rm inner}\times \frac{2 + \Delta R/R}{2 - \Delta R/R}. \label{eq:32outer95disk95radius}\tag{5}\]

The debris disk comprises 2000 massless debris particles. Each particle has its semimajor axis \(a_{\rm deb}\) drawn such that the disk surface-density profile goes as \(a_{\rm deb}^{-3/2}\), like the Minimum-Mass Solar Nebula ([48], [49]). Each particle has an individual eccentricity \(e_{\rm deb}\) uniformly drawn in the range \({0 \leq e_{\rm deb}< 10^{-4}}\), and an inclination \(i_{\rm deb}\) relative to the planet’s orbital plane uniformly drawn in the range \({0 \leq i_{\rm deb}\leq 10^{-4} {\rm \; radians}}\). Each particle’s longitude of ascending node, argument of pericentre, and mean anomaly are individually drawn from uniform distributions between \(0\) and \({2\pi {\rm \; radians}}\).

The above setup describes simulations without additional, massive projectiles, such that the planet is the only body perturbing the debris disk. For these simulations, the setup is almost complete. We set the whfast fixed timestep to be \({1\; {\rm per \; cent}}\) of the planet’s period, and the minimum value of the ias15 variable timestep to be \({10^{-6}}\) times the planet’s period. During the integration, any body is removed from the simulation if it reaches a threshold distance of \({10^4\; {\rm au}}\).

We run each simulation until 10 diffusion timescales is reached, where the diffusion timescale \(t_{\rm diff}\) characterises the scattering timescale. This ensures that scattering is essentially complete by the end of the simulation, as we will later demonstrate in Section 4.3.1. Formally, one diffusion timescale quantifies how long it takes for the energy of a particle undergoing scattering to change by order of itself, and is given by

\[t_{\rm diff}\approx 0.01 T_{\rm plt}\left(\frac{a_{\rm plt}}{a}\right)^{\frac{1}{2}} \left(\frac{m_{\rm plt}}{m_*}\right)^{\rm -2}, \label{eq:32diffusion95time}\tag{6}\]

where \(T_{\rm plt}\) is the planet’s orbital period, and \(a\) is the semimajor axis of the scattered body ([50]; Equation 18 in [51]). To determine the end time of our simulations, we evaluate \(t_{\rm diff}\) at the planet’s location, i.e. at \({a=a_{\rm plt}}\).

3.1.3 Projectile simulations↩︎

Half of our simulations are set up as in Section 3.1.2. For the other half, we also include massive debris bodies (projectiles). These bodies typically represent dwarf planets, and are initialised such that they are very likely to get scattered onto high-eccentricity orbits by the planet. In the projectile simulations, the star, planet, massless disk, and integrator parameters are again set up exactly as in Section 3.1.2. We then initialise 10 projectiles, each with the same mass \(m_{\rm proj}\), where for each simulation \(m_{\rm proj}\) is drawn from a uniform-logarithmic distribution between a Pluto mass (\(2.2\times10^{-3}\; {\rm M_\oplus}\)) and \({10\; {\rm per \; cent}}\) of the planet’s mass. Each projectile’s semimajor axis \(a_{\rm proj}\) is individually drawn between the planet’s semimajor axis and 3 planetary Hill radii beyond this, such that the projectiles are initially located close enough to the planet to undergo scattering. Like massless debris, the projectile’s semimajor axes are drawn from a \(a_{\rm proj}^{-3/2}\) profile, and their individual eccentricities and inclinations are uniformly drawn between 0 to \({10^{-4}}\) and 0 to \({10^{-4} {\rm \; radians}}\) respectively. Also, each projectile’s longitude of ascending node, argument of pericentre, and mean anomaly are individually drawn from uniform distributions between \(0\) and \({2\pi {\rm \; radians}}\). Simulations are again run to 10 diffusion timescales, as calculated from the initial orbit of the planet.

3.2 Simulation strategy↩︎

For each system setup we run two simulations; one with projectiles, and one without. Both simulations use the same seed, such that the star, planet, and massless disk are identical between the two. We run 160 pairs of randomized simulations. We omit any simulations where projectiles significantly excite the planet eccentricity (above 0.01), because this results in the eccentric planet secularly stirring the disk, which is not the scenario we investigate.

After each N-body simulation finishes, we run the numerical stirring analysis from Section 2.2 to assess the state of the disk. The code classes each debris particle as either “scattered” (i.e. it was ejected, or its semimajor axis changed by at least \({20\; {\rm per \; cent}}\)), or “stirred” or “unstirred” (from the assessment described in Section 2.2 and Appendix 8.3). We bin the debris particles by their initial semimajor axis, to assess the degree of scattering and stirring at multiple locations across the disk.

3.3 Example simulation results↩︎

Figure 2 shows an example pair of simulations; the left panels show the result of the planet-only simulation, and the right panels the equivalent simulation with the addition of 10 massive projectiles.

In both simulations, the planet clears the majority of debris initially within 3 Hill radii of its orbit, with most debris beyond this surviving until the end of the simulation. In the planet-only simulation (left plots), the planet excites debris at the disk inner edge, as well as near the nominal location of several MMRs, but the majority of the surviving disk remains unstirred. Conversely, the addition of projectiles greatly increases the stirring level (right plots); debris eccentricities are excited to a much higher degree, and the surviving disk is almost completely stirred. Nine of the ten projectiles were ejected during the simulation.

The above results are for one specific system setup. Across all of our setups, we generally find that the addition of projectiles increases the excitation level of the disk, but the degree of stirring varies across systems. We also find that some planet-only simulations are able to stir very broad regions through MMRs alone. In Sections 4 and 5 we assess projectile and resonant stirring across all of our simulations, including predicting when these mechanisms are important.

4 Projectile stirring↩︎

Over half of our projectile simulations show a significant increase in debris-stirring level, like the simulation on the right panels of Figure 2. In this section we analyse the mechanism in greater detail. In Section 4.1 we determine the dynamical interaction by which projectiles stir debris, in Section 4.2 we examine how the number and mass of projectiles affects stirring, and in Section 4.3 we derive an analytic criterion for when we expect projectile stirring to occur, and verify it across our simulations.

4.1 Do projectiles excite debris via secular or scattering interactions?↩︎

We first determine the dynamical effect by which projectiles excite debris eccentricities. There are two plausible interactions: secular or scattering. In this section we analyse both, and demonstrate that secular interactions are the means by which projectiles stir debris. We will do this by comparing the timescales associated with each interaction.

4.1.1 Secular interactions between debris and projectiles↩︎

The first dynamical effect we consider is the secular effect of a high-eccentricity, disk-crossing projectile on other debris particles. Secular effects are long-term interactions, occurring on timescales much longer than orbital periods. They cause orbits to periodically oscillate in eccentricity, inclination, and orientation, whilst semimajor axes remain constant. Secular oscillations will be further discussed in Section 5.2, but for this section it is sufficient to note that secular interactions cause debris eccentricities to oscillate.

The period of these oscillations is the secular timescale, \(t_{\rm sec}\):

\[t_{\rm sec}\approx \begin{cases} 4 T_{\rm plt}\left( \frac{m_{\rm plt}}{m_*}\right)^{-1} \left(\frac{a_{\rm plt}}{a_{\rm deb}}\right)^{-\frac{5}{2}} \left[ b_{3/2}^{(1)} \left(\frac{a_{\rm plt}}{a_{\rm deb}}\right)\right]^{-1}, & \text{if } a_{\rm plt}< a_{\rm deb};\\ 4 T_{\rm plt}\left( \frac{m_{\rm plt}}{m_*}\right)^{-1} \left(\frac{a_{\rm deb}}{a_{\rm plt}}\right)^{-\frac{1}{2}} \left[ b_{3/2}^{(1)} \left(\frac{a_{\rm deb}}{a_{\rm plt}}\right)\right]^{-1}, & \text{else}, \end{cases} \label{eq:32secular95time}\tag{7}\]

where \(b_{\rm 3/2}^{(\rm 1)} (\alpha)\) is a Laplace coefficient:

\[b_{\rm s}^{(\rm j)} (\alpha) \equiv \frac{1}{\pi} \int_{\rm 0}^{\rm 2\pi} \frac{\cos({j\psi})}{(1 - 2\alpha\cos{\psi} + \alpha^{\rm 2})^{\rm s}} d\psi \label{eq:32laplace95coeff}\tag{8}\]

[31]. Since the secular time is dependent on particle semimajor axis, debris particles farther away from the massive body will take longer to undergo one full period of the secular interaction.

The above secular framework is based on Laplace–Lagrange theory, which is strictly only valid for low-eccentricity, non-crossing orbits. However, the rough principles still hold when applied to high-eccentricity, crossing orbits [51]–[55]. We will therefore use the Laplace–Lagrange secular timescale to quantify the secular effect of high-eccentricity projectiles on debris.

4.1.2 Scattering of debris by projectiles↩︎

The second dynamical effect we consider is scattering. Scattering occurs when two bodies make a close approach to each other; in the context of projectile stirring, we investigate whether projectiles on highly eccentric, disk-crossing orbits excite debris through close approaches to individual debris. Scattering is a short-term interaction, occurring on timescales much shorter than the orbital periods.

In Section 3.1.2 we argued that scattering is quantified by the diffusion timescale (Equation 6 ). Specifically, we argued that this timescale is related to how long it takes a larger body to excite or eject smaller bodies via scattering. Like the secular timescale (Equation 7 ), the diffusion timescale was originally derived assuming the massive body to be on a low-eccentricity orbit [50]; however, it is reasonably accurate for high-eccentricity orbits as well [51], [54]. We therefore use the diffusion timescale to quantify the scattering of debris by a high-eccentricity projectile.

4.1.3 Secular or scattering?↩︎

To assess whether projectiles would excite debris via secular or scattering interactions, we compare the secular and scattering timescales from Equations 6 and 7 . Figure 3 shows the ratio of these timescales for an interaction between a massive body and a test particle; for an eccentric, massive body, secular interactions will dominate if the secular timescale is shorter than the diffusion timescale, and vice-versa.

Figure 3 is general for any interaction between an eccentric, massive body and a test particle. We now apply this to projectile stirring; specifically, to determine the dynamical effect of a projectile (i.e. the massive body on Figure 3) on massless debris particles. The projectiles we consider are typically large planetesimals or dwarf planets, so the mass ratio \(m_{\rm proj}/m_*\) is always very small (\(\lesssim 10^{-3}\)). For this regime, Figure 3 shows that in almost all cases the diffusion timescale will be longer than the secular timescale, i.e. in the shaded region of the plot. The conclusion is that, in almost all cases, the secular effect of the projectile will have greater impact on debris than scattering; this means that projectiles stir debris via secular interactions, rather than by scattering.

4.2 Projectile-stirring efficiency as a function of mass and number of projectiles↩︎

In Section 3.3 we presented an example pair of simulations, which are good representations of the majority of our simulations. These showed that stirring is significantly increased when projectiles are included. In that example there were 10 equal-mass projectiles, with a total mass of \({0.88\; {\rm M_\oplus}}\). We now examine how the stirring level changes if we vary the mass and number of projectiles; the quantitative results will be specific to that example system, but they qualitatively hold for all of our simulations.

To better visualize the results of the increase in stirring done by projectiles for a particular system (instead of generalizing across any system), we can also plot the stirring percentage increase against the total-projectile mass to determine a critical mass for when stirring is likely to occur.

To quantify projectile stirring, we first define a metric. This will provide a single number that quantifies how much the stirring level increases due to the presence of projectiles. Since the planet on its own is able to stir some debris through MMRs (Section 5), we must account for this planet-only stirring when assessing the effect of projectiles. To proceed for a given system, we run our numerical stirring analysis on the planet-only simulation, and then separately on the projectile simulation. We then quantify the degree of projectile stirring as

\[\begin{align} &\text{projectile-stirring efficiency} = \\ &\frac{\text{particles~stirred~in~proj.~sim. - particles~stirred~in~planet~sim.}}{\text{particles~unstirred~in~planet~sim.}}, \end{align} \label{eq:32stir95percent95increase}\tag{9}\]

i.e. the increase in the number of debris particles stirred in the projectile simulation relative to the planet-only simulation, divided by the number that were not stirred in the planet-only simulation. For the purposes of this evaluation, we also class scattered and ejected particles as “stirred”, since they have been excited. According to this definition, if projectiles do not affect the stirring outcome then the projectile-stirring efficiency is \({0\; {\rm per \; cent}}\); conversely, if projectiles stir all bodies that were not stirred by the planet alone, then the efficiency is \({100\; {\rm per \; cent}}\). We apply this metric to quantify how the stirring level changes as a function of the mass and number of projectiles in the following sections.

4.2.1 Projectile mass↩︎

To assess the effect of projectile mass on the stirring level, we re-run the projectile simulation from Figure 2 with different projectile masses. In each simulation the 10 projectiles always have equal masses. For each projectile mass we test, we run multiple simulations; the planet, star, and disk parameters are kept constant, whilst the initial projectile orbits are randomised.

Figure 4 shows the results, where we quantify the projectile-stirring efficiency using Equation 9 . As could be expected, the higher the projectile mass, the more efficiently the disk is stirred. For this specific system, projectile stirring has a significant effect if the total-projectile mass is greater than \({\sim0.05\; {\rm M_\oplus}}\). Later, in Section 4.3.3, we will use an analytic argument to predict the minimum total-projectile mass required for projectile stirring to occur; that analysis predicts a minimum total-projectile mass of \({0.040\; {\rm M_\oplus}}\), which agrees with these simulations.

4.2.2 Number of projectiles↩︎

We arbitrarily considered 10 projectiles in our simulations, but we need to determine whether changing this number has any effect. We therefore re-run the example from Figure 2, but this time varying the number of projectiles. In each case the total mass of projectiles is kept at the original value, but is distributed between different numbers of equal-mass bodies. As in Section 4.2.1, we run multiple simulations for each number of projectiles tested; in each case the planet, star and disk parameters are kept constant, whilst the initial projectile orbits are randomised.

Figure 5 shows the results. The plot shows that, provided the number of projectiles is not too small, the number has minimal effect on the simulation outcome; for this specific setup, the outcome is similar regardless of the number of projectiles, provided it is at least three. The outcome is much more variable if just one projectile is initialised, because in that case the stirring level would be strongly linked to the stochastic evolution of the projectile as it is repeatedly scattered. Comparing Figure 5 to Figure 4, it appears that the total mass of projectiles is more important than the number of projectiles in setting the stirring outcome.

4.3 Analytical prediction of projectile stirring↩︎

We now aim to find a general analytical prediction for whether projectile stirring would be important in a debris-disk system.

Since projectile stirring is a secular effect (Section 4.1), we expect projectile stirring will occur if the projectiles can secularly affect debris before the projectiles are ejected by the planet. We first define the relevant ejection timescale in Section 4.3.1, then the relevant secular timescale in Section 4.3.2, before combining these to produce an analytic criterion in Section 4.3.3.

4.3.1 Scattering of projectiles by the planet↩︎

We first confirm that the diffusion timescale is the best means to quantify how long it takes projectiles to be scattered and ejected by the planet. We have so far argued that particles initialized within 3 Hill radii of a planet are likely to be scattered and ultimately ejected, and that this is characterised by the diffusion timescale; here we test both of these assumptions across all of our simulated systems.

To gain general insight into scattering, we consider our planet-only simulations. For each simulation, we examine all of the massless debris particles that were initialised within 3 Hill radii of the planet. We then plot the fraction of those particles that are ejected or scattered to high eccentricity, as a function of time. The results for all simulations are shown in Figure 6, where the time in each simulation is expressed in terms of the diffusion timescale (Equation 6 , evaluated at \({a=a_{\rm plt}}\)).

Figure 6 shows that the diffusion timescale is a good metric to quantify the scattering timescale. After 1 diffusion timescale, roughly \({70\; {\rm per \; cent}}\) of particles initialised on low-eccentricity orbits within 3 Hill radii of a planet have been scattered up to eccentricities of at least 0.6, and roughly \({50\; {\rm per \; cent}}\) have reached eccentricities of at least 0.9. Both of these percentages include particles ejected from the system. The scattered fractions are much smaller after just 0.1 diffusion timescales, being roughly 40 and \({20\; {\rm per \; cent}}\) respectively, whilst after 10 diffusion timescales they are both about \({90\; {\rm per \; cent}}\). Note that the ejection fraction may not necessarily reach \({100\; {\rm per \; cent}}\), because some particles initialised near the planet can occupy stable resonant orbits (e.g. particles in the 1:1 MMR in Figure 2). This exercise demonstrates that the diffusion time well-quantifies the timescale for scattering by one planet over a broad range of system parameters, and so should be a good estimate for projectile lifetimes.

4.3.2 Defining the projectile secular timescale↩︎

We now define the relevant secular timescale for our prediction of projectile-stirring efficiency. The timescale over which a projectile would secularly perturb a debris particle is given by Equation 7 (replacing the planet parameters with those of the projectile); however, evaluating this timescale is not straightforward, because it strongly depends on the semimajor axes of the debris and projectile. The projectile semimajor axis varies greatly throughout a simulation as it is repeatedly scattered by the planet, and the debris spans a range of semimajor axes. To proceed, we must choose “typical” values for these to use in our analytic prediction.

First, we choose to consider debris at the outer edge of the initial disk, i.e. setting \({a_{\rm deb}\equiv a_{\rm outer}}\) in Equation 7 (where \(a_{\rm outer}\) is the semimajor axis of the outermost debris body). This is because we are interested in stirring the entire disk and, since processes at the outer edge are typically slower than at the inner edge, if the outer edge is stirred then the inner edge would probably be stirred too. Second, for the projectile semimajor axis we choose a value between the disk inner and outer edges, because the projectiles start at the inner edge and are then scattered outwards. We arbitrarily choose to evaluate the secular time assuming a projectile semimajor axis of 95% of the distance from the initial inner edge to the outer edge (i.e. a projectile semimajor axis of \({a^\prime \equiv 0.05a_{\rm plt}+ 0.95a_{\rm outer}}\)); we will later show that this value lets us accurately predict when projectile stirring will occur. Finally, we must decide what mass to use when calculating the secular timescale; we choose to use the total mass of all projectiles (rather than an individual projectile) when evaluating Equation 7 , because we find that the stirring level depends more on the total-projectile mass than on the number of projectiles (Figures 4 and 5). So to summarise, when evaluating the secular timescale (Equation 7 ) for our analytic prediction of projectile stirring, we set \(m_{\rm plt}\) to be the total-projectile mass, \({a_{\rm deb}\equiv a_{\rm outer}}\), and consider a projectile semimajor axis of \({a^\prime \equiv 0.05a_{\rm plt}+ 0.95a_{\rm outer}}\).

4.3.3 Analytical criterion for projectile stirring↩︎

Using the results of the previous two sections, we now produce an analytic criterion for when projectile stirring should be important. We argued that projectiles could stir a disk if the secular timescale over which projectiles excite debris, \(t_{\rm sec,proj}\), is shorter than the diffusion timescale over which the planet would eject projectiles, \(t_{\rm diff,plt}\). In Figure 7 we plot the projectile-stirring efficiency measured from each of our \(N\)-body simulations, as a function of the ratio \({t_{\rm sec,proj}/t_{\rm diff,plt}}\) (evaluated for each simulation using the assumptions in Sections 4.3.1 and 4.3.2). The figure shows that the theoretical criterion \({t_{\rm sec,proj}< t_{\rm diff,plt}}\) is a good indicator of whether projectile stirring will occur.

We can substitute the expressions for the two timescales in the criterion \({t_{\rm sec,proj}< t_{\rm diff,plt}}\). The result is that projectile stirring is predicted to occur if the total mass of projectiles satisfies

\[\begin{align} & \left(\frac{m_{\rm proj, tot}}{{\; {\rm M_\oplus}}}\right) \gtrsim 120 \left(\frac{a^\prime}{\rm au}\right)^{-1} \left(\frac{a_{\rm plt}}{\rm au}\right)^{-3/2} \left(\frac{a_{\rm outer}}{\rm au}\right)^{5/2} \\ & \times \left(\frac{m_{\rm plt}}{\; {\rm M_{Jup}}}\right)^{2} \left(\frac{m_*}{{\; {\rm M_\odot}}}\right)^{-1} \left[ b_{\rm 3/2}^{(1)} \left(\frac{a^\prime}{a_{\rm outer}}\right) \right]^{-1}, \end{align} \label{eq:32theoMinProjMassToStir}\tag{10}\]

where \({a^\prime \equiv 0.05~a_{\rm plt}+ 0.95~a_{\rm outer}}\), and the Laplace coefficient \(b_{\rm 3/2}^{(1)}(\alpha)\) is given by Equation 8 (a code to evaluate Laplace coefficients is available online⁶).

For the simulation in the right of Figure 2, Equation 10 predicts a total-projectile mass of at least \({0.040\; {\rm M_\oplus}}\) is needed to stir the disk. This is consistent with the simulation; the total-projectile mass in that simulation is much higher (\({0.88\; {\rm M_\oplus}}\)), and indeed the simulated disk is stirred. This theoretical prediction is also shown in Figure 4, and it agrees with projectile mass at which those simulations transition from projectile stirring being inefficient to efficient.

To verify Equation 10 across all our simulations, in Figure 8 we plot our measured projectile-stirring efficiencies, against the ratio of the actual projectile mass to the Equation 10 prediction. The plot shows that, in simulations where the total-projectile mass is less than the theoretical requirement from Equation 10 , projectile stirring is absent or minimal; conversely, projectile stirring is significant in simulations where the total-projectile mass exceeds the theoretical requirement. Hence Equation 10 is a good indicator of the total-projectile mass required to stir a disk, and can be used to predict the outcome of projectile stirring.

Equation 10 lets us identify the properties of systems where projectile stirring is more likely to occur. Projectile mass is a key factor; larger projectiles are better able to stir debris. Projectile stirring is also more likely if the mass of the scattering planet is small, because smaller planets take longer to eject projectiles, and hence allow them more time to stir. However, if the planet mass is too small, then it may take too long to scatter projectiles in the first place; specifically, if the planet’s diffusion time (Equation 6 ) is much longer than the system age, then it is unlikely that the planet will have yet scattered projectiles onto high-eccentricity orbits. Projectile stirring is also more efficient around high-mass stars, because this increases the ratio of the diffusion and secular timescales. Finally, projectile stirring is most efficient in small, narrow debris disks.

5 Resonant stirring↩︎

In the previous sections, we analysed projectile stirring by comparing the debris-excitation level in N-body simulations that contained a planet, projectiles and debris to equivalent simulations without projectiles. However, we noticed that many simulations without projectiles were stirred to a much greater degree than would be expected from secular-planet-stirring, despite the planet being the only possible driver of stirring. In these simulations, it transpired that stirring was being performed by the planet’s mean-motion resonances, which excite material across a much larger region of the disk than expected. In this section we discuss “resonant stirring” as a means to increase planet-stirring efficiency. We plan to examine this effect in greater detail in a future publication, including analytic predictions of resonant-stirring efficiency; for the current paper, we demonstrate the effect and discuss its prevalence in our N-body simulations.

In Section 5.1 we present example N-body simulations of resonant stirring, and in Section 5.2 we demonstrate that mean-motion resonances are the stirring mechanism. In Section 5.3 we analyse resonant stirring across all of our simulations, and demonstrate that it seems viable for planet masses above \({\sim0.5\; {\rm M_{Jup}}}\).

5.1 Example resonant-stirring simulations↩︎

Figure 9 shows an example N-body simulation where resonant stirring occurs. The system comprises a \({1.3\; {\rm M_\odot}}\) star, a \({1.9\; {\rm M_{Jup}}}\) planet at \({60\; {\rm au}}\), and a debris disk initially extending from the planet out to \({110\; {\rm au}}\). By the end of the simulation (10 diffusion timescales), the planet has ejected non-resonant debris originating within 3 Hill radii of its orbit, but debris beyond this is stirred across almost the entire disk, despite no projectiles being present. Specifically, the majority of debris particles are excited to eccentricities above \(10^{-2}\) to \(10^{-1}\), compared to their initial eccentricities of less than \(10^{-4}\). The stirring level is linked to the planet’s mass; Figure 10 shows a simulation with an identical setup but with the planet mass reduced by a factor of 10 (to \({0.19\; {\rm M_{Jup}}}\)), and the stirring level is greatly reduced.

5.2 Identifying mean-motion resonances as the stirring mechanism↩︎

We now demonstrate that MMRs are the stirring mechanism in the planet-only simulation in Figure 9. In such simulations comprising a star, a planet and massless debris, there are only three interactions that could excite debris: scattering interactions, secular interactions, and mean-motion resonances. We assess each of these here.

Planet-debris scattering can be quickly discounted as the stirring mechanism. This is because scattering exchanges both energy and angular momentum, so it would change debris semimajor axes as well as eccentricities; however, the semimajor axes of surviving debris in the Figure 9 simulation remain essentially unchanged (as denoted by particle colours). In addition, the outer edge of the debris disk is at \({110\; {\rm au}}\), which is 9.7 Hill radii beyond the planet; this is far outside the expected scattering region.

The second potential stirring mechanism is a secular interaction between the planet and debris; this is the “traditional” planet-stirring mechanism considered in the literature [12]. Unlike scattering, secular interactions conserve energy, so would be consistent with the debris semimajor axes remaining unchanged in the simulation. However, secular stirring can also be discounted, because the planet is not eccentric enough. For debris with initially low eccentricity, an eccentric planet would cause the debris eccentricity to oscillate between roughly its initial value and twice the forcing eccentricity \(e_{\rm forced}\), where

\[2 e_{\rm forced}(a_{\rm deb}) = 2\times \frac{5}{4}\frac{a_{\rm plt}}{a_{\rm deb}}e_{\rm plt}. \label{eq:32forced95ecc}\tag{11}\]

The forcing eccentricity, and hence the maximum eccentricity attained by secular debris, is set by the planet eccentricity and the planet-debris-semimajor-axis ratio. For the simulation in Figure 9, the the maximum debris eccentricity expected at the disk outer edge from secular interactions is therefore

\[ 2 e_{\rm forced}(110\; {\rm au}) = 2 \times \frac{5}{4} \frac{60\; {\rm au}}{110\; {\rm au}} \times 10^{-4} = 1.4\times10^{-4}, \label{eq:32forced95ecc95calculated}\tag{12}\]

which is far smaller than the eccentricities of \({4\times10^{-3}}\) at the outer edge of the simulated disk. Even nearer to the planet, where the forcing eccentricity would be higher, debris is still much more eccentric than would be expected from secular interactions; the dot-dash line on Figure 9 shows \({2e_{\rm forced}}\) across the disk, which is far smaller than the actual debris-excitation level. Hence secular interactions can also be discounted as the stirring mechanism.

Mean-motion resonances are therefore the only possible mechanism that could have stirred the disk. MMR interactions cause debris to oscillate in eccentricity and semimajor axis over a narrow range, so are again consistent with debris semimajor axes being essentially unchanged in the simulation. From Figure 9 it is clear that MMRs are operating; debris has well-defined eccentricity peaks just outside the nominal 3:2, 5:3 and 2:1 MMRs, as expected from resonant interactions. Importantly, these peaks are not infinitely narrow in semimajor-axis space; they have width, and the maximum debris eccentricity smoothly decreases at semimajor axes on either side of the nominal MMR locations. This width is a fundamental property of MMRs, and the shapes of the debris profiles in eccentricity-semimajor-axis space around e.g. the nominal 2:1 MMR on Figure 9 well match those expected from MMRs [31]. We therefore conclude that mean-motion resonances are responsible for stirring the disk in this case, and importantly, that their non-zero widths mean that MMRs can stir debris over a much broader region than just their nominal locations.

5.3 Resonant stirring across our simulations↩︎

Resonant stirring occurs in at least one third of our planet-only simulations (i.e. those without projectiles). By this, we mean that the debris exhibits the characteristic eccentricity-semimajor-axis profile seen in Figure 9, and a substantial region of the disk is stirred (not only at nominal MMR locations, as in Figure 10, but between them too). In this section we analyse resonant stirring across our planet-only simulations; this effect is probably also present in our projectile simulations, but we omit those from this analysis to avoid confusing the effects of resonant and projectile stirring. Since MMR theory is more complicated than the relatively simple analyses we performed for projectile stirring (Section 4), we intend to perform a detailed resonant-stirring analysis in a dedicated future paper; in this section, we instead aim to identify the key parameters that set the importance of resonant stirring.

For massless, resonant debris, the main parameters predicted to set the debris eccentricity and MMR width are the planet-to-star-mass ratio \(m_{\rm plt}/m_*\), and the terms associated with the specific MMR (i.e. \(p\) and \(q\) in the MMR notation \({p+q:p}\); [31], [37], [56]). Typically, higher planet-to-star-mass ratios and lower-order resonances result in higher debris eccentricities and broader MMRs. We therefore expect eccentricity excitation to be highest for high-mass planets, and for disks near the strong 3:2 and/or 2:1 MMRs. This is consistent with our simulations; Figures 9 and 10 show that debris is most excited around these two MMRs, and the degree of excitation is higher for the higher-mass planet.

However, \(p\), \(q\) and \(m_{\rm plt}/m_*\) alone are not sufficient to determine whether resonant stirring would occur. The reason is that the debris eccentricity required for stirring also depends on other system parameters; specifically, Equations 1 to 3 show that the minimum eccentricity for stirring depends on the debris material, disk location and star mass. Hence the resonant-stirring efficiency may differ in systems with identical \(m_{\rm plt}/m_*\) and relative debris locations, but different star masses and/or planet locations. We therefore hypothesise that the main parameters setting the efficiency of resonant stirring are the star mass, planet mass, locations of the disk edges, and debris composition.

We now test resonant-stirring efficiency against these parameters, to establish which are most important. We first define a rough proxy for resonant-stirring efficiency in our planet-only simulations. Since we define resonant stirring to have occurred if debris is stirred both at and between nominal MMR locations, and since the 3:2 and 2:1 MMRs are typically the strongest, we choose to take the stirring level at the midpoint between the nominal 3:2 and 2:1 locations as a proxy for the degree of resonant stirring in a planet-only simulation. For this analysis we exclude any simulations with massive projectiles, or where the midpoint between the 3:2 and 2:1 MMRs lies beyond the disk outer edge, or where the midpoint lies within 5 Hill radii of the planet (to exclude cases where debris is excited by scattering rather than MMRs). For each of the included simulations, we run our numerical stirring analysis (Section 2.2), then bin each debris particle into one of 20 equal-width bins depending on its initial semimajor axis. Following this process, if there are more than 10 debris particles in the bin coinciding with the midpoint between the 3:2 and 2:1 MMRs, then we use the fraction of stirred particles in this bin as our proxy for the degree of projectile stirring in the simulation.

Figure 11 shows the results of this analysis as a function of planet mass. Of the parameters listed above, planet mass appears to be the most important in setting the efficiency of resonant stirring; for our tested parameter space (Table 1), it appears that resonant stirring is significant for planet masses above \({\sim0.5\; {\rm M_{Jup}}}\). This critical planet mass appears to be relatively independent of star mass and disk-edge locations. Other parameters do not seem to affect resonant stirring as clearly; we show the resonant-stirring proxy as a function of star mass, planet-to-star-mass ratio, and disk-edge locations in Figure 16, but these individual parameters do not have as tight relations to resonant-stirring efficiency as planet mass does.

In summary, broad mean-motion resonances from an internal planet can efficiently stir debris disks, and we show that this effect can become significant for planet masses above \({\sim0.5\; {\rm M_{Jup}}}\). We identify this mechanism as a potential means to significantly increase planet-stirring efficiency, especially for low-eccentricity planets. However, we do not conduct a detailed numerical or theoretical exploration of resonant stirring here; that analysis is planned for a future paper.

6 Discussion↩︎

We demonstrated that the efficiency of planet-stirring can be increased through two additional mechanisms: projectile stirring and resonant stirring. In this section we discuss how projectile and resonant stirring fit in with the other stirring mechanisms (Section 6.1), the viability of projectile stirring (Section 6.2), the limitations of our analyses (Section 6.3), and what future studies could focus on (Section 6.4).

6.1 Implications of projectile and resonant stirring↩︎

The nature of debris stirring is unknown. Several mechanisms have been proposed, but these have problems; secular-planet-stirring requires the planet to have high eccentricity ([12]), which is not necessarily consistent with known Solar-System and extrasolar planets. It may also run into difficulty when disk mass is considered, since disk self-gravity may resist secular-planet-stirring (Löhne et al., in prep.). The most widely considered alternative, self-stirring due to debris-disk self gravity, also has significant problems; it cannot be occurring in some debris disks unless their size distributions are very different to conventional debris theory ([18]–[20]).

We showed that planet-stirring could be significantly more efficient than previously thought, due to the additional effects of projectiles and mean-motion resonances. In particular, these effects would occur even for non-eccentric planets. They add to the arsenal of mechanisms by which planets could stir debris; as well as secular stirring by eccentric planets [12], these include stirring by MMR sweeping as planets migrate [47], and interplays between planet-debris interactions and disk self gravity (‘mixed stirring’; [38], [39]). These mechanisms open the potential for planet-stirring to be more plausible than previously thought, with the possibility that planets are the main debris-stirring mechanism in nature.

Planet-stirring could potentially be tested soon; the James Webb Space Telescope (JWST) is searching for planets near the inner edges of debris disks, and detections (or non detections) would help assess the viability and prevalence of planet-stirring. Our paper constrains the planet masses required for projectile and/or resonant stirring to occur, which could be used to interpret JWST planet searches. If no sufficiently massive planets were found, then this could point towards the prevalence of less-considered stirring mechanisms; examples include pre-stirring and flyby stirring.

Finally, there is a potential implication of projectile stirring for the Solar System. Applying Equation 10 to the Kuiper Belt and Neptune, we see that projectile stirring of the Belt is possible if projectiles had a total mass of at least \({4 \times 10^{-4}\; {\rm M_\oplus}}\). This is just \({10 \; {\rm per \; cent}}\) of the mass of the dwarf planet Eris which, as noted in Section 1, has a highly eccentric, Kuiper-Belt crossing orbit, consistent with historical scattering by Neptune. This raises the possibility that the Kuiper Belt may have been dynamically excited by projectile stirring. Neptune’s diffusion timescale is \({0.6\; {\rm Gyr}}\) (Equation 6 ), which is also an upper limit on the secular timescale of projectiles with a total mass of at least \({4 \times 10^{-4}\; {\rm M_\oplus}}\). Hence Neptune has had enough time to scatter projectiles, and those projectiles have had enough time to stir the Kuiper Belt, within the age of the Solar System. Similar ideas have been suggested before; for example, secular perturbations by an Earth-mass rogue planet could have excited Kuiper-Belt Objects (KBOs) and generated the detached population [33], [57]. Our results suggest that Kuiper-Belt stirring could have been performed by smaller dwarf planets, like Eris, which were scattered by Neptune. However, we are cautious with this conclusion, because we have not performed specific Solar-System modelling; in particular, we do not simulate Neptune’s migration, nor include the Kuiper-Belt mass, nor attempt to reproduce the specific populations of resonant and detached KBOs. Nonetheless, the existence of high-eccentricity Eris, combined with the dynamical excitation of many KBOs, could imply that some form of projectile stirring has occurred in our own Solar System.

6.2 Viability of projectile and resonant stirring↩︎

There are several requirements of a planetary system in order for projectile stirring to occur. First, a planet must be present, with a supply of nearby, massive bodies to scatter. This appears plausible, based on our knowledge of the Solar System and extrasolar systems; Neptune scattered a large quantity of massive debris from the inner edge of the Kuiper Belt [32], and many extrasolar debris disks have sharp inner edges and/or extended radial profiles indicative of planetary sculpting and/or scattered populations [34]–[37], [58]. This implies that, in at least some cases, planets are (or were) within a few Hill radii of debris populations, which they then scattered. For our simulations, we imposed a scattered population by initializing the planet at the disk inner edge, but it is not a requirement of projectile stirring that planets formed in such locations; for example, a planet forming away from a debris disk and then migrating towards it could also result in projectile stirring. Since we have evidence of such planetary migration in both the Solar System [59]–[61] and possibly extrasolar debris-disk systems [47], [62], post-formation migration is a plausible means to initiate projectile stirring. Hence the first requirement for projectile stirring, a supply of massive bodies close to a planet, appears plausible.

The second requirement is that the total mass of projectiles must be sufficiently high. Equation 10 shows that, for a debris disk centred on \({100\; {\rm au}}\) with width \({50\; {\rm au}}\), the required combined mass of projectiles is \({\sim 0.16 (m_{\rm plt}/\; {\rm M_{Jup}})^2 (m_*/\; {\rm M_\odot})^{-1}}\). So if, for example, a Jupiter-mass planet were present at the disk inner edge, and the star were A-type (\({\sim 2 \; {\rm M_\odot}}\)), then the required total mass of projectiles would be \({\sim30\; {\rm M_\oplus}}\). This mass is plausible; the early Kuiper Belt is thought to have had a mass of \({\sim10\; {\rm M_\oplus}}\) [63], [64], and extrasolar debris disks could have masses up to \(100-1000~{\; {\rm M_\oplus}}\) [18]. Hence sufficiently massive projectile reservoirs should be available. Taken together, the two main requirements for projectile stirring seem likely to be satisfied in at least some systems; based on this, we argue that projectile stirring is a viable means to stir debris disks.

For resonant stirring, the requirements should be even easier to satisfy. The minimum-required planetary mass is just \({\sim0.5\; {\rm M_{Jup}}}\), and such planets are known to exist in our Solar System and beyond. Furthermore, the planet need not even be located close to the debris disk; in the example on Figure 9, debris is stirred all the way out to the nominal location of the 2:1 MMR. In such a system, debris with semimajor axis \(a\) would be stirred if \({a \leq 1.6a_{\rm plt}}\), so even well-separated planets could stir debris disks. Whilst Jupiter and Saturn may hence be too distant to stir the Kuiper Belt, there are known extrasolar systems with massive planets sufficiently close to debris-disk inner edges (e.g. \({{\rm HR} \;8799}\); [58], [65]). Hence the requirement of a sufficiently massive planet close enough to a debris disk seems plausible from currently known planets, and there may also be additional planets even closer to disks that are yet to be discovered. In addition, were a planet to migrate, then MMR sweeping would increase both the trapping efficiency and eccentricity of debris [47], [62], so migration would further enhance resonant stirring. In summary, both resonant and projectile stirring appear to be viable mechanisms to stir debris disks, with the required conditions thought to be satisfied in at least some systems.

6.3 Limitations of our analyses↩︎

We now discuss the limitations of our numerical and analytic analyses, and their potential impact on our results.

6.3.1 Debris-disk masses↩︎

We modelled most debris as massless; only a few projectiles were assigned masses in our N-body simulations. This meant we neglected disk self-gravity, particularly in simulations without massive projectiles, or in the outer regions of the disk (away from projectiles). There are two main effects of this omission. First, our disks are unable to undergo self-stirring, so the stirring level may be underestimated in some cases (especially for broad disks in simulations without projectiles). Second, the disks are unable to resist perturbations from the planet and projectiles; Löhne et al. (in prep.) show that massive debris disks can resist secular-planet-stirring, and it is therefore possible that they would resist projectile and resonant stirring too. Hence the stirring level could be overestimated in some of our simulations.

However, we justify omitting debris-disk masses for several reasons. First, debris-disk masses are unknown, and the plausible, theoretically justified estimates of disk mass span several orders of magnitude [18]. Therefore, if we were to include disk mass, it is unclear what mass should actually be assumed. Second, we sought to study the effects of projectile and resonant stirring, and this was done by isolating those mechanisms; including disk mass would have led to self-stirring as well, which would have made the effects of projectile and resonant stirring much harder to disentangle. Third, in projectile simulations, it could be argued that disk mass is actually included; debris-disk masses are thought to be dominated by the largest bodies, so the massive projectiles could be thought of as the largest bodies in the disk. This means that we assume the disk mass is initially concentrated at the disk inner edge (where the projectiles are initialised), but this may be a valid approach; since planetesimal-formation efficiency probably decreases with stellar distance (e.g. [14]), and large bodies may also migrate inwards through planetesimal scattering [47], it could be expected that the mass of a broad debris disk is concentrated at its inner edge. For these reasons, we argue that we were justified in omitting disk mass. Doing so significantly sped up our N-body simulations, letting us run many more simulations and explore a much larger region of parameter space.

6.3.2 Assumption of an initially unstable debris disk↩︎

Our N-body simulations were initialised in an unstable configuration, with debris extending all the way to the planet’s orbit. This was deliberate, because we sought to generate a scattered-projectile population. The main alternative for investigating projectile stirring would be to assume a more-distant disk, and insert projectiles on eccentric orbits that cross both the planet and disk. However, the latter approach would be hard to justify, for two reasons. First, the projectile orbits would be artificial and less realistic; the distribution of projectile orbits in semimajor-axis, eccentricity, inclination and orientation space would almost certainly have been biased had we initialised them manually, a bias we avoided by letting the planet naturally generate this scattered population. Second, if we had initialised projectiles on high-eccentricity, disk-crossing orbits, then we would have omitted any system evolution that occurred as the projectiles’ eccentricities increased from the low-eccentricity orbits they presumably formed on. By initialising projectiles near the planet, we captured the evolution of the disk and planet as projectile eccentricities grew.

To some degree, our setup is physically justified. As noted in Sections 1 and 6.2, we have considerable evidence that debris scattering occurs in nature. Whilst the specific system architectures and histories that lead to this scattering are debatable, the outcome is the same; somehow, significant quantities of debris have come within a few Hill radii of a planet, and been scattered. We choose to model this in the simplest-possible way, by initialising debris near the planet. In reality, it may be that planets form away from debris disks and then migrate towards them (like Neptune in the Nice Model; [59]–[61]), but we chose not to model migration because it requires significant assumptions that would have biased our results. Regardless, if migration were rapid (like Neptune in the Nice Model), then the planet would quickly encounter debris, and the result would roughly resemble our simulation setup. Alternatively, debris could form in the disk and diffuse through self-scattering towards the planet, again leading to scattering. In this case, the projectile-stirring outcome would be similar to our simulations; once massive debris encountered the planet, it would be scattered onto eccentric orbits and projectile stirring would begin.

By initialising particles close to the planet, we created a scattered population of massless debris. Before being ejected, this population acquired high-eccentricity orbits that crossed the unperturbed outer regions of the disk. However, this was accounted for in our stirring analyses, so as not to overestimate the stirring level in the outer regions. As explained in Section 2.2 and Appendix 8.3, our numerical stirring analysis only compared pairs of debris particles with similar semimajor axes, and whose semimajor axes remained close to their initial values. This ensured that high-eccentricity, massless debris scattered from the inner edge did not contaminate the measured stirring level further out.

6.3.3 Omission of inclination in the numerical stirring analyses↩︎

With our numerical stirring analysis, we aim to properly account for different debris eccentricities and orientations when assessing stirring levels. This produces a more-accurate picture than the simple, minimum-eccentricity value from Equations 1 3 , and works by assessing collision speeds at orbit-intersection points. However, we implemented the numerical analysis in just two dimensions; mutual inclinations are not considered when calculating whether debris orbits intersect. We omitted inclinations because we would otherwise need additional, non-physical variables in the code; for example, to account for finite numbers of debris particles in our N-body simulations, we would have had to include numerical impact parameters around each orbit when calculating three-dimensional intersection points. The measured stirring levels could be very sensitive to such artificial parameters, so we chose to omit them by not considering mutual inclinations.

This decision should be valid, because the disks are initialised to be very thin, and the non-scattered regions (where stirring is measured) stay relatively thin throughout the simulations. For example, debris at the outer edge of the disk on the right panels of Figure 2 gets excited to a maximum inclination of just \({0.6^\circ}\). For discs comprising a fixed number of bodies, an increased inclination spread would reduce the collision probability, because that probability is proportional to the geometric thickness of the disk. This effect could lead to our stirring percentages being slightly overestimated. However, a counterpoint is that mutual inclinations would also increase collision speeds, meaning that our stirring levels could alternatively be slightly underestimated. Either way, for the thin discs we consider, the omission of inclination is unlikely to strongly influence our stirring analyses. Also, since the code only considers pairs of debris bodies with small mutual inclinations, our stirring results are not affected by highly inclined debris that do not intersect the disk.

6.3.4 Assumed projectile semimajor axis in secular calculations↩︎

In Sections 4.3.2 and 4.3.3 we predicted the outcome of projectile stirring by considering the secular timescale for projectiles to stir a disk. This analytic calculation required a value for the projectile semimajor axis. However, a projectile’s semimajor axis constantly varies throughout an N-body simulation, as it is repeatedly scattered by the planet. We therefore had to assume an “effective” semimajor axis in our calculations, which accounted for all projectiles at all times. Our chosen value of \({0.05a_{\rm plt}+ 0.95a_{\rm outer}}\) is justified because the resulting equations well reproduce the simulations outcomes (Figures 7 and 8), but assuming a single value across all of our simulations could mask dynamical subtleties. However, since our chosen value appears to result in accurate predictions across our broad parameter space, we do not believe this to be a significant issue.

6.4 Future directions↩︎

We conducted a preliminary analysis of projectile and resonant stirring. In this section we outline additional analyses that should be performed in the future, to better understand these mechanisms, their prevalence in nature, and their effect on planet-stirring efficiency.

First, whilst we produced an analytic prediction for projectile-stirring efficiency (Equation 10 ), we did not produce an equivalent for resonant stirring. This is because the required MMR theory is more complicated, and such an analysis lies beyond the scope of this project. We plan to produce a detailed theoretical analysis of resonant stirring in a future paper; specifically, we intend to take analytic prescriptions for the shapes of external MMRs in eccentricity-semimajor-axis space, and combine these with theoretical stirring requirements (e.g. Equation 1 ) to derive analytic predictions of resonant-stirring efficiency. These predictions would then be verified against a broad suit of N-body simulations, and could be applied to observed debris-disk systems.

Second, future works should analyse the interplay between projectile, resonant, and self-stirring. This would involve similar analyses to those conducted here, only with massive debris disks. This would let us assess whether disk self gravity increases or reduces stirring efficiency. Similar studies have been conducted by [38], [39]; comparing results such as these to dynamical predictions for projectile and resonant stirring would help us understand potential connections between different stirring mechanisms.

Future investigations could also include planet migration, which has the potential to both increase and decrease the effectiveness of projectile and resonant stirring. For example, migration could enhance projectile stirring by causing the planet to encounter and scatter more projectiles as it moves through a disk, but conversely it could also hinder stirring because potential projectiles could get trapped in sweeping mean-motion resonances and hence become protected from scattering ([47]). Since planet migration is a plausible means to bring planets into contact with projectile reservoirs, the effect of migration on these stirring mechanisms should be investigated.

Finally, projectile and resonant stirring should be applied to real debris disks, to assess their viability in such systems, and to make specific predictions for unseen stirring planets. This could be especially useful for debris disks where the current, traditional stirring models do not fit well; in such cases, the new mechanisms of projectile and resonant stirring could offer insights into how debris disks are stirred.

7 Conclusions↩︎

In this paper, we examined two new mechanisms that could increase the effect of planet-stirring. The first mechanism is projectile stirring, where massive debris bodies that have been scattered onto eccentric orbits by a planet can stir the remaining debris disk. The second is resonant stirring, where a planet (without any massive projectiles) stirs a debris disk through broad mean-motion resonances alone.

To determine the viability of these stirring mechanisms, we ran N-body simulations of debris-disk systems across a broad parameter space. We then analysed these simulations using a bespoke python code, which numerically quantifies the degree of stirring in debris-disk systems simulated with rebound. We have made this stirring-analysis code publicly available.

We demonstrated that both projectile and resonant stirring have the potential to significantly increase planet-stirring efficiency. In particular, neither requires the planet’s orbit to be eccentric, in contrast to the traditional, secular model of planet-stirring.

For projectile stirring, we identified an analytic condition (Equation 10 ) to predict when this type of stirring is likely to occur. For resonant stirring, we demonstrated that a planet’s broad mean-motion resonances can substantially excite debris eccentricity across a wide region of a disk, rather than only near the nominal locations of MMRs. We showed that this stirring mechanism is viable for planet masses above \({\sim0.5\; {\rm M_{Jup}}}\).

We conclude that both massive projectiles scattered by a planet, and broad mean-motion resonances, can significantly increase debris eccentricity and hence greatly increase planet-stirring efficiency. Previous models of stirring, while applicable in some cases, had difficulty describing how some observed disks came to be stirred; the two new mechanisms we investigate could help explain how stirring is achieved in debris disks.

Acknowledgements↩︎

We would like to thank Marco Muñoz-Gutiérrez for his thorough review, which substantially improved the quality of this paper. We also thank him, Jonathan Marshall, and Antonio Peimbert for sharing their manuscript before its publication, and for useful discussions regarding alternate forms of stirring. TDP and AVK were supported by Deutsche Forschungsgemeinschaft (DFG) grants Kr 2164/14-2 and Kr 2164/15-2. TDP is also supported by a Warwick Prize Fellowship, made possible by a generous philanthropic donation.

Data availability↩︎

The data underlying this article will be shared upon reasonable request to the corresponding author.

8 Stirring analyses↩︎

This appendix provides more details of our analytic and numerical stirring analyses.

8.1 Catastrophic collisions↩︎

We consider dust released through catastrophic collisions between planetesimals. For a collision to be catastrophic, its kinetic energy needs to be greater than the critical fragmentation energy, \(Q_{\rm D}^*\), such that the largest remnant after the collision has half the original mass of the planetesimal that broke apart ([11]). Approximating the debris particles as basalt, \(Q_{\rm D}^*\) is given by

\[Q_{\rm D}^*(s,v_{\rm coll}) = 9.13~{\rm J~kg^{-1}} \left(\frac{v_{\rm coll}}{\rm m~s^{-1}}\right)^{1/2} \left[ \left(\frac{s}{\rm m}\right)^{\rm -0.36} + \left(\frac{s}{\rm km}\right)^{\rm 1.4}\right], \label{eq:32qdstar}\tag{13}\]

where \(v_{\rm coll}\) is the collisional speed and \(s\) is the debris-body radius (see Equation 1 in [66]).

This \(Q_{\rm D}^*\) prescription results in two regimes in the shape of a “v”. For basalt, the minimum \(Q_{\rm D}^*\) is for debris with size 110 m ([19]), meaning that bodies of this size are easiest to break apart. Smaller debris is considered to be in the strength regime; in this regime, \(Q_{\rm D}^*\) is approximately equal to the shattering strength ([11]). Bodies larger than 110 m are considered to be in the gravity regime, where additional energy is needed to overcome the planetesimal’s gravity and break it apart.

Catastrophic collisions will occur at the onset of the collisional cascade between particles of similar sizes ([19]). From this, we can define the fragmentation speed (dependant on debris size), under the assumption that the material is basalt (see [67]):

\[v_{\rm frag}(s) = 17.5\left[\left(\frac{s}{\rm m}\right)^{-0.36} + \left(\frac{s}{\rm km}\right)^{1.4}\right]^{2/3}~\rm {m/s}. \label{eq:32fragmentation95velocity95app}\tag{14}\]

Note that the equation for fragmentation speed may vary for a given size when considering materials other than basalt.

8.2 Simple analytic stirring analysis↩︎

In this section we detail the derivation of Equations 1 to 3 , which provide a simple approximation of the minimum debris eccentricity required for stirring.

Figure 12 depicts two debris particles with crossing orbits, with the assumptions of co-planar particles of equal size, semimajor axis and eccentricity. We aim to solve for the collisional speeds at the points where the orbits cross, in terms of the rotational angle \(\omega\). Since the orbits (defined with subscripts 1 and 2 for the two bodies in subsequent equations) are rotated from one another, they will have differing true anomalies, f, at the intersection points. At the intersecting points, the orbits will have the same \(x\) and \(y\) positions; hence we can equate \({x_1=x_2}\) and \({y_1=y_2}\) at these locations. The distances \(x\) and \(y\) are given by equations \(x = r\cos{f}\) and \(y = r\sin{f}\), with \(r = a(1-e^{2})(1+e \cos{f})^{-1}\) (Equations 2.21 and 2.20 in [31], respectively). This yields

\[ \cos{f_{\rm 1}} + e\cos{f_{\rm 1}}\cos{f_{\rm 2}} = \cos{(f_{\rm 2}+\omega)} + e\cos{f_{\rm 1}}\cos{(f_{\rm 2}+\omega)}, \label{eq:32initial95vel95simplification95partial951}\tag{15}\]

and

\[ \sin{f_{\rm 1}} + e\sin{f_{\rm 1}}\cos{f_{\rm 2}} = \sin{(f_{\rm 2}+\omega)} + e\cos{f_{\rm 1}}\sin{(f_{\rm 2}+\omega)}. \label{eq:32initial95vel95simplification95partial952}\tag{16}\]

The two solutions for the two intersection points are

\[f_{\rm 1} = \frac{\omega}{2}, \pi + \frac{\omega}{2}, \label{eq:32vel95solution95f951}\tag{17}\]

and

\[ f_{\rm 2} = -\frac{\omega}{2}, \pi - \frac{\omega}{2}, \label{eq:32vel95solution95f952}\tag{18}\]

where all angles are in radians.

With the true anomalies at the intersection points calculated, it is possible to solve for the collisional speed of these debris particles. The collisional speed is

\[v_{\rm coll} = \sqrt{(v_{\rm x,1} - v_{\rm x,2})^2 + (v_{\rm y,1} - v_{\rm y,2})^2}, \label{eq:32initial95collision95vel}\tag{19}\]

where \(v_{\rm x,i}\) and \(v_{\rm y,i}\) are the velocities in the x and y directions for each particle. They are

\[\begin{align} v_{\rm x,1} &= -A \sin{f_{\rm 1}}\\ v_{\rm y,1} &= A(e+\cos{f_{\rm 1}}), \label{eq:32velocity95eq95particle951} \end{align}\tag{20}\]

and

\[\begin{align} v_{\rm x,2} &= -A \sin{f_{\rm 2}}\cos{\omega} - A(e+\cos{f_{\rm 2}})\sin{\omega}\\ v_{\rm y,2} &= -A\sin{f_{\rm 2}}\sin{\omega} + A(e+\cos{f_{\rm 2}})\cos{\omega}, \label{eq:32velocity95eq95particle952} \end{align}\tag{21}\]

where the velocities for particle 2 have been rotated by angle \(\omega\), and the constant A is given by \(A = \left(Gm_*\right)^{1/2}\left[{a(1-e^2)}\right]^{-1/2}\). These velocities are given by Equation 2.36 in [31], with the pre-factor A being the same but in a different form.

To determine whether these particles are considered stirred, we solve for the collisional speed at the locations where the orbits cross. Substituting Equations 20 and 21 into Equation 19 results in the collisional speed equation:

\[v_{\rm coll,simple} = e \sqrt{\frac{2 Gm_*(1-\cos{\omega})}{a(1-e^2)}}. \label{eq:32final95collision95vel}\tag{22}\]

Note that with the assumption that both debris particles have the same semimajor axis and eccentricity, and are offset by angle \(\omega\), both solutions for \(f_1\) and \(f_2\) result in the same simplified collisional speed equation. Therefore, the solution of \({f_1 = \omega /2, f_2 = -\omega/2}\) is used.

To consider the eccentricity below which debris particles will certainly not be stirred, Equations 22 and 14 are compared to find where the collisional speed is less than the fragmentation speed. We are specifically looking for the boundary below which it is certain the debris particles will be unstirred (for the given setup), so to minimize the eccentricity the particles’ orbits are considered to be anti-aligned, i.e. \(\omega = 180°\) (see Figure 13). Doing so, the unstirred eccentricity boundary can be defined as:

\[e_{\rm unstirred} < \frac{v_{\rm frag}}{2} \sqrt{\frac{a_{\rm deb}}{G m_*}}. \label{eq:32unstirred95ecc95appendix}\tag{23}\]

We now plot the collisional speed against semimajor axis, eccentricity, grain size, and orbital rotation \(\omega\) to see how the fragmentation and collisional speeds vary.

Figure 13 shows how the collisional speed changes with orbit-rotational offset \(\omega\). The collisional speed increases the more anti-aligned the orbits are, reaching a maximum at 180. As seen in the figure, higher eccentricities result in a higher collisional speed. It can also be seen that it is possible for the collisional speed to be larger than the fragmentation speed even if the orbits are not fully anti-aligned, as long as the eccentricity of the debris particles is sufficiently high.

Figure 14 shows how the collisional speed varies compared to the semimajor axis for given eccentricities. The collisional speed decreases with increasing semimajor axis, and larger eccentricities provide a greater collisional speed.

The collisional speed can also be compared to the fragmentation speed by considering the grain size, as on Figure 15. The collisional speed is constant with grain size for a given eccentricity, however the fragmentation speed varies with grain size. For basalt, the weakest debris particles would have a radius of 110 m, which can be seen on Figure 15. The fragmentation speed increases in either direction from the grain size of 110 m. As also seen on the plot, there is a minimum particle eccentricity needed so that the collisional speed can be greater than the fragmentation speed at the weakest point. This minimum eccentricity increases when considering debris particles smaller or larger than 110 m.

8.3 Numerical stirring analysis↩︎

In Appendix 8.2 we derived a rough criterion for the minimum eccentricity required for debris to be stirred. However, this is not sufficient to define whether debris is actually stirred, because it assumes that two bodies are on identical, optimally misaligned orbits. To properly assess stirring in our N-body simulations, we create a numerical code to assess the actual collision speeds. The workings of the code are described in this section.

8.3.1 Generalised collision speed↩︎

For the case where two coplanar particles have different orbits with any orientation, we can once again solve for the collisional speed by considering the points where the two orbits intersect. The solutions to the angles of the collision points are given by Equation 4 in [68]:

\[\cos{\theta} = \frac{-AB \pm C\sqrt{C^2 + B^2 - A^2}}{B^2 + C^2}, \label{eq:32crossing95orbit95angle95solution}\tag{24}\]

where

\[\begin{align} A = & \; p_{\rm 2} - p_{\rm 1} \\ B = & \; e_{\rm 1}p_{\rm 2}\cos(\Delta \varpi) - e_{\rm 2}p_{\rm 1} \\ C = & \; e_{\rm 1}p_{\rm 2}\sin(\Delta \varpi), \label{eq:32crossing95orbit95constants} \end{align}\tag{25}\]

where the subscripts again correspond to debris particles 1 and 2. In these equations, \(p_{\rm i}\) is the semilatus rectum of the particle, and \(\Delta \varpi\) is the difference in the two particles’ longitudes of pericentre. The orbits will cross if there are real solutions to Equation 24 , which occurs when the discriminant is positive, i.e. if

\[C^2 + B^2 - A^2 \ge 0. \label{eq:32doOrbitsIntersect}\tag{26}\]

If there are real solutions, then the relative velocities at the intersection points can be calculated as

\[\begin{align} v_{\rm coll} =&\left\{\left[e_{\rm 1}\sin{\left(\theta_{\pm}-\Delta \varpi\right)} \left(\frac{p_{\rm 1}}{\rm au}\right)^{-1/2} - e_{\rm 2}\sin{\left(\theta_{\pm}\right)}\left(\frac{p_{\rm 2}}{\rm au}\right)^{-1/2}\right]^2 \right. \\ &\left.+ \left[\left(\frac{p_{\rm 1}}{\rm au}\right)^{1/2}\left(\frac{r_{\rm \pm}}{\rm au}\right)^{-1} - \left(\frac{p_{\rm 2}}{\rm au}\right)^{1/2}\left(\frac{r_{\rm \pm}}{\rm au}\right)^{-1}\right]^{2}\right\}^{1/2} \\ &\times~29.78~{\rm km/s} \left( \frac{m_*}{\; {\rm M_\odot}}\right)^{1/2}, \end{align} \label{eq:32relative95velocities}\tag{27}\]

which is modified from Equation 7 in [68] to include the orbital speed mentioned in their text. In this equation \(r_{\rm \pm}\) and \(\theta_{\rm \pm}\) define the locations of the intersection points, where

\[r_{\rm \pm} = \frac{p_{\rm 2}}{1+e_{\rm 2} \cos{\theta_{\rm \pm}}}. \label{eq:32r95pm}\tag{28}\]

8.3.2 Numerical implementation↩︎

We want to calculate the level of stirring in a rebound N-body simulation. To do so, we develop a bespoke python script capable of performing this calculation. The code considers pairs of debris particles, and checks whether their orbits intersect using Equation 26 . If they intersect, then the code calculates their collision speed using Equation 27 ; if their collision speed is greater than the fragmentation speed for centimetre-sized debris (Equation 14 , where \({v_{\rm frag}(1~{\rm cm}) \approx 52.8\; {\rm m \; s^{-1}}}\)), then the two particles are considered stirred. In this section we describe the numerical implementation, which we apply to each of our rebound simulations once those simulation have finished.

To check the amount of stirring in an N-body simulation, each particle in the simulated debris disk needs to be checked against other disk particles, to determine whether the conditions mentioned above are fulfilled. This is done at various times throughout the disk’s evolution, to account debris orbits changing over time. To avoid unrealistic scenarios, such as high-eccentricity particles that get scattered from the inner edge of the disk colliding with those at the outer edge before colliding with a nearby particle, or one debris particle colliding with many other particles, some criteria are implemented when comparing pairs of debris particles:

• If two debris particles collide and are considered stirred, they cannot collide with any other particles, either at this time or at any other time in the simulation (i.e. once a particle is stirred, it is assumed to have been collisionally destroyed).

• Debris-particle orbits must be bound (i.e. \(e < 1\)).

• Debris particles’ semimajor axes must have been within a certain range of each other at the start of the simulation (to avoid scattered particles moving through the disk at high eccentricities, and colliding with everything).

• Debris particles’ semimajor axes must be within a certain range of each other at the current time.

• Debris particles’ semimajor axes at the current time must be within a certain range of their initial values (particles whose semimajor axes have changed considerably have probably been scattered, and are therefore not considered).

• Debris particles’ inclinations must be within a certain range of each other at the current time (orbits with high mutual inclination may satisfy the two-dimensional intersection criterion from Equation 26 , but would not intersect in reality).

Only pairs of debris particles satisfying all of the above criteria are considered in the collision-speed analysis.

Different values are tested to determine an appropriate set of numerical values for the above conditions. These parameters are “Particle Near Range” (PNR), which defines how close the semimajor axes of a particle pair should be (at both the initial and current times), “Particle Initial Range” (PIR), which defines how close each particle’s semimajor axis should be to its initial value, and “Inclination Range” (IR), which defines the maximum mutual inclination between particle pairs. We test multiple values for each of these numerical parameters, across multiple simulations, to ensure that they provide realistic and expected results. Table 2 shows the results of one such test, demonstrating how the measured stirring level in the projectile simulation on Figure 2 changes as the numerical parameters are varied. Based on these tests, we select values of \({{\rm PNR} = 1\; {\rm per \; cent}}\), \({{\rm PIR} = 1\; {\rm per \; cent}}\), and \({{\rm IR} = 1°}\). These values are chosen as a compromise between accuracy and speed.

The other main numerical parameter is the time interval between successive checks. Each of our N-body simulations outputs 1000 equally spaced snapshots, where the first snapshot is at time zero, and the final snapshot is at 10 diffusion timescales. Table 3 shows how the measured stirring level for the projectile simulation on Figure 2 changes as a function of the interval between checked snapshots. Based on these results, we set the numerical stirring analysis to check every 20 snapshots.

Using these analyses, the code classifies each debris body as either stirred, scattered or unstirred. A stirred particle is defined as one that satisfies the above conditions, and a scattered body is defined as one whose semimajor axis changes by more than \({20\; {\rm per \; cent}}\) throughout the simulation; all other debris bodies are classed as unstirred.

9 Resonant-stirring efficiency versus other system parameters↩︎

Figure 16 shows the efficiency of resonant stirring as a function of several system parameters.

References↩︎

---

1. E-mail: tysoncosta10@gmail.com↩︎

2. E-mail: tim.pearce@warwick.ac.uk↩︎

3. tdpearce.uk/public-code↩︎

4. A non-zero eccentricity is required to ensure proper behaviour in rebound, but the planet’s initial eccentricity is always much smaller than that required to secularly stir the disk.↩︎

5. A planet will typically scatter non-resonant material originating within a little over 3 Hill radii of its orbit [37], [43]–[47].↩︎

6. tdpearce.uk/public-code↩︎