Origin and Evolution of Short-period Comets

Our numerical integrations consist of tracking the orbits of the four giant planets (Jupiter to Neptune) and a large number of small bodies representing the outer planetesimal disk. The terrestrial planets are not included. To set up an integration, Jupiter and Saturn are placed on their current orbits. Uranus and Neptune are placed inside their current orbits and are migrated outward. The initial semimajor axis ${a}_{{\rm{N}},0}$, eccentricity ${e}_{{\rm{N}},0}$, and inclination ${i}_{{\rm{N}},0}$ define Neptune's orbit before the main stage of migration/instability. Here we use ${a}_{{\rm{N}},0}=22\,\mathrm{au}$, ${e}_{{\rm{N}},0}=0$, and ${i}_{{\rm{N}},0}=0$.

The swift_rmvs4 code, part of the Swift N-body integration package (Levison & Duncan 1994), is used to follow the orbital evolution of all bodies. The code was modified to include artificial forces that mimic the radial migration and damping of planetary orbits. These forces are parametrized by the exponential e-folding timescales, ${\tau }_{a}$, ${\tau }_{e}$, and ${\tau }_{i}$, where ${\tau }_{a}$ controls the radial migration rate, and ${\tau }_{e}$ and ${\tau }_{i}$ control the damping rates of e and i (NV16). We set ${\tau }_{a}={\tau }_{e}={\tau }_{i}$ because such roughly comparable timescales were found in NM12.

The numerical integration is divided into two stages with migration/damping timescales ${\tau }_{1}$ and ${\tau }_{2}$ (NV16). The first migration stage is stopped when Neptune reaches ${a}_{{\rm{N}},1}\,\simeq 27.7\,\mathrm{au}$. Then, to approximate the effect of planetary encounters during the dynamical instability in NM12, we apply a discontinuous change of Neptune's semimajor axis and eccentricity, ${\rm{\Delta }}{a}_{{\rm{N}}}$ and ${\rm{\Delta }}{e}_{{\rm{N}}}$. Motivated by the results of NM12 and Nesvorný (2015b), we set ${\rm{\Delta }}{a}_{{\rm{N}}}=0.5\,\mathrm{au}$ and ${\rm{\Delta }}{e}_{{\rm{N}}}=0.1$. The second migration stage starts with Neptune having the semimajor axis ${a}_{{\rm{N}},2}={a}_{{\rm{N}},1}+{\rm{\Delta }}{a}_{{\rm{N}}}$. We use the swift_rmvs4 code and migrate the semimajor axis (and damp the eccentricity) on an e-folding timescale ${\tau }_{2}$. The migration amplitude was adjusted such that the planetary orbits obtained at the end of the simulations were nearly identical to the real orbits.⁹

We found from NM12 that the orbital behavior of Neptune during the first and second migration stages can be approximated by ${\tau }_{1}\simeq 10\,\mathrm{Myr}$ and ${\tau }_{2}\simeq 30\,\mathrm{Myr}$ for a disk mass ${M}_{\mathrm{disk}}=20$ ${M}_{\oplus }$, and ${\tau }_{1}\simeq 30\,\mathrm{Myr}$ and ${\tau }_{2}\simeq 50\,\mathrm{Myr}$ for ${M}_{\mathrm{disk}}=15$ ${M}_{\oplus }$. The real migration slows down, relative to a simple exponential, at late stages. Here we therefore set ${\tau }_{1}=10\mbox{--}30$ Myr and ${\tau }_{2}=30\mbox{--}100$ Myr. All migration simulations were run to 0.5 Gyr. They were extended to 4.5 Gyr with the standard swift_rmvs4 code (i.e., without migration/damping after 0.5 Gyr).

We developed an approximate method to represent the jitter that Neptune's orbit experiences due to close encounters with massive planetesimals. The method has the flexibility to use any smooth migration history of Neptune as an input, include any number of massive planetesimals in the original disk, and generate a new migration history where the random element of encounters with the massive planetesimals is taken into account. This approach is useful, because we can easily control how grainy the migration is while preserving the global orbital evolution of planets from the smooth simulations. See NV16 for a detailed description of the method. Here we set the mass of massive planetesimals to be equal to that of Pluto. This choice is motivated by the fact that two Pluto-class objects are known in the Kuiper Belt today (Pluto and Eris).

The migration graininess is included in the present integrations, because NV16 showed that grainy migration of Neptune is required to get the right proportion of resonant populations in the classical belt. The migration graininess may not be important for cometary reservoirs, but we include it in the present work for completeness.

The planetesimal disk was divided into two parts. The part from just outside Neptune's initial orbit ($\simeq 22$ au) to ${r}_{\mathrm{edge}}$ represents the massive inner part of the disk (NM12). We used ${r}_{\mathrm{edge}}=28\mbox{--}30$ au, because our previous simulations in NM12 showed that the massive disk's edge must be at 28–30 au for Neptune to stop at ≃30 au (Gomes et al. 2004). The estimated mass of the planetesimal disk below 30 au is ${M}_{\mathrm{disk}}\simeq 15\mbox{--}20\ {M}_{\oplus }$ (NM12). The massive disk has a crucial importance here, because it was the main source from which cometary reservoirs formed (Dones et al. 2015). The planetesimal disk had an outer, low-mass extension reaching from ${r}_{\mathrm{edge}}$ to at least $\simeq 45\,\mathrm{au}$. The disk extension is needed to explain why the Cold Classicals (hereafter CCs) have several unique physical and orbital properties, but it should not substantially contribute elsewhere, because of the small original mass of the extension.

Specifically, Fraser et al. (2014) estimated that the current CC population at 42–47 au represents only ∼0.0003 Earth masses. Assuming that CCs did not lose much of its original population during planetary migration (e.g., Nesvorný 2015b found that the primordial population of CCs was reduced by a factor of ∼2 during early stages), this shows that the surface density of solids must have dropped substantially from 30 au to 42 au. The profile of the surface density at 30–42 au is not well constrained, but given that Neptune stopped at 30 au, it is reasonable to assume that the density decreased immediately beyond 30 au. Here we therefore choose to ignore the outer extension of the disk at >30 au. Including it would probably not change our results substantially. Presumably, the dynamics of objects starting at 30–35 au during Neptune's migration would be similar to those starting at <30 au. They would become scattered by Neptune (Levison et al. 2008) with a small fraction ending in the scattered disk and Oort cloud, but this contribution should presumably be minor compared to that of the inner, more massive disk. A detailed investigation of this subject is beyond the scope of the present work.

Each of our simulations included one million disk bodies distributed from outside Neptune's initial orbit to ${r}_{\mathrm{edge}}$. The radial profile was set such that the disk surface density ${\rm{\Sigma }}\propto 1/r$, where r is the heliocentric distance. The initial eccentricities and initial inclinations of disk bodies in our simulations were distributed according to the Rayleigh distribution with ${\sigma }_{e}=0.05$ and ${\sigma }_{i}=2^\circ$, where σ is the usual scale parameter of the Rayleigh distribution (the mean of the Rayleigh distribution is equal to $\sqrt{\pi /2}\sigma$). The disk bodies were assumed to be massless, such that their gravity did not interfere with the migration/damping routines.

The gravitational effects of the hypothetical fifth giant planet (Nesvorný 2011; Batygin et al. 2012, NM12) on the disk planetesimals were ignored. The fifth giant planet orbit probably crossed into the trans-Neptunian reservoirs only briefly, during some ∼10⁵ yr, before it was ejected by Jupiter into interstellar space. It likely did not cause major perturbations of orbits in the trans-Neptunian region, although this may depend on how exactly planets evolved during the instability (Batygin et al. 2012). In any case, it is difficult to account for the fifth giant planet with the numerical scheme used here, because its orbit evolves chaotically during the instability and cannot be easily parametrized. To include the fifth planet in a simulation, the orbital histories of planets would need to be taken directly from the self-consistent simulations of planetary instability/migration (e.g., Nesvorný et al. 2013).

P9 is hypothesized, via its gravitational shepherding effects, to produce the non-uniform distribution of orbital angles of known extreme KBOs ($a\gt 150\,\mathrm{au}$, $q\gt 35\,\mathrm{au};$ Trujillo & Sheppard 2014; Batygin & Brown 2016a). Its existence could also explain the tilt of the plane of the solar system with respect to the solar equator (Bailey et al. 2016; Lai 2016; Gomes et al. 2017), and the high inclinations of large semimajor axis Centaurs (Gomes et al. 2015, Batygin & Brown 2016b). The predicted parameters of P9 are ${M}_{9}\gtrsim 10$ ${M}_{\oplus }$, $400\lesssim {a}_{9}\lesssim 900\,\mathrm{au}$, $0.4\lesssim {e}_{9}\lesssim 0.8$, and ${i}_{9}\lesssim 30^\circ$. In this work, we performed several simulations with P9 to investigate the influence of P9 on the dynamical structure of distant cometary reservoirs and on comet delivery. The goal was to diagnose properties of P9 from orbital characteristics of the cometary populations.

To construct an adequate model with P9, we first need to know how and when P9 reached its current orbit. We considered several possibilities (e.g., Kenyon & Bromley 2016; Li & Adams 2016). It seems to us, for example, that P9 cannot be related to the hypothesized fifth giant planet (Nesvorný 2011; Batygin et al. 2012, NM12). This is especially clear if the planetary instability happened late, because in that case it is hard to imagine any plausible mechanism that could have stabilized the fifth planet on a wide orbit. Instead, we find it more plausible that P9 reached its wide orbit well before the epoch of planetary instability (i.e., before Neptune dispersed the massive planetesimal disk below 30 au; e.g., Izidoro et al. 2015). This would mean that the dynamical origin of the trans-Neptunian populations, including distant cometary reservoirs, postdates the chain of events that ended with P9 on its wide orbit. Working under this assumption, here we performed several simulations where P9 was included as a massive perturber since t = 0.

We assumed that the Galaxy is axisymmetric and the Sun follows a circular orbit of radius R₀ in the Galactic midplane. The Sun's angular speed about the Galactic center is ${{\rm{\Omega }}}_{0}$. Then, from Levison et al. (2001; see also Heisler & Tremaine 1986; Wiegert & Tremaine 1999), the Galactic tidal acceleration is given by

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{\mathrm{tide}}={{\rm{\Omega }}}_{0}^{2}\,\left[(1-2\delta )x{{\boldsymbol{e}}}_{x}-y{{\boldsymbol{e}}}_{y}-\left(\displaystyle \frac{4\pi G{\rho }_{0}}{{{\rm{\Omega }}}_{0}^{2}}-2\delta \right)z{{\boldsymbol{e}}}_{z}\right],\end{eqnarray} \tag{ 2 }$$

where $\delta =-(A+B)/(A-B)$, A and B are the Oort constants, G is the gravitational constant, and ${\rho }_{0}$ is the mass density in the solar neighborhood. Here, the coordinate system $({{\boldsymbol{e}}}_{x},{{\boldsymbol{e}}}_{y},{{\boldsymbol{e}}}_{z})$ was chosen such that ${{\boldsymbol{e}}}_{x}$ points away from the Galactic center, ${{\boldsymbol{e}}}_{y}$ points in the direction of the Galactic rotation, and ${{\boldsymbol{e}}}_{z}$ points toward the south Galactic pole. Rotations were applied to move between $({{\boldsymbol{e}}}_{x},{{\boldsymbol{e}}}_{y},{{\boldsymbol{e}}}_{z})$ and the reference system of our integrations, which has the Z-axis pointing along the initial angular momentum vector of the planets.

Numerically, we used ${\rho }_{0}=0.1$ M_⊙ pc⁻³ (M_⊙ is the solar mass) in most of our simulations, and tested ${\rho }_{0}=0.2$ M_⊙ pc⁻³ in one case as well. The Oort constants were set to be A = 14.82 km s⁻¹ kpc⁻¹ and $B=-12.37$ km s⁻¹ kpc⁻¹, giving $\delta =-0.09$. Also, ${{\rm{\Omega }}}_{0}=A-B=2.78\times {10}^{-8}$ yr⁻¹.

The effect of stellar encounters was modeled in the N-body code by adding a star at the beginning of its encounter with the Sun and removing it after the encounter was over. The stars were released and removed at the heliocentric distance of 1 pc (206,000 au). We used the model of Heisler et al. (1987) to generate stellar encounters. The model accounts for 20 species of main-sequence stars and white dwarfs. The stellar mass and number density of different stellar species were computed following the procedure outlined in Heisler et al. (1987). For each species, the velocity distribution was approximated by an isotropic Maxwellian with one-dimensional variance. The number of stellar encounters below perihelion distance q therefore followed $N(\lt \ q)\propto {q}^{2}$. The dynamical effects of passing molecular clouds were ignored.

The early stages, when the Sun presumably interacted with other stars in an embedded globular cluster (Adams 2010), are not considered here. On one hand, the effect of stellar encounters during these stages may be needed to explain the detached orbits of some extreme KBOs (e.g., Sedna and VP113; Levison et al. 2004). On the other hand, cometary-size disk planetesimals have been strongly affected by aerodynamic gas drag during the early stages (before the nebular gas was removed by photoevaporation). Instead of being ejected to large heliocentric distances, the orbits of small bodies were probably circularized by gas drag inside and outside of planetary orbits (e.g., Brasser et al. 2007). If so, these early stages would not substantially contribute to the formation of cometary reservoirs.

We performed 14 simulations in total (Table 1). Two reference simulations were performed without P9, the Galactic tide, or stellar encounters. They differed in the timescale of Neptune's migration: ${\tau }_{1}=30\,\mathrm{Myr}$ and ${\tau }_{2}=100\,\mathrm{Myr}$ (CASE 1 or C1 for short) and ${\tau }_{1}=10\,\mathrm{Myr}$ and ${\tau }_{2}=30\,\mathrm{Myr}$ (CASE 2 or C2). Other simulations used the same migration parameters as C1 or C2, but also included some combination of P9, Galactic tide, and/or stellar encounters. Three simulations were performed in C1 with no P9. In one job, we used ${\rho }_{0}=0.1$ M_⊙ pc⁻³ (C1G1) and no stellar encounters. The remaining two jobs were done with ${\rho }_{0}=0.1$ M_⊙ pc⁻³ and ${\rho }_{0}=0.2$ M_⊙ pc⁻³, and stellar encounters (C1G1S and C1G2S, respectively).

Table 1.  A Summary of the Numerical Integrations Performed in this Work

	${\tau }_{1}$	${\tau }_{2}$	M9	a9	e9	i9	${\rho }_{0}$	SE
	(Myr)	(Myr)	(${M}_{\oplus }$)	(au)		(° )	(${M}_{\odot }\,{\mathrm{pc}}^{-3}$)
C1	30	100	0	⋯	⋯	⋯	0	no
C1G1	30	100	0	⋯	⋯	⋯	0.1	no
C1G1S	30	100	0	⋯	⋯	⋯	0.1	yes
C1G2S	30	100	0	⋯	⋯	⋯	0.2	yes
C1M15	30	100	15	700	0.6	30	0	no
C1M20a	30	100	20	500	0.5	15	0	no
C1M20b	30	100	20	700	0.6	30	0	no
C1M20c	30	100	20	900	0.78	30	0	no
C1I0	30	100	20	700	0.6	0	0	no
C1ALL	30	100	15	700	0.6	30	0.1	yes
C2	10	30	0	⋯	⋯	⋯	0	no
C2M10	10	30	10	700	0.6	30	0	no
C2M20	10	30	20	700	0.6	30	0	no
C2M30	10	30	30	700	0.6	30	0	no

Note. Column SE indicates whether stellar encounters were included in each job.

Download table as:  ASCIITypeset image

In addition, we performed nine simulations with different masses and orbits of P9. We used ${M}_{9}=10$, 15, 20, and 30 M_⊕; ${a}_{9}=500$, 700, or 900 au; and i₉ = 0°, 15°, or 30°. The eccentricity of P9 was set in each case from the solar obliquity constraint (Bailey et al. 2016; Lai 2016; Gomes et al. 2017), except for one case with ${i}_{9}=0$, where we used ${M}_{9}=20$ M_⊕, ${a}_{9}\,=700\,\mathrm{au}$, e₉ = 0.6, and C1 migration parameters. None of these simulations, except one, included effects of the Galactic tide or stellar encounters. Our most complete job with C1 migration parameters included P9 with ${M}_{9}=15$ M_⊕, ${a}_{9}=700\,\mathrm{au}$, e₉ = 0.6, ${i}_{9}=30^\circ$, Galactic tide with ${\rho }_{0}=0.1$ M_⊙ pc⁻³, and stellar encounters.

The simulations were performed on NASA's Pleiades supercomputer (10 jobs in C1, one in C2) and Prague cluster Tiger (three jobs in C2). One C1 simulation over 4.5 Gyr required about 600 hr on 25 Ivy Bridge nodes (20 cores each) of Pleiades, totaling over 34 CPU-years per simulation.

The last integration segment, between $t=3.5\,\mathrm{Gyr}$ and $t=4.5\,\mathrm{Gyr}$ (${\rm{\Delta }}T=1\,\mathrm{Gyr};$ time t is defined such that t = 0 at the start of our integrations about 4.5 Gyr ago and t runs forward in time to t = 4.5 Gyr at the current epoch), was performed with a code specialized for the analysis of cometary orbits. First, we used cloning to improve the statistics of orbits reaching below Saturn's orbit. This was done by monitoring the heliocentric distance, r, of each body at each time step (0.5 yr). If $r\lt 9\,\mathrm{au}$ for the first time, the body was cloned 100 times producing 100 new (cloned) orbits. The cloned orbits were generated by small random perturbation of the velocity vector of the original orbit. A similar method was used in BM13.

Second, in addition to the normal output of Swift, we modified the code to output the orbital elements of comets in 100 yr intervals. The information was written in separate output files if $a\lt 35\,\mathrm{au}$, corresponding to $P\lesssim 200\,\mathrm{yr}$, and $q\lt 5.2\,\mathrm{au}$. SPCs with perihelion distances beyond the orbit of Jupiter were not recorded in the file, but some statistics for them can be obtained from the standard Swift output. The orbital elements written in the output file were rotated to the reference plane defined by the instantaneous angular momentum vector of the four outer planets (Jupiter to Neptune). This is required because P9, included in some simulations, acts to tilt the angular momentum vector of the Jovian planets by several degrees over 4.5 Gyr.

Third, a detailed output was implemented for LPCs. This was done by monitoring the heliocentric distance of each body in the simulation, including clones. If a body reached $r\lt 5.2\,\mathrm{au}$, we recorded the body's and planetary state vectors into a special "LPC" file. A separate output was written in the LPC file for each perihelion passage with $q\lt 5.2\,\mathrm{au}$. After the whole simulation was over, in a subsequent set of simulations, we used the LPC file to set up backward integrations such that we can determine the orbital elements of each comet before it entered into the planetary region. The orbital elements were calculated near orbital aphelion, if the aphelion distance $Q\lt 200\,\mathrm{au}$, or near r = 200 au, if $Q\gt 200$ au (and for hyperbolic orbits).

Nongravitational forces on comets were ignored. We used a relatively long time step, 0.5 yr, in all main integrations, which is roughly 1/20 of Jupiter's orbital period, and verified that using a shorter time step (we tried 0.25, 0.1, and 0.05 yr) does not significantly affect the results. The simulation results were used to build a steady-state model of comets. Initially, we used the full length of the last integration segment (${\rm{\Delta }}T=1$ Gyr) to obtain the best statistics. Subsequently, we also tested how the results depend on the length of the time segment used for the analysis. Using a short time segment ${\rm{\Delta }}T$ near the current epoch should more closely reflect the present population of comets. On the other hand, the statistics become inadequate if ${\rm{\Delta }}T$ is too short, especially for HTCs and if a short physical lifetime is assumed (see the next section). We did not find any significant differences in the results and used ${\rm{\Delta }}T=1\,\mathrm{Gyr}$, which has the best statistics, in the rest of this work. Note that the transfer time of comets from $q\lt 9\,\mathrm{au}$ to $q\lt 2.5\,\mathrm{au}$ is short (∼6 Myr) compared to ${\rm{\Delta }}T=1\,\mathrm{Gyr}$. Thus, not cloning bodies that reached $q\lt 9\,\mathrm{au}$ just before $t=3.5\,\mathrm{Gyr}$ is an adequate approximation.

The steady-state model of SPCs is compared to observations in Section 4. In addition to the distribution of orbital elements a, q, and i, we also consider the distribution of ${T}_{{\rm{J}}}$. To do this comparison correctly, as pointed out in LD97, we must account for the physical lifetime of active comets (i.e., how long comets remain active). We considered three different parametrizations of the physical lifetime:

• Number of perihelion passages with $q\lt 2.5\,\mathrm{au}$. In our simplest parametrization of the physical lifetime, we count the number of perihelion passages with $q\lt 2.5\,\mathrm{au}$ ${N}_{{\rm{p}}}(2.5)$, and assume that a comet becomes inactive if ${N}_{{\rm{p}}}(2.5)$ exceeds some threshold. The threshold is determined by orbital fits to observations.
• Time spent with $r\lt 2.5\,\mathrm{au}$. We determine the time spent by each body with $r\lt 2.5\,\mathrm{au}$, $T(2.5)$, and assume that a comet becomes inactive if $T(2.5)$ exceeds some threshold. Relative to the ${N}_{{\rm{p}}}(2.5)$ criterion, $T(2.5)$ penalizes orbits with low q and/or low a values, because bodies with these orbits spend more time below 2.5 au.
• Heliocentric distance weighted effective erosion time. Comets reaching low heliocentric distances are heated by solar irradiation and are expected to erode faster than more distant comets. The nature of the relationship between the effective¹⁰ erosion rate and heliocentric distance is uncertain. Here we assume, motivated by the heliocentric distance dependence of the water ice sublimation rate (Marsden et al. 1973), that the erosion rate is proportional to ${r}^{-2}$ if $r\lt 2.5\,\mathrm{au}$. The time spent at each r is then weighted by ${r}^{-2}$ and accumulated in ${T}_{{\rm{e}}}(2.5)$. Relative to $T(2.5)$, ${T}_{{\rm{e}}}(2.5)$ penalizes bodies reaching low heliocentric distances.

The parametrizations described above are a compromise between complexity and realism. More complex models, such as the splitting model of Di Sisto et al. (2009), are not considered. These models may be more realistic but have more parameters and are therefore difficult to constrain. We do not use LD97's parametrization of the physical lifetime (see also BM13), because their parametrization was developed for ECs, and is not applicable to HTCs or LPCs, which have much longer orbital periods. Two possibilities exist for a comet to become inactive: it either becomes dormant or it disrupts and disappears. We do not distinguish between these different possibilities in this work and attempt to constrain our model from observations of active comets.

Figures 4–6 show the orbital structure of the trans-Neptunian region at the end of our simulations (t = 4.5 Gyr). The results of the simulations without P9 and with different assumptions about external perturbations from the Galaxy and passing stars are compared in Figure 4. There are two notable structures in Figure 4. The first one is the scattered disk that extends along the Neptune-crossing line to ∼1000 au. The scattered disk is created when bodies encounter Neptune and are scattered along the $q\sim 30\,\mathrm{au}$ line to very large heliocentric distances. The detailed orbital structure of the inner part of the scattered disk, and how it depends on Neptune's migration, was recently discussed in Kaib & Sheppard (2016) and Nesvorný et al. (2016).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We find that the scattered disk, here defined as orbits with $50\lt a\lt 1000\,\mathrm{au}$, contains ∼3000 bodies at $t=4.5\,\mathrm{Gyr}$, of which ≃80% have $50\lt a\lt 200\,\mathrm{au}$ (hereafter the inner SDOs). This means that the population of inner SDOs is much larger, roughly four times larger, than the population of outer SDOs ($200\lt a\lt 1000$ au). The total number of SDOs with $50\lt a\lt 1000\,\mathrm{au}$ represents a fraction $\simeq 3\times {10}^{-3}$ of the original 10⁶ disk bodies at t = 0, or, in terms of mass, $\simeq 0.06$ ${M}_{\oplus }$ for ${M}_{\mathrm{disk}}=20$ ${M}_{\oplus }$. This estimate is consistent with the results reported in Nesvorný et al. (2016). BM13 found from their simulations that SDOs should represent a fraction ≃6 × 10⁻³ of the original disk, which is a value $\simeq 2$ times larger than what we found here. The difference can be attributed to the different orbital evolutions of planets (BM13 adopted orbital histories of planets from the original Nice model and Levison et al. 2008).

In our simulations, $\simeq 1500$ bodies ended on stable orbits in the classical KB with $a\lt 50$ au (this includes hot classicals and resonant populations), corresponding to the fraction $\simeq 1.5\times {10}^{-3}$ of the original disk, or $\simeq 0.03$ ${M}_{\oplus }$. According to these results, the scattered disk should presently be ∼2 times more populous/massive than the classical KB. For comparison, Trujillo et al. (2001) estimated from observations that the mass of the scattered disk is ∼0.03 ${M}_{\oplus }$, which is a value two times lower than the one found here, but their 1σ uncertainty admits masses as high as ∼0.06 ${M}_{\oplus }$, which would be in good agreement with our work. Trujillo et al. (2001) also suggested that the mass of the scattered disk is similar to that of the classical KB with $a\lt 50\,\mathrm{au}$, while Fraser et al. (2014) found instead that the mass of the classical KB should only be ∼0.01 ${M}_{\oplus }$, which is three times lower than Trujillo's estimate for the scattered disk. Thus, while there is general agreement to within a factor of a few among different works, a better characterization of the trans-Neptunian population from observations will be needed to test our model in detail.

The orbital structure of the inner scattered disk is very different from that of the outer scattered disk. Most inner SDOs (≃80%) are fossilized, meaning that their (barycentric) semimajor axis did not change by more than 1.5 au over the last 1 Gyr. This includes objects that interacted with Neptune's orbital resonances in the past and subsequently decoupled from Neptune through various dynamical processes (Kaib & Sheppard 2016; Nesvorný et al. 2016). The remaining ≃20% of inner SDOs are being actively scattered by Neptune (hereafter the scattering SDOs; Gladman et al. 2008). Nearly all outer SDOs, on the other hand, are on scattering orbits (i.e., their semimajor axis changes by more than 1.5 au in the last 1 Gyr). Thus, even though the inner scattered disk is more massive than the outer one, the number of scattering objects in each population is roughly the same.

The second notable feature in Figure 4 is the Oort cloud. In Figures 4(a) and (b), where the stellar encounters were ignored, the Oort cloud has a well-defined structure with inner ($a\lt {\rm{20,000}}$ au) and outer parts ($a\gt {\rm{20,000}}$ au). The outer part of the Oort cloud forms first and is already present in our simulations in the first 10 Myr. By checking on the orbital histories of outer Oort-cloud bodies, we found that most of them reached $a\gt {\rm{20,000}}\,\mathrm{au}$ after having encounters with Saturn (and typically without having encounters with Jupiter; Dones et al. 2004). In addition, a significant fraction of outer Oort-cloud bodies reached their distant orbits by being scattered by Uranus or Neptune (and without having encounters with Jupiter or Saturn). The inner Oort cloud formed as a "wave front" of orbits in our simulations that was moving from the outside in as time advanced. Most bodies that ended up in the inner Oort cloud were scattered to $a\gt 1000\,\mathrm{au}$ by Neptune (some after having encounters with Uranus, but rarely with Jupiter/Saturn).

The gap between the inner and outer parts of the Oort cloud at $a\simeq {\rm{20,000}}\,\mathrm{au}$ also formed gradually in our simulations as orbits were slowly removed from this region. The removal process is controlled by the period of Kozai cycles produced by the Galactic tide. For $a\lt {\rm{10,000}}\,\mathrm{au}$, the Kozai period is longer than the age of the solar system (Higuchi et al. 2007), and orbits that become decoupled from the Jovian planets by the Galactic tide do not have time to complete one Kozai cycle. These orbits persist to the end of the simulations. The orbits with $a\gt {\rm{10,000}}\,\mathrm{au}$, on the other hand, have shorter Kozai periods and can complete one or more Kozai cycles. Once the semimajor axis is above $a\sim {\rm{20,000}}\,\mathrm{au}$, however, the time for a comet to cycle from $q\gt 30$ au to $q\lt 30\,\mathrm{au}$ to back to $q\gt 30\,\mathrm{au}$ is substantially shorter than its orbital period. Thus, even if q drops below 30 au, the comet may never make a passage near the planets, and the planets are thus less efficient at influencing orbits with $a\gt {\rm{20,000}}\,\mathrm{au}$.

The inner Oort cloud has an anisotropic distribution of inclinations. Two features can be noted in Figure 4(b): (1) the retrograde orbits with $a\lt {\rm{20,000}}\,\mathrm{au}$ generally do not have $i\gt 150^\circ$, and (2) there is a concentration of prograde orbits with $i\sim 30^\circ$. Issue (1) is related to the fact that the Galactic tide below 20,000 au can be closely approximated by the quadrupole term. In the quadrupole approximation, the orbits that start with inclinations $i^{\prime} \lt 90^\circ$ with respect to the Galactic plane cannot swap to $i^{\prime} \gt 90^\circ$ (e.g., Naoz et al. 2013). The highest inclination that the inner Oort cloud orbits can reach with respect to the solar system plane is thus $\simeq 90^\circ +60^\circ$, where $\simeq 60^\circ$ is the angle between the Galactic and solar system planes, or $\simeq 150^\circ$ in total. Issue (2) is also related to Kozai dynamics. The concentration for $i\sim 30^\circ$ appears because the inner Oort cloud orbits spend the most time with $i^{\prime} \simeq 90^\circ$ (e.g., Higuchi et al. 2007), meaning that they are perpendicular to the Galactic plane.

The orbital features discussed above are smeared when it is accounted for stellar encounters, but they are still visible in Figures 4(c)–(f). With ${\rho }_{0}=0.2$ ${M}_{\odot }$ pc⁻³ (simulation C1G2S; panels e and f), the inner Oort cloud extends to slightly lower semimajor axes, but its overall structure remains the same. We find that ≃6.5% of the original disk bodies starting at 22–30 au at t = 0 end up in the Oort cloud ($a\gt 1000$ au) at $t=4.5\,\mathrm{Gyr}$. This is similar to the results of BM13 and somewhat higher than estimates obtained in previous dynamical models (3%–5% e.g., Dones et al. 2004; Kaib & Quinn 2008; Brasser et al. 2010; Kaib et al. 2011). With ${M}_{\mathrm{disk}}=20$ ${M}_{\oplus }$, we therefore find that the total mass of today's Oort cloud should be ∼1.3 ${M}_{\oplus }$. Of this, roughly 60% should be in the inner Oort cloud ($1000\lt a\lt {\rm{20,000}}$ au) and ≃40% in the outer Oort cloud ($a\gt {\rm{20,000}}$ au). The Oort-cloud-to-scattered-disk ratio obtained in our simulations is found to be ≃20, while BM13 reported $\simeq 12$ from their simulations.

The orbital structure of the trans-Neptunian region dramatically changes when P9 is included in the model. To start with, we first discuss our P9 models without the Galactic tide or stellar encounters (Figure 5). The dynamical effects of P9, mainly the Kozai resonance, act to decouple SDOs from Neptune and produce a nearly isotropic cloud of bodies roughly centered at P9's semimajor axis location. In the following, we call this hypothetical structure the P9 cloud.

Figure 5 shows how the orbital structure of the P9 cloud depends on the orbital parameters of P9. We find that $\simeq 1.7\times {10}^{4}$ bodies end up in the P9 cloud ($200\lt a\lt 1000$ au) at $t=4.5\,\mathrm{Gyr}$, corresponding to 1.7% of the original 10⁶ disk bodies at t = 0, or $\simeq 0.34$ ${M}_{\oplus }$ for ${M}_{\mathrm{disk}}=20$ ${M}_{\oplus }$. If real, the P9 cloud would represent a population ∼5 times larger than the classical KB and scattered disk below 200 au combined. The number of inner SDOs with $50\lt a\lt 200\,\mathrm{au}$, and the number and orbital structure of the classical KB are not affected by P9.

Figure 6 summarizes the structure of the trans-Neptunian region in different models. Without any external perturbations, only the scattered disk is present (panels a and b), and the orbits of extreme SDOs, such as Sedna and 2012 VP113, are not obtained in the model. With P9 (panels c–f), the P9 cloud forms with an estimated mass of ≃0.3–0.4 ${M}_{\oplus }$. This would provide an explanation for the high-q orbits of Sedna and 2012 VP11. Shankman et al. (2017), however, claimed that the detection of known extreme KBOs would imply a very massive P9 cloud (tens of ${M}_{\oplus }$). This exceeds, by roughly two orders of magnitude, the P9-cloud masses inferred from our dynamical modeling. From the simulations with P9, we find that orbits similar to Sedna and 2012 VP11 have perihelion longitudes ϖ concentrated near $\varpi -{\varpi }_{9}=180^\circ$ (Batygin & Brown 2016a). This concentration, however, is not strong enough, at least for the P9 parameters investigated here, to explain the current observations. We fail to identify any anisotropy in the distribution of nodal longitudes, Ω, and perihelion arguments, ω.

The orbital distribution obtained in the C1ALL model with P9 (${M}_{9}=15$ ${M}_{\oplus }$), Galactic tide (${\rho }_{0}=0.1$ M_⊙ pc⁻³), and stellar encounters is shown in Figures 6(e) and (f). Both the P9 and Oort clouds form in this model with an approximate division between them at $a\simeq 3000\,\mathrm{au}$. We find that the populations of the classical Kuiper Belt ($a\lt 50$ au), inner scattered disk ($50\lt a\lt 200$ au), and P9 cloud ($200\lt a\lt 3000$ au) represent fractions of $1.3\times {10}^{-3}$, $3.2\times {10}^{-3}$, and 0.017, respectively, which is very similar to the fractions reported for the other models above. The Oort cloud population ($a\gt 3000$ au) in C1ALL is somewhat smaller than in the models without P9, representing a fraction of 0.044 of the original disk (while we found a larger fraction of 0.060 for $a\gt 3000\,\mathrm{au}$ in the C1G1S model). This probably means that the presence of P9 makes it somewhat more difficult for bodies to reach the Oort cloud.

Using the methods described in Sections 3.7 and 3.8, we determined the orbital distribution of ECs in our models. Here we first discuss the results obtained without P9. Figure 7 shows the best result from the C1G1S model (Table 1). To limit the effect of observational biases discussed in Section 2, here we considered comets with $P\lt 20\,\mathrm{yr}$, $2\lt {T}_{{\rm{J}}}\lt 3$, $q\lt 2.5\,\mathrm{au}$, and ${H}_{{\rm{T}}}\lt 10$. The best fit shown in Figure 7 was obtained with ${N}_{{\rm{p}}}(2.5)=500$ (Section 3.8). It turns out that similarly good fits can be obtained with other parametrizations of the physical lifetime described in Section 3.8 (e.g., $T(2.5)=400\,\mathrm{yr}$ or ${T}_{{\rm{e}}}(2.5)=100$ yr). A realistic range of ${N}_{{\rm{p}}}(2.5)$, as determined by the Kolmogorov–Smirnov (K–S) test, is 300–800. The models with ${N}_{{\rm{p}}}(2.5)\lt 300$ or ${N}_{{\rm{p}}}(2.5)\gt 800$ do not fit the observed inclination distribution of ECs. The model distributions are narrower for ${N}_{{\rm{p}}}(2.5)\lt 300$ and broader for ${N}_{{\rm{p}}}(2.5)\gt 800$ than the observed distribution.

Zoom In Zoom Out Reset image size

The physical lifetime of comets is consistent with the results of Di Sisto et al. (2009) who found ${N}_{{\rm{p}}}(2.5)=300\mbox{--}450$ (and ${N}_{{\rm{p}}}(1.5)=170\mbox{--}200$). Our results are also consistent with LD97, where the physical lifetime of ECs was parametrized by the time, ${T}_{\mathrm{act}}$, during which a comet remained active after first becoming visible (i.e., after first reaching $q\lt 2.5$ au). Specifically, LD97 found that ${T}_{\mathrm{act}}\simeq 12$,000 yr was required to fit the inclination distribution of ECs (see also BM13). Since, according to Figure 7(a), the median orbital period of ECs is $\simeq 8\,\mathrm{yr}$, ${N}_{{\rm{p}}}(2.5)=500$ implies ${T}_{\mathrm{act}}\simeq 4000$ yr. This is a factor of $\simeq 3$ shorter than LD97's best estimate of ${T}_{\mathrm{act}}$, mainly because the source of ECs in our model is the scattered disk with a wide inclination distribution (while LD97 considered the classical KB with $i\lt 5^\circ$). The new ECs, reaching orbits with $q\lt 2.5\,\mathrm{au}$ for the first time, thus have a slightly wider inclination distribution in our model than in LD97. This implies shorter ${T}_{\mathrm{act}}$.¹¹

The best-fit result shown in Figure 7 is very good. We do not need to invoke observational biases to obtain a good fit. This is satisfactory, because it leaves ${N}_{{\rm{p}}}(2.5)$ (or, equivalently, $T(2.5)$ or ${T}_{{\rm{e}}}(2.5)$) as the only significant free parameter that needs to be adjusted in the model. The orbital distribution of ECs is independent of the timescale of Neptune's migration (models C1 and C2 produce similar results) and of whether Galactic tide or stellar encounters are included in the model (models C1, C1G1, and C1G1S produce the same result). Our model is thus identified as the simplest physical/dynamical model that is capable of matching the orbital distribution of active ECs. Other, more elaborate physical models have been developed in the past (e.g., Di Sisto et al. 2009; Rickman et al. 2017), but these models have more parameters and are more difficult to constrain.

Note that our physical model must be, to some degree, unrealistic, because many known active ECs have $q\gt 2.5\,\mathrm{au}$. It is therefore not true that ECs can be active only when $q\lt 2.5\,\mathrm{au}$. When we consider ECs with $q\gt 2.5\,\mathrm{au}$, however, we immediately run into a problem with observational biases. First, many ECs with $q\gt 2.5\,\mathrm{au}$ are probably undetected, even if they become active, because they appear too faint for a terrestrial observer. Second, only a fraction of ECs probably become active when reaching, say, $2.5\lt q\lt 5\,\mathrm{au}$. Accounting for the observational incompleteness is tricky, and we do not feel confident that expanding the model in this direction would produce meaningful results. Still, for the sake of argument, we attempted to match the observed distribution of active ECs with $q\lt 5\,\mathrm{au}$. We found that acceptable fits can be found, for example, with ${N}_{{\rm{p}}}(2.5)=500$, and assuming that all comets with $q\lt 2.5\,\mathrm{au}$ are detected, while only a fraction ${(q/2.5)}^{-\gamma }$ of those with $2.5\lt q\lt 5\,\mathrm{au}$ are detected, where $\gamma \sim 5$. This would indicate that only ∼3% of ECs are detected when they reach $q\sim 5\,\mathrm{au}$, and this fraction becomes ∼100% for $q\lesssim 2.5\,\mathrm{au}$.

Figure 8 compares the inclination distribution of ECs to those of SDOs and Centaurs. For SDOs, we used the the C1G1S simulation results at $t=4.5\,\mathrm{Gyr}$ and selected orbits with $50\lt a\lt 200\,\mathrm{au}$ and $q\lt 35\,\mathrm{au}$. This approximates the source region of ECs in our model (see Section 4.7). For Centaurs, we plotted the inclination distribution of known Centaurs from the JPL Small-Body Database Search Engine. Figure 8 shows that both the SDO and Centaur inclination distributions are significantly wider than the inclination distribution of ECs. This may seem surprising because the dynamical processes that mediate the delivery of ECs from the scattered disk, mainly the scattering encounters with planets, should act to increase the orbital inclinations and not to decrease them. By testing this, we found that the handover of bodies from the Neptune-crossing orbits toward Jupiter favors orbits with the Tisserand parameter with respect to Neptune, ${T}_{{\rm{N}}}$, somewhat smaller, but not much smaller, than 3. This naturally selects the low-inclination orbits (see discussion in LD97). In addition, $2\lt {T}_{{\rm{J}}}\lt 3$, used here to define ECs, also favors the low-inclination orbits. Both these effects therefore contribute to create an unfamiliar situation, where the inclination distribution of the target population (ECs) is narrower than that of the source (SDOs).

Zoom In Zoom Out Reset image size

Figure 9 shows the orbital distribution of ECs obtained in the model with P9. This is the best result that we were able to obtain with P9 in the C1M15 simulation. Other simulations with P9 produced similar results. The fit is not as good as the one in Figure 7, because in this case the inclination distribution of model ECs is somewhat broader than the inclination distribution of real ECs.¹² We applied the K–S test to understand how significant the difference is. Because, as we explained in Section 2, the inclination distribution of known ECs is not sensitive to the ${H}_{{\rm{T}}}$ cutoff, here we did not use any ${H}_{{\rm{T}}}$ cutoff to maximize the statistics. We found that the K–S probability with P9 is ${p}_{{\rm{K}}-{\rm{S}}}=0.008$, which is to be compared to ${p}_{{\rm{K}}-{\rm{S}}}=0.91$ obtained in the model without P9.

Zoom In Zoom Out Reset image size

The difference therefore seems to be significant, indicating that ECs could be used to provide a useful constraint on P9. Related to that, we note that the difference of the inclination distributions in Figure 9 is somewhat diminished by selecting orbits with $2\lt {T}_{{\rm{J}}}\lt 3$. If, instead, we compare cometary populations with $0\lt {T}_{{\rm{J}}}\lt 3$ (and $P\lt 20\,\mathrm{yr}$ and $q\lt 2.5\,\mathrm{au}$, as usual), both the real and model inclination distributions become broader, but the discrepancy becomes more significant, because the model distribution with P9 has many high-i orbits with $P\lt 20\,\mathrm{yr}$. This issue cannot be resolved by considering ${N}_{{\rm{p}}}(2.5)\lt 300$. This is because, with ${N}_{{\rm{p}}}(2.5)\lt 300$, the model distribution with P9 has a different profile from the observed distribution (Figure 10). We found that ${p}_{{\rm{K}}-{\rm{S}}}=6\times {10}^{-7}$ for the model with P9, $0\lt {T}_{{\rm{J}}}\lt 3$, and ${N}_{{\rm{p}}}(2.5)=100$, while ${p}_{{\rm{K}}-{\rm{S}}}=0.89$ for the model without P9, $0\lt {T}_{{\rm{J}}}\lt 3$, and ${N}_{{\rm{p}}}(2.5)=500$. Adding to that, with ${N}_{{\rm{p}}}(2.5)\sim 100$, there are not enough active/visible ECs in the model to explain the observed population (Section 4.6).

Zoom In Zoom Out Reset image size

The problems with P9 discussed above are related to the fact that the scattering disk at $50\lt a\lt 200\,\mathrm{au}$, which is the main source reservoir of ECs (Section 4.7), has significantly larger inclinations than in the model without P9. This happens because many SDOs interact with P9, with their orbital inclination being excited, and then return into the scattering disk with $a\lt 200\,\mathrm{au}$. Specifically, we find that all our simulations without P9 show similar inclination distributions of inner SDOs ($50\lt a\lt 200$ au) with 20%–25% of scattering orbits having $i\gt 30^\circ$, and only 5%–7% of scattering orbits having $i\gt 40^\circ$. Thus, the scattering disk without P9 is relatively flat. With P9, instead, roughly 60% of scattering orbits with $50\lt a\lt 200\,\mathrm{au}$ have $i\gt 30^\circ$ and roughly 50% of scattering objects with $50\lt a\lt 200\,\mathrm{au}$ have $i\gt 40^\circ$. The scattering disk with P9 is thus apparently puffed up by the dynamical effects of P9 (see Figure 5). This is reflected by the broad inclination distribution of ECs obtained with P9.

There are several possible solutions to this problem, some of which we were able to rule out. For example, we tested P9 with zero orbital inclination with respect to the invariant plane of the solar system, and found that the inclination distribution of ECs obtained in this model (C1I0; Table 1) is practically the same as in models with ${i}_{9}\gt 0$. We also verified that the same results were obtained when we used a shorter integration time step. Another possibility would be to consider P9 on an orbit with ${q}_{9}\gt 300\,\mathrm{au}$, such that the effect of P9 on inner SDOs is diminished ($200\lt {q}_{9}\lt 300\,\mathrm{au}$ in all models investigated here), and/or P9 with a lower mass. It is not clear, however, whether these cases could match other constraints such as the orbital alignment of extreme KBOs, solar obliquity, etc. A detailed investigation of this is beyond the scope of this work. Here, we found that cases with ${M}_{9}\lt 15$ ${M}_{\oplus }$ (and $200\lt {q}_{9}\lt 300$ au) did not produce a sufficiently strong orbital alignment of extreme KBOs, and cases with ${M}_{9}\gtrsim 15$ ${M}_{\oplus }$ produced a plausible alignment in ϖ, but not in Ω or ω.

Figure 11 shows the orbital distribution of HTCs obtained in the C1G1S model (Galactic tide and stellar encounters included, no P9). HTCs are produced from the Oort cloud in this model. They are a low orbital period extension of the returning OCCs (e.g., Nurmi et al. 2002). The range of orbital parameters in Figure 11 is restricted to $q\lt 2\,\mathrm{au}$ and $a\lt 20\,\mathrm{au}$, such that we avoid issues with the observational incompleteness of known HTCs with $q\gt 2\,\mathrm{au}$ and/or $a\gt 20\,\mathrm{au}$ (Section 2). In addition, here we only consider the HTC orbits with $a\gt 10\,\mathrm{au}$ to limit the potential pollution of the sample from ECs (both HTCs and ECs are considered in the next section). While the restricted range of orbital elements is probably the best-characterized part of the HTC population, the orbital distributions within this range may still be affected by observational biases. Thus, as a word of caution, we note that the comparison of model results with observations in this section is subject to some uncertainty.

Zoom In Zoom Out Reset image size

With these precautions in mind, Figure 11 appears to show a relatively good agreement. HTCs produced from the Oort cloud have a nearly isotropic inclination distribution with a slight preference for prograde orbits. The model distributions of a, q, and ${T}_{{\rm{J}}}$ look good as well. To obtain these results, we assumed ${N}_{{\rm{p}}}(2.5)=3000$, which is a factor of $\simeq 6$ higher than what was needed to fit the inclination distribution of ECs. For HTCs, the orbital distributions are not very sensitive to ${N}_{{\rm{p}}}(2.5)$, and ${N}_{{\rm{p}}}(2.5)=500$ gives qualitatively similar results to those shown in Figure 11. The value of ${N}_{{\rm{p}}}(2.5)$ is thus not very well constrained by the fit to the observed orbital distribution of HTCs. Instead, ${N}_{{\rm{p}}}(2.5)\gt 1000$ is driven by the requirement to produce a number of active HTCs that is consistent with observations (to be discussed in Section 4.6).

We tested different parametrizations of the physical lifetime of comets described in Section 3.8 and found that they do not help solve this problem. There are several other possibilities. (1) For some reason, ${N}_{{\rm{p}}}(2.5)\simeq 500$ derived from the fit to the inclination distribution of ECs is too low. For example, our scattered disk (in the simulations without P9) may be too excited in inclinations, which could then drive ${N}_{{\rm{p}}}(2.5)$ to low values (because the new ECs would already have a broad inclination distribution). It is hard to imagine that this might be the case, because the scattered disk is gradually excited by encounters of SDOs with Neptune. The excitation therefore does not depend on some simulation detail. (2) The number of HTCs obtained in our model is too low, by a factor of several. This possibility is discussed in Section 5 together with our preferred resolution of this problem, where the physical lifetime of a comet is a function of the comet's size.

Figure 12 shows the orbital distributions of HTCs obtained in the C1M15 model (P9 with ${M}_{9}=15$ ${M}_{\oplus }$, and no Galactic tide or stellar encounters; Table 1). The results obtained with other parameters of P9 were similar. In this model, HTCs are produced from the P9 cloud that is roughly centered at the semimajor axis of P9 (Figure 5). The delivery of HTCs from the P9 cloud is a two-step process. First, the secular effects of P9 act to decrease the perihelion distance of an object in the P9 cloud. Subsequently, when $q\lt 30\,\mathrm{au}$, the orbital period of SDOs can be shortened by the gravitational effects of the Jovian planets. Since the period of secular cycles in the P9 cloud is >100 Myr, the first step is a very slow process. This is an important difference with respect to the delivery of HTCs from the Oort cloud, where bodies can be placed on orbits with very low perihelia in one orbital period.

Zoom In Zoom Out Reset image size

The HTC population obtained from the P9 cloud model shows similarities to the observed population (Figure 12), but it does not fit the orbital distribution as well as the Oort cloud model (Figure 11). The inclination distribution of model HTCs in Figure 12(b) is very nearly isotropic with a small preference for retrograde orbits (median inclination ≃100° ). While this preference does not seem to be reflected in the existing observational data, it cannot be ruled out either. The perihelion distance distribution of model HTCs does not fit the data too well, showing a convex profile with the number of orbits below q proportional to q². The observed distribution is flatter. This may be a consequence of larger observational incompleteness for orbits with higher perihelion distances.

So far we have discussed the population of HTCs with $10\lt a\lt 20\,\mathrm{au}$ and $q\lt 2\,\mathrm{au}$. This is because this part of orbital space should be best characterized from observations (Section 2). In a recent paper, Fernández et al. (2016) opted to use the full range $7.4\lt a\lt 34.2\,\mathrm{au}$ (corresponding to $20\lt P\lt 200$ yr) and $q\lt 1.3\,\mathrm{au}$ (16 known comets), and compared the orbital distribution of HTCs to those obtained from LPCs and Centaurs. They argued that HTCs should have had at least one perihelion passage in modern history and should thus have a good chance of being detected. They found that the distribution of orbital energy (or equivalently, of semimajor axis) obtained from LPCs does not fit the distribution of HTCs. Instead, they argued that the immediate source of HTCs are Centaurs. Here we confirm that HTC orbits evolving from the Oort cloud have a cumulative semimajor axis distribution $N(\lt \ a)\propto {a}^{2}$, while HTCs from the P9 cloud would have a flatter distribution (that better fits the known HTCs with $20\lt P\lt 200\,\mathrm{yr}$ and $q\lt 1.3$ au). A careful characterization of the HTC population with $a\gt 20\,\mathrm{au}$ will be needed before these arguments can be placed on firmer ground.

In summary, we find that both the Oort and P9 clouds are potential sources of HTCs, but the Oort cloud model fits the existing orbital data of HTCs better than the P9 cloud model. As we will discuss in Section 4.6, the Oort cloud is a more prolific source of HTCs than the P9 cloud (by a factor of ∼2–3). This shows that P9 is not required from the considerations based on HTCs. On the other hand, inclusion of P9 in a model does not harm the orbital distribution of HTCs.

Here we remove the distinction between ECs and HTCs and attempt to fit the orbits of all SPCs together. Figure 13 shows the result for the C1G1S simulation (no P9) and ${N}_{{\rm{p}}}(2.5)=500$, which we established in Section 4.1 to be the best-fit value for ECs. The model distributions in Figure 13 have profiles similar to the observed distributions but do not fit them very well. The distribution in panel (a) can be interpreted as evidence that we have too many HTCs, relative to the EC population, in our model. This would be surprising, however, because with ${N}_{{\rm{p}}}(2.5)=500$, the number of HTCs is severely reduced (Section 4.6). If ${N}_{{\rm{p}}}(2.5)=3000$ instead, which is the preferred value to obtain the right number of HTCs, the problem in Figure 13(a) would appear to be much worse.

Zoom In Zoom Out Reset image size

We believe that this problem is related to the observational incompleteness of comets with long orbital periods. We find that the semimajor axis distributions in Figure 13 would match when it is assumed that ∼70% of HTCs with $8\lt a\lt 20\,\mathrm{au}$ have been discovered so far (while the population of ECs is assumed to be nearly complete). Note that, however, because of the issues with ${N}_{{\rm{p}}}$ discussed above, the intrinsic population of HTCs may in reality be larger than shown in Figure 13. If so, a larger incompleteness of HTCs would need to be invoked to bring the model into agreement with observations (but see discussion in Section 5, where we argue that ${N}_{{\rm{p}}}$ is a function of comet size and ${N}_{{\rm{p}}}(2.5)\sim 500$ should apply to kilometer-sized comets in general).

The simulations presented here allow us to link the number of ECs and HTCs to the number of planetesimals in the original trans-Neptunian disk below 30 au. This has not been done before, except by BM13, at least not in a self-consistent dynamical model that was also shown to reproduce many properties of other small-body populations in the solar system. In previous works, ECs and HTCs were considered separately and different schemes were developed to deal with different comet categories. In some cases, the number of ECs was linked, through a chain of multiplicative factors, to the number of objects in the present scattered disk. In other cases, the number of HTCs was calibrated by the number of observed new LPCs (e.g., Rickman et al. 2017). While all of these works have their own merits, here we prefer to emphasize the link to the original planetesimal disk, which presumably is the common source of ECs and HTCs. This is done as follows.

NV16 calibrated the number of bodies in the original disk. For that, they assumed that the size distribution of disk planetesimals followed the size distribution of today's Jupiter Trojans, which is well-characterized down to at least $D\simeq 5$ km (Wong & Brown 2015; Yoshida & Terai 2017). This assumption is based on previous modeling efforts, which showed that Jupiter Trojans were implanted from the original planetesimal disk (Morbidelli et al. 2005; Nesvorný et al. 2013) and that the collisional evolution of Jupiter Trojans after their implantation was not strong enough to substantially modify their size distribution (e.g., Wong & Brown 2015). To absolutely calibrate the number of original disk planetesimals, NV16 used the estimate of the implantation probability determined in Nesvorný et al. (2013), who showed that a fraction of ≃7 × 10⁻⁷ of the disk planetesimals becomes implanted on stable Trojan orbits.

The uncertainty of this estimate is not well established. In the three simulations presented in Nesvorný et al. (2013), the implantation probability was found to vary by only ≃15%. On the other hand, our new simulations with very slow planetary migration rates reveal how the survival rate of Jupiter Trojans depends on the migration timescale. Some of the results with the longest migration timescales show probabilities as low as ≃3 × 10⁻⁷. Here we therefore choose to adopt the implantation probability 5 × 10⁻⁷, which is in the middle of the values discussed above, and use this value to calibrate the number of planetesimals in the original disk.

Additional constraints on the size distribution come from the mass of the original disk needed to generate the plausible dynamical evolution of the planetary system (${M}_{\mathrm{disk}}\simeq 20$ ${M}_{\oplus };$ NM12; Deienno et al. 2017), the expected number of Pluto-class objects in the original disk (${N}_{\mathrm{Pluto}}\,=1000\mbox{--}4000$; NV16), and various Kuiper Belt constraints (see NV16). Figure 14 shows the reconstructed cumulative size distribution of disk planetesimals. This figure indicates that there were approximately $6\times {10}^{9}$ disk planetesimals with $D\gt 10$ km. This estimate is uncertain by a factor of ∼2, mainly due to the uncertainty in the implantation probability of Jupiter Trojans and its dependence on planetary migration.

Zoom In Zoom Out Reset image size

The number of comets expected in a steady state, ${N}_{\mathrm{com}}(\gt D)$, can be computed from

$$\begin{eqnarray}&&{N}_{\mathrm{com}}(\gt D)={N}_{\mathrm{rec}}\displaystyle \frac{{N}_{\mathrm{disk}}(\gt D)}{{N}_{\mathrm{sim}}}\displaystyle \frac{{\rm{\Delta }}t}{{\rm{\Delta }}T},\end{eqnarray} \tag{ 3 }$$

where ${N}_{\mathrm{rec}}$ is the number of cometary orbits recorded in ${\rm{\Delta }}T$, ${N}_{\mathrm{disk}}(\gt D)$ is the number of original disk planetesimals larger than D, ${N}_{\mathrm{sim}}={10}^{8}$ stands for the number of bodies used in our simulations (10⁶ original bodies times 100 for cloning), ${\rm{\Delta }}t=100\,\mathrm{yr}$ is the sampling interval, and ${\rm{\Delta }}T$ is the time interval used to accumulate good statistics. Here we use ${\rm{\Delta }}T=1\,\mathrm{Gyr}$ (${\rm{\Delta }}T\lt 1\,\mathrm{Gyr}$ leads to similar results but the statistics are worse). Note that ${N}_{\mathrm{rec}}$ depends on the assumed physical lifetime of comets.

In the following text, we compare ${N}_{\mathrm{com}}$ with the number of known comets with $D\gt 10$ km. There are several reasons behind this choice, perhaps the most important being that comets with small nuclei are probably more difficult to detect than those with $D\gt 10$ km. The known population of small comets is therefore incomplete and biased in uncertain ways. In previous work, the total absolute magnitude ${H}_{{\rm{T}}}$ was often taken as a proxy for the nuclear size of a comet, but we do not find any correlation when the nuclear size of the comets determined in Fernández et al. (2013) is compared to their ${H}_{{\rm{T}}}$ (or their nuclear magnitude) reported at the JPL site. We therefore believe that attempts to relate ${H}_{{\rm{T}}}$ to the nuclear size may be misguided.

We searched various catalogs to determine the number of known ECs and HTCs with $D\gt 10$ km. There are two large ECs listed in the JPL database: 10P/Tempel 2 (D = 10.6 km) and 28P/Neujmin 1 (D = 21.4 km) (see Lamy et al. 2004 for a discussion). In addition, Fernández et al. (2013) reported four additional ECs with $D\gt 10$ km, two of which have $q\lt 2.5\,\mathrm{au}$ at the present time. These are 162P/Siding Spring (D = 14.1 km) and 315P/LONEOS (D = 10.8 km).¹³ So, together, there appear to be four known ECs with $D\gt 10$ km and $q\lt 2.5\,\mathrm{au}$.

It is not known how solid this estimate is. On one hand, some size estimates were obtained from an assumed visual albedo (often taken as low as ∼2%). These determinations are less reliable than those derived from IR observations and thermal modeling. In addition, it is often assumed that the absolute nuclear magnitude can be determined from observations, either because the contribution of the coma is thought to be negligible (e.g., observations close to the orbital aphelion), or because the nucleus appears to be resolved. On the other hand, the known sample of ECs with $D\gt 10$ km and $q\lt 2.5\,\mathrm{au}$ may still be incomplete. Indeed, both 162P and 315P were only discovered in 2004.

As for HTCs, there are 1P/Halley (D = 11 km), C/1991 L3 Levy (D = 11.6 km), and C/2001 OG108 LONEOS (D = 13.6 km). The diameter of C/1991 L3 Levy was not measured in the thermal IR and is not reliable, while the (effective) diameters of 1P/Halley and C/2001 OG108 LONEOS are well established. 109P/Swift-Tuttle with D = 26 km has semimajor axis $a=26.1\,\mathrm{au}$ and is outside the range considered here ($10\lt a\lt 20$ au). We conclude that there are ≃2–3 known HTCs with $D\gt 10$ km, $q\lt 2\,\mathrm{au}$, and $a\lt 20\,\mathrm{au}$. Again, it is not clear how complete this sample is, but it should probably be more incomplete than ECs. The large HTCs may therefore be as common as large ECs, if not even more common.

Using Equation (3), we find from our simulations without P9 that ${N}_{\mathrm{EC}}=1\mbox{--}2$ for $D\gt 10$ km, $2\lt {T}_{{\rm{J}}}\lt 3$, $P\lt 20\,\mathrm{yr}$, $q\lt 2.5\,\mathrm{au}$, and $300\lt {N}_{{\rm{p}}}(2.5)\lt 800$ (i.e., the range of ${N}_{{\rm{p}}}(2.5)$ required to fit the inclination distribution; Figure 15(a) shows how ${N}_{\mathrm{EC}}$ depends on ${N}_{{\rm{p}}}(2.5)$).¹⁴ This is somewhat lower than the number of large ECs discussed above, thus indicating that our model may be anemic, by a factor of ∼2–4, when compared to observations. This factor would be larger if the observational incompleteness of large ECs is significant. With P9 and $300\lt {N}_{{\rm{p}}}(2.5)\lt 800$, we obtain ${N}_{\mathrm{EC}}=0.7\mbox{--}1$ for $D\gt 10$ km, $2\lt {T}_{{\rm{J}}}\lt 3$, $P\lt 20\,\mathrm{yr}$, and $q\lt 2.5\,\mathrm{au}$. The number of ECs in a model with P9 is thus somewhat lower than in a model without P9. This is probably related to a larger excitation of the orbital inclinations of SDOs when P9 is present (see discussion in Section 4.3).

Zoom In Zoom Out Reset image size

As for the HTCs, we find that the Oort cloud should produce $\simeq 2.3$ HTCs in a steady state with $D\gt 10$ km, $10\lt a\lt 20\,\mathrm{au}$, ${T}_{{\rm{J}}}\lt 2$, and $q\lt 2\,\mathrm{au}$. This is right in the middle of the range inferred from existing observations above, if the incompleteness of the existing sample could be ignored. This estimate was obtained for ${N}_{{\rm{p}}}(2.5)=3000$ and model C1G1S. The number of large HTCs obtained in our other Oort-cloud models is similar (2.7 in C1G1 and 1.7 in C1G2S). If ${N}_{{\rm{p}}}(2.5)\lesssim 1000$ is assumed instead, then ${N}_{\mathrm{HTC}}\lesssim 1$ (Figure 15(b)), at least $\simeq 2$ times below the value indicated by observations.

The P9 cloud is less efficient in producing HTCs than the Oort cloud. In particular, we find that ${N}_{\mathrm{HTC}}=0.9$ with $D\gt 10$ km, $10\lt a\lt 20\,\mathrm{au}$, ${T}_{{\rm{J}}}\lt 2$, and $q\lt 2\,\mathrm{au}$ in the C1M15 model and ${N}_{{\rm{p}}}(2.5)=3000$. The population estimates obtained in other models with P9 are similar. This is a factor of $\simeq 2.5$ smaller than the number of HTCs obtained from the Oort cloud. The comparison of different models therefore shows that most HTCs should be coming from the Oort cloud, and the contribution of P9, if real, should be relatively minor. In the C1ALL model, where both P9 and the Oort cloud contribute to the population of HTCs, ${N}_{\mathrm{HTC}}=2.1$ for ${N}_{{\rm{p}}}(2.5)=3000$, which is very similar to the estimates obtained without P9. There are fewer HTCs coming from the Oort cloud in the C1ALL model, because the Oort cloud population is smaller (Section 4.1), but the contribution from the P9 cloud nearly compensates for the difference.

We used our simulation results to characterize the source reservoirs of SPCs. For that, we selected ECs with $2\lt {T}_{{\rm{J}}}\lt 3$, $P\lt 20\,\mathrm{yr}$, $q\lt 2.5\,\mathrm{au}$, and ${N}_{{\rm{p}}}(2.5)=500$ and HTCs with ${T}_{{\rm{J}}}\lt 2$, $10\lt a\lt 20\,\mathrm{au}$, $q\lt 2\,\mathrm{au}$, and ${N}_{{\rm{p}}}(2.5)=3000$. The source orbits of selected ECs and HTCs in the C1G1S simulation are shown in Figure 16. For ECs, the orbits are shown at $t=1.5\,\mathrm{Gyr}$ after the start of the C1G1S simulation, or roughly 3 Gyr ago. The migration phase of Neptune has ended at this point. For HTCs, we prefer to plot their orbits at $t=3.5\,\mathrm{Gyr}$, or roughly 1 Gyr ago. This is because the orbital structure of the Oort cloud, which is the main source reservoir of HTCs, continues to evolve over gigayears.

Zoom In Zoom Out Reset image size

Most ECs (≃75%) were produced from the scattered disk with $50\lt a\lt 200\,\mathrm{au}$ (Figure 17). About 20% of ECs started with $a\lt 50\,\mathrm{au}$. Of these, most bodies had stable orbits that remained with $a\lt 50\,\mathrm{au}$ from $t=1.5\,\mathrm{Gyr}$ to $t=3.5\,\mathrm{Gyr}$. The classical KB, including various resonant populations below 50 au (about 4% of ECs evolved from the Plutino population in 3:2 resonance with Neptune), is therefore a relatively important source of ECs. Interestingly, ≃3% of our model ECs started in the Oort cloud (see also Emel'yanenko et al. 2013). The orbital evolution of these comets is similar to returning LPCs or HTCs, except that they were able to reach orbits with very low orbital periods and low inclinations. The median semimajor axis of source EC orbits is $\simeq 60\,\mathrm{au}$. The median inclination of source EC orbits 1 Gyr ago was ≃20°.

Zoom In Zoom Out Reset image size

Figures 16(c) and (d) show that a great majority (≃95%) of HTCs in C1G1S come from the Oort cloud, and only ≃5% from the $a\lt 100\,\mathrm{au}$ region. The inclination distribution of source orbits is slightly anisotropic with the median inclination $\simeq 70^\circ$ (Figure 18(c)). This is similar to the median inclination of new HTCs in our simulations. The inner and outer parts of the Oort cloud, as defined in Section 4.1 ($1000\lt a\lt {\rm{20,000}}\,\mathrm{au}$ inner, $a\gt {\rm{20,000}}\,\mathrm{au}$ outer), contribute in nearly equal proportions to the HTC population. We see some exchange of orbits between the inner and outer Oort clouds in our simulations. The partition of source orbits therefore depends on the time when the source orbits are extracted. For example, if the orbits are extracted at $t=0.5\,\mathrm{Gyr}$, or roughly 4 Gyr ago, then we find that ≃70% of current-day HTCs started in the inner Oort cloud (see also Kaib & Quinn 2009).

Zoom In Zoom Out Reset image size

We can now estimate the number of bodies in the source reservoirs. Because of the uncertain nature of P9, we limit the following discussion to our models without P9. Summarizing the findings discussed in the previous sections, we found that the inner SD ($50\lt a\lt 200$ au) and Oort cloud ($a\gt {\rm{10,000}}$ au) represent the fractions $\simeq 2.5\times {10}^{-3}$ and $\simeq 0.065$ of the original planetesimal disk. The numbers of current-day ECs and HTCs in steady state are the fractions $\sim 2.5\times {10}^{-10}$ and $\sim 4.2\times {10}^{-10}$ of the original disk, respectively. The quoted fractions apply to active ECs (${N}_{{\rm{p}}}(2.5)=500$) on orbits with $P\lt 20\,\mathrm{yr}$, $2\lt {T}_{{\rm{J}}}\lt 3$, and $q\lt 2.5\,\mathrm{au}$, and to active HTCs (${N}_{{\rm{p}}}(2.5)=3000$) on orbits with $10\lt a\lt 20\,\mathrm{au}$, ${T}_{{\rm{J}}}\lt 2$, and $q\lt 2\,\mathrm{au}$. If ${N}_{{\rm{p}}}(2.5)=3000$ is assumed for large ECs instead (see discussion in Section 5), the fraction becomes $\sim 6.7\times {10}^{-10}$.

From these, we estimate that the ratio of active ECs with $q\lt 2.5\,\mathrm{au}$ to inner SDOs is $\sim 2.5\times {10}^{-10}/2.5\times {10}^{-3}\,=1.0\times {10}^{-7}$ for ${N}_{{\rm{p}}}(2.5)=500$ or $\sim 2.7\times {10}^{-7}$ for ${N}_{{\rm{p}}}(2.5)\,=3000$. The former value is more similar to the fraction $6.7\times {10}^{-8}$ obtained in BM13 for ${T}_{\mathrm{act}}=12$,000 yr. The ratio of active HTCs ($q\lt 2\,\mathrm{au}$, $10\lt a\lt 20$ au) to Oort cloud bodies is $\sim 4.2\times {10}^{-10}/0.065=6.5\times {10}^{-9}$ for ${N}_{{\rm{p}}}(2.5)\,=3000$. Since, as we discussed in Section 4.6, there are some four known active ECs with $D\gt 10$ km, there should be $\sim 4/1.0\times {10}^{-7}=4.0\times {10}^{7}$ $D\gt 10$ km bodies in the inner SD if ${N}_{{\rm{p}}}(2.5)=500$ or $\sim 1.5\times {10}^{7}$ $D\gt 10$ km in inner SDOs if ${N}_{{\rm{p}}}(2.5)=3000$ (most of these have detached orbits; Section 4.1). We prefer the latter estimate for reasons that will be explained in Section 5. Also, from two to three HTCs with $D\gt 10$ km (Section 4.6), we estimate that the Oort cloud should contain $\sim 2.5/6.5\times {10}^{-9}=3.8\times {10}^{8}$ $D\gt 10$ km comets.

According our work, the ratio of the Oort cloud to scattered disk should be $\mathrm{OC}/\mathrm{SD}\sim 20$ (Section 4.1). BM13 obtained $\mathrm{OC}/\mathrm{SD}=12\pm 1$ from their simulations based on the original Nice model, and inferred $\mathrm{OC}/\mathrm{SD}\sim 44$ from observations. The latter estimate has a large uncertainty mainly due to the uncertain size and number of new LPCs. For example, BM13 pointed out that the flux of new LPCs may be lower than assumed before, because some LPCs thought previously to be new are actually returning LPCs (Królikowska & Dybczyński 2010). They ended up giving preference to OC/SD ∼23, roughly in the middle of the values quoted above. Our new estimate, $\mathrm{OC}/\mathrm{SD}\sim 20$, is spot-on their preferred value.

In previous works, briefly discussed in Section 1, various estimates were given for the number of $D\gt 2$ km or $D\gt 2.3$ km bodies in the scattered disk and Oort cloud. To be able to compare with these works, we use the distribution shown in Figure 14, where the ratio of $D\gt 10$ km to $D\gt 2$ km bodies is ≃29, and the ratio of $D\gt 10$ km to $D\gt 2.3$ km bodies is ≃22. From this we find that there should be $\sim 4.4\times {10}^{8}$ $D\gt 2$ km bodies in the inner scattered disk, and $\sim 1.1\times {10}^{10}$ $D\gt 2$ km bodies in the Oort cloud. The former estimate is a factor of ∼2 lower than the one reported in Rickman et al. (2017), who found, by combining several factors from LD97 and other works, that the capture rate of JFCs requires $\sim {10}^{9}$ $D\gt 2$ km bodies in the scattered disk. Duncan & Levison (1997) reported $\sim 6\times {10}^{8}$ SDOs from their modeling of ECs (the size range to which this estimate applies was not specified), which is only slightly larger than our estimate for $D\gt 2$ km.

BM13 and Brasser & Wang (2015) obtained somewhat higher estimates: ∼2 × 10⁹ and ∼6 × 10⁹ SDOs with $D\gt 2.3$ km, respectively. These estimates are a factor of ∼6–18 higher than ours. In addition, BM13 found that the observed flux of new LPCs implies that there are $\sim 4\times {10}^{10}$ to ∼10¹¹ $D\gt 2.3$ km comets in the Oort cloud. These estimates are a factor of ∼5–12 higher than ours. Thus, while we agree with BM13 on the OC/SD ratio, for some reason, our best estimates are at least a factor of ∼5 lower.

Some of the differences quoted above can be explained by the uncertain relationship between total absolute magnitude and nuclear size. As we already mentioned, we do not find any correlation when we compare ${H}_{{\rm{T}}}$ from the JPL database with the nuclear diameter estimates from Fernández et al. (2013). This could mean that ${H}_{{\rm{T}}}$ expresses comet activity rather than the nuclear size, perhaps because only a small part of a comet's surface is typically active (e.g., Sosa & Fernández 2011).

If that is the case, it may be incorrect to use ${H}_{{\rm{T}}}$ as a proxy for size. Brasser & Wang (2015), for example, assumed that ${H}_{{\rm{T}}}=10.8$ corresponds to D = 2.3 km and estimated that there are $\simeq 300$ JFCs with $D\gt 2.3$ km and $q\lt 2.5\,\mathrm{au}$.¹⁵ Adopting these numbers and assuming that there are ∼4 JFCs with $D\gt 10$ km and $q\lt 2.5\,\mathrm{au}$ (Section 4.6), the cumulative power-law slope at $2\lt D\lt 10$ km would be $\sim -3$, which is much steeper than the $\sim -2$ typically found from observations (e.g., see Lamy et al. 2004 for a review). We therefore believe that the number of small JFCs is significantly lower, possibly $\gtrsim 4$ times lower, than the one estimated in Brasser & Wang (2015).

Here we calibrated the number of large $D\gt 10$ km SDOs and Oort cloud bodies from Jupiter Trojans and large ECs/HTCs for which the nuclear size is relatively well known from observations (e.g., thermal IR). These two calibrations are consistent in that they lead to the same population estimates. We then used the size distribution of Jupiter Trojans to extrapolate our estimates for $D\gt 10$ km to $D\gt 2$ km and $D\gt 2.3$ km. This method should provide more robust results than the previous works, because it circumvents the problems with the uncertain relationship between ${H}_{{\rm{T}}}$ and nuclear size.