© The Authors 2023

1 Introduction

The Oort Cloud remains the last region of the Solar System where no objects have ever been observed at the area they are believed to occupy (heliocentric distance r ~ 10^4–10⁵ au). Most of our knowledge about this area comes from analyses of the orbital element distributions of long-period comets observed from Earth, as done in the pioneering work of Oort (1950). Later, a more detailed analysis by Delsemme (1987) provided evidence for the role of Galactic tides in the injection of comets from the Oort Cloud toward observable regions. Since then, most studies devoted to the Oort Cloud and its relation to the observed flux of long-period comets have focused on estimating the number of comets in the cloud (Emel’yanenko et al. 2013; Kaib & Quinn 2009; Brasser & Morbidelli 2013). In the present study, we considered distributions of the orbital elements of observable comets that might carry valuable information about the shape of the Oort Cloud in the present and the past.

We had already conducted this line of study using synthetic models (Fouchard et al. 2020). In Fouchard et al. (2020) we show that the flux of long-period comets can yield features related to the initial shape of the Oort cloud. We drew conclusions from the numerical orbit propagation of two samples of objects. One is a swarm of comets with small ecliptic inclinations and perihelia at less than 45 au from the Sun, the so-called disk model. The other is a fully thermalized (Hills 1981) sample of objects, named the isotropic model.

The present work is a natural extension of Fouchard et al. (2020), who sought information on the shape of the Oort Cloud that the known long-period comets may contain. In Sect. 2, we recap the main results of Fouchard et al. (2020) and the relevant theoretical analysis in Higuchi (2020). In Sect. 3, the databases employed in this work are presented, together with the methodology used to produce the useful distributions. The results of our exploration of these databases are presented in Sect. 4. Section 5 presents our conclusions.

2 Choice of the features of interest

2.1 Synthetic samples

Fouchard et al. (2020) investigated whether the memory of the initial shape of the Oort Cloud could be detected in the distribution of new long-period comets observable after roughly 4.5 Gyr.

This preliminary investigation was based on the numerical orbit propagation of two large samples of synthetic comets over 4.5 Gyr, taking the effects of the Galactic tides, passing stars, and planetary perturbations into consideration. The Galactic parameters – a density of the local Galactic disk ρ₀ = 0.1 M_⊙ · pc⁻³ and an angular velocity of the Sun around the Galactic center Ω₀ = −26km s⁻¹ kpc⁻¹ – were taken from Levison et al. (2001), the stellar encounters sequence was computed randomly from the density of stars in the solar neighborhood (Rickman et al. 2004, 2008), and the planetary perturbations were computed taking the four giant planets on circular and coplanar orbits into account and using their present day semimajor axis.

The two samples are representative of two extreme cases of a proto-Oort Cloud with the following characteristics. The characteristics of the disk model include: uniform distributions for the orbital energy z = −1/a with a semimajor axis a between 500 au and 20000 au, the perihelion distance, q, between 3 and 50 au, the ecliptic-inclination, i_E, between 0° and 20°, and for the ecliptic argument of perihelion, ω_E, longitude of the ascending node, Ω_E, and mean anomaly, M, all between 0° and 360^°. The total number of initial comets for this set was 13 428 570.

Characteristics of the isotropic model are uniform distributions: for the orbital energy, z, with semimajor axis a between 500 au and 50 000 au; for cos i_G between −1 and 1, where i_G is the Galactic inclination; and for the Galactic argument of perihelion, ω_G, Galactic longitude of the ascending node, Ω_G, and mean anomaly, M, all between 0° and 360°. For the eccentricity, e, we used the distribution dn/de ∝ e, where n is the number of comets between e and de. This distribution corresponds to a fully thermalized Oort Cloud (Hills 1981). The initial perihelion distance, q, was greater than 15 au. For this set, the total number of objects was 20 307 700.

From each simulation, we constructed a sample of new observable comets. To realize this, five different stopping times between 3.5 and 4.5 Gyr separated by 250 Myr were considered. From each stopping time, the sample of remaining comets in the Oort Cloud was propagated forward for an additional 30 Myr; however, all passing stars that could produce a cometary shower were removed. After 30 Myr, each comet was propagated until the first perihelion passage. The comet was considered observable if this perihelion was at less than 5 au from the Sun. In addition, if its original¹ orbital energy, z = −1/a, was greater than 10⁻⁴au⁻¹, then the comet was considered a “new” observable long-period comet. The advantage of the five stopping times was that they worked as if we had considered an initial sample that was five times larger. A detailed description and justification of this strategy can be found in Fouchard et al. (2018, 2020).

It is important to note that the different observable comets in the models do not contribute similarly to the flux of new long-period comets per million years. Indeed, in the models a comet is propagated until its perihelion passage to determine whether it is observable. Thus, when a comet is observable, because of the orbital period of the comet, it must be converted into a flux of comet per million years. As such, the larger the orbital period of an observable comet, the smaller its contribution is to the flux per million year of new comets.

2.2 Transport and detection of the memory of a disk-like proto-Oort Cloud

As reported by Heisler & Tremaine (1986), the Galactic tides induce the von Zeipel-Lidov-Kozaï mechanism (von Zeipel 1910; Lidov 1961; Kozai 1962; Ito & Ohtsuka 2019) in the dynamics of Oort Cloud comets with large oscillation of the perihelion distance, related to periodic motion of the eccentricity, Galactic inclination, argument of perihelion and longitude of the ascending node. For comets whose semimajor axis remains below 15 000 au during the age of the Solar System, these long-term periodic motions could be largely conserved despite stellar perturbations. For such comets, detecting the memory of its initial position is possible, even when it becomes observable after 4.5 Gyr.

In Fouchard et al. (2020), three distinct categories of long-term dynamical evolution typical of the disk model were able to carry such a memory: (i) B1 group comets, whose orbits remain largely unchanged during the entire propagation except at the very end, when a rapid increase in the orbital energy induces a rapid decrease in the perihelion distance under the action of the Galactic tides; (ii) B2 group comets, with a single large perihelion distance oscillation caused by Galactic tides, with the other elements remaining nearly constant; and (iii) B3 group comets, for which the perihelion distance and argument of perihelion performed a single period, whereas the longitude of the ascending node performed half a period. Such comets returned to the planetary region with an orbital plane close to the empty ecliptic (EE) plane (Higuchi 2020), which was obtained by rotating the ecliptic plane around the Galactic pole by 180°. An example of each behavior is shown in Fig. 1.

The memory of their behavior and thus of their initial condition, was revealed through a detailed investigation of their orbital distributions when the comets became observable. In particular, it was shown that one had to consider the original orbital elements at the previous perihelion passage, that is, the one preceding the passage where the comet was observed. Thus, memory was detectable in the distribution of the ecliptic inclination and Galactic longitude of the ascending node. To highlight the features related to this memory, the sample of new long-period comets was divided in subsamples according to the previous perihelion distance, q_p. Indeed, the B3 group was mainly found for q_p < 7 au, whereas the B1 and B2 group mainly found for 13 < q_p < 31 au.

Details are provided here on the dynamic aspect of each group. From Higuchi (2020) it appears that when the perihelion distance q is making a full period together with the Galactic inclination i_G and the Galactic argument of perihelion, ω_G, the Galactic ascending node Ω_G performs half a period; second when the perihelion distance is going through its minimum, Ω_G and i_G change very quickly. In addition, for an orbit starting with small inclination with respect to the ecliptic, Ω_G precesses under the action of the Galactic tides

Therefore, groups B2 and B3 appear to share the same long-term dynamics. The difference lies in the variation in the speed at which the precession of Ω_G is performed: for B2, the speed is nearly zero all along the trajectory, except when the perihelion is going through a minimum of its Galactic oscillation. At this point, a sharp decrease in Ω_G was observed. For B3, precession occurred at a nearly constant speed all along the trajectory.

For a typical B2 comet, at the perihelion preceding the perihelion passage where the comet is observable, Ω_G has not yet experienced the sharp decrease, but it has weakly precessed from its initial value (which is close to 186°, i.e., the Galactic ascending node of the ecliptic). Consequently, always at this preceding perihelion, the orbital plane of the comet was close to the ecliptic plane with the perihelion moving back to the planetary region. These conditions favor the scattering of the comet’s semimajor axis by a giant planet. For our comet, this scattering resulted in an increase in its semimajor axis, which strengthened the Galactic tides. Hence, all evolutions are accelerated, and during the last orbital period, the perihelion distance of the comet “jumps” in the region where the comet is observable, while Ω_G undergoes a sharp decrease.

This helpful scattering by the giant planets that occurred at the previous perihelion passage is known as the Kaib and Quinn (KQ) effect (Kaib & Quinn 2008). Fouchard et al. (2017b) showed that such effects are mainly due to Neptune and Uranus; the scatter by Jupiter and Saturn is strong, which generally expels the comets from the Oort Cloud.

For a typical B3 comet, at the perihelion preceding the perihelion passage where the comet is observable, nearly half the period has been achieved by Ω_G. Meanwhile, the Galactic inclination and argument of perihelion performed a full period, such that when the perihelion is back to the planetary region, the orbital plane of the comet is close to the EE plane. One can find a continuity of behavior between B2 and B3, but what is important is the type of imprint they give to the final distributions.

Group B1 contains comets with a rather small semimajor axis (approximately 1000 au, against 7000 au for a typical comet of the B2 or B3 group) during the entire time span. Such a small semimajor axis implies that the typical timescale of the Galactic tides effect is much longer than the age of the Solar System. Hence, these comets are frozen, except for the scattering of the orbital energy caused by the giant planets. In the disk model, these frozen orbits correspond to orbits close to the ecliptic plane and perihelia close to the planetary region. Such comets are good candidates for supporting the KQ effect at the previous perihelion passage, which made them observable.

Fig. 1 Long-term evolution of observable comets. The time is given in millions of years before the time of observation. The open circles correspond to the initial position, the filled circles correspond to the original orbit at the perihelion preceding the observation, and the stars correspond to the original orbit when the comet is observable. The evolution is given by one point every 250 Myr. For each panel, the upper diagrams show the evolution versus time of the semimajor axis (dashed line) and perihelion distance (full line); the two lower diagrams are, from left to right, the motion of the comet in the (ΩG, cos iE) plane and in the (lG, sin bG) plane. The black region in the (ΩG, cos iE) plane is a forbidden region. In the (lG, sin bG) plane, the full black curve and the dashed black curve are, respectively, the trace of the ecliptic plane and the EE plane. The boldface letter in the upper-right corner of each panel refers to the long-term behavior described in the main text.

2.3 Chosen features

The present study aimed to detect some features typical of a disklike proto-Oort Cloud in the distributions of well-chosen parameters of long-period comets. Let and be the Galactic ascending node and the ecliptic inclination, respectively, both of which are computed at the perihelion passage preceding the perihelion where the comets are observable. In Fouchard et al. (2020), it has been shown that the KQ process supported by B1 and B2 comets, which required an orbital plane close to the ecliptic and a perihelion in the region of Uranus of Neptune, yields an excess of comets with in [90,180°], and excess of cos close to one, for previous perihelion distance q_p > 13 au.

Based on their different dynamic histories, should be closer to 180^° for B1 comets than for B2 comets. For the later could also be smaller than 90°. However, the number of comets in the databases did not allow us to distinguish between these two behaviors.

For comets in the B3 group, the orbital plane was close to the EE plane, which implied an excess of close to 0° and an excess of cos close to −0.5. These features should be detected in comets with q_p < 7 au.

These features are related to previous orbital elements. This choice was necessary, as explained by Fouchard et al. (2020) because drastic changes in the orbital parameters occurred during the last orbital period, canceling out the features of interest. Unfortunately, this is problematic when known comets are considered instead of synthetic comets coming from numerical simulations. This will be discussed in Sect. 3.3.

Nevertheless, Higuchi (2020) showed that when the perihelion is going through its minimal value of its Galactic oscillation, the direction of the perihelion is well conserved, to the contrary of Ω_G and i_G. Consequently, one can use this direction considering its latitude with respect to the ecliptic plane 6 and to the EE plane δ^*. We have the following relations: (1) (2)

where η = 60.19° is the inclination between the Galactic and the ecliptic plane, Ω₀ = 186.38° is the Galactic longitude of the ascending node of the ecliptic², lG and bG are the Galactic longitude and latitude of the perihelion of the comet.

Then the KQ effect for comets in the B1 or B2 group will be related to small values of | sin δ|, whereas comets in the B3 group that are observed close to the EE plane will be related to small values of | sin δ^*|. In addition, owing to the KQ effect, which increases the orbital energy, features related to the B1 and B2 groups are observed for larger semimajor axis than the features related to the B3 group. In Sect. 3.3 we describe our splitting of our comet sample.

Thus, we have two sets of parameters that will be considered: and . These two sets are complementary and do not correspond to the same osculating orbit. could be sufficient, but their large variation close to the perihelion where the comets are observable can degrade their reliability significantly, whereas are more stable but less complete. One notes that | sin δ| and | sin δ*| have the nice property of having an uniform distribution in [0, 1] when an isotropic distribution of the perihelion is considered.

We can then split the features of interest in two sets. The first contains the KQ features that consist in an excess of in [0,180°], a maximum of the cos distribution close to one, and a maximum of the | sin δ| distribution close to zero. The second are the EE features that consist in a maximum of the distribution close to 0°, a maximum of the cos distribution close to -0.5, and a maximum of the | sin δ*| distribution close to zero. In addition, a concentration of perihelia close to the ecliptic yields a maximum of the | sin δ*| distribution at sin ŋ ≈ 0.87, and similarly a concentration of perihelion close to the EE plane yields a maximum of the | sin δ| distribution at 0.87.

3 Databases of cometary orbital elements

In this study we consider only the comets that satisfy the following conditions: (1) comets whose perihelion distance q in the present apparition is smaller than 5 au, and (2) comets whose original semimajor axis a_orig is greater than 10⁴ au. Candidate comets were selected from two databases. We note that the cometary orbital elements stored in these databases do not differ from one another at a significant level (Ito & Higuchi 2020; Ito 2021).

3.1 Warsaw catalog (CODE)

Królikowska published a series of cometary catalogs called the Warsaw catalog of cometary orbits, also known as the CODE catalog. They are available in the VizieR database³ and on their own websites⁴. Details regarding the theoretical background and the implementation of their method are available in their series of publications: Królikowska (2014, 2020), Królikowska et al. (2014), Królikowska & Dybczynski (2017, 2018, 2020). As of 2023 January, the original orbital elements of 278 comets (along with their previous, current, next, and future orbits) are available with 1σ uncertainties. Królikowska’s online comet catalog stores a maximum of 12 orbital solutions for a single comet. In this study we consider only their “preferred” orbit. We note that one of the main specificities of this catalog is the strong attention given to the modeling of non-gravitational forces (e.g., Marsden et al. 1973). We discuss this point in Sect. 4. A total of 163 in this database comets fulfilled our criteria.

3.2 MPC ephemeris

The Minor Planet Center (MPC) has an online ephemeris on their website⁵. This ephemeris includes 2984 C/ comets as of December 1, 2022. Among all, 1500/2984 are non-SOHO comets, among the non-SOHO comets, 986/1500 are given the reciprocal of their original semimajor axis , and among these comets, 878/ 986 have their uncertainty values, . For some comets, several solutions for the best present orbit are proposed. We chose the first solution for every comet, as it always has the best quality available.

Unfortunately, the MPC ephemeris yields only the original values of the semimajor axes, but not of other orbital elements as the perihelion distances, ecliptic inclinations, arguments of perihelion and longitudes of the ascending node. However, it is well known that the planetary scattering in quasi-parabolic orbits affect mainly the semimajor axis (Šteins & Kronkalne 1964; Yabushita 1972). Therefore, in the present study we assume that the present orbital elements (except for semimajor axis) have the same values as their original counterparts. A total of 231 comets fulfilled the inclusion criteria in the CODE database.

3.3 Data processing

The following procedure was developed to account for uncertainties in the original orbital energy, which are crucial for the origin of long-period comets. In both databases, the original, or assimilated, orbital elements were provided for each comet. These elements are called the nominal orbital elements. Each comet was then cloned 199 times; all clones had the same orbital elements as the nominal one except for the original orbital energy, which was computed randomly using a Gaussian distribution with the mean corresponding to the nominal orbital energy and the standard deviation equal to the uncertainty. Using the comet clones in the database, 200 different comet samples were built. Each sample contained the same comets as the database, but with different original orbital energies. Hereafter, a sample from a database is called a clone of the database.

A comet in a clone of a database is a “new” long-period comet when its original perihelion distance is smaller than 5 au, and its original orbital semimajor axis is larger than 10⁴ au or the original orbit is hyperbolic. This implies that the number of new long-period comets in the different clones from the same database may not be always the same.

As it was discussed in Sect. 2.3, some features of interest can be detected considering the distributions of and cos These elements were computed at the perihelion passage preceding the perihelion where the comet was observable. These data were systematically stored in the numerical experiments conducted by Fouchard et al. (2020). This is clearly not the case for real comets, mainly because it is not possible to reconstruct precisely the effect of passing stars during the last orbital period. In the present study, the orbital elements at the previous perihelion, in particular the perihelion distance q_p ecliptic inclination and Galactic longitude of the ascending node , were obtained using a backward integration from the observed perihelion until the first passage at perihelion considering only the Galactic tides. Królikowska & Dybczyński (2020) proposed such orbital elements; however, they used a different model for the backward integration. To ensure consistency with the strategy used for the synthetic models, we decided to compute our own previous orbital elements.

Neglecting stellar perturbations during this one orbital period backward propagation could be a drawback for our study. However, using the data from the synthetic models, the relative error made on the previous perihelion distance was smaller than 20% for half of the comets with an original semimajor axis smaller than 40 000 au. In addition, Fouchard et al. (2020) observed that our features of interest on the and cos distributions were still present even if passing stars were non including for the computation of the previous orbital elements. Consequently, our strategy to investigate the presence or not of some features of interest in the and cos distributions is reliable. The direction of perihelion used to compute δ and δ^* was that of the original orbit, that is, it was directly obtained from the databases.

Each clone of a database yields distributions for cos , , | sin δ| and | sin δ^*|. Hence, these 200 clones allowed the computation of the means and standard deviations for each distribution.

3.4 Splitting the samples

As discussed in Sect. 2, in order to highlight the different features of interest, the sample of comets has to be split according to either q_p for the and cos distributions or the original orbital energy for the | sin δ| and | sin δ^*| distributions.

Because the number of comets in the databases was much lower than the number of synthetic comets obtained from the disk or isotropic models, only three zones were defined for q_p: (i) q zone #1: q_p < q_c; (ii) q zone #2: q_c < q_p < 31 au; and (iii) q zone #3: q_p > 3l au. Here, we used the value q_c = 7 au. Therefore, according to Fouchard et al. (2020), comets with 7 < q_p < 13 au were merged with comets with 13 < q_p < 31 au. This choice allows for the best isolation of the EE features. That is, the EE features should be observed in distributions computed from comets in zone #1, and the KQ features in the distributions computed from the comets in zone #2.

The three regions, hereafter referred as q zones, required the computation of the previous perihelion distance q_p. It was natural to split our comet sample to investigate the distributions of and cos . It was not possible to define the q zone for a comet on an original hyperbolic orbits because the previous perihelion passage was not defined. These comets were then added to the q zone #3.

The q zones are not adequate for the other parameters δ and δ^*, since they are obtained directly from the original orbits. However, we still needed to split the sample to identify the features of interest. To achieve this, we defined three different regions according to the original semimajor axis. So that these regions are as close as possible to the q zones, the border of the region are defined as follows. We considered the orbital energy z = −1/a and computed the mean values, , of the nominal orbital energy of comets in q zone #i (hyperbolic comets are excluded for q zone #3). Then the medium value between and , z_m, and the medium value between and , z_M, were computed. Thus, the z zones are defined as follows: (i) z zone #1 contains comets with original orbital energy between −10⁻⁴ au⁻¹ and z_m, (ii) z zone #2 contains comets with original orbital energy between z_m and z_M, and (iii) z zone #3 contains comet with original orbital energy greater than z_M. As for the q zones, the EE-features should be observed in the distributions obtained from z zone #1 and KQ-features from comets in z zone #2.

A final comment on data processing: the only difference between the clones of the same database was the original orbital energy of the comets. This has a direct consequence on q_p, , and cos , but it does not affect the value of | sin δ| or | sin δ^*|. However, the distributions of the latter two elements will still be affected because the comets in the different z zones will be different.

With regard to our synthetic samples coming from the disk and the isotropic models, the previous orbital elements are computed just as they are for the comets in the database, that is, neglecting passing stars. In addition, we decided not to apply any statistical process to generate the distributions from the models. The number of comets used in our simulations was huge (tens of millions) and the final sample size was rather large. In addition, the original orbital elements were not subject to any uncertainty, except for a small numerical one. We could have introduced artificial errors in the orbital energy; however, such errors are strongly dependant on the condition of observation of a comet, which is not a dynamical problem. Therefore, we preferred not to introduced such an artificial error in the synthetic data.

4 Results

4.1 Preliminary investigations

Our first exploration of the data involved counting the new observable comets. The number of comets in the different q and z zones was also computed. In addition, when origin of long-period comets is considered, the original orbital energy is the key parameter; thus, the mean value of the orbital energy was computed for the total group and for the comets in the different zones. The original hyperbolic orbits were not considered when computing the mean orbital energy. All such comets were added to zones #3.

The results obtained for the two synthetic models and the two database are listed in Table 1. We recall that the number of observable comets, or data, in the synthetic models must be converted to a flux per million years in order to take into account the orbital period of the comets. That is the contribution of a comet to the flux scales as P⁻¹ where P is the orbital period of the comet. A comparison of the fluxes was not be performed here. Such a study was already performed in Fouchard et al. (2017a).

It should be noted that the semimajor axes corresponding to the mean orbital energies for the q zones #1 and #2 were very similar for all sample of comets (databases or synthetics). This is a direct consequence of the fact that q_p depends mainly on the orbital energy. It also depends on the Galactic tides model used; however, the same model was used for our synthetic comets as for the comets in the databases for the computation of q_p. For the q zone #3, the mean orbital energy was higher for the models than for the databases. This will be discussed later in this paper.

It is important to note that for the databases the computed values of q_p may differ from the real perihelion distances of the comets at their previous perihelion passage. This is not only because passing stars were neglected, but also because the Galactic tides model we used to compute the previous orbital elements might not correspond to the real ones, particularly the density of the Galactic disk in the solar neighborhood.

Regarding the distribution in the different zones, it should be noted that the original orbital energy of the comets in the MPC database was globally larger than of the CODE catalogs. One reason for this could be that in the CODE catalog non gravitational forces are modeled with careful attention, and introduced more frequently than in the MPC database. Indeed, the best solutions for the osculating orbit in CODE were obtained by introducing non-gravitational forces (NGFs) for 50% of the comets. In contrast, NGFs were introduced for only 6% of the best solutions in MPC database. Because non-gravitational forces generally increase the orbital energy, the original semimajor axis of a non-gravitational solution is generally smaller than that of a purely gravitational solution (Królikowska 2006). A singular consequence of this is that the number of hyperbolic comets was much larger for the MPC database than for the CODE catalog.

Another important difference between the databases and the models is that the outer zones are less densely populated for the synthetic models than for the databases. This may be attributed to several factors. For instance, the real outer bound of the Oort Cloud may be more stable than that of the synthetic models. Hence, stellar perturbations were overestimated in our synthetic models (on gigayear timescales). Or, on the contrary, the stellar perturbations in the models may have been underestimated, which means that they are less efficient into injecting comets in the regions from which the Galactic tides can make them observable. Another reason could be that the original semimajor axes stored in the database were still overestimated despite the introduction of non-gravitational forces. This last point could also explain why the number of hyperbolic comets was higher in the databases than in the synthetic models and also why the median orbital energy for q zone #3 is higher for the databases than for the models. So far, we have not been able to determine which of these reasons are more likely. The answer resides in part in a better understanding of the neighborhood of our Solar System on a long timescale, in particular by astronomical observations.

The proportion of comets in q zone #1 showed the largest fluctuation among the sample considered. The differences between the CODE and the MPC databases may arise again from the orbital energy, which is generally larger for the MPC database than for the CODE database. Thus, the variations in the perihelion distance in one orbital period are larger for MPC comets, which implies that fewer comets are found in q zone #1. For the models, considering again that the tides are less powerful for comets in q zone #1, comets in this zone have the perihelion that may stay near the planetary region of the Solar System for many periods. A large ecliptic inclination is preferable for such comets to survive planetary scattering. Clearly, the isotropic model contains a larger number of such comets than the disk model; hence, more comets are found in q zone #1 in the isotropic model than in the disk model. Comparison between models and databases is difficult because the justifications would differ according to which database–model pairs are considered.

4.2 cos and distributions

We first analyzed the pair of variables cos and . Figure 2 shows the position of the observable comets in the (, cos ) plane for the three q zones. The main features observed in the two models were already discussed by Fouchard et al. (2020). As shown in Fig. 2, the comets were concentrated along the i_G = 90° curve. Indeed, it has already been shown that this is a consequence of Galactic tides injection, and is almost independent of the route taken by the comets toward the observable region (see Higuchi et al. 2007, for instance). The nonuniform distribution of comets along this line is related to this route, and it seems that some nonuniformity is also present in the data stored in the databases.

We then considered the distribution of and cos . Using the strategy outlined in Sect. 3.3, we considered 200 clones for each database. The means and standard deviations for each distribution were computed. As mentioned in Sect. 3.3, no clone were computed for the synthetic models; thus, the distributions were directly obtained from the data. The results are shown in Fig. 3.

First, we clearly define the dynamical features of interest. To achieve this, the distribution of the models was considered. For a feature to be meaningful, it must be clearly visible in the distribution obtained from the disk model and not present in the distributions obtained from the isotropic model. With this constraint, the following features were identified: (i) , an EE feature that corresponds to a maximum of cos around –0.5 in q zone 1, and (ii) E₁ and E₂, two KQ features that correspond respectively to a maximum of in the interval [90°, 180°] and a maximum of cos at 1 in q zone 2. These features are denoted by blue and orange patches, respectively, in Fig. 3.

This raises the most important question in the present study, of whether the above three features are present in the orbital distributions of the comets stored in the databases. As shown in Fig. 3, it is not possible to provide a clear answer. However, the following observations were made.

For , the typical EE feature is rather clear on the MPC database but not on the CODE database. However, Fouchard et al. (2017b) show that a comet has a greater chance of surviving planetary scattering when its ecliptic inclination at the perihelion passage preceding the observed perihelion is close to 180° (Fouchard et al. 2017b). Hence, for an isotropic source of comets, a maximum of the cos distribution should located in −1. This is the case for the isotropic model. However, none of the database cos distribution shows such a maximum; thus, the comets in the databases in q zone #1 are not coming from an isotropic reservoir.

For E₁, in the distributions, the KQ feature is clearly present in the CODE catalog. In addition, this feature is present in the MPC catalog but is less striking. Finally, the KQ feature E₂ is observed for the CODE and MPC databases.

We found several interesting characteristics in the databases we used. First, in q zone #1, whereas for the two models a lack of comets with small ecliptic inclination was found, such orbits were common in the two databases. This is striking: as already discussed in Sect. 4.1, comets in this zone should have large ecliptic inclinations in order to survive planetary scattering. An excess of comets with a moderate semimajor axis (smaller than 20 000 au), perihelion in the region of the planetary region, and a small ecliptic inclination may be explained by either observational biases or the models failing to reproduce the dynamical routes that lead to such final orbits.

In fact, orbital data from observations must contain comets that are also coming from the scattering disk, and such long-period comets would likely end in q zone #1 with small ecliptic inclination⁶. This route is not present in any of the models. This is evident for the isotropic model, which does not contain any disk in its initial stage. This is also the case in the disk model, which initially contains orbits with much larger eccentricities than in the typical scattered disk. Hence, even our disk model could not reproduce the contribution of the scattered disk to the population of new long-period comets.

Fig. 2 Position of the comets obtained through the two databases and the two models (CODE, MPC, disk, and isotropic, from left to right) in the (, cos ) plane for the three q zones from top to bottom (see the main text for details). The double-hashed area is a forbidden region, and the solid black line corresponds to the Galactic inclination iG = 90°.

Fig. 3 Distributions of (left column) and cos (, cos ) (right column) for the three different q zones from top to bottom (see the main text for details). Red lines (with error bars) correspond to the comets stored in the CODE databases and blue lines to the MPC database. The black and gray lines (without error bars) correspond to the disk and isotropic models, respectively. The patches indicate the KQ (blue) and EE (orange) features of interest.

4.3 Perihelion distributions

We next considered the direction of the perihelion. As discussed in Sect. 2.3, its direction is well conserved when the perihelion distance reaches its minimum during its Galactic oscillation. Figure 4 shows the position of the perihelion direction of comets on the (l_G, b_G) plane where l_G and b_G are the Galactic longitude and latitude, respectively. The panels are categorized into the three z zones. The ecliptic and EE planes are also plotted. At first glance, it seems that for the two databases and the two models, the directions of perihelion are denser in the ecliptic and EE neighborhood.

In order to see if some features of interest are present, the distributions of | sin δ| and | sin δ^*| are considered. The same methodology as for the and cos distributions was used. Thus, the difference between two clones of a database consists only in the populations of the z zones because the direction of the perihelion is not affected.

Figure 5 shows the distributions of | sin δ| and | sin δ^*| for the two databases and the two models for the three different z zones. As shown in Fig. 3, patches indicate the feature of interest where the disk distributions highlight a typical and distinguishable EE or KQ features. These features are, for the KQ feature, a maximum of the |δ| distribution at zero for z zone #2. For the EE features, they are a maximum of the |δ^*| distributions toward zero for z zones #1 and #2 and a maximum of the | sin δ| distributions in one in the same two zones. These last two maxima are related to the concentration of perihelion along the EE as discussed in Sect. 2.3.

Notably, EE features were also present in z zone #2. Indeed, it is more difficult to separate the two types of features using the original semimajor axis rather than q_p. All the above features are present in the two databases’ distributions.

4.4 Comparison to a complete sample of observed comets

It is not within the scope of the present study to discuss observational biases in detail. However, such bias may have been present in our samples. A sample with no biases is a complete sample, containing all the existing objects that satisfy the constraints defining the sample. Fouchard et al. (2017a) observed that since 1998, the beginning of the LINEAR survey⁷, the known flux of comets with a total absolute magnitude smaller than 11 and a perihelion distance smaller than 4 au has not increased despite the advent of new surveys such as the Catalina Sky Survey⁸ and Pan-STARRS⁹. In view of this, the sample of known long-period comets with total magnitude smaller than 11 and perihelion distance smaller than 4 au was considered complete.

In the present study, the same strategy was applied: we picked the comets from the databases that were discovered after January 1, 1998, with a total magnitude smaller than 11, and which passed less than 5 au from the Sun (that is, 1 au higher than in the 2017 study). The magnitudes were taken from the MPC. With the above filters, the numbers of comets is 53 (2 of which are hyperbolic) in the CODE database and it 69 (7 of which are hyperbolics) for the MPC database.

In q zone #1 and #2, the number of comets is now 36 in the CODE database, and 39 for the MPC database, compared to 112 and 115 for the complete database, respectively; thus, the filtered samples contain only approximately one-third the comets contained in the full ones. This clearly have affected the statistical reliability of the results.

The distributions obtained using the filtered database are shown in Fig. 6. First, we consider the comet population in q zone #1 with small a ecliptic inclination. A shown in Fig. 3, the unfiltered databases contained an excess of such comets. By using the filtered databases, this excess was significantly reduced. Thus, some observational bias should exist in the complete sample, meaning that the number of comets along the ecliptic plane is overestimated. However, this does not rule out the possibility that the databases must contain comets coming from the scattered disk, and this flux was not present in our models.

Regarding the features of interest, except for the EE feature, all the others are distinguishable on the distributions obtained using the filtered MPC database. Compared to Figs. 3 and 5, the features are all weakly altered, but they are still distinguishable.

The situation is different for the CODE database, in which many of the features are invisible when filtered data are adopted. This is particularly true for EE features. However, as shown in Fig. 6, although the features are unclear, the distributions using the filtered CODE database were clearly more similar to the distributions obtained with the disk model than with the isotropic model. This alteration in the results could be explained by considering the quality of the cometary orbits. The orbit quality index for the CODE catalog (Królikowska & Dybczyński 2013) is slightly different from that used for the MPC (Marsden et al. 1978). Specifically, the best 1A class of the MPC database corresponds roughly to the 1a+, 1a and 1b class of the CODE database. The proportions of 1a+, 1a, and 1b quality in the full CODE database were approximately 60%, whereas the proportion of the 1A quality in the MPC database was 48%. We cannot directly compare these proportions because they depend on the manner in which the data are used. However, when we consider the filtered databases, the best quality orbits for the CODE catalogs reached 43%, whereas they were 65% for the MPC database. This indicates that, the filtered data for the CODE catalogs are of lower quality than those for the whole catalog, whereas the opposite is true for the MPC database. Notably, 48 comets were included in both filtered databases. That is, almost all the comets in the CODE database were also in the MPC database¹⁰. Thus, the orbit quality is intrinsic to each database.

Fig. 4 Position of the comets for the two databases and the two models (CODE, MPC, disk, and isotropic, from left to right) in the (, sin ) plane for the three z zones from top to bottom (see the main text for details). The solid black line corresponds to the ecliptic plane, and the dotted black line corresponds to the EE plane.

4.5 Statistical tests

To obtain a comprehensive picture of our results, the following experiments were performed. For each feature of interest, E₁₋₃ for the KQ ones and for the EE ones, the number of distributions among the 200 ones (the nominal one plus the 199 clones) for which the maximum of the distribution coincides with the feature is counted. The proportion obtained here can be compared to the proportion that is expected from a uniform distribution, 0.25. This value comes from the fact that for any distributions four bins are considered. The results obtained for each database, filtered and not, are shown in Table 2.

We note, however, that with this method a distribution cannot have a maximum that coincides with both E3 and . And indeed, one notes that E₃ is almost always detected, whereas is not. However, this does not imply the absence of this feature, just that the feature E3 dominates.

From Table 2 we can see that the CODE database matches all three KQ features for both the full database and the fll-tered database. Three EE features were matched using the full database but only one using the filtered one.

The MPC database had slightly different results for all three KQ features using the full database but only two with the filtered database. Three EE features were detected in the full database and only two in the filtered database.

Fig. 5 Distributions of | sin δ| (left column) and | sin δ*| (right column) for the three different z zones from top to bottom (see the main text for details). The red lines (with error bars) correspond to the CODE database and blue lines to the MPC database. The black and gray lines (without error bars) correspond to the disk and isotropic models, respectively. The patches indicate the KQ (blue) and EE (orange) features of interest.

5 Conclusions

Fouchard et al. (2020) show that a memory of Oort Cloud’s initial structure can be detected in the distribution of observable comets after four billion years. Their study was performed considering two models: a disk model, in which the initial synthetic orbits had a small inclination and a perihelion below 50 au from the Sun, and an isotropic model, in which the initial orbits were fully thermalized.

In the present study, we attempted to detect such traces in the observational databases of cometary orbits. We used two databases of long-period comets that were constructed independently. From each database, we selected new long-period comets that had an original semimajor axis greater than 10⁴ au or were on their original hyperbolic orbit, and with perihelion distances smaller than 5 au. In addition, to cancel out the effect of observational biases, each database was filtered in order to obtain a sample of comets whose completeness is much higher than that of the entire database.

Thus, we attempted to retrieve the typical dynamical features that the disk model would yield in the observated sample of long-period comets. Two sets of features were defined: one linked to the KQ process, and the other linked to the orbital plane close to the EE plane. This plane corresponds to a 180° rotation of the ecliptic plane around the Galactic poles.

The KQ process consists of planetary scattering at the perihelion that precedes the observations, which helps the Galactic tides make comets observable. Two dynamical routes typical of the disk model are linked to this process. The first one, called B1, consists of orbits in the initial disk that are nearly frozen during the entire propagation and only at the very end does the scattering of the orbital energy by Neptune lead to the final KQ process. The second one, called B2, also consists of orbits that start in the initial disk; however, during propagation the eccentricity, Galactic inclination, and argument of perihelion make a full Galactic period, whereas the Galactic longitude of the ascending node slowly decreases until the perihelion where the KQ process occurs.

A third dynamical route typical of the disk model was also identified to generate the EE features. This is very similar to route B2 except that the Galactic longitude of the ascending node is decreasing smoothly such that when the perihelion is back to the planetary region, the longitude of the ascending node has accomplished half a period. In this case, the final orbital plane is close to the EE plane.

Our results show that both features (KQ and EE) were detectable in the distributions obtained from the cometary orbital databases. Such a detection is of great importance for understanding the shape of the current Oort Cloud. Indeed, the detection of the KQ feature indicates that the B1 and B2 dynamical routes might be common for comets in the databases, and the detection of EE features indicates that the B3 dynamical route might also be common. Hence, a disk should be present in the current and initial Oort Cloud at the end of its formation process. This result confirms what was already proposed by Fernández (2002), who observed an excess of prograde comets among “new” long-period comets. The size of the disk can be estimated using the results obtained by Fouchard et al. (2020). Indeed, the B1 route is related to a disk of objects with a semimajor axis at about 1000 au, whereas the B2 and B3 routes are related to a disk of objects with a semimajor axis at about 7000 to 10000 au.

It is important to note that these limits come from the disk model. As such, they depend on the Galactic tide parameters used in the model, in particular on the density of the Galactic disk, ρ_o. We took ρ_o.1 M_⊙ • pc⁻³, but using a Galactic disk density of 0.15 M_⊙ • pc⁻³, as in Nesvorný et al. (2017), would imply a reduction of the disk from 10 000 to 5–7 000 au.

In any case, the presence of such a disk in the Oort Cloud places constraints not only on the long-term dynamics, but also on the formation of the initial Oort Cloud. Indeed, to build the initial Oort Cloud, two main processes are considered: scattering by the giant planets, eventually including a migration of the planets (Brasser & Morbidelli 2013; Vokrouhlický et al. 2019), and scattering caused by the stellar cluster where the Sun was born (Fernández & Brunini 2000; Kaib & Quinn 2008; Brasser et al. 2012; Wiśniowski et al. 2022). The capture of comets from other stars may have also contributed to the formation of the Oort Cloud (Levison et al. 2010). The estimated shape of the initial Oort Cloud differs according to the role assigned to each process in its formation. For instance, Rickman et al. (2023) show that even the orbit of the giant planets could be destabilized by scattering from the birth stellar cluster of the Sun. This highlights the importance of the presence of a disk in the Oort Cloud with regard to the process at work during its formation.

Although some KQ and EE features were detected in the distributions obtained from the databases, the fit with the disk model was far from perfect. Clearly, our disk model is a rough simplification of reality. No formation Oort Cloud process leads to the perfect initial disk considered in the model. In addition, the billions of years of dynamical evolution after the formation is important: the Galactic tides, history of passing stars, and encounters with giant molecular clouds (Jakubík & Neslušan 2008), all of which have shaped the present Oort Cloud.

Our results should be completed by other investigation on populations of objects coming from the region beyond Neptune that may contain some fingerprint of the structure of the outer Solar System. For instance, Saillenfest et al. (2019) show that a transition between the domination of the effect of the Galactic tides and the domination of the giant planets on the long-term dynamics of trans-Neptunian objects is observed at approximately 1000 au. This could have some implication for the distribution of other families of small bodies (e.g., “old” long-period comets, Halley-type or Jupiter family comets, and Centaurs).

Two other aspects of our study should also be improved. The first regards the Galactic tides environment on a long timescale, from which the Galactic tides model is built. This is crucial because the features observed are directly related to the model. With the help of Gaia releases (Gaia Collaboration 2021), this is now within reach. The second aspect concerns the completeness of the cometary databases. Indeed, some bias was present when the full databases were used. The filtered databases were more complete in the sense that they were less subject to bias. Unfortunately, the quality of the results was altered by the number of data points used and the quality of the orbits, which were globally worse than that of the full sample for the CODE database. To solve this problem, a detailed investigation of the biases present in the main surveys, which are the main discovers of comets, should be conducted. However, here again, the future will be our best friend, with more reliable observations anticipated, mainly from the Simonyi Survey Telescope at the Vera C. Rubin Obser-vatory¹¹, whose full survey operations are expected to begin in 2024.

Fig. 6 Same as Figs. 3 and 5 but with the databases restricted to supposedly complete subsets of comets.

Acknowledgements

We are grateful to Melaine Saillenfest, Vacheslav Emel’Yanenko and Giovanni Valsecchi for helpful discussions that allowed to improve the quality of the present paper. We are also grateful to Hans Rickman and a second anonymous referee for their helpful comments. We acknowledge NASA’s Astrophysics Data System Bibliographic Services. T.I. acknowledges the research support from the JSPS Kakenhi Grant (JP18K03730/2018–2021). The cometary orbital elements retrieved from various databases were processed and reorganized for analysis at the Center for Computational Astrophysics (CfCA), National Astronomical Observatory of Japan. This work was supported by the Programme de Planétologie (PNP) of CNRS/INSU, co-funded by CNES.

References

All Tables

All Figures

Fig. 1 Long-term evolution of observable comets. The time is given in millions of years before the time of observation. The open circles correspond to the initial position, the filled circles correspond to the original orbit at the perihelion preceding the observation, and the stars correspond to the original orbit when the comet is observable. The evolution is given by one point every 250 Myr. For each panel, the upper diagrams show the evolution versus time of the semimajor axis (dashed line) and perihelion distance (full line); the two lower diagrams are, from left to right, the motion of the comet in the (ΩG, cos iE) plane and in the (lG, sin bG) plane. The black region in the (ΩG, cos iE) plane is a forbidden region. In the (lG, sin bG) plane, the full black curve and the dashed black curve are, respectively, the trace of the ecliptic plane and the EE plane. The boldface letter in the upper-right corner of each panel refers to the long-term behavior described in the main text.

In the text

	Fig. 2 Position of the comets obtained through the two databases and the two models (CODE, MPC, disk, and isotropic, from left to right) in the (, cos ) plane for the three q zones from top to bottom (see the main text for details). The double-hashed area is a forbidden region, and the solid black line corresponds to the Galactic inclination iG = 90°.
In the text

Fig. 3 Distributions of (left column) and cos (, cos ) (right column) for the three different q zones from top to bottom (see the main text for details). Red lines (with error bars) correspond to the comets stored in the CODE databases and blue lines to the MPC database. The black and gray lines (without error bars) correspond to the disk and isotropic models, respectively. The patches indicate the KQ (blue) and EE (orange) features of interest.

In the text

	Fig. 4 Position of the comets for the two databases and the two models (CODE, MPC, disk, and isotropic, from left to right) in the (, sin ) plane for the three z zones from top to bottom (see the main text for details). The solid black line corresponds to the ecliptic plane, and the dotted black line corresponds to the EE plane.
In the text

Fig. 5 Distributions of | sin δ| (left column) and | sin δ*| (right column) for the three different z zones from top to bottom (see the main text for details). The red lines (with error bars) correspond to the CODE database and blue lines to the MPC database. The black and gray lines (without error bars) correspond to the disk and isotropic models, respectively. The patches indicate the KQ (blue) and EE (orange) features of interest.

In the text

	Fig. 6 Same as Figs. 3 and 5 but with the databases restricted to supposedly complete subsets of comets.
In the text