© ESO, 2014

1. Introduction

The solar system is host to a large population of comets, which tend to be concentrated in two sources, i.e. the Oort Cloud and scattered disc (Duncan & Levison 1997). When examining the orbital data of new comets entering the inner solar system, Oort (1950) noted that the distribution of reciprocal semi-major axis showed a sharp increase when 1/a < 5 × 10^-4 AU^-1. The observed semi-major axis distribution led Oort to suggest that the Sun is surrounded by a cloud of comets in the region between 20 000 AU to 150 000 AU, and that it contains approximately 10¹¹ comets with isotropic inclination and random perihelia. Oort showed that the action of passing stars is able to account for the isotropic inclination distribution, random perihelia, the continuous injection of new comets on trajectories that lead to perihelia close enough to the Sun to become visible, and the fact that most new comets have near-parabolic orbits. This hypothesized cloud of comets surrounding the Sun is now called the Oort Cloud (OC).

The only method with which we can infer the number of comets in the OC is by determining the flux of long-period comets (LPCs) and Oort Cloud comets (OCCs). The former are comets with a period longer than 200 yr that have had multiple passages through the solar system while the latter are comets that are having their first perihelion passage in the inner solar system.

For a new comet to enter the planetary region from the OC, it must have sufficiently low angular momentum to allow it to penetrate to small perihelion distances. For a given comet, its angular momentum per unit mass is given by v_tr, where v_t is its transverse velocity and r is the distance from the Sun. In the OC, where r is large, v_t must be very low to allow the comet to eventually come close to the Sun. Oort (1950) recognized this and added that the maximum allowable v_t for a given distance r would be a function of r. He gave this the name “loss cone” because comets within the velocity cone could enter the planetary region and be ejected or captured to a shorter period orbit by the perturbations of the giant planets, principally Jupiter and Saturn (which would require a perihelion distance q < 15 AU). In other words, the comets would be lost from the OC. These lost comets remain in the loss cone until ejected or destroyed. Thus, the presence of Jupiter and Saturn presents an obstacle for new comets to overcome in order to enter the inner solar system. This is often called the Jupiter-Saturn barrier.

There are two agents that perturb the comets in the cloud onto orbits that enter the inner solar system: passing stars (Weissman 1980; Hills 1981) and the Galactic tide (Heisler & Tremaine 1986; Levison et al. 2001). The passing stars usually cause small random deviations in the orbital energy and other orbital elements. The Galactic tide on the other hand systematically modifies the angular momentum of the comets at constant orbital energy. Heisler & Tremaine (1986) discovered that if the semi-major axis of the comet is large enough then the comet’s change in perihelion can exceed 10 AU in a single orbit and jump across the orbits of Jupiter and Saturn and thus not suffer their perturbations before becoming visible; in this particular case the comet is considered a new comet. However, the Heisler & Tremaine (1986) approximation of the Galactic tide only used the vertical component, which is an order of magnitude stronger than the radial components. In this approximation the z-component of the comet’s orbital angular momentum in the Galactic plane is conserved and the comets follow closed trajectories in the q − ω plane (with q being the perihelion distance and ω being the argument of perihelion). Including the radial tides breaks this conservation and the flux of comets to the inner solar system from the OC is increased (Levison et al. 2006). The trajectories that lead comets into the inner solar system would be quickly depleted, were it not for the passing stars to refill them (Rickman et al. 2008). The synergy between these two perturbing agents ensures there is a roughly steady supply of comets entering the inner solar system. The story is slightly more complicated because not all comets need to jump over the Jupiter-Saturn barrier: a small subset is able to sneak through and undergo small enough energy changes that they still formally reside in the cloud (Królikowska & Dybczyński 2010; Fouchard et al. 2013), but the main argument still holds.

The OCCs form just a subset of comets that are visible from the Oort Cloud. The OCCs are generally thought to be on their first passage through the inner solar system, but may have had a few passages beyond ~5 AU. However, we also have observational evidence of comets that are dynamically evolved. These are the LPCs: they have undergone many passages through the realm of the giant planets and/or the inner solar system and most of them are dynamically decoupled from the Oort Cloud. A subset of these are the Halley-type comets, or HTCs, named after its progenitor, comet 1P/Halley.

The second reservoir of cometary objects is the scattered disc (Duncan & Levison 1997). This source population resides on dynamically unstable orbits beyond the orbit of Neptune. These bodies dynamically interact with Neptune and reside generally on low-inclination orbits. The scattered disc is the most likely source of the so-called Jupiter-family comets (JFCs; Duncan & Levison 1997; Volk & Malhotra 2008; Brasser & Morbidelli 2013).

Historically, visible comets tend to be categorised into several groups: the JFCs, HTCs, LPCs, and OCCs. Here we adhere to the definitions of Levison (1996) with a few modifications. These are:

• JFCs: T_J ∈ (2,3) and a < 7.35 AU (period P < 20 yr);

• HTCs: T_J < 2, q > 0.01 AU and a < 34.2 AU (period P < 200 yr);

• LPCs: T_J < 2, q > 0.01 AU and 34.2 < a < 10 000 AU (period P > 200 yr);

• OCCs: q > 0.01 AU and a > 10 000 AU.

Here a is the semi-major axis, P is the period of a comet, and T_J is the Tisserand parameter with respect to Jupiter. The criterion q > 0.01 AU was applied in order to exclude the Sun-grazing Kreutz comets. Dormant comets were also excluded in our work. We realize that the above populations are not all inclusive: we did not include the Encke type and the Chiron type. There is no definition for comets with T_J > 2 and P > 20 yr, although these are generally included in the Centaur population.

The origin of the HTCs has long been controversial. The best attempts to determine the source of the HTCs were done by Fernández & Gallardo (1994), Levison et al. (2001), Nurmi et al. (2002), and Levison et al. (2006).

Fernández & Gallardo (1994) investigated the dynamical behaviour of LPCs coming from the Oort Cloud and the role of perturbations by the planets. They used a combination of a Monte Carlo random walk method and Öpik’s equations (Öpik 1951) and discovered that an initially isotropic distribution may result in a flattened distribution of evolved comets because of the inclination dependence on planetary energy kicks and the tendency of near-perpendicular comets to become retrograde. They find that the typical comet undergoes 400 revolutions before becoming extinct, and combined with their capture probability they suggest there are about 300 active HTCs.

Levison et al. (2001) studied an Oort Cloud origin for the HTCs, but concluded that the non-isotropy of the HTCs could only be reproduced with an Oort Cloud that had a non-isotropic structure: it had to be flattened in the inner parts and isotropic in its outer parts. Levison et al. (2001) included cometary fading in their simulations, but they found that this alone was not enough to match the HTC inclination distribution. Nurmi et al. (2002) reached a different conclusion based on results of their Monte Carlo simulations: the outer Oort Cloud could be the source of the HTCs and their number and inclination distribution could be fairly well reproduced when cometary fading and disruption was taken into account. However, when adding the inner Oort Cloud they produced far too many HTCs, which Levison et al. (2001) suggested was needed in their model. Unfortunately Nurmi et al. (2002) did not try to match the semi-major axis distribution, something that Levison et al. (2001) did try to do. In addition, the initial conditions of both of these works were artificial because they started the comets on planet-crossing orbits.

These inconsistencies led Levison et al. (2006) to suggest that the HTCs originate from the outer edge of the scattered disc. This reservoir has a flattened inclination distribution (Levison & Duncan 1997) and the effects of the Galactic tide at the outer edge is able to supply the retrograde HTCs. Their simulations match the inclination and semi-major axis distribution of the HTCs reasonably well. The problem with this scenario is that it requires the scattered disc to be more than an order of magnitude more massive than most estimates thus far (Gomes et al. 2008; Brasser & Morbidelli 2013). This leads us back to the Oort Cloud being the most likely possible source.

Here we aim to demonstrate that the Oort Cloud is the most likely source of the HTCs and that a single dynamical model with cometary fading is sufficient to reproduce their current inclination and semi-major axis distributions.

This paper is organised as follows. In the next section we present the observational dataset that we use for comparison. Section 3 details our numerical methods. In Sect. 4 we present the results of our simulations and comparison with the observed dataset followed by a short discussion in Sect. 5. We present a summary and conclusions in the last section.

2. Observational dataset

In this study we determine whether or not the Oort Cloud is the major source of the HTCs. Even though the scattered disc and Kuiper belt do generate some HTCs (Levison & Duncan 1997), the production probability is low and their orbital distribution is inconsistent with the observed distribution. Thus, in this work we shall not consider the scattered disc and Kuiper belt as a source.

To compare our simulation with the observational data we need to have an observational catalogue of comets that is as complete as possible. There are two easily accessible sources of cometary catalogues: the Catalogue of Cometary Orbits 2008 (Marsden & Williams 2008) and the JPL Small-Body Search Engine¹.

The main difference between these two catalogues is that Marsden & Williams (2008) did a backward integration and evaluated the orbital elements when the comets were far away from the planetary region. They thus obtained the original orbital elements. Instead, JPL evaluates the orbital elements when the comets are still in the planetary region. At present the Catalogue of Cometary Orbits is no longer being updated (G. Williams, pers. comm.) and therefore it does not contain many recent cometary discoveries from surveys such as WISE, Pan-STARRS and Catalina. These new comets are listed in the JPL catalogue, however. Therefore, in our work we decided to make use of the JPL catalogue. As of July 2013, based on the criteria mentioned in previous section, the JPL catalogue contains 427 LPCs and 60 HTCs. All HTCs based on our definition are listed in Table1. One thing to note in the table is that all these comets were discovered at various times, some even dating back to the 18th century.

One can imagine that, because of instrumental limitations, comets discovered even as late as the early 20th century needed to be either very bright and/or very close to Earth, and hence were most likely observationally biased (Everhart 1967). To support our claim, the perihelion distance versus inclination for LPCs (top) and HTCs (bottom) is shown in Fig. 1. Blue circles are comets discovered before 1950, and red circles are comets discovered after 1950. For both LPCs and HTCs the blue circles were concentrated near low perihelia. Thus, it appears that comets discovered before 1950 had a stronger observational bias toward low perihelia because of the instrumental limitations and limited sky coverage at the time. Comets discovered after 1950 are mostly detected by surveys that have fewer biases such as sweeping larger portions of the sky, and higher limiting magnitude. In addition, the HTC sample contains more blue dots at low inclination than at high inclination, hinting at an inclination bias as well.

Fig. 1 Perihelion versus inclination of comets from the JPL Small-Body Search Engine. Top panel (plots only comets with q < 5 AU, 92% of all JPL LPCs) is the plot of LPCs, and the bottom panel (plots only comets with q < 5 AU, 100% of all JPL HTCs) is for HTCs. Red circles are comets discovered after 1950 and blue circles are comets discovered before 1950. As can be seen, comets discovered before 1950 were all concentrated at low perihelion distances.

It is well known that the HTCs show a predominantly prograde inclination distribution (Fernández & Gallardo 1994; Levison et al. 2001). It is possible that observational biases play a substantial role in the determination of the cumulative inclination distribution of HTCs and this could be the reason why the median inclination value has thus far been substantially lower than the one expected from an isotropic distribution (Levison et al. 2001, 2006). Levison et al. (2001) and Levison et al. (2006) tried to achieve observational completeness by concentrating only on HTCs with q < 1.3 AU or 2.5 AU, following Duncan et al. (1988) and Quinn et al. (1990). However, it is not clear whether this prograde excess of HTCs is a result of observational biases or their definition. Indeed, the relation T_J < 2 implies $\sqrt{\mathit{q}}\cos \mathit{i}$ is approximately constant if a ≫ q, which holds for most HTCs. Thus using an HTC catalogue with a maximum q of 1.3 AU may result in a lower median inclination than when imposing a maximum q of 2 AU. This phenomenon is clearly visible in Table 2 for the observational data, which lists the number of discovered HTCs in 50-year intervals, together with their median perihelion distance and median inclination. Dates closer to the current epoch detect HTCs with higher perihelion distance because of their increased limiting magnitude and sky coverage, but the median inclination also increases and approaches an isotropic distribution. However, data from our numerical simulations, discussed below, contain no such trend and instead have a median inclination close to 90° for all values of q, supporting the idea that the prograde excess is caused by early observational biases.

At the same time we cannot rule out that the definition of HTCs is solely responsible for the prograde excess. Indeed, none of the previous studies made a cut to the dataset based on the epoch of observations, potentially still leading to a biased sample. Here we make use of the JPL catalogue from Fig. 1 to verify this point. Apparently, if one only keeps all HTCs with q < 1.3 AU by arguing that this yields an observationally complete sample for the entire dataset, then many new discoveries due to advanced technology after 1950 will be discarded. Not knowing the observational conditions for comets discovered before 1950, attempting to de-bias the old discoveries would be difficult. Fortunately, we have many new discoveries in the last 20 years (see Table 2) and therefore we decided that the easiest way to de-bias the old discoveries is to remove comets discovered earlier than 1950 that did not pass perihelion between 1950 and 2013.

Fig. 2 Histogram of perihelion distance q of comets with TJ < 2 from the JPL catalogue. It shows a sharp drop of the number of comets with q > 2 AU.

The last point that warrants discussion is the maximum imposed perihelion distance of the observed sample. Duncan et al. (1988) and Levison et al. (2001) imposed different threshold values of the maximum perihelion distances below which the comet sample was thought to be observationally complete. In other words, it was assumed that below a maximum perihelion distance q_m there was a 100% detection efficiency, up to a certain total limiting magnitude, for all comets. For this reason we chose to keep the maximum perihelion distance a free parameter below 2 AU. Our choice for the upper limit is motivated by the distribution of perihelion distance of comets with T_j < 2 in the JPL catalogue, which shows a sharp decrease beyond 2 AU (see Fig. 2).

We ended up with a total of 158 LPCs and 40 HTCs. We now have a fairly unbiased sample of LPCs and HTCs from the JPL comet database for comparison between theory and observations. Thus, we move on to describe our numerical methods in the next section.

3. Numerical simulation and initial conditions

We determined whether or not the Oort Cloud is a feasible mechanism to generate the HTCs by running numerical simulations of the evolution of Oort Cloud comets under the gravitational perturbations of the Sun, giant planets, Galactic tide, passing field stars and by comparing the ensuing orbital element distributions against those of the observed sample.

We generated a fictitious Oort Cloud by using the method described in Kaib et al. (2009) from the Oort Cloud formation simulations presented in Brasser & Morbidelli (2013). We divided each Oort Cloud into six equal sections of semi-major axis and computed the cumulative distributions of the argument of perihelion (ω), longitude of the ascending node (Ω), specific orbital angular momentum (L), and the z-component of the specific orbital angular momentum (L_z) for each section. The cumulative distribution of the semi-major axis was computed for the whole Oort Cloud. We generated a comet by picking a uniform random number between 0 and 1 and selecting the corresponding semi-major axis from the semi-major axis distribution. The value of the semi-major axis identified the section of the cloud the comet would reside in. We then picked another four uniform random numbers to determine the values of ω, Ω, L, and L_z. The eccentricity was then computed from L and the inclination from L_z. The mean anomaly of the comet was picked at random uniformly between 0 and 360°. We repeated this process 32 768 times for each simulation and generated 40 separate simulations, resulting in a total of 1.3 million test particles.

The motion of the comets was subsequently integrated using the SWIFT RMVS3 integrator (Levison & Duncan 1994) without the giant planets present, but with the perturbations from the Galactic tide and the passing stars. The perturbations from the Galactic tide and stars were modelled after Levison et al. (2001) for the tides and Rickman et al. (2008) for the stars. The comets were removed from the simulation once they got closer than 38 AU to the Sun or farther than 1 pc from the Sun. These simulations were run for 4 Gyr with a time step of 50 yr. Different stellar encounter files were used for each simulation.

The second part of our simulations consisted of re-integrating only those comets that were removed from the first set of simulations because they came closer than 38 AU to the Sun. We also integrated a single clone of these objects to increase the total particle count per simulation. On average only 18% of all Oort Cloud objects penetrated the planetary region; the rest remained in the Oort Cloud or were removed by other means. The clone had the same original elements apart from a random mean anomaly. This typically resulted in 12 000 test particles per simulation and a total of 40 simulations. The difference with the first set of simulations is that now we included the giant planets. This procedure of running the simulations in two parts increases the number of comets we can integrate through the realm of the giant planets.

We ran these second set of simulations with SCATR (Symplectically-Corrected Adaptive Timestepping Routine; Kaib et al. 2011). It is based on SWIFT’s RMVS3, but it has a speed advantage over SWIFT’s RMVS3 or MERCURY (Chambers 1999) for objects far away from both the Sun and the planets where the time step is increased. We set the boundary between the regions with short and long time steps at 300 AU from the Sun (Kaib et al. 2011). Closer than this distance the computations are performed in the heliocentric frame, like SWIFT’s RMVS3, with a time step of 0.4 yr. Farther than 300 AU, the calculations are performed in the barycentric frame and we increased the time step to 50 yr. The error in the energy and angular momentum that is incurred every time an object crosses the boundary at 300 AU is significantly reduced through the use of symplectic correctors (Wisdom et al. 1996). For the parameters we consider, the cumulative error in energy and angular momentum incurred over the age of the solar system is of the same order or smaller than that of SWIFT’s RMVS3. The same Galactic and stellar parameters as in the first simulation were used. Comets were removed once they were farther than 1 pc from the Sun, or if they collided with the Sun or a planet.

To determine how the comets fade we modified SCATR to keep track of the number of perihelion passages that each comet made. Levison et al. (2001) suggest that comets fade the fastest when their perihelion distance q < 2.5 AU, so we only counted the perihelion passages when the perihelion was closer than 2.5 AU. The number of passages together with time, osculating barycentric semi-major axis, eccentricity, inclination, perihelion distance, and mean anomaly were stored in a separate file with an output interval of 100 yr; but only when the comet was closer than 300 AU to the Sun. The format of the output made for easy comparison with the JPL catalogue.

The sublimation of the cometary material due to the solar radiation will make the comets become fainter and fainter after each perihelion passage. Depending on the nature of the cometary material, comets will eventually become too faint after a high number of perihelion passages and turn into dormant or extinct comets that easily escape detection. This is the so-called fading problem first suggested by Oort (1950) to explain the deficiency of observed comets. A more extended discussion of different fading functions can be found in Wiegert & Tremaine (1999). In this work, we applied a post-processing simple power-law fading of the form ${\mathrm{\Phi }}_{\mathit{m}}\mathrm{=}{\mathit{m}}^{\mathrm{-}\mathit{k}}\mathit{,}$(1)where Φ_m is the remaining visibility function introduced by Wiegert & Tremaine (1999). A comet that has its mth perihelion passage with q < 2.5 AU will have its remaining visibility be Φ_m, where k is a positive constant, and Φ₁ = 1 for the first out-going perihelion passage of a comet. Without fading, every comet in our simulation has the same weighting (Φ_m = 1, m = 1,2,3...,n_p) in constructing the cumulative distribution of semi-major axis or inclination of active comets. When we applied the fading effect to comets, the remaining visibility Φ_m is considered to be a weighting factor. The larger the number of perihelion passages, the less each comet contributes to the cumulative inclination or semi-major axis distribution of the active comets.

For example, in constructing the inclination distribution with fading, we first sort the comets by inclination in increasing order and make the cumulative fraction based on their remaining visibilities. In this way we can construct the cumulative distribution of orbital elements with fading by changing the power-law index k in the remaining visibility. Whipple (1962) and Wiegert & Tremaine (1999) found k ~ 0.6 has the best match between simulation and observation and thus we expect to obtain a similar value here. Our database has a baseline of 63 years, which is shorter than the period of some HTCs. Therefore, before comparing our simulation to the observation, the cumulative distributions of inclination and semi-major axis from simulation have to be weighted by the ratio of the period of comets and the observed window. We roughly corrected for the bias with a weighting factor to the remaining visibility of ${\mathrm{\Phi }}_{\mathit{m}}^{\mathrm{\prime }}\mathrm{=}{\mathit{m}}^{\mathrm{-}\mathit{k}}(\frac{\mathrm{63}\hspace{0.17em}\mathrm{yr}}{{\mathit{a}}^{\mathrm{3}\mathit{/}\mathrm{2}}})\mathrm{;}\mathit{a}\mathit{>}\mathrm{15.8}\hspace{0.17em}\mathrm{AU}\mathit{.}$(2)For every discovered comet with high semi-major axis (a > 15.8 AU), there are actually a^3/2/63 more comets with a similar semi-major axis comets. We ran all of our simulations on the CP computer cluster at ASIAA. The total time taken per simulation was about one month.

4. Simulation results

In this section we present the results of our numerical simulations and the comparison with observational data. Our simulations yielded 7140 individual LPCs and 1159 individual HTCs that pass the same criteria, with each category containing a few million entries.

4.1. Oort Cloud model and fading

We first need to establish whether our OC model is compatible with the observed OOCs and LPCs. We followed the analysis in Wiegert & Tremaine (1999) to determine whether or not our simulations could simultaneously match different parts of the observed LPC population from the JPL catalogue. Wiegert & Tremaine (1999) introduce three quantities that measure the relative contributions to the flux of comet passing perihelion: Ψ₁ is the fraction of all comets with a > 10 kAU; Ψ₂ is the fraction of all comets that have 34.2 < a < 69.8 AU, i.e. the tail end of the LPC distribution and Ψ₃ is the fraction of all comets that are prograde. In addition to the Ψ values, we also follow Wiegert & Tremaine (1999) and define the simulated to observed ratio ${\mathit{X}}_{\mathit{i}}\mathrm{=}{\mathrm{\Psi }}_{\mathit{i}}^{\mathrm{\prime }}\mathit{/}{\mathrm{\Psi }}_{\mathit{i}}$, where ${\mathrm{\Psi }}_{\mathit{i}}^{\mathrm{\prime }}$ is the simulated value after applying fading. A good model for the Oort Cloud should yield X_i = 1 simultaneously at certain fading index k.

In order to compare the flux of comets in our simulations with those of Wiegert & Tremaine (1999) we need to impose a maximum perihelion distance for the LPCs, as opposed to the subsequent HTC analysis where we consider them to always be detected. For our purposes we decided to use q_m = 2.0 AU which equals the median value in the current LPC sample. From the JPL catalogue for comets with q < 2.0 AU that were discovered or observed after 1950 we have

• Ψ₁ = 0.044 ± 0.020,

• Ψ₂ = 0.152 ± 0.043,

• Ψ₃ = 0.519 ± 0.100,

where the errors are assumed to be Poisson statistics. We did not include hyperbolic comets in the determination of Ψ₁. The Ψ values should be compared to those in Wiegert & Tremaine (1999), who reported

• Ψ₁ = 0.380 ± 0.043,

• Ψ₂ = 0.066 ± 0.015,

• Ψ₃ = 0.505 ± 0.051.

The largest discrepancy occurs in Ψ₁, with our value approximately an order of magnitude lower than theirs. We attribute this to the fact that Wiegert & Tremaine (1999) used the Catalogue of Cometary Orbits of 1993 (Marsden & Williams 1993) rather than JPL’s, and they calculated Ψ₁ from the original semi-major axis of the comets. In addition, Wiegert & Tremaine (1999) did not apply a cut in q when calculating Ψ₁ and used all the LPCs in the catalogue, which could also increase the value of Ψ₁.

Here we decided to use the JPL catalogue, which publishes the semi-major axis of comets when they are still in the realm of the giant planets. Doing this skews the semi-major axis distribution with the result that the contribution from the Oort spike is substantially reduced. This reduction in the Oort spike is the result of planetary perturbations which cause a change in binding energy that is typically greater than the original binding energy (Fernández 1981). Thus, the use of the JPL catalogue and evaluating the orbits in the realm of the giant planets may cause errors in classifying a comet as LPC or OCC. This means we need to determine whether or not our low value of Ψ₁ is somehow consistent with that of Wiegert & Tremaine (1999) by correcting for the different methods to compute the semi-major axes.

Even though the initial conditions of Wiegert & Tremaine (1999) were very different from ours, we can use the results of their numerical simulations to estimate our value of Ψ₁ when planetary perturbations are removed. Without cometary fading Wiegert & Tremaine (1999) find X₁ = 0.075 ± 0.011, and because ${\mathrm{\Psi }}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{=}{\mathrm{\Psi }}_{\mathrm{1}}{\mathit{X}}_{\mathrm{1}}$ our corrected value is Ψ₁ ~ 0.508 ± 0.185. This is consistent with Wiegert & Tremaine (1999) within the error bars.

Wiegert & Tremaine (1999) report that numerical simulations without fading yield values of X_i that are very different from unity, with X₁ being too low and X₂ being too high. For this reason they introduce several fading laws and compute the X_i values for various fading parameters. They argue that introducing a fading law should result in all X_i values being near unity, with some laws better than others. They obtain the best result from the power law fading Eq. (1) above with k ~ 0.6. We applied the same fading law to our simulations to test our initial conditions and numerical methods. In Fig. 3 we depict the X values as a function of fading index k. The red line plots X₁, the green line is for X₂ and the blue line depicts X₃. One can tell from the plot that after applying the fading, our simulations are a good match to the observed distribution when the fading index is in the range 0.6 < k < 0.8, with a slight preference towards 0.6 rather than 0.7. This is close to the same value found by Whipple (1962) and Wiegert & Tremaine (1999) and suggests that our Oort Cloud model can be made consistent with observations if the comets fade via this law.

Fig. 3 Values of Xi as a function of fading index k. The three curves cross each other near the region 0.6 < k < 0.8, indicating that after applying fading to the simulations, the best match should fall in this region.

The result depicted in Fig. 3 begs the following question: can the same fading law be used to reproduce the observed orbital distribution of HTCs? We attempt to answer this in the next subsection.

4.2. The Halley-type comets

In order to find out how well our simulations match with the observed HTCs, we computed the cumulative inclination and semi-major axis distributions of the observed and simulated comets for a range of maximum perihelia, q_m, and fading parameter, k. Once these distributions were generated we performed a Kolmogorov-Smirnov (K-S) test (Press et al. 1992), which searches for the maximum absolute deviation D_max between the observed and simulated cumulative distributions. The K-S test assumes that the entries in the distributions are statistically independent. The probability of a match, P_D, as a function of D_max, can be calculated to claim whether these two populations were from the same parent distribution. However, in our simulations, a single comet would be included many times in the final distribution as long as the comet met our criteria for being an HTC during each of its perihelion passages. Including its dynamical evolution in this manner would cause many of the entries in the final distribution to become statistically dependent and the K-S test would be inapplicable. In the meantime, there are only between 28 and 40 HTCs in the JPL catalogue that passed our criteria, depending on q_m, which means our results may suffer because of the small data set. These limitations prevent us from using the analytical formula described in Press et al. (1992) to calculate the K-S probability.

We solved this issue by applying a Monte Carlo method to perform the K-S test as described in Levison et al. (2006). Once we have the inclination and semi-major axis distributions from the simulation, we then randomly selected 10 000 fictitious samples from the simulation. Each sample has the same number of data points as the HTCs from the JPL catalogue. The Monte Carlo K-S probability is then the fraction of cases that have their D values between the fictitious samples and the real HTC samples larger than the D_max found from the real HTC samples and cumulative distributions from simulation.

Before making the fictitious datasets, we needed to find the empirical probability density functions of inclination and semi-major axis from which we then sampled the fictitious HTCs. Here we generated these from a normalized histogram. One crucial point in making the histogram is that we weighed each entry by its remaining visibility (Eq. (1)) and its orbital period (Eq. (2)). The fictitious comets were then sampled from the distributions with the von Neumann rejection technique (von Neumann 1951). This sampling method relies on generating two uniform random numbers on a grid. An entry is accepted when both numbers fall under the probability density curve.

Fig. 4 Contour plots of the K-S probability as a function of the fading parameter and maximum perihelion distance. The top panel shows the fit for the inclination distribution. The bottom panel pertains to the semi-major axis. The parameters that best fit both distributions are k = 1.05 and q = 1.77 AU.

Fig. 5 Inclination and semi-major axis distributions from observation and simulation for HTCs that gave the best match from Fig. 4.

The results of applying the fading to the simulated data in the HTC region are shown in Fig. 4. Here we plot contours of the K-S probability as a function of k and q_m for the inclination distribution (top) and semi-major axis distribution (bottom). The highest K-S probability for the inclination occurs at low q_m ~ 1.5 AU and relatively high k ~ 1.4, while for the semi-major axis distribution the maximum K-S probability occurs when q_m ~ 1.8 AU and k ~ 1. Combining these results we find that the combination of q_m and k that best fits both the semi-major axis and inclination distribution is q = 1.77 AU and k = 1.05, with K-S probabilities of 84.3% for the semi-major axis and 31% for the inclination. We show the resulting fits to the data in Fig. 5 below. The top panel shows the inclination distribution, the bottom panel depicts the semi-major axis distribution. The median simulated inclination is 80.4° (observed 74.2°) and the median simulated semi-major axis is 14.3 AU (observed 16.3 AU). Levison et al. (2001) apply a different fading law than the simple one we proposed here, but they conclude that their model may work when k ~ 0.8. This is in the same ballpark as ours, but they state that their value is very uncertain.

Summarizing, it appears that an Oort Cloud source combined with a simple fading law can match the observed HTC inclination and semi-major axis distribution rather well when imposing a maximum perihelion distance close to 1.8 AU and taking the fading parameter to be ~ 1. The imposed upper value of the perihelion distance is below the maximum of 2 AU beyond which the catalogue shows a steep decrease in the number of comets, and is close to the median value obtained from observations in the last 13 years. This suggests that the current catalogue of HTCs with q < 1.8 AU is nearing completion.

Unfortunately, the optimum value of k is inconsistent with that of the LPCs found in the previous section. It appears that the HTCs fade more quickly than the LPCs. It is now well recognized that cometary nuclei develop non-volatile, lag-deposit crusts that reduce the fraction of the nucleus surface available for sublimation (Brin & Mendis 1979; Fanale & Salvail 1984). For most Jupiter-family comets (JFCs), the active fraction, which is the fraction of the nucleus surface area that must be active to explain the comet’s water production rate, is typically only a few per cent, or even a fraction of a per cent (Fernández et al. 1999). Most likely for HTCs the active surface fraction is similar to that of JFCs. Activity could be re-ignited if the comet undergoes a substantial decrease in perihelion (Rickman et al. 2008). The age of the comet and a non-pristine surface could account for the fading index being higher than for new comets that first enter the inner solar system.

4.3. Expected HTC population

We can use the results of the numerical simulations above to constrain the expected number of HTCs larger than a given size using the method described in Levison & Duncan (1997) and Brasser & Morbidelli (2013). These works focused on the active lifetimes of Jupiter-family comets (JFCs). The method consists of finding the best match of simulated comets to the observed semi-major axis and inclination distributions as a function of active lifetime τ_vHTC. In other words one computes the K-S probability for matching the simulated inclination or semi-major axis distribution to the observed one for each value of the active lifetime. The maximum K-S probability corresponds to the best estimate of τ_vHTC. However, both of the above studies assumed the comets retained their original activity for a time τ_vHTC and would fade completely immediately afterwards. Our matching procedure here must rely on a different approach because we fade the comets with time. Thus, the choice of τ_vHTC should depend on the fading index k and the number of perihelion passages after which we consider the comet to have faded completely. There is also a q dependence, but in what follows we only consider HTCs with q < 1.8 AU.

The formula relating the number of active HTCs (N_vHTC) to the number of Oort Cloud comets (N_OC) is ${\mathit{N}}_{\mathrm{vHTC}}\mathrm{=}{\mathit{N}}_{\mathrm{OC}}{\mathit{\tau }}_{\mathrm{vHTC}}\mathrm{\left|}{\mathit{r}}_{\mathrm{OC}}\mathrm{\right|}{\mathit{f}}_{\mathrm{vHTC}}\mathit{,}$(3)where r_OC is the fractional decay rate of the Oort Cloud at the current time, f_vHTC is the fraction of the comets escaping from the Oort Cloud that penetrate into the visible HTC region with q < 1.8 AU and N_vHTC and N_OC are the number of visible HTCs and Oort Cloud comets. We evaluate these three quantities below.

Defining by f_OC(t) the fraction of Oort Cloud objects surviving in the Oort Cloud at time t, the fractional decay rate of the Oort Cloud population is defined as r_OC(t) = (df_OC(t)/dt)/f_OC(t). We measured r_OC from the simulations in Brasser & Morbidelli (2013) over the last 3 Gyr and obtained ⟨ r_OC ⟩ = −1.48 × 10^-10 per year.

The fraction of Oort Cloud objects that became HTCs with q < 1.8 AU, f_vHTC, is tied to the active lifetime, τ_vHTC. When the HTCs enter the inner solar system for the first time with q < 1.8 AU many of them will have m ≫ 1 and so some of these will already have faded and will not be detected. If we were to assume that the average time any HTC spends with q < 1.8 AU while active at each output entry is equal to the entire 100 yr between outputs, we would obtain τ_vHTC = 310 kyr, which is longer than the median dynamical lifetime of 112 kyr (cf. Levison et al. 2006; and Nurmi et al. 2002; report 69 kyr and 100 kyr). Thus a different approach is needed.

Brasser & Morbidelli (2013) report an active JFC lifetime of τ_vJFC = 12 kyr. This corresponds to roughly 2000 revolutions, but only ~400 of these are spent with q < 2.5 AU. Imposing the latter as an upper limit on the number of perihelion passages for the HTCs with q < 1.8 AU, n_p, results in a remaining visibility of the order of 0.25% when k ~ 1. By this time the comet has faded substantially: the change in absolute magnitude is ΔH = 2.5klog n_p = 6.5. We take this as a good upper limit to consider the comet to have faded for the following reasons.

Sosa & Fernández (2011) showed that an LPC with total absolute magnitude H_T = 6.5 has a nuclear magnitude of H_N = 17.3, corresponding to a diameter of 2.3 km assuming an albedo of 4%. A JFC with the same size has a typical total absolute magnitude H_T = 9.3 (Fernández et al. 1999; Brasser & Morbidelli 2013). Thus, on average, one may consider the JFCs to have faded by 2.8 mag compared to LPCs of the same size. Assuming that the HTCs follow the same trend, our earlier limit of 6.5 mag for complete fading suggested above seems reasonable.

Assuming once again that the comet spends the full 100 yr in the region of interest between outputs, and also assuming that the JFCs and HTCs spend a similar amount of time with q > 2.5 AU we obtain τ_vHTC = 23.2 kyr.

We then compute f_vHTC by identifying all unique HTCs in our simulations that have n_q ≤ 2000 when they first reach q < 1.8 AU and divide by the total number of comets. This results in f_vHTC = 4.2 × 10^-4. Our value is approximately a factor of 1.5 lower than that found by Levison et al. (2001).

The last ingredient of the formula is the number of comets in the Oort Cloud. Brasser & Morbidelli (2013) report there are 7.6 × 10¹⁰ comets in the Oort Cloud with diameter D > 2.3 km. Combining these together results in N_vHTC ~ 100 with D > 2.3 km and q < 1.8 AU. This estimate is quite uncertain, with the largest uncertainty being in τ_vHTC, which can easily vary by a factor of a few. Incidentally the number of active HTCs is comparable to the number of active JFCs of the same size reported in Levison & Duncan (1997) and Brasser & Morbidelli (2013). The ratio of dormant to active HTCs is approximately 4, given by dividing the dynamical lifetime by the active lifetime.

5. Discussion

In the previous section we demonstrated the ability to reproduce the currently observed inclination and semi-major axis distribution of the HTCs from an Oort Cloud source using a simple cometary fading law. The result is simple, works very well, and thus begs further discussion.

We want to emphasize here that we have been careful to minimise possible sources of error and contamination. First, we have justified our rationale for using the JPL catalogue and for comparing our simulations to the observational data. The initial conditions from our simulations come from earlier Oort Cloud formation simulations from Brasser & Morbidelli (2013) using a robust method for obtaining a highly populated Oort Cloud without needing to resort to many Oort Cloud formation simulations (Kaib et al. 2009). We used a computer code that has been tested extensively and used the full models for the Galactic tides, and passing stars (Levison et al. 2001; Rickman et al. 2008). Thus, we see no obvious substantial difference in the numerical methodology with earlier works that attempted to reproduce the orbit distribution of HTCs.

There are two points where we differ from earlier attempts by Levison et al. (2001) and Levison et al. (2006). First, we included more HTCs in our observational data. Both Levison et al. (2001) and Levison et al. (2006) used the observed HTCs with q < 1.3 AU, arguing that HTCs with larger perihelion distance were observationally incomplete. Their sample remained small, however, containing just ~20 comets. However, they also included comets discovered before 1950 which we discarded if they did not pass perihelion after 1950 for reasons discussed in Sect. 2. Second, the last decade has seen a vast increase in the number of observed comets, including new HTCs. Compared to Levison et al. (2001) and Levison et al. (2006), we have different inclination and semi-major axis distributions because of the cut in the epoch of observations, more new observations, and higher maximum perihelion distances. In other words, in our work we were comparing our simulated comets to a different inclination and semi-major axis distributions than Levison et al. (2001, 2006) because of the increased number of comets in the last decade. Specifically, the median inclination in the sample of Levison et al. (2006) is 55° and in our best-fit sample it is 75°, inching towards a more isotropic distribution that would be expected from the models. Similar arguments apply to the work of Nurmi et al. (2002), with the added emphasis that we had very different initial conditions, in particular the initial perihelion distribution.

A second point that warrants discussion is our method of calculating the fading of the comets. We based the fading on the methods of Wiegert & Tremaine (1999) because of their surprisingly good result for power-law fading when k ~ 0.6. Unfortunately, Wiegert & Tremaine (1999) did not describe how they applied their fading in great detail, something we did here. It is possible that our methods differ from those of Wiegert & Tremaine (1999), but the fact that we obtain similar results suggests otherwise.

6. Summary and conclusions

In this work we have used numerical simulations to determine whether the distribution of the HTCs can be reproduced using the Oort Cloud as their source distribution and when considering cometary fading. We generated a large number of Oort Cloud comets that were evolved numerically under the gravitational influence of the Sun, giant planets, Galactic tide and passing stars using well established methods. Guided by the work of Wiegert & Tremaine (1999) we post-processed the data with a power-law fading mechanism that was applied to all comets that passed perihelion closer than 2.5 AU. We compared the numerical output to the JPL catalogue, keeping only those comets discovered and observed after 1950 to limit observational biases.

We find that the Oort Cloud with cometary fading is the most likely source of the HTCs. Both of their semi-major axis and inclination distributions can be accurately reproduced when the fading index k ~ 1 and when one considers only those comets for which q ≲ 1.8 AU. The HTC fading index is somewhat higher than the value of 0.6 we find for LPCs. This increase in fading index could be the result of physical processes on the cometary surface. With these results we estimate there are about 100 active HTCs with diameters D > 2.3 km and q < 1.8 AU at any time.

Acknowledgments

We thank Julio A. Fernández and Nathan Kaib for helpful reviews that greatly improved the quality of this paper.

References

All Tables

All Figures

Fig. 1 Perihelion versus inclination of comets from the JPL Small-Body Search Engine. Top panel (plots only comets with q < 5 AU, 92% of all JPL LPCs) is the plot of LPCs, and the bottom panel (plots only comets with q < 5 AU, 100% of all JPL HTCs) is for HTCs. Red circles are comets discovered after 1950 and blue circles are comets discovered before 1950. As can be seen, comets discovered before 1950 were all concentrated at low perihelion distances.

In the text

	Fig. 2 Histogram of perihelion distance q of comets with TJ < 2 from the JPL catalogue. It shows a sharp drop of the number of comets with q > 2 AU.
In the text

	Fig. 3 Values of Xi as a function of fading index k. The three curves cross each other near the region 0.6 < k < 0.8, indicating that after applying fading to the simulations, the best match should fall in this region.
In the text

	Fig. 4 Contour plots of the K-S probability as a function of the fading parameter and maximum perihelion distance. The top panel shows the fit for the inclination distribution. The bottom panel pertains to the semi-major axis. The parameters that best fit both distributions are k = 1.05 and q = 1.77 AU.
In the text

	Fig. 5 Inclination and semi-major axis distributions from observation and simulation for HTCs that gave the best match from Fig. 4.
In the text