An Oort Cloud origin of the Halleytype comets
Institute for Astronomy and Astrophysics, Academia Sinica; 11F AS/NTU
building, 1 Roosevelt Rd., Sec. 4, 10617 Taipei, Taiwan
email:
brasser_astro@yahoo.com
Received: 20 August 2013
Accepted: 3 February 2014
The origin of the Halleytype comets (HTCs) is one of the last mysteries of the dynamical evolution of the solar system. Prior investigation into their origin has focused on two source regions: the Oort Cloud and the scattered disc. From the former it has been difficult to reproduce the nonisotropic, prograde skew in the inclination distribution of the observed HTCs without invoking a multicomponent Oort Cloud model and specific fading of the comets. The scattered disc origin fares better, but needs an order of magnitude more mass than is consistent with theory and observations. Here we revisit the Oort Cloud origin and include cometary fading. Our observational sample stems from the JPL catalogue. We only keep comets discovered and observed after 1950, but place no a priori restriction on the maximum perihelion distance of observational completeness. We then numerically evolve half a million comets from the Oort Cloud through the realm of the giant planets and keep track of their number of perihelion passages with perihelion distance q < 2.5 AU, below which the activity is supposed to increase considerably. We can simultaneously fit the HTC inclination and semimajor axis distribution very well with a powerlaw fading function of the form m^{−k}, where m is the number of perihelion passages with q < 2.5 AU and k is the fading index. We match both the inclination and semimajor axis distributions when k ~ 1 and the maximum imposed perihelion distance of the observed sample is q_{m} ~ 1.8 AU. The value of k is higher than the one obtained for the longperiod comets (LPCs), for which typically k ~ 0.7. This increase in k is most likely the result of cometary surface processes. We argue the HTC sample is now most likely complete for q_{m} < 1.8 AU. We calculate that the steadystate number of active HTCs with diameter D > 2.3 km and q < 1.8 AU is on the order of 100.
Key words: comets: general / Oort Cloud
© 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 semimajor axis showed a sharp increase when 1/a < 5 × 10^{4} AU^{1}. The observed semimajor 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^{11} 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 nearparabolic 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 longperiod 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_{t}r, 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 JupiterSaturn 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 semimajor 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 zcomponent 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 JupiterSaturn 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 Halleytype 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 lowinclination orbits. The scattered disc is the most likely source of the socalled Jupiterfamily 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 semimajor 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 Sungrazing 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 nearperpendicular 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 nonisotropy of the HTCs could only be reproduced with an Oort Cloud that had a nonisotropic 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 semimajor 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 planetcrossing 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 semimajor 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 semimajor 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 SmallBody Search Engine^{1}.
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, PanSTARRS 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.
A total of 60 HTCs from the JPL SmallBody Search Engine meet our HTC definition.
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 SmallBody 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. 

Open with DEXTER 
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 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 50year 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 debias 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 debias the old discoveries is to remove comets discovered earlier than 1950 that did not pass perihelion between 1950 and 2013.
Number of HTCs discovered every 50 years starting from 1850 to 2013.
Fig. 2 Histogram of perihelion distance q of comets with T_{J} < 2 from the JPL catalogue. It shows a sharp drop of the number of comets with q > 2 AU. 

Open with DEXTER 
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 semimajor axis and computed the cumulative distributions of the argument of perihelion (ω), longitude of the ascending node (Ω), specific orbital angular momentum (L), and the zcomponent of the specific orbital angular momentum (L_{z}) for each section. The cumulative distribution of the semimajor 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 semimajor axis from the semimajor axis distribution. The value of the semimajor 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 reintegrating 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 (SymplecticallyCorrected 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 semimajor 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 socalled 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 postprocessing simple powerlaw fading of the form (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} = 1 for the first outgoing 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 semimajor 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 semimajor 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 powerlaw 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 semimajor 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 (2)For every discovered comet with high semimajor axis (a > 15.8 AU), there are actually a^{3/2}/63 more comets with a similar semimajor 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: Ψ_{1} is the fraction of all comets with a > 10 kAU; Ψ_{2} is the fraction of all comets that have 34.2 < a < 69.8 AU, i.e. the tail end of the LPC distribution and Ψ_{3} 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 , where 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

Ψ_{1} = 0.044 ± 0.020,

Ψ_{2} = 0.152 ± 0.043,

Ψ_{3} = 0.519 ± 0.100,
where the errors are assumed to be Poisson statistics. We did not include hyperbolic comets in the determination of Ψ_{1}. The Ψ values should be compared to those in Wiegert & Tremaine (1999), who reported

Ψ_{1} = 0.380 ± 0.043,

Ψ_{2} = 0.066 ± 0.015,

Ψ_{3} = 0.505 ± 0.051.
The largest discrepancy occurs in Ψ_{1}, 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 Ψ_{1} from the original semimajor axis of the comets. In addition, Wiegert & Tremaine (1999) did not apply a cut in q when calculating Ψ_{1} and used all the LPCs in the catalogue, which could also increase the value of Ψ_{1}.
Here we decided to use the JPL catalogue, which publishes the semimajor axis of comets when they are still in the realm of the giant planets. Doing this skews the semimajor 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 Ψ_{1} is somehow consistent with that of Wiegert & Tremaine (1999) by correcting for the different methods to compute the semimajor 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 Ψ_{1} when planetary perturbations are removed. Without cometary fading Wiegert & Tremaine (1999) find X_{1} = 0.075 ± 0.011, and because our corrected value is Ψ_{1} ~ 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_{1} being too low and X_{2} 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_{1}, the green line is for X_{2} and the blue line depicts X_{3}. 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 X_{i} 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. 

Open with DEXTER 
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 Halleytype comets
In order to find out how well our simulations match with the observed HTCs, we computed the cumulative inclination and semimajor 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 KolmogorovSmirnov (KS) test (Press et al. 1992), which searches for the maximum absolute deviation D_{max} between the observed and simulated cumulative distributions. The KS 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 KS 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 KS probability.
We solved this issue by applying a Monte Carlo method to perform the KS test as described in Levison et al. (2006). Once we have the inclination and semimajor 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 KS 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 semimajor 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 KS 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 semimajor axis. The parameters that best fit both distributions are k = 1.05 and q = 1.77 AU. 

Open with DEXTER 
Fig. 5 Inclination and semimajor axis distributions from observation and simulation for HTCs that gave the best match from Fig. 4. 

Open with DEXTER 
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 KS probability as a function of k and q_{m} for the inclination distribution (top) and semimajor axis distribution (bottom). The highest KS probability for the inclination occurs at low q_{m} ~ 1.5 AU and relatively high k ~ 1.4, while for the semimajor axis distribution the maximum KS 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 semimajor axis and inclination distribution is q = 1.77 AU and k = 1.05, with KS probabilities of 84.3% for the semimajor 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 semimajor axis distribution. The median simulated inclination is 80.4° (observed 74.2°) and the median simulated semimajor 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 semimajor 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 nonvolatile, lagdeposit crusts that reduce the fraction of the nucleus surface available for sublimation (Brin & Mendis 1979; Fanale & Salvail 1984). For most Jupiterfamily 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 reignited if the comet undergoes a substantial decrease in perihelion (Rickman et al. 2008). The age of the comet and a nonpristine 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 Jupiterfamily comets (JFCs). The method consists of finding the best match of simulated comets to the observed semimajor axis and inclination distributions as a function of active lifetime τ_{vHTC}. In other words one computes the KS probability for matching the simulated inclination or semimajor axis distribution to the observed one for each value of the active lifetime. The maximum KS 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 (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^{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 semimajor 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 semimajor 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 semimajor 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 bestfit 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 powerlaw 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 postprocessed the data with a powerlaw 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 semimajor 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
 Brasser, R., & Morbidelli, A. 2013, Icarus, 225, 40 [NASA ADS] [CrossRef] [Google Scholar]
 Brin, G. D., & Mendis, D. A. 1979, ApJ, 229, 402 [NASA ADS] [CrossRef] [Google Scholar]
 Chambers, J. E. 1999, MNRAS, 304, 793 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Duncan, M. J., & Levison, H. F. 1997, Science, 276, 1670 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Duncan, M., Quinn, T., & Tremaine, S. 1988, ApJ, 328, L69 [NASA ADS] [CrossRef] [Google Scholar]
 Dybczyński, P. A., & Królikowska, M. 2011, MNRAS, 416, 51 [NASA ADS] [Google Scholar]
 Everhart, E. 1967, AJ, 72, 716 [NASA ADS] [CrossRef] [Google Scholar]
 Fanale, F. P., & Salvail, J. R. 1984, Icarus, 60, 476 [NASA ADS] [CrossRef] [Google Scholar]
 Fernández, J. A. 1981, A&A, 96, 26 [NASA ADS] [Google Scholar]
 Fernández, J. A., & Gallardo, T. 1994, A&A, 281, 911 [NASA ADS] [Google Scholar]
 Fernández, J. A., Tancredi, G., Rickman, H., & Licandro, J. 1999, A&A, 352, 327 [NASA ADS] [Google Scholar]
 Fouchard, M., Rickman, H., Froeschlé, C., & Valsecchi, G. B. 2013, Icarus, 222, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Gomes, R. S., Fernández, J. A., Gallardo, T., & Brunini, A. 2008, in The Scattered Disk: Origins, Dynamics, and End States, eds. M. A. Barucci, H. Boehnhardt, D. P. Cruikshank, A. Morbidelli, & R. Dotson (Tucson: University of Arizona press), 259 [Google Scholar]
 Heisler, J., & Tremaine, S. 1986, Icarus, 65, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Hills, J. G. 1981, AJ, 86, 1730 [NASA ADS] [CrossRef] [Google Scholar]
 Kaib, N. A., Becker, A. C., Jones, R. L., et al. 2009, ApJ, 695, 268 [NASA ADS] [CrossRef] [Google Scholar]
 Kaib, N. A., Quinn, T., & Brasser, R. 2011, AJ, 141, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Królikowska, M., & Dybczyński, P. A. 2010, MNRAS, 404, 1886 [NASA ADS] [Google Scholar]
 Levison, H. F. 1996, in Completing the Inventory of the Solar System, eds. T. Rettig, & J. M. Hahn, ASP Conf. Ser., 107, 173 [Google Scholar]
 Levison, H. F., & Duncan, M. J. 1994, Icarus, 108, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Levison, H. F., & Duncan, M. J. 1997, Icarus, 127, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Levison, H. F., Dones, L., & Duncan, M. J. 2001, AJ, 121, 2253 [NASA ADS] [CrossRef] [Google Scholar]
 Levison, H. F., Duncan, M. J., Dones, L., & Gladman, B. J. 2006, Icarus, 184, 619 [NASA ADS] [CrossRef] [Google Scholar]
 Marsden, B. G., & Williams, G. V. 1993, Catalogue of Cometary Orbits 1993, Eighth edn. (Cambridge, MA: IAU) [Google Scholar]
 Marsden, B. G., & Williams, G. V. 2008, Catalogue of Cometary Orbits 2008 [Google Scholar]
 Nurmi, P., Valtonen, M. J., Zheng, J. Q., & Rickman, H. 2002, MNRAS, 333, 835 [NASA ADS] [CrossRef] [Google Scholar]
 Oort, J. H. 1950, Bull. Astron. Inst. Netherlands, 11, 91 [NASA ADS] [Google Scholar]
 Öpik, E. J. 1951, Proc. R. Irish Acad. Sect. A, 54, 165 [Google Scholar]
 Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical recipes in C, The art of scientific computing (Cambridge: University press) [Google Scholar]
 Quinn, T., Tremaine, S., & Duncan, M. 1990, ApJ, 355, 667 [NASA ADS] [CrossRef] [Google Scholar]
 Rickman, H., Fouchard, M., Froeschlé, C., & Valsecchi, G. B. 2008, Celest. Mech. Dyn. Astron., 102, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Sosa, A., & Fernández, J. A. 2011, MNRAS, 416, 767 [NASA ADS] [Google Scholar]
 Tisserand, F. 1889, Bulletin Astronomique, Serie I, 6, 289 [Google Scholar]
 Volk, K., & Malhotra, R. 2008, ApJ, 687, 714 [NASA ADS] [CrossRef] [Google Scholar]
 von Neumann, J. V. 1950, Nat. Bureau Standards 12, 36 [Google Scholar]
 Weissman, P. R. 1980, A&A, 85, 191 [NASA ADS] [Google Scholar]
 Whipple, F. L. 1962, AJ, 67, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Wiegert, P., & Tremaine, S. 1999, Icarus, 137, 84 [NASA ADS] [CrossRef] [Google Scholar]
 Wisdom, J., Holman, M., & Touma, J. 1996, Fields Institute Communications, 10, 217 [Google Scholar]
All Tables
A total of 60 HTCs from the JPL SmallBody Search Engine meet our HTC definition.
All Figures
Fig. 1 Perihelion versus inclination of comets from the JPL SmallBody 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. 

Open with DEXTER  
In the text 
Fig. 2 Histogram of perihelion distance q of comets with T_{J} < 2 from the JPL catalogue. It shows a sharp drop of the number of comets with q > 2 AU. 

Open with DEXTER  
In the text 
Fig. 3 Values of X_{i} 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. 

Open with DEXTER  
In the text 
Fig. 4 Contour plots of the KS 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 semimajor axis. The parameters that best fit both distributions are k = 1.05 and q = 1.77 AU. 

Open with DEXTER  
In the text 
Fig. 5 Inclination and semimajor axis distributions from observation and simulation for HTCs that gave the best match from Fig. 4. 

Open with DEXTER  
In the text 