# What long-period comets tell us about the Oort Cloud

Marc Fouchard, Arika Higuchi, Takashi Ito

► To cite this version:

Marc Fouchard, Arika Higuchi, Takashi Ito. What long-period comets tell us about the Oort Cloud. Astronomy & Astrophysics - A&A, 2023, 676, 10.1051/0004-6361/202243728. insu-04849157

HAL Id: insu-04849157

https://insu.hal.science/insu-04849157v1

Submitted on 19 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# What long-period comets tell us about the Oort Cloud

Marc Fouchard<fnref n="1" />, Arika Higuchi<fnref n="2" />, and Takashi Ito<fnref n="3" />,<fnref n="4" />

Received 7 April 2022 / Accepted 14 June 2023

## ABSTRACT

*Context.* The Oort Cloud is located in the farthest outskirts of the Solar System, extending to a heliocentric distance of several tens of thousands of au, and remains the last region of the Solar System where no object has been detected. Thus, all our knowledge of the Oort Cloud has been deduced from the observed long-period comets that are thought to originate from it.

*Aims.* We aimed to retrieve valuable information that might be hidden in the orbital distributions of the observed long-period comets. Such information will allow us to impose constraints not only on the present shape of the Oort Cloud but also on its initial shape 4.5 Gyr ago. This has direct implications for the scenario proposed for its formation.

*Methods.* We used two different databases of long-period comets. First, we calculated the distribution of orbital elements that might carry valuable information about the shape of the Oort Cloud. Then, we compared the distribution with that obtained from two synthetic samples of observable comets. These samples correspond to two considerably different initial configurations: one is a disk model, where we consider a swarm of comets with orbits aligned to the ecliptic plane and with a cometary perihelion close to the giant planets. The other is an isotropic model, where we consider a fully isotropic and thermalized initial distribution of comets.

*Results.* The comparison revealed that the databases contained several features that were in better agreement with the disk model than with the isotropic model. The Oort Cloud contained an initial disk of objects with perihelia close to the planetary region of the Solar System and aphelia extending out to roughly 20 000 au. Some parts of this disk likely remain in the present Solar System. However, the fit to the disk model is poor. The discrepancy between the observational and synthetic results indicates that some dynamical processes in the current Oort Cloud were not included in either model.

*Conclusions.* This initial shape of the Oort Cloud implies that planetary scattering was crucial during its formation. In addition, the fact that some dynamical features are still detectable 4.5 Gyr after the cloud formation imposes constraints on the role of exosolar effects, such as giant molecular clouds, Galactic tides, and the stellar cluster surrounding the Solar System at the time of its formation.

**Key words.** celestial mechanics – astronomical databases: miscellaneous – methods: data analysis

## 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 \sim 10^4-10^5$ 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.

<footnote n="1">¹ IMCCE, Observatoire de Paris, 77 av. Denfert-Rochereau, 75014 Paris; Observatoire de Lille, Université de Lille, 1 imp. de l'Observatoire, 59000 Lille, France
e-mail: marc.fouchard@univ-lille.fr</footnote>

<footnote n="2">² Kyoto Sangyo University, Motoyama, Kamigamo, Kita-ku, Kyoto 603-8555, Japan</footnote>

<footnote n="3">³ Center for Computational Astrophysics, National Astronomical Observatory of Japan, Osawa 2-21-1, Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote n="4">⁴ Planetary Exploration Research Center, Chiba Institute of Technology, 2-17-1 Tsudanuma, Narashino, 275-0016 Chiba, Japan</footnote>
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 param-
eters – a density of the local Galactic disk ρ₀ = 0.1 M☉ · pc⁻³
and an angular velocity of the Sun around the Galactic center
Ω₀ = −26 km s⁻<fnref n="1" /> kpc⁻<fnref n="1" /> – 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 char-
acteristics of the disk model include: uniform distributions for
the orbital energy z = -1/a with a semimajor axis a between
500 au and 20 000 au, the perihelion distance, q, between 3
and 50 au, the ecliptic-inclination, iE, 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 distribu-
tions: for the orbital energy, z, with semimajor axis a between
500 au and 50 000 au; for cos i<sub>G</sub> between -1 and 1, where i<sub>G</sub> is
the Galactic inclination; and for the Galactic argument of peri-
helion, ω<sub>G</sub>, Galactic longitude of the ascending node, Ω<sub>G</sub>, and
mean anomaly, M, all between 0° and 360°. For the eccentric-
ity, 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 consid-
ered. 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 <fnref n="1" /> orbital energy, z = -1/a, was
greater than 10⁻⁴ au⁻<fnref n="1" />, then the comet was considered a “new”
observable long-period comet. The advantage of the five stop-
ping 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 eccentric-
ity, 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 perihe-
lion 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 per-
formed a single period, whereas the longitude of the ascending
node performed half a period. Such comets returned to the plane-
tary region with an orbital plane close to the empty ecliptic (EE)
plane (Higuchi 2020), which was obtained by rotating the eclip-
tic 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 con-
dition, 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 origi-
nal 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 eclip-
tic 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<sub>p</sub>. Indeed, the B3 group was
mainly found for q<sub>p</sub> < 7 au, whereas the B1 and B2 group mainly
found for 13 < q<sub>p</sub> < 31 au.

Details are provided here on the dynamic aspect of each
group. From Higuchi (2020) it appears that when the perihe-
lion distance *q* is making a full period together with the Galactic
inclination *i*<sub>G</sub> and the Galactic argument of perihelion, ω<sub>G</sub>, the
Galactic ascending node Ω<sub>G</sub> performs half a period; second
when the perihelion distance is going through its minimum, Ω<sub>G</sub>
and i<sub>G</sub> change very quickly. In addition, for an orbit starting with
small inclination with respect to the ecliptic, Ω<sub>G</sub> 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 per-
ihelion 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

<footnote n="1">¹ The original orbital elements for a given perihelion passage of a comet in the planetary region of the Solar System correspond to the barycentric orbital elements of the comet at 250 au from the Sun before the perihelion passage.</footnote>
<figure><img src="image_2.png" /><figcaption>**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 ($\Omega_G$, cos $i_E$) plane and in the ($l_G$, sin $b_G$) plane. The black region in the ($\Omega_G$, cos $i_E$) plane is a forbidden region. In the ($l_G$, sin $b_G$) 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.</figcaption></figure>

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 $\Omega_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 $\Omega_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.

## 2.3. Chosen features

The present study aimed to detect some features typical of a disk-like proto-Oort Cloud in the distributions of well-chosen parameters of long-period comets. Let $\Omega_{G_p}$ and $i_{E_p}$ 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 $\Omega_{G_p}$ in [90, 180°], and excess of cos $i_{E_p}$ close to one, for previous perihelion distance $q_p > 13$ au.

Based on their different dynamic histories, $\Omega_{G_p}$ should be closer to 180° for B1 comets than for B2 comets. For the later $\Omega_{G_p}$ 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 $\Omega_{G_p}$ close to 0° and an excess of cos $i_{E_p}$ 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 $\Omega_G$ and $i_G$. Consequently, one can use this direction considering its latitude with respect to the ecliptic plane $\delta$ and to the EE plane $\delta^*$. We have the following relations:

$$ \sin \delta = \cos \eta \sin b_G - \sin \eta \cos b_G \sin(l_G - \Omega_0), \quad (1) $$

$$ \sin \delta^* = \cos \eta \sin b_G - \sin \eta \cos b_G \sin(l_G - \Omega_0 - \pi), \quad (2) $$

where $\eta = 60.19°$ is the inclination between the Galactic and the ecliptic plane, $\Omega_0 = 186.38°$ is the Galactic longitude of the ascending node of the ecliptic$^2$, $l_G$ and $b_G$ 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 \delta|$, whereas comets in the B3 group that are observed close to the EE plane will be related to small values of $|\sin \delta^*|$. In addition, owing to the KQ effect, which

<footnote>$^2$ These value comes from the resolution of the XI IAU General assembly in 1961 where the north Galactic pole and the origin of the Galactic longitude were defined. See also Leinert et al. (1998).</footnote>
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: ($\Omega_{G_p}$, $\cos i_{E_p}$) and ($|\sin \delta|, |\sin \delta^*|$). These two sets are complementary and do not correspond to the same osculating orbit. ($\Omega_{G_p}$, $\cos i_{E_p}$) could be sufficient, but their large variation close to the perihelion where the comets are observable can degrade their reliability significantly, whereas ($|\sin \delta|, |\sin \delta^*|$) are more stable but less complete. One notes that $|\sin \delta|$ and $|\sin \delta^*|$ have the nice property of having a 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 $\Omega_{G_p}$ in [0, 180°], a maximum of the $\cos i_{E_p}$ distribution close to one, and a maximum of the $|\sin \delta|$ distribution close to zero. The second are the EE features that consist in a maximum of the $\Omega_{G_p}$ distribution close to 0°, a maximum of the $\cos i_{E_p}$ distribution close to -0.5, and a maximum of the $|\sin \delta^*|$ distribution close to zero. In addition, a concentration of perihelia close to the ecliptic yields a maximum of the $|\sin \delta^*|$ distribution at $\sin \eta \approx 0.87$, and similarly a concentration of perihelion close to the EE plane yields a maximum of the $|\sin \delta|$ 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^4$ 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<fnref n="3" /> and on their own websites<fnref n="4" />. 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 & Dybczyński (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\sigma$ 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<fnref n="5" />. 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 ($a_{orig}^{-1}$), and among these comets, 878/986 have their uncertainty values, $\delta(a_{orig}^{-1})$. 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^4$ 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 $\Omega_{G_p}$ and $\cos i_{E_p}$. 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 $i_{E_p}$ and Galactic longitude of the ascending node $\Omega_{G_p}$, 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 $\Omega_{G_p}$ and $\cos i_{E_p}$ distributions were still present even if passing stars were non

<footnote n="3">³ https://vizier.u-strasbg.fr</footnote>

<footnote n="4">⁴ https://pad2.astro.amu.edu.pl/comets/</footnote>

<footnote n="5">⁵ https://minorplanetcenter.net/</footnote>
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 $\Omega_{G_p}$ and $\cos i_{E_p}$ distributions is reliable. The direction of perihelion used to compute $\delta$ and $\delta^*$ was that of the original orbit, that is, it was directly obtained from the databases.

Each clone of a database yields distributions for $\cos i_{E_p}$, $\Omega_{G_p}$, $|\sin \delta|$ and $|\sin \delta^*|$. 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 $\Omega_{G_p}$ and $\cos i_{E_p}$ distributions or the original orbital energy for the $|\sin \delta|$ and $|\sin \delta^*|$ 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 > 31$ 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 $\Omega_{G_p}$ and $\cos i_{E_p}$. 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 $\delta$ and $\delta^*$, 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, $\bar{z}_i$, of the nominal orbital energy of comets in $q$ zone #i (hyperbolic comets are excluded for $q$ zone #3). Then the medium value between $\bar{z}_1$ and $\bar{z}_2$, $z_m$, and the medium value between $\bar{z}_2$ and $\bar{z}_3$, $z_M$, were computed. Thus, the $z$ zones are defined as follows: (i) $z$ zone #1 contains comets with original orbital energy between $-10^{-4}$ au$^{-1}$ 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$, $\Omega_{G_p}$, and $\cos i_{E_p}$, but it does not affect the value of $|\sin \delta|$ or $|\sin \delta^*|$. 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 dependent on the condition of observation of a comet, which is not a dynamical problem. Therefore, we preferred not to introduce 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^{-1}$ 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.
<table><caption>**Table 1.** General statistics.</caption><thead><tr><th></th><th></th><th></th><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td>Disk</td><td colspan="6">n<sub>tot</sub> = 252.0 com · Myr<sup>-1</sup> (1 464 + 4 hyp.) ¯a = 25 722 au</td></tr><tr><td>q zones (au):</td><td>0</td><td>&lt; q<sub>p</sub> &lt;</td><td>13</td><td>&lt; q<sub>p</sub> &lt;</td><td>31</td><td>&lt; q<sub>p</sub> &lt;</td><td>+∞</td></tr><tr><td>n, f</td><td colspan="3">70.9 (199) 28 %</td><td colspan="3">124.4 (626) 49%</td><td>56.7 (635) 23%</td><td>+hyp.</td></tr><tr><td>ā (au)</td><td colspan="3">17 278</td><td colspan="3">28 200</td><td>44 235</td><td></td></tr><tr><td>z zones (au):</td><td>10<sup>4</sup></td><td>&lt; -1/z &lt;</td><td>22 739</td><td>&lt; -1/z &lt;</td><td>36 217</td><td>&lt; -1/z &lt;</td><td>+∞</td></tr><tr><td>n, f</td><td colspan="3">61.4 (137) 24%</td><td colspan="3">132.3 (655) 53%</td><td>58.3 (668) 23%</td><td>+hyp.</td></tr><tr><td>Isot</td><td colspan="6">n<sub>tot</sub> = 105.1 com · Myr<sup>-1</sup> (537 + 1 hyp.) ¯a = 23 969 au</td></tr><tr><td>q zones (au):</td><td>0</td><td>&lt; q<sub>p</sub> &lt;</td><td>13</td><td>&lt; q<sub>p</sub> &lt;</td><td>31</td><td>&lt; q<sub>p</sub> &lt;</td><td>+∞</td></tr><tr><td>n, f</td><td colspan="3">41.2 (128) 39%</td><td colspan="3">49.6 (245) 47%</td><td>14.3 (163) 14%</td><td>+hyp.</td></tr><tr><td>ā (au)</td><td colspan="3">18 259</td><td colspan="3">27 454</td><td>44 292</td><td></td></tr><tr><td>z zones (au):</td><td>10<sup>4</sup></td><td>&lt; -1/z &lt;</td><td>22 856</td><td>&lt; -1/z &lt;</td><td>35 873</td><td>&lt; -1/z &lt;</td><td>+∞</td></tr><tr><td>n, f</td><td colspan="3">35.2 (82) 34%</td><td colspan="3">52.8 (260) 50%</td><td>17.2 (194) 16%</td><td>+hyp.</td></tr><tr><td>CODE</td><td colspan="6">n<sub>tot</sub> = 151 + 12 ¯a = 24 118 au</td></tr><tr><td>q zones (au):</td><td>0</td><td>&lt; q<sub>p</sub> &lt;</td><td>13</td><td>&lt; q<sub>p</sub> &lt;</td><td>31</td><td>&lt; q<sub>p</sub> &lt;</td><td>+∞</td></tr><tr><td>n, f</td><td colspan="3">80 49%</td><td colspan="3">33 20%</td><td>51 31%</td><td>+hyp.</td></tr><tr><td>ā (au)</td><td colspan="3">17 731</td><td colspan="3">29 126</td><td>59 327</td><td></td></tr><tr><td>z zones (au):</td><td>10<sup>4</sup></td><td>&lt; -1/z &lt;</td><td>22 043</td><td>&lt; -1/z &lt;</td><td>39 071</td><td>&lt; -1/z &lt;</td><td>+∞</td></tr><tr><td>n, f</td><td colspan="3">55 34%</td><td colspan="3">61 37%</td><td>48 29%</td><td>+hyp.</td></tr><tr><td>MPC</td><td colspan="6">n<sub>tot</sub> = 178 + 53 ¯a = 26 577 au</td></tr><tr><td>q zones (au):</td><td>0</td><td>&lt; q<sub>p</sub> &lt;</td><td>13</td><td>&lt; q<sub>p</sub> &lt;</td><td>31</td><td>&lt; q<sub>p</sub> &lt;</td><td>+∞</td></tr><tr><td>n, f</td><td colspan="3">71 31%</td><td colspan="3">44 19%</td><td>116 50%</td><td>+hyp.</td></tr><tr><td>ā (au)</td><td colspan="3">17 710</td><td colspan="3">28 697</td><td>54 535</td><td></td></tr><tr><td>z zones (au):</td><td>10<sup>4</sup></td><td>&lt; -1/z &lt;</td><td>21 903</td><td>&lt; -1/z &lt;</td><td>37 605</td><td>&lt; -1/z &lt;</td><td>+∞</td></tr><tr><td>n, f</td><td colspan="3">51 22%</td><td colspan="3">66 29%</td><td>124 49%</td><td>+hyp.</td></tr></tbody></table>

**Notes.** For each database, the first row gives the total number of comets on elliptic orbits or the flux per Myr for the models + the number of comets on hyperbolic orbits, followed by the semimajor axis corresponding to the mean orbital energy of all non-hyperbolic comets $ā$. In the whole table, numbers in brackets for the models correspond to the number of data; see the main text for details. For the three *q* zones, the first line gives the border of each zone, the second line gives the number and the proportion of comets in each zone. The semimajor axis corresponding to the mean orbital energy of the comets in each *q* zone are on the third line. For the three *z* zones, the first line gives the border of each zone and the second line the number and proportion of comets in each zone.

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 $i_{E_p}$ and $\Omega_{G_p}$ distributions

We first analyzed the pair of variables $cos i_{E_p}$ and $\Omega_{G_p}$. Figure 2 shows the position of the observable comets in the ($\Omega_{G_p}$, $cos i_{E_p}$) 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^\circ$ 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
<figure><img src="image_3.png" /><figcaption>**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 ($\Omega_{G_p}$, cos $i_{E_p}$) 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 $i_G = 90^\circ$.</figcaption></figure>

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 $\Omega_{G_p}$ and cos $i_{E_p}$. 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) $E_1^*$, an EE feature that corresponds to a maximum of cos $i_{E_p}$ around -0.5 in *q* zone 1, and (ii) $E_1$ and $E_2$, two KQ features that correspond respectively to a maximum of $\Omega_{G_p}$ in the interval [90°, 180°] and a maximum of cos $i_{E_p}$ 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 $E_1^*$, 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 $i_{E_p}$ distribution should located in -1. This is the case for the isotropic model. However, none of the database cos $i_{E_p}$ distribution shows such a maximum; thus, the comets in the databases in *q* zone #1 are not coming from an isotropic reservoir.

For $E_1$, in the $\Omega_{G_p}$ 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_2$ 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 semi-major 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.
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Distributions of $\Omega_{Gp}$ (left column) and $\cos i_{E_p}$ (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.</figcaption></figure>

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<fnref n="6" />. 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.

### 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 \delta|$ and $|\sin \delta^*|$ are considered. The same methodology as for the $\Omega_{Gp}$ and $\cos i_{E_p}$ 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 \delta|$ and $|\sin \delta^*|$ 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 $|\delta|$ distribution at zero for $z$ zone #2. For the EE features, they are a maximum of the $|\delta^*|$ distributions toward zero for $z$ zones #1 and #2 and a maximum of the $|\sin \delta|$ 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<fnref n="7" />, 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<fnref n="8" /> and Pan-STARRS<fnref n="9" />. 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

<footnote n="6">⁶ The conservation of the Tisserand parameter during a close encounter with a giant planet for prograde and very eccentric orbits implies that an increase in the orbital energy is related to a decrease in the ecliptic inclination.</footnote>

<footnote n="7">⁷ http://www.ll.mit.edu/impact/watch-potentially-hazardous-asteroids</footnote>

<footnote n="8">⁸ https://catalina.lpl.arizona.edu/</footnote>

<footnote n="9">⁹ https://neo.ifa.hawaii.edu/</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Distributions of $|\sin \delta|$ (left column) and $|\sin \delta^*|$ (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.</figcaption></figure>

indeed, one notes that $E_3$ is almost always detected, whereas $E_4^*$ is not. However, this does not imply the absence of this feature, just that the feature $E_3$ dominates.

From Table 2 we can see that the CODE database matches all three KQ features for both the full database and the filtered 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.

## 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^4$ au or were on their original hyperbolic orbit, and with perihelion distances

<table><caption>**Table 2.** Tests.</caption><thead><tr><th rowspan="2">#</th><th colspan="4">KQ features</th></tr><tr><th colspan="2">CODE</th><th colspan="2">MPC</th></tr></thead><tbody><tr><td>E<sub>1</sub></td><td>0.93</td><td>0.28</td><td>0.62</td><td>0.59</td></tr><tr><td>E<sub>2</sub></td><td>0.89</td><td>0.79</td><td>0.32</td><td>0.03</td></tr><tr><td>E<sub>3</sub></td><td>0.96</td><td>1.00</td><td>0.91</td><td>1.00</td></tr></tbody><thead><tr><th rowspan="2">#</th><th colspan="4">EE features</th></tr><tr><th colspan="2">CODE</th><th colspan="2">MPC</th></tr></thead><tbody><tr><td>E<sub>1</sub><sup>*</sup></td><td>0.21</td><td>0.02</td><td>0.55</td><td>0.10</td></tr><tr><td>E<sub>2</sub><sup>*</sup></td><td>0.76</td><td>0.20</td><td>0.94</td><td>0.92</td></tr><tr><td>E<sub>3</sub><sup>*</sup></td><td>0.66</td><td>0.00</td><td>0.99</td><td>0.96</td></tr><tr><td>E<sub>4</sub><sup>*</sup></td><td>0.05</td><td>0.00</td><td>0.09</td><td>0.00</td></tr><tr><td>E<sub>5</sub><sup>*</sup></td><td>0.74</td><td>0.83</td><td>0.24</td><td>0.03</td></tr></tbody></table>

**Notes.** For each feature of interest (first column), gives the proportion of cases over the 200 clones of each database where the maximum of the distribution coincides with the feature. For each feature, the proportion is first given using the entire database and the second using the supposed complete subset of comets.

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 observed 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
<figure><img src="image_6.png" /><figcaption>**Fig. 6.** Same as Figs. 3 and 5 but with the databases restricted to supposedly complete subsets of comets.</figcaption></figure>

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 10 000 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, $\rho_o$. We took $\rho_o = 0.1 M_\odot \cdot pc^{-3}$, but using a Galactic disk density of $0.15 M_\odot \cdot pc^{-3}$, 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 Observatory<fnref n="11" />, whose full survey operations are expected to begin in 2024.

*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.

<footnote n="11">¹¹ https://www.lsst.org/</footnote>
References

Brasser, R., & Morbidelli, A. 2013, Icarus, 225, 40

Brasser, R., Duncan, M. J., Levison, H. F., Schwamb, M. E., & Brown, M. E. 2012, Icarus, 217, 1

Delsemme, A. H. 1987, A&A, 187, 913

Emel'yanenko, V. V., Asher, D. J., & Bailey, M. E. 2013, Earth Moon and Planets, 110, 105

Fernández, J. A. 2002, in ESA Special Publication, 500, 303

Fernández, J. A., & Brunini, A. 2000, Icarus, 145, 580

Fouchard, M., Rickman, H., Froeschlé, C., & Valsecchi, G. B. 2017a, A&A, 604, A24

Fouchard, M., Rickman, H., Froeschlé, C., & Valsecchi, G. B. 2017b, Icarus, 292, 218

Fouchard, M., Higuchi, A., Ito, T., & Maquet, L. 2018, A&A, 620, A45

Fouchard, M., Emel'yanenko, V., & Higuchi, A. 2020, Celest. Mech. Dyn. Astron., 132, 43

Gaia Collaboration (Brown, A. G. A., et al.) 2021, A&A, 649, A1

Heisler, J., & Tremaine, S. 1986, Icarus, 65, 13

Higuchi, A. 2020, in European Planetary Science Congress, EPSC2020-265

Higuchi, A., Kokubo, E., Kinoshita, H., & Mukai, T. 2007, AJ, 134, 1693

Hills, J. G. 1981, AJ, 86, 1730

Ito, T. 2021, Earth and Space Science Open Archive, 32, presented at JpGU Meeting 2021, May 30–June 6, 2021, PPS04-P01

Ito, T., & Higuchi, A. 2020, Earth and Space Science Open Archive, 7, presented at JpGU–AGU Joint Meeting, 2020 July 12–16, 2020, PPS07-P22

Ito, T., & Ohtsuka, K. 2019, Monogr. Environ. Earth Planets, 7, 1

Jakubík, M., & Neslušan, L. 2008, Contrib. Astron. Observatory Skalnate Pleso, 38, 33

Kaib, N. A., & Quinn, T. 2008, Icarus, 197, 221

Kaib, N. A., & Quinn, T. 2009, Science, 325, 1234

Kozai, Y. 1962, AJ, 67, 591

Królikowska, M. 2006, Acta Astronomica, 56, 385

Królikowska, M. 2014, A&A, 567, A126

Królikowska, M. 2020, A&A, 633, A80

Królikowska, M., & Dybczyński, P. A. 2013, MNRAS, 435, 440

Królikowska, M., & Dybczyński, P. A. 2017, MNRAS, 472, 4634

Królikowska, M., & Dybczyński, P. A. 2018, MNRAS, 477, 2393

Królikowska, M., & Dybczyński, P. A. 2020, A&A, 640, A97

Królikowska, M., Sitarski, G., Pittich, E. M., et al. 2014, A&A, 571, A63

Leinert, C., Bowyer, S., Haikala, L. K., et al. 1998, A&AS, 127, 1

Levison, H. F., Dones, L., & Duncan, M. J. 2001, AJ, 121, 2253

Levison, H., Duncan, M., Brasser, R., & Kaufmann, D. 2010, Science, 329, 187

Lidov, M. L. 1961, Artif. Earth Satellite, 8, 5

Marsden, B. G., Sekanina, Z., & Yeomans, D. K. 1973, AJ, 78, 211

Marsden, B. G., Sekanina, Z., & Everhart, E. 1978, AJ, 83, 64

Nesvorný, D., Vokrouhlický, D., Dones, L., et al. 2017, ApJ, 845, 27

Oort, J. H. 1950, Bull. Astron. Inst. Neth., 11, 91

Rickman, H., Froeschlé, C., Froeschlé, C., & Valsecchi, G. B. 2004, A&A, 428, 673

Rickman, H., Fouchard, M., Froeschlé, C., & Valsecchi, G. B. 2008, Celest. Mech. Dyn. Astron., 102, 111

Rickman, H., Wajer, P., Przyłuski, R., et al. 2023, MNRAS, 520, 637

Šaillenfest, M., Fouchard, M., Ito, T., & Higuchi, A. 2019, A&A, 629, A95

Šteins, K., & Kronkalne, S. 1964, Acta Astron., 14, 311

Vokrouhlický, D., Nesvorný, D., & Dones, L. 2019, AJ, 157, 181

von Zeipel, H. 1910, Astron. Nachr., 183, 345

Wiśniowski, T., Rickman, H., & Wajer, P. 2022, MNRAS, 510, 1470

Yabushita, S. 1972, A&A, 16, 471