The Contribution of Dwarf Planets to the Origin of Jupiter Family Comets

Table 1 list the physical parameters of the 34 largest TNOs in the solar system observed as of 2018 February (the objects in the Minor Planet Center database, or MPC, with absolute magnitude H_V ≤ 4); this set includes the four TNOs currently classified as DPs by the IAU (Pluto, Eris, Haumea, and Makemake); however, for simplicity, we call all of the large TNOs used in this study DPs, regardless of their actual shapes or official classifications.

Table 1.  The Sample of the 34 Largest TNOs in the Solar System Used in This Work

Name	R (km)	ρ (gr cm−3)	M ($\times {10}^{-3}$ M⊕)
Eris	1200.0	2.3	2.7956
Pluto	1188.3	1.854	2.4467
Makemake	717.0	2.14	0.5531
Haumea	806.8	1.821	0.6706
2007 OR10	767.2	1.9273	0.6110
Orcus	450.0	1.676	0.1073
Quaoar	551.7	1.99	0.2343
Varda	371.5	1.24	0.0446
2007 UK126	324.0	1.74	0.0415
2002 TX300	143.0	0.933	0.0018
2002 UX25	332.0	0.82	0.0209
Varuna	334.0	0.992	0.0259
2003 AZ84	385.5	0.87	0.0349
Sedna	497.5	1.4532	0.1254
2002 AW197	384.0	1.5277	0.0606
2014 UZ224	317.5	1.1865	0.0266
2005 UQ513	249.0	1.1206	0.0121
Ixion	308.5	1.2811	0.0263
2002 MS4	467.0	1.9176	0.1369
2005 QU182	208.0	1.3957	0.0088
2015 RR245	335.0	1.2577	0.0331
2005 RN43	339.5	0.9326	0.0255
2002 TC302	292.0	1.3725	0.0239
2015 KH162	400.0	1.3236	0.0594
2010 EK139	235.0	0.6002	0.0054
2004 GV9	340.0	0.7910	0.0218
2010 JO179	375.0	1.7194	0.0635
2013 FY27	362.5	1.1722	0.0391
2014 EZ51	443.0	1.6585	0.1012
2010 RF43	305.0	1.1458	0.0227
2003 OP32	313.0	1.3422	0.0288
2012 VP113	250.0	1.0785	0.0118
2010 KZ39	260.0	1.3606	0.0167
2014 WK509	218.5	1.4378	0.0105

Download table as:  ASCIITypeset image

Note that only the 13 firsts objects in Table 1 have reliable measurements of their masses and densities (from Eris to 2003 AZ₈₄). The next 14 objects have only good estimations of their sizes (from Sedna to 2010 JO₁₇₉), while the last seven count only with a measurement of their absolute magnitudes (from 2013 FY₂₇ to 2014 WK₅₀₉). Our procedure to assign masses and densities to the 21 objects without known data is described in detail in the Appendix.

We consider these 34 objects our DP population, similarly to what we have done in previous works when exploring the secular effect of DPs over cold debris disks (Muñoz-Gutiérrez et al. 2015, 2017, 2018). The initial conditions in heliocentric Cartesian elements (positions and velocities) for all of the massive objects in our simulations were retrieved from the JPL Horizons system for the Julian day 2458176.5, corresponding to 2018 February 27.

Figure 1 shows the distribution of mass versus semimajor axis, a, and perihelion distance, q, of the 34 DPs; the values for the mass of Pluto, Ceres, and Mimas are shown for reference. Note that, although 2002 TX₃₀₀ is much less massive than Mimas (blue circle below the Mimas line), we choose to include it in our set despite its small size, because it has a measured density.

Zoom In Zoom Out Reset image size

For the massless disk particles that integrate the population of cometary nuclei in the Kuiper Belt, we used the L7 synthetic model described by Petit et al. (2011) and Gladman et al. (2012), which was made publicly available by the Canada France Ecliptic Plane Survey (CFEPS) team.⁶ This is a model based on the CFEPS survey, aiming to represent the true, unbiased orbital distribution of Kuiper Belt objects down to an absolute magnitude of H_g ≤ 8.5. The model comprises: 50,975 classical particles (representing the cold, hot, and kernel populations of the classical KB); 13,450 resonant particles populating 11 MMRs with Neptune, up to fourth order; and 1612 particles that represent the scattering population.

In Figure 2, we plot the distribution, in the a versus e and a versus i planes, of all the particles of the L7 synthetic model. Note that we do not make a distinction between classical and detached particles. This model suits the needs of our simulations, given that it represents a plausible distribution of objects in the Kuiper Belt, while at the same time it provides a large enough population from which to draw statistically significant conclusions.

Zoom In Zoom Out Reset image size

The implementation and use in our simulations of the classical and scattering populations in the L7 model is straightforward, as the initial conditions in osculating orbital elements are given by the CFEPS team. However, the resonant population requires a more careful treatment: we intend to use a population of particles well-trapped in resonances. This is done to avoid any contamination introduced by the short-term stability of particles that are not protected by the resonant mechanism, despite having a semimajor axis similar to the one of the MMR; such particles, if used, would contaminate the rate of escapees from the Kuiper Belt toward the Neptune crossing region that we are interested in measuring, falsely increasing the rate at which the secular perturbations from DPs are able to destabilize the objects of the KB.

In order to make use only of the actually resonant particles, we ran 10 Myr long simulations including only the four giant planets and the resonant population of the L7 model. In Section 2.4, we summarize the results from this simulation, providing some statistics regarding the resonant particles. For the full-term simulations, we make use of a subpopulation of 8 371 particles from the L7 model with libration amplitudes below 175° in any of the 11 MMRs.

We performed six simulations using the symplectic integrator included in the MERCURY package of Chambers (1999). The initial time-step was set to 400 days, with an accuracy parameter of 10⁻¹⁰ for the Bulirsch–Stöer integrator, which goes into action when particles come closer than 4 Hill radii around each massive particle (giant planets and DPs). All of the simulations included a central star of 1+ M_⊙ mass, where is the summed mass of the interior planets plus the Moon and Ceres. Three of the simulations included only the four giant planets as massive objects, and the other three included the four giant planets and the 34 DPs as massive bodies, as well as the test particles in the Kuiper Belt, with separate simulations for the classical (50,975 particles), resonant (8371 particles), and scattering (1612 particles) populations.

The simulations were 1 Gyr each. We have used this integration time as a compromise between giving enough time for the secular effects to take hold while maintaining sensible CPU time requirements. Also, these simulations are designed to represent (more or less) the current behavior of the solar system: during its formation, the inner solar system must have been a much more violent place; in the far future, less cometary nuclei are expected to remain, and thus fewer new comets are expected to be injected.

Note that one additional simulation used to characterize the resonant population is only 10 Myr long and includes only the four giant planets as massive objects, with the entire set of 13,450 resonant particles from the L7 model.

Depending on their angular elements, particles with semimajor axes with values very close to the nominal location of Neptune's MMRs may or may not be librating in resonance with the planet. Particles actually trapped in MMRs with Neptune should be librating with amplitudes below 180°. Those particles are protected from getting close enough to the giant as to be scattered off by it; thus, they could, in principle, be long-term stable. Of course, this is not always the case, as chaotic perturbations can result in the stirring of eccentricities and in an increasing of the librating amplitude, which finally leads particles into encountering the giant. In any case, it is preferable to only make use of initially librating particles in the long-term simulations, in order to avoid a source of contamination coming from nonresonant particles that may be short-lived, even when they are located close to the resonance (but are not protected by it).

In Figure 2, the initial distribution of particles in the resonant population of the L7 synthetic model is shown by red dots. The resonant particles populate 11 MMRs with Neptune, from zero to fourth order, which are labeled at the top of the figure. Those resonances are the following: 1:1, 2:1, 3:2, 4:3, 5:4, 3:1, 5:3, 5:2, 7:4, 5:1, and 7:3. After 10 Myr of dynamical evolution under the influence of the four giant planets, out of the 13,450 particles initially in the 11 MMrs, 24 were ejected from the system (when their a grew larger than 1000 au) and four collided with the Sun or a planet. From the remainder, we found that 8371 particles librate during more than half of the integration with an amplitude below 175°, while the remaining 5051 particles were nonresonant. We therefore define our resonant population, for the long-term integrations of the next section, as the subset of 8371 librating particles.

Table 2 shows the number of particles initially in each resonance, as well as the number of particles that, after 10 Myr, were librating or turned out to be nonresonant. Note that, for our long-term simulations, the 3:2 MMR contains the largest population, followed by the 5:2 MMR.

To analyze the evolution of particles, originally located in any of the families of the trans-Neptunian region (as given in the L7 model), we begin by defining three stages which any cometary nuclei, that effectively becomes a comet, is most likely to go through. These categories are as follows:

• 1.  
  Crossers: particles whose perihelion, q, is smaller than ${a}_{\mathrm{Nep}}+3\sqrt{2}{R}_{\mathrm{HN}}\approx 33\,\mathrm{au}$, where ${a}_{\mathrm{Nep}}$ and ${R}_{\mathrm{HN}}$ are the semimajor axis and Hill radius of Neptune, respectively (i.e., particles that cross into Neptune's region of influence). This limit is based on previous studies (Gladman & Duncan 1990; Bonsor & Wyatt 2012; Muñoz-Gutiérrez et al. 2018), which found that particles inside this region have a high probability of experiencing a close interaction with the giant planet, thereby jeopardizing their orbital stability.
• 2.  
  Nearly Interactive Particles (NIPs): particles that get closer than one Hill radius to Neptune. These particles experience a close encounter with the giant and are strongly affected in their orbital motion. They can move inward, toward Uranus, or outward, away from the solar system.
• 3.  
  Jupiter Family Comets (JFCs). These are the ones that went inward after interacting with Neptune, and then were handed down through the rest of the giant planets until interacting with Jupiter; each new step has the probability of passing the particle inward or expelling it from the solar system. Once the cometary nuclei is close to Jupiter, the particle is defined as a JFC if its Tisserand parameter with respect to Jupiter, ${T}_{{\rm{J}}}$, remains between 2 and 3 (e.g., Levison & Duncan 1997).

The Tisserand parameter is given by

$$\begin{eqnarray}&&{T}_{{\rm{J}}}=\displaystyle \frac{{a}_{{\rm{J}}}}{a}+2\sqrt{(1-{e}^{2})\displaystyle \frac{a}{{a}_{{\rm{J}}}}}\cos i,\end{eqnarray} \tag{ 1 }$$

where ${a}_{{\rm{J}}}$ is Jupiter's semimajor axis, and a, e, and i are the semimajor axis, eccentricity, and inclination of the particle. It has been shown that ${T}_{{\rm{J}}}$ is a suitable parameter to differentiate between the Jupiter Family and Halley-type comets, both of which are considered SPCs (historically, those with periods below 200 yr; the JFCs having periods below 20 yr). However, by only considering the period as a classifying criteria, one ends up with a large overlap between both families, and thus a dynamically motivated criterion such as ${T}_{{\rm{J}}}$ turns out to be a better choice in defining families among the observed comets, with Halley types having a ${T}_{{\rm{J}}}\lt 2$ (Kresák 1972). In Figure 3, we show the distribution of all the comets known as of October 2018 (obtained from the JPL Small-Body Database Search Engine).

Zoom In Zoom Out Reset image size

We consider the perihelion, q, in order to count a particle in our simulations as a comet, assuming that the cometary nuclei would become visible the first time q becomes less than 2.5 au (i.e., it would display cometary activity in the form of a coma and a tail); this is a common cut used in different works (e.g., Levison & Duncan 1997; Nesvorný et al. 2017) when attempting to compare the results from numerical simulations with the observed population of comets in the solar system.

From the simulations, we saved the information relative to close encounters of all particles below four Hill radii of any major body; from those files, we can extract the time at which one particle first becomes a member of any of the three categories defined above. To measure the effect of DPs over the evolution of the Kuiper Belt's cometary nuclei, we compare the cumulative fraction of particles for simulations with and without DPs; we do this separately for each population of the L7 model, i.e., Classical, Resonant, and Scattering. Please note that objects that start in the Classical (or Resonant) population that become crossers, NIPs, or JFCs will most likely be part of the Scattered population for a period of time, but we will consider them nonetheless as members of the Classical (or Resonant) population, respectively, for the remaining of the paper.

The crosser stage represents an "initial" step in the evolution toward the inner planetary system: after being perturbed inside the non–planet-crossing region of the Kuiper Belt, a particle enters a region where it can be strongly perturbed and its orbit significantly modified by Neptune. If a particle gets close enough to Neptune (an NIP), the interchange of energy and angular momentum with the giant will alter the orbit severely, sending the particle inward or outward. If inward, the same interchange process can occur with the rest of the giant planets until the particle is stabilized in a cometary orbit by Jupiter, with a nearly constant Tisserand parameter between 2 and 3; the particle will remain in this region for at least a few thousands of years, being observable as a comet while it remains active.

In what follows, we consider a JFC as a particle that satisfies the two above criteria ($2\lt {T}_{{\rm{J}}}\lt 3$ and q < 2.5 au) for the first time during its orbital evolution. It is noteworthy that the ability of our simulations to precisely follow orbits after they first become JFCs is limited (this being due to our integration timesteps, our output cadence, the absence of inner planets, and the lack of a model for the outgassing of the comet). Therefore, we prefer to restrict ourselves up to this instant in the history of the cometary nuclei only, when accounting for the effect of DPs on the injection rate of such comets.

The separated results for each of the L7 model populations are presented below.

We performed two simulations (1 Gyr long) with the 8371 resonant particles identified in Section 2.4, this population represents 13.7% of all test particles used in this work. The first simulation included only the four giant planets as massive objects, while the second also included the 34 DPs of Tables 4–6.

The cumulative fraction of particles that become crossers, NIPs, and JFCs are shown in Figure 4. The black lines are the cumulative fractions from the simulation without DPs, while the blue lines are from the simulation with DPs.

Zoom In Zoom Out Reset image size

We observe a significant increase in the number of JFCs in the simulation with DPs, confirming that the interactions with the DPs have a significant secular effect in the evolution of the resonant populations on Gyr timescales. Indeed, the three families of particles experience an increase in the total number of their respective members. This effect has a delay of ∼150 Myr in all cases, showing the secular nature of this phenomenology. After this time, the blue curves remain above the black curves, revealing that the perturbations produced by DPs are enough to increase the population of unstable resonant particles in the Kuiper Belt.

The numbers of crossers, NIPs, and JFCs at the end of the simulation without DPs are 3147, 2585, and 799, respectively;for the simulation with DPs, those numbers are 3571, 2885, and 895. The effective increases of ∼13.5%, 11.6%, and 12.0%, respectively, are due to DPs. We can also observe that, in both simulations, around 80% of the crosser particles become NIPs, and from those, around 30% become JFCs (25% of the crossers become JFCs). Those last fractions are independent of the presence of DPs, as one would expect: after the particles reach the dominion of the giant planets, their outcome will mostly depend on the interchanges of energy and angular momentum with Neptune and the other giant planets, compared to which the perturbations from DPs are negligible. Previous works have also found similar fractions for the particles that become visible JFCs after first encountering Neptune (Levison & Duncan 1997; Volk 2013).

3.2.1. The Contribution of Individual Resonances to the Population of JFCs

In Section 2.4, we define the resonant population on the basis of their libration in each of the 11 MMRs described in the L7 model. We found that only a little more than 60% of the original particles—of those defined as resonant in the L7 model—actually librate in each of the MMRs with Neptune, while a little less than 40% are nonresonant. The MMRs are not homogeneously populated nor do they have a homogeneous librating fraction; in fact, the 3:2 and 5:2 populations alone comprise almost 50% of the original population in the L7 model (see Table 2), while after our characterization, they comprise around 60% of our librating sample. Therefore, it would be expected that the majority of new JFCs that escaped from the resonant regions would actually come from the 3:2 and 5:2 MMRs.

Figure 5 shows the number of new visible JFCs proceeding from each MMR with Neptune for the simulations with and without DPs. The black columns show the number of comets from each MMR for the simulation without DPs, while the red columns stand for the simulation including DPs. A great deal of new JFCs originate from only two resonances: the 3:2 (43.7% of all JFCs without DPs, and 45.0% when DPs are present) and the 5:3 (20.1% without DPs, 14.7% with DPs).

Zoom In Zoom Out Reset image size

Interestingly, the behavior of the 5:3 MMR (and to a lesser extent, the 7:3 and the 4:3 MMRs) is contrary to the rest of the MMRs: when the DPs are included in the simulations, the number of JFCs produced by most MMRs increases, whereas this number decreases by 18% for the 5:3 MMR. Overall, the total number of JFCs produced by MMRs increases by 12% when DPs are present.

Given its contribution of up to 45% of all the JFCs coming from MMRs, the Plutino population deserves an independent analysis. We observe an increment of around 15% in the number of comets produced by the 3:2 resonance as a result of the inclusion of DPs in the simulation; this is the second-largest percentile increment for any of the MMRs analyzed (second only to the 56% increase in the 5:2 MMR) but the largest in actual numbers (54 additional comets). We recall that our population of DPs includes four Plutinos: Pluto, Orcus, 2003 AZ₈₄, and Ixion. To understand the effect of these objects on the diffusion of the whole Plutino population, we perform two additional 1 Gyr simulations including only the 3022 massless particles of the 3:2 MMR, with the giant planets and Pluto as large objects in the first case, and the giant planets and Pluto, Orcus, 2003 AZ₈₄, and Ixion as large objects in the second case.

For the two above simulations, Figure 5 presents the number of JFCs produced in each case. With only Pluto, the number of JFCs produced is 384, an effective increase of 10% with respect to the 349 JFCs obtained when no DPs are included in the simulations. On the other hand, when the four massive Plutinos are included, there are 399 JFCs, nearly identical to the 403 JFCs obtained in the simulation that included all the 34 DPs.

In a previous work, Tiscareno & Malhotra (2009) found that Pluto has only a modest effect, decreasing the mean particle lifetime in the 3:2 resonance by 3%, i.e., an expected decrease of 2% in the number of particles remaining in resonance after 1 Gyr. Although our estimation of the significance of Pluto when producing JFCs may appear significantly larger than the Tiscareno & Malhotra (2009) estimate, in reality the fraction of particles remaining in the resonance is only reduced by around 8%, due to the presence of Pluto. This percentage increases slightly in the simulation that included the four massive Plutinos, as well as in the one with the 34 DPs.

The classical population is initially formed by 50,975 particles with a < 500 au and q > 33 au. This population represents about 83.6% of all the particles simulated in this work, so it is expected to be the main source of JFCs. In Figure 6, we show the evolution of the fraction of particles that become crossers, NIPs, and JFCs in the simulations of the classical population, with DPs (blue lines) and without DPs (black lines).

Zoom In Zoom Out Reset image size

In the same way as for the resonant population, we can see an increment in the fractional number of particles that become crossers, NIPs, and JFCs, due to the presence of DPs; this reaffirms their significance, this time for the broader population of the classical Kuiper Belt.

The numbers of crossers, NIPs, and JFCs obtained in the simulation without DPs are 8586, 5741, and 1332, respectively, while in the simulations with DPs, these numbers become 9211, 6018, and 1561. The observable effect of DPs over the fraction of particles in each stage shows a delay of approximately 400 Myr for the crossers and almost 700 Myr for the NIPs, but only 300 Myr for the JFCs. The respective fractions of particles induced by DPs increased by 7.3%, 4.8%, and 17.2% for each stage. For the classical population, without DPs, 66.8% of crossers become NIPs, and from these, 23.2% become JFCs (15.5% of crossers become JFCs); on the other hand, when DPs are included, these percentages change slightly, with 65.3% of the crossers becoming NIPs and 25.9% of those reaching the JFCs stage (in this case, 16.9% of crossers become JFCs). These numbers show that, for the classical population, the added particles reaching the stage of NIPs due to perturbations from DPs have a slightly larger probability of becoming JFCs.

The scattering population in the L7 model is formed by only 1612 particles, which account for 2.7% of all the particles simulated in this work. In Figure 7, we show the evolution of the fractional number of particles for the crossers, NIPs, and JFCs categories in the simulations with DPs (in blue lines) and without DPs (in black lines). For this population, there is not an appreciable increase in the number of crossers and NIPs due to the presence of DPs, and only a slight increase in the production of JFCs is obtained when DPs are included.

Zoom In Zoom Out Reset image size

The number of crossers, NIPs, and JFCs in the simulation without DPs are 1342, 1001, and 235, respectively; with DPs, on the other hand, these numbers are 1336, 1003, and 242, respectively. In this case, it is not straightforward to link the slight increase of comets to the effect of DPs, given that such small variations could result from statistical fluctuations expected from the two different simulations (yet the presence of DPs results in a 3% increase in the total number of JFCs).

Most of the scattering objects are destined to evolve in this track. Note that more than 80% of the total population of scattering objects become crossers, while more than 60% of the whole population become NIPs. In fact, for both simulations, around 75% of crossers become NIPs—and from those, around 23% become JFCs (around 17% of crossers become JFCs).

We can further confirm that the increment in the observed number of JFCs in our simulations is a direct product of the dynamical perturbations produced by DPs. To do so, we have compared the degree of overall dynamical excitation on the trans-Neptunian cometary particles produced by the giant planets alone with the excitation produced by the combined effect of the giant planets and the 34 DPs.

In general, we expect that the inclusion of massive perturbers in the trans-Neptunian region would result in an increase of dynamical excitation of disk particles. From our simulations, we observe that such an increase, although small, is significant when compared with the dynamical excitation produced by the giant planets alone. In Figure 8, we show the root-mean-square (rms) eccentricity and inclination of the L7 model disk particles as a function of time. Note that we include all the particles that began the simulation, including the ones that were ejected from the system (and would not be observable in the solar system).

Zoom In Zoom Out Reset image size

In the upper panel of Figure 8, we can see how the rms eccentricities of the three simulated populations evolve to larger values when DPs are present: this is clear for the Scattering and the Resonant populations, although for the Classical population the difference is marginal but noticeable. The latter is due mainly to the orbital distribution of DPs. Given that only 15 DPs are located beyond 50 au, they are unable to exert a more significant perturbation on a large fraction of particles in a vast region of space (44% of the Classical population have initial a > 50 au). Note that periods of at least 100 Myr for the Scattering population, and around 300 Myr for the Resonant population, are required for the effect of DPs to become noticeable. Finally, we can conclude that, as expected, the DPs contribute to globally stirring the eccentricities of the objects in the trans-Neptunian region.

In the lower panel of Figure 8, we show the evolution of the rms inclinations of the three simulated populations as a function of time. As for the rms eccentricities, the rms inclinations of the Scattering and Resonant populations are larger at the end of the simulations, and the final values are larger with the presence of DPs. For the Classical population, a reverse behavior is observed, as the final value of the rms inclination of the population is slightly smaller than at the beginning of the simulation; this could signify that the L7 model's inclination is slightly overestimated, as we would expect the Classical population to be nearly in equilibrium; i.e., this could be due to an excess of large inclination particles in the L7 Classical population.

Figure 9 is similar to Figure 8, but in this case we exclude the particles that are ejected from the simulation (i.e., we only use the particles that would be observed in the solar system). As with the case in which all particles are counted, here we see that the rms eccentricities for the three populations increase as a function of time, with the final values being clearly larger when DPs are present for the Scattering and Resonant populations—and again, only marginally larger for the Classical population.

Zoom In Zoom Out Reset image size

The rms inclinations of the three populations as a function of time are shown in the lower panel of Figure 9. Here, we observe the opposite behavior to the one shown by the rms eccentricities (upper panel), with the final rms inclination values being smaller, for the Classical and Resonant populations, than their initial value at the beginning of the simulations—and such values are even lower when DPs are present for the Resonant population. In the case of the Scattering population, after an initial increase of the rms inclination (lasting up to 300 Myr), the subsequent trend shows a clear decrement of this value, which is more pronounced with the presence of DPs. Although the resonant population shows an overall decrement in the final inclination, there is also a small increase in the first 100–150 Myr. The difference in behavior between Figure 9 and Figure 8 is due to the selection effect caused by the escaping cometary nuclei, i.e., escaping cometary nuclei will have, on average, larger inclinations and eccentricities than the typical disk particle.

Figure 10 shows the remaining fraction of particles in each population of the L7 model, with and without DPs, as a function of time, using the same color code as that of the previous figures. It is interesting to note how almost immediately the number of particles of the Scattering population begins to decrease, while at the same time the rms eccentricity and inclination began their increasing trends. Because the scattering particles are the most easily perturbed, due to their proximity to Neptune, many of them are lost quickly, either by being ejected from the solar system or by becoming JFCs.

Zoom In Zoom Out Reset image size