Anisotropy of Long-period Comets Explained by Their Formation Process

The tidal force from the Galactic disk is the dominant external force in the evolution of the bodies in the Oort cloud (e.g., Dones et al. 2004). The vertical component (i.e., perpendicular to the Galactic plane) of the Galactic tide plays the most important role in the formation of the Oort cloud and the production of long-period comets from it (e.g., Harrington 1985; Byl 1986; Heisler & Tremaine 1986). The Galactic potential is often approximated as axisymmetric by neglecting the radial component.

The vertical component of the tide acts on comets in the Oort cloud like a secular perturbation from a planet does on asteroids, and it drives the von Zeipel–Lidov–Kozai mechanism (von Zeipel 1910; Kozai 1962; Lidov 1962; Ito & Ohtsuka 2019) in the orbital evolution, as first shown in Heisler & Tremaine (1986). The time-averaged disturbing function that arises from the vertical component of the Galactic tide is obtained by averaging the Galactic potential over one orbital period of the comet (Heisler & Tremaine 1986). By substituting the time-averaged disturbing function into the variational equations of orbital elements (e.g., Murray & Dermott 1999), time variations of the secular orbital elements are obtained (e.g., Matese & Whitman 1989, 1992; Brasser 2001; Breiter & Ratajczak 2005; Brasser et al. 2006; Higuchi et al. 2007). Higuchi et al. (2007) applied the solutions to examine the formation of the Oort cloud from planetesimals with large semimajor axes initially on the ecliptic plane.

The analytical solutions to the orbital elements presented by the above authors are useful for understanding the evolution of the Oort cloud and the overall behavior of the distribution of comets generated by the Galactic tide. However, the solutions are not so useful for the discussion about observed long-period comets returning to the planetary region for the following two reasons. First, the time variation of the longitude of the ascending node in the Galactic coordinates, dΩ/dt, becomes large as the eccentricity e approaches 1 (Higuchi et al. 2007). This means that Ω of a long-period comet and the inclination with respect to the ecliptic plane i_E, which is a function of Ω and the inclination with respect to the Galactic plane i, are drastically changing with the perihelion distance q when it is in the observable region (i.e., e ≃ 1). Consequently, there is no firm relation between the initial orbital elements and the observed orbital elements in the planetary region. Second, the angular momentum of a comet with e ≃ 1 is quite small and it is easily changed by perturbations from passing stars and/or the radial component of the Galactic tide, both of which are neglected in the derivation of the analytical solutions (e.g., Matese et al. 1999; Higuchi et al. 2007). The conservation of the vertical component of the angular momentum, which is defined as $j=\sqrt{1-{e}^{2}}\cos i$, is crucial in order to derive the analytical solution to the inclination at small q. Therefore, the accurate prediction of i at e ≃ 1 is difficult. For the above two reasons, the analytical solutions to the orbital elements, especially to Ω, i, and i_E, are not so useful for describing the orbits of long-period comets.

Besides Ω, i, and the argument of perihelion, ω, in the Galactic coordinates, the Galactic longitude and latitude of the direction of aphelion, L and B (or those of perihelion, l and b), are also used to evaluate the distribution of observed long-period comets. Many authors have pointed out anomalies in the distributions of L and B. For example, Luest (1984) and Delsemme (1987) found depletions around b = 0 and b = ±90°. Delsemme (1987) explained that the depletions are the result of the strength of the Galactic tide, which is minimum for b = 0 and b = ±90°. Matese & Whitmire (1996) evaluated the effect of the radial component of the Galactic tide in the distributions of l and b of long-period comets. Biermann et al. (1983), Luest (1984), and Delsemme (1986) investigated aphelion clustering on the L–B plane and Matese et al. (1999) identified an anomalously overpopulated "great circle" as two peaks centered on L = 135° and 315°. These concentrations of the aphelia were explained by introducing a hypothetical perturber that encountered the solar system.

The above investigations of the distribution of aphelia are made on the assumption that the Oort cloud, which stores the long-period comets, has an isotropic distribution of comets. However, based on the standard formation scenario, the Oort cloud comets are planetesimals formed in the protoplanetary disk and initially on the ecliptic plane with perihelion distances near the giant planets (e.g., Dones et al. 2004). The role of stars in the evolution of the Oort cloud has been examined by many authors (e.g., Dybczyński 2002; Fouchard et al. 2011). They showed that passing stars act like random noise on the distribution of comets in the Oort cloud. As long as the Oort cloud is not completely destroyed by close stellar encounters, the memory of the initial distribution can be found as anisotropies in the present distribution, which can be explained without assuming any hypothetical perturber.

In this paper, we investigate the evolution of the aphelia of comets initially on the ecliptic plane under the axisymmetric approximation of the Galactic tide with the same procedure as in Higuchi et al. (2007). Using the analytical solutions, we predict the distribution of long-period comets on the L–B plane. The solutions to L, B, and other orbital elements are derived in Section 2. In Section 3, the analytic solutions are evaluated by comparisons with numerical integrations of the equation of motion that take into account not only the vertical but also the radial component of the Galactic tide. In Section 4, we approximate the analytical solutions for the special case of Oort cloud comets and propose the concentration of comets on the ecliptic plane and the "empty ecliptic" plane, which is defined as a plane formed by a rotation of the ecliptic around the Galactic pole by 180°. In Section 5, the distribution of observed small bodies is briefly examined to find the concentrations on the ecliptic and the empty ecliptic in the L–sin B plane. Section 6 is devoted to a summary and discussion.

In this section, we derive the Galactic longitude and latitude of the direction of the aphelion, L and B, respectively, and their time variations and solutions. Time variations and solutions for the eccentricity e and the longitude of the ascending node Ω are also shown but we use slightly different expressions from Higuchi et al. (2007) for the purpose of this paper. The orbital elements are given in Galactic coordinates except for the ecliptic inclination i_E.

2.1. The Galactic Longitude L and Latitude B

Using the orbital elements, the unit vector of the direction of aphelion in the Galactic coordinates r_Q is written as

$$\begin{eqnarray}&&\begin{array}{c}{{\boldsymbol{r}}}_{Q}=\left(\begin{array}{c}{Q}_{x}\\ {Q}_{y}\\ {Q}_{z}\end{array}\right)\,=\,\left(\begin{array}{l}-\cos \omega \cos {\rm{\Omega }}+\sin \omega \sin {\rm{\Omega }}\cos i\\ -\cos \omega \sin {\rm{\Omega }}-\sin \omega \cos {\rm{\Omega }}\cos i\\ -\sin \omega \sin i\end{array}\right).\end{array}\end{eqnarray} \tag{ 1 }$$

Then L and sin B are written as

$$\begin{eqnarray}&&L={\rm{\arctan }}\left(\displaystyle \frac{{Q}_{y}}{{Q}_{x}}\right)={\rm{\Omega }}+\theta ,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\sin B={Q}_{z}=-\sin \omega \sin i,\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray}&&\begin{array}{c}\theta =\left\{\begin{array}{c}{\rm{a}}{\rm{r}}{\rm{c}}{\rm{t}}{\rm{a}}{\rm{n}}\left(\displaystyle \frac{\sin \omega \cos i}{\cos \omega }\right)\,{\rm{f}}{\rm{o}}{\rm{r}}\,\cos \omega \lt 0\\ \pi +{\rm{a}}{\rm{r}}{\rm{c}}{\rm{t}}{\rm{a}}{\rm{n}}\left(\displaystyle \frac{\sin \omega \cos i}{\cos \omega }\right)\,{\rm{f}}{\rm{o}}{\rm{r}}\,\cos \omega \gt 0.\end{array}\right.\end{array}\end{eqnarray} \tag{ 4 }$$

2.2. Conserved Quantities

Assume that the Galactic tide is much smaller than the solar gravity. The time-averaged Hamiltonian of a body moving under the approximated Galactic potential is given as

$$\begin{eqnarray}&&\langle H\rangle =-\displaystyle \frac{{{GM}}_{\odot }}{2a}-R,\end{eqnarray} \tag{ 5 }$$

where G is the gravitational constant, M_⊙ is the solar mass, a is the semimajor axis of the body, and R is the disturbing function

$$\begin{eqnarray}&&R=-\displaystyle \frac{{\nu }_{0}^{2}}{4}{a}^{2}{\sin }^{2}i\left(1-{e}^{2}+5{e}^{2}{\sin }^{2}\omega \right),\end{eqnarray} \tag{ 6 }$$

where ${\nu }_{0}=\sqrt{4\pi G\rho }$ is the vertical frequency and ρ is the total density in the solar neighborhood (e.g., Heisler & Tremaine 1986). From Equation (6) and Lagrange's planetary equation for da/dt, we know a is constant. Then we introduce a new simplified Hamiltonian:

$$\begin{eqnarray}&&c={\sin }^{2}i\left(1-{e}^{2}+5{e}^{2}{\sin }^{2}\omega \right).\end{eqnarray} \tag{ 7 }$$

The simplified z-component of the angular momentum, which is a conserved quantity under the axisymmetric approximation of the potential, is written as

$$\begin{eqnarray}&&j=\sqrt{1-{e}^{2}}\cos i.\end{eqnarray} \tag{ 8 }$$

Substituting Equation (8) into Equation (7), one can draw equi-Hamiltonian curves on the ω−e plane for given c and j using Equation (7). From the Hamiltonian curves, we can learn the overall behavior of ω and e without solving the equation of motion. For some cases, equi-Hamiltonian curves circulate with ω, and for other cases they librate around ω = 90° or 270°. This libration is essentially the von Zeipel–Lidov–Kozai mechanism (von Zeipel 1910; Kozai 1962; Lidov 1962; Ito & Ohtsuka 2019). The condition for circulation is to have a solution to e for ω = 0, i.e.,

$$\begin{eqnarray}&&\begin{array}{c}\begin{array}{l}c+{j}^{2}\lt 1\,({\rm{circulation}})\\ c+{j}^{2}\gt 1\,({\rm{libration}})\end{array}\end{array}\end{eqnarray} \tag{ 9 }$$

and c + j² = 1 gives the separatrix (e.g., Higuchi et al. 2007). This leads to the necessary condition on i for libration, $\sin i\gt \sqrt{1/5}$.

Substituting Equations (3) and (8) into Equation (7), the Hamiltonian is given with B instead of i and ω,

$$\begin{eqnarray}&&c=1-{e}^{2}-{j}^{2}+5{e}^{2}{\sin }^{2}B.\end{eqnarray} \tag{ 10 }$$

The separatrix with Equation (10) is written as

$$\begin{eqnarray}&&c+{j}^{2}=1-{e}^{2}\left(1-5{\sin }^{2}B\right)=1.\end{eqnarray} \tag{ 11 }$$

For e² > 0, the sufficient condition on B for libration is given as

$$\begin{eqnarray}&&\sin | B| \gt \sqrt{\displaystyle 1/5}.\end{eqnarray} \tag{ 12 }$$

Matese & Whitman (1989) defined the value $B\,=\arcsin \sqrt{1/5}\simeq \pm 26\buildrel{\circ}\over{.} 6$ as a barrier that the latitude of perihelion cannot migrate across.

2.3. Time Variations

Substituting Equation (6) into Lagrange's planetary equations (Murray & Dermott 1999), we obtain the time variations of e, i, Ω, and ϖ as

$$\begin{eqnarray}&&\displaystyle \frac{{de}}{{dt}}=\displaystyle \frac{5{\nu }_{0}^{2}}{2n}e\sqrt{1-{e}^{2}}{\sin }^{2}i\sin \omega \cos \omega ,\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{di}}{{dt}}=-\displaystyle \frac{5{\nu }_{0}^{2}}{2n}\displaystyle \frac{{e}^{2}}{\sqrt{1-{e}^{2}}}\sin i\cos i\sin \omega \cos \omega ,\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{\rm{\Omega }}}{{dt}}=-\displaystyle \frac{{\nu }_{0}^{2}}{2n}\displaystyle \frac{\cos i}{\sqrt{1-{e}^{2}}}(1-{e}^{2}+5{e}^{2}{\sin }^{2}\omega ),\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d\varpi }{{dt}} & = & \displaystyle \frac{{\nu }_{0}^{2}}{2n}\displaystyle \frac{1}{\sqrt{1-{e}^{2}}}\left(1-{e}^{2}+5{e}^{2}{\sin }^{2}\omega \right.\\ & & -\ \left.5{\sin }^{2}\omega {\sin }^{2}i\right)+\displaystyle \frac{d{\rm{\Omega }}}{{dt}},\end{array}\end{eqnarray} \tag{ 16 }$$

where n = a^−3/2 is the mean motion and ϖ = ω + Ω is the longitude of the pericenter. Note that the short-period terms arising from the variation of the mean longitude are dropped. From the definition ϖ = ω + Ω,

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d\omega }{{dt}} & = & \displaystyle \frac{d\varpi }{{dt}}-\displaystyle \frac{d{\rm{\Omega }}}{{dt}}\\ & = & \ \displaystyle \frac{{\nu }_{0}^{2}}{2n}\displaystyle \frac{1}{\sqrt{1-{e}^{2}}}\left(1-{e}^{2}+5{e}^{2}{\sin }^{2}\omega \right.\\ & & -\ \left.5{\sin }^{2}\omega {\sin }^{2}i\right).\end{array}\end{eqnarray} \tag{ 17 }$$

From Equation (2), we have

$$\begin{eqnarray}&&\displaystyle \frac{{dL}}{{dt}}=\displaystyle \frac{d{\rm{\Omega }}}{{dt}}+\displaystyle \frac{d\theta }{{dt}}.\end{eqnarray} \tag{ 18 }$$

From Equation (4), we have

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d\theta }{{dt}} & = & {\left(1+{\tan }^{2}\theta \right)}^{-1}\left({\cos }^{-2}\omega \displaystyle \frac{d\omega }{{dt}}\cos i\right.\\ & & -\ \left.\tan \omega \sin i\displaystyle \frac{{di}}{{dt}}\right).\end{array}\end{eqnarray} \tag{ 19 }$$

Substituting Equations (4), (3), (14), and (17) into Equation (19), we have

$$\begin{eqnarray}&&\displaystyle \frac{d\theta }{{dt}}=-\displaystyle \frac{d{\rm{\Omega }}}{{dt}}-\displaystyle \frac{2{\nu }_{0}^{2}}{n}\sqrt{1-{e}^{2}}\cos i{\tan }^{2}B.\end{eqnarray} \tag{ 20 }$$

Substituting Equation (20) into Equation (18), we have

$$\begin{eqnarray}&&\displaystyle \frac{{dL}}{{dt}}=-\displaystyle \frac{2{\nu }_{0}^{2}}{n}\sqrt{1-{e}^{2}}\cos i{\tan }^{2}B.\end{eqnarray} \tag{ 21 }$$

2.4. Solutions

2.4.1. Solutions to e, i, ω, and B

Introducing χ = 1 − e², Equation (13) can be rewritten as

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{d\chi }{{dt}}\,=\,\displaystyle \frac{d\chi }{{de}}\displaystyle \frac{{de}}{{dt}}\\ \,=\,-\displaystyle \frac{5{\nu }_{0}^{2}}{n}{\chi }^{-1/2}(1-\chi )(\chi -{j}^{2})\sin \omega \cos \omega .\end{array}\end{eqnarray} \tag{ 22 }$$

Using Equations (7) and (8), ω can be rewritten as

$$\begin{eqnarray}&&{\sin }^{2}\omega =\displaystyle \frac{\chi (c+{j}^{2}-\chi )}{5(1-\chi )\left(\chi -{j}^{2}\right)}.\end{eqnarray} \tag{ 23 }$$

Substituting $\sin \omega$ and $\cos \omega$ from Equation (23) into Equation (22),

$$\begin{eqnarray}&&\displaystyle \frac{d\chi }{{dt}}=-{A}_{1}\sqrt{({\chi }_{0}^{* }-\chi )({\chi }_{2}^{* }-\chi )(\chi -{\chi }_{1}^{* })},\end{eqnarray} \tag{ 24 }$$

where

$$\begin{eqnarray}&&{A}_{1}=\displaystyle \frac{2{\nu }_{0}^{2}}{n},\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{\chi }_{0}^{* }=c+{j}^{2},\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}\begin{array}{rcl}Q(\chi ) & = & -4{\chi }^{2}+\left(5+4{j}^{2}-c\right)\chi -5{j}^{2}\\ & \equiv & \ 4(\chi -{\chi }_{1}^{* })({\chi }_{2}^{* }-\chi )\,({\chi }_{1}^{* }\lt {\chi }_{2}^{* }).\end{array}\end{eqnarray} \tag{ 27 }$$

The solution to Equation (24) is expressed using a Jacobian elliptic function, sn (Byrd & Friedman 1971),

$$\begin{eqnarray}&&\chi =({\alpha }_{1}-{\alpha }_{0}){\mathrm{sn}}^{2}u+{\alpha }_{0},\end{eqnarray} \tag{ 28 }$$

where we define ${\alpha }_{0}=\min \{{\chi }_{0}^{* }$, ${\chi }_{1}^{* }$, ${\chi }_{2}^{* }\}$, ${\alpha }_{1}=\mathrm{med}\{{\chi }_{0}^{* }$, ${\chi }_{1}^{* }$, ${\chi }_{2}^{* }\}$, and ${\alpha }_{2}=\max \{{\chi }_{0}^{* }$, ${\chi }_{1}^{* }$, ${\chi }_{2}^{* }\}$,

$$\begin{eqnarray}&&u=\displaystyle \frac{{A}_{1}}{g}t+{u}_{0},\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&g=\displaystyle \frac{2}{\sqrt{{\alpha }_{2}-{\alpha }_{0}}},\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{u}_{0}=\pm F({\phi }_{0},k),\end{eqnarray} \tag{ 31 }$$

where F is a normal elliptic integral of the first kind with modulus k and amplitude ϕ₀,

$$\begin{eqnarray}&&{k}^{2}=\displaystyle \frac{{\alpha }_{1}-{\alpha }_{0}}{{\alpha }_{2}-{\alpha }_{0}}\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&\sin {\phi }_{0}=\sqrt{\displaystyle \frac{{\chi }_{0}^{* }-{\alpha }_{0}}{{\alpha }_{1}-{\alpha }_{0}}},\end{eqnarray} \tag{ 33 }$$

The sign of F(ϕ₀, k) in Equation (31) depends on ω_i, the value of ω at t = 0; it is negative for $\sin (2{\omega }_{i})\gt 0$ and positive for $\sin (2{\omega }_{i})\lt 0$.

From Equation (28), we learn that χ oscillates between the maximum (α₁) and minimum (α₀) values according to the parameter u, with the period u = 2K, where K = K(k) is a complete elliptic integral of the first kind with modulus k. The period given in time is

$$\begin{eqnarray}&&{P}_{\chi }=2K\displaystyle \frac{g}{{A}_{1}}.\end{eqnarray} \tag{ 34 }$$

Figure 1 shows periods of evolution against B_i, the value of B at t = 0. The light blue curve in Figure 1 is P_χ for a = 2 × 10⁴ au, q_i = 10 au, and i_i = 60° obtained from Equation (34), where q_i and i_i are the values of q and i at t = 0, respectively. It varies with B_i and becomes a maximum at the separatrix. However, the dependence on B_i except around the separatrix is not as strong as that on a, which is proportional to a^−3/2. The typical value of P_χ for bodies whose eccentricity at t = 0 is e_i ≃ 1 is derived in Section 4 as a function of a (Equation (89)). The period of the oscillation of i and that of ω in the case of libration are the same as P_χ. For ω in the case of circulation, the period of circulation is 2P_χ. The period of oscillation of B is the same as that for ω, i.e., P_χ or 2P_χ.

Zoom In Zoom Out Reset image size

Using Equation (28), one can calculate e at time t from e_i, i_i, and ω_i. Once e is calculated, i is obtained from Equation (8) and then sin²ω and $\sin B$ are obtained from Equations (23) and (3), respectively. To calculate the real ω from ${\sin }^{2}\omega$, one needs to know if ω circulates or librates from Equation (9) and at which phase of the evolution the time t is located using the period given by Equation (34). One can also use Equation (10) instead of (7) to obtain $| B|$.

2.4.2. Solution to Ω

Substituting Equations (7) and (8) into Equation (15) and replacing 1 − e² with χ, we have

$$\begin{eqnarray}&&\displaystyle \frac{d{\rm{\Omega }}}{{dt}}=-\displaystyle \frac{{A}_{2}}{1-{w}_{1}^{2}{\mathrm{sn}}^{2}(u,k)},\end{eqnarray} \tag{ 35 }$$

where

$$\begin{eqnarray}&&{A}_{2}=\displaystyle \frac{{A}_{1}}{4}j\displaystyle \frac{{\chi }_{0}^{* }-{j}^{2}}{{\alpha }_{0}-{j}^{2}},\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&{w}_{1}^{2}=\displaystyle \frac{{\alpha }_{1}-{\alpha }_{0}}{{j}^{2}-{\alpha }_{0}}.\end{eqnarray} \tag{ 37 }$$

The integration of (35) with t is rewritten as

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Omega }} & = & {{\rm{\Omega }}}_{i}-{A}_{2}{\displaystyle \int }_{0}^{t^{\prime} }\displaystyle \frac{{dt}}{1-{w}_{1}^{2}{\mathrm{sn}}^{2}(u,k)}\\ & = & \ {{\rm{\Omega }}}_{i}-{A}_{3}{\displaystyle \int }_{{u}_{0}}^{u^{\prime} }\displaystyle \frac{{du}}{1-{w}_{1}^{2}{\mathrm{sn}}^{2}(u,k)},\end{array}\end{eqnarray} \tag{ 38 }$$

where Ω_i is the value of Ω at t = 0 and

$$\begin{eqnarray}&&u=\displaystyle \frac{{A}_{1}}{g}t+{u}_{0},\quad {dt}=\displaystyle \frac{g}{{A}_{1}}{du},\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}&&{A}_{3}=g\displaystyle \frac{{A}_{2}}{{A}_{1}}=\displaystyle \frac{1}{2}\displaystyle \frac{j({j}^{2}-{\chi }_{0}^{* })}{{j}^{2}-{\alpha }_{0}}\displaystyle \frac{1}{\sqrt{{\alpha }_{2}-{\alpha }_{0}}}.\end{eqnarray} \tag{ 40 }$$

Since sn²(u, k) oscillates with a period of u = 2K, we split u' as

$$\begin{eqnarray}&&u^{\prime} ={u}_{r}+2{mK},\end{eqnarray} \tag{ 41 }$$

where m is an integer that gives 0 ≤ u_r < 2K. Then, Equation (38) is written and integrated using an elliptic integral of the third kind Π (Byrd & Friedman 1971) as

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Omega }} & = & {{\rm{\Omega }}}_{i}-{A}_{3}\left\{{\displaystyle \int }_{0}^{{u}_{r}}\displaystyle \frac{{du}}{1-{w}_{1}^{2}{\mathrm{sn}}^{2}(u,k)}+2m\right.\\ & & \times \ \left.{\displaystyle \int }_{0}^{K}\displaystyle \frac{{du}}{1-{w}_{1}^{2}{\mathrm{sn}}^{2}(u,k)}-{\displaystyle \int }_{0}^{{u}_{0}}\displaystyle \frac{{du}}{1-{w}_{1}^{2}{\mathrm{sn}}^{2}(u,k)}\right\}\\ & = & \ {{\rm{\Omega }}}_{i}-{A}_{3}\left[{{\rm{\Pi }}}_{{\rm{\Omega }}}^{{\prime} }+2m{\rm{\Pi }}(K,{w}_{1}^{2},k)-{\rm{\Pi }}({u}_{0},{w}_{1}^{2},k)\right],\end{array}\end{eqnarray} \tag{ 42 }$$

where ${\rm{\Pi }}(K,{w}_{1}^{2},k)$ is a complete elliptic integral of the third kind and

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Pi }}}_{{\rm{\Omega }}}^{{\prime} }\\ =\ \left\{\begin{array}{lll}{\rm{\Pi }}({u}_{r},{w}_{1}^{2},k) & \mathrm{for} & {u}_{r}\lt K\\ 2{\rm{\Pi }}(K,{w}_{1}^{2},k)-{\rm{\Pi }}(2K-{u}_{r},{w}_{1}^{2},k) & \mathrm{for} & {u}_{r}\gt K.\end{array}\right.\end{array}\end{eqnarray} \tag{ 43 }$$

The period of Ω can be estimated by linear approximation as

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Omega }}}^{* } & = & {{\rm{\Omega }}}_{i}-\displaystyle \frac{{A}_{3}{\rm{\Pi }}[K,{w}_{1}^{2},k]}{K}(u-{u}_{0})\\ & = & \ {{\rm{\Omega }}}_{i}-{n}_{{{\rm{\Omega }}}^{* }}t,\end{array}\end{eqnarray} \tag{ 44 }$$

where

$$\begin{eqnarray}&&{n}_{{{\rm{\Omega }}}^{* }}=\displaystyle \frac{{A}_{2}}{K}{\rm{\Pi }}[K,{w}_{1}^{2},k].\end{eqnarray} \tag{ 45 }$$

The period of Ω* is obtained as

$$\begin{eqnarray}&&{P}_{{{\rm{\Omega }}}^{* }}=\displaystyle \frac{2\pi K}{{A}_{2}{\rm{\Pi }}[K,{w}_{1}^{2},k]}.\end{eqnarray} \tag{ 46 }$$

Figure 1 shows ${P}_{{{\rm{\Omega }}}^{* }}$ in orange against B_i for a = 2 × 10⁴ au, q_i = 10 au, and i_i = 60°. The behavior of ${P}_{{{\rm{\Omega }}}^{* }}$ is quite similar to that of P_χ. The approximate relation between ${P}_{{{\rm{\Omega }}}^{* }}$ and P_χ is shown in Section 4.

2.4.3. Solution to L

From Equations (3), (8), and (23), ${\sin }^{2}B$ is expressed with c, j, and χ as

$$\begin{eqnarray}&&{\sin }^{2}B\,=\,{\sin }^{2}\omega {\sin }^{2}i=\displaystyle \frac{c+{j}^{2}-\chi }{5(1-\chi )}.\end{eqnarray} \tag{ 47 }$$

Then, ${\tan }^{2}B$ is written as

$$\begin{eqnarray}&&{\tan }^{2}B=\displaystyle \frac{{\sin }^{2}B}{1-{\sin }^{2}B}=\displaystyle \frac{1}{4}\left[1+\displaystyle \frac{5(c+{j}^{2}-1)}{5-c-{j}^{2}-4\chi }\right].\end{eqnarray} \tag{ 48 }$$

Substituting Equation (28) into Equation (48), ${\tan }^{2}B$ is expressed as

$$\begin{eqnarray}&&{\tan }^{2}B=\displaystyle \frac{1}{4}\left[1+\displaystyle \frac{{A}_{4}}{1-{w}_{2}^{2}{\mathrm{sn}}^{2}(u,k)}\right],\end{eqnarray} \tag{ 49 }$$

where

$$\begin{eqnarray}&&{A}_{4}=\displaystyle \frac{5(c+{j}^{2}-1)}{5-c-{j}^{2}-4{\alpha }_{0}},\end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray}&&{w}_{2}^{2}=\displaystyle \frac{4({\alpha }_{1}-{\alpha }_{0})}{5-c-{j}^{2}-4{\alpha }_{0}}.\end{eqnarray} \tag{ 51 }$$

Substituting Equation (48) into Equation (18), we have

$$\begin{eqnarray}&&\displaystyle \frac{{dL}}{{dt}}=-\displaystyle \frac{{A}_{1}}{4}j-\displaystyle \frac{{A}_{5}}{1-{w}_{2}^{2}{\mathrm{sn}}^{2}(u,k)},\end{eqnarray} \tag{ 52 }$$

where

$$\begin{eqnarray}&&{A}_{5}=\displaystyle \frac{1}{4}{A}_{1}{A}_{4}j={A}_{1}\displaystyle \frac{5j(c+{j}^{2}-1)}{4(5-c-{j}^{2}-4{\alpha }_{0})}.\end{eqnarray} \tag{ 53 }$$

As well as Ω, the integration of (52) with t is expressed using an elliptic integral of the third kind as

$$\begin{eqnarray}&&L={L}_{i}-{A}_{7}u-{A}_{6}\left[{{\rm{\Pi }}}_{L}^{{\prime} }+2m{\rm{\Pi }}(K,{w}_{2}^{2},k)-{\rm{\Pi }}({u}_{0},{w}_{2}^{2},k)\right],\end{eqnarray} \tag{ 54 }$$

where L_i is the value of L at $t=0$ and

$$\begin{eqnarray}&&{A}_{6}=g\displaystyle \frac{{A}_{5}}{{A}_{1}}=\displaystyle \frac{5j(c+{j}^{2}-1)}{2\sqrt{{\alpha }_{2}-{\alpha }_{0}}(5-c-{j}^{2}-4{\alpha }_{0})},\end{eqnarray} \tag{ 55 }$$

$$\begin{eqnarray}&&{A}_{7}=\displaystyle \frac{j}{4}g=\displaystyle \frac{j}{2\sqrt{{\alpha }_{2}-{\alpha }_{0}}},\end{eqnarray} \tag{ 56 }$$

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Pi }}}_{L}^{{\prime} }\\ \ =\ \left\{\begin{array}{l}{\rm{\Pi }}({u}_{r},{w}_{2}^{2},k)\ \mathrm{for}\ {u}_{r}\lt K\\ 2{\rm{\Pi }}(K,{w}_{2}^{2},k)-{\rm{\Pi }}(2K-{u}_{r},{w}_{2}^{2},k)\ \mathrm{for}\ {u}_{r}\gt K.\end{array}\right.\end{array}\end{eqnarray} \tag{ 57 }$$

The period of L is also estimated in the same manner as Ω using the linear approximation,

$$\begin{eqnarray}\begin{array}{rcl}{L}^{* } & = & {L}_{i}-\displaystyle \frac{{A}_{1}}{4}{jt}-\displaystyle \frac{{A}_{6}{\rm{\Pi }}[K,{w}_{2}^{2},k]}{K}(u-{u}_{0}),\\ & = & \ {L}_{i}-{n}_{{L}^{* }}t,\end{array}\end{eqnarray} \tag{ 58 }$$

where

$$\begin{eqnarray}&&{n}_{{L}^{* }}=\displaystyle \frac{{A}_{1}}{4}j\left(1+\displaystyle \frac{{A}_{4}}{K}{\rm{\Pi }}[K,{w}_{2}^{2},k]\right).\end{eqnarray} \tag{ 59 }$$

The period of ${L}^{* }$ is obtained as

$$\begin{eqnarray}&&{P}_{{L}^{* }}=\displaystyle \frac{8\pi }{{A}_{1}j}{\left(1+\displaystyle \frac{{A}_{4}}{K}{\rm{\Pi }}[K,{w}_{1}^{2},k]\right)}^{-1}\end{eqnarray} \tag{ 60 }$$

The black dashed curve in Figure 1 shows ${P}_{{L}^{* }}$ plotted against B_i for a = 2 × 10⁴ au, q_i = 10 au, and i_i = 60°. The behavior of ${P}_{{L}^{* }}$ looks identical to that of ${P}_{{{\rm{\Omega }}}^{* }}$ for B_i > 266. In contrast, for B_i < 266, ${P}_{{L}^{* }}$ suddenly becomes 10^2–5 times larger than that for B_i > 266. For this example, ${P}_{{L}^{* }}$ for c < 1 is much longer than the age of the solar system.

We test the analytical solutions to the orbital elements, L, and B by comparing with the orbital evolution obtained by numerical integrations. We are especially interested in checking two approximations that we have made in the derivation of the analytical solutions: the axisymmetric approximation of the Galactic tide by neglecting the radial component and the time-averaging of the Hamiltonian assuming that the Galactic tide is much smaller than the solar gravity. Brasser (2001) examined both approximations using numerical orbital integrations for comets mainly with e ≪ 1. Higuchi et al. (2007) considered comets with e_i ∼ 1 and showed that the time-averaging is plausible for comets with orbital periods T_K ≲ 10P_χ; however, they neglected the radial component of the Galactic tide in their numerical integrations. In this paper, we focus on comets with e ≃ 1 and examine how the analytical solutions are useful for the discussion about long-period comets in the observable region.

3.1. Equation of Motion

Under the epicyclic approximation (e.g., Binney & Tremaine 1987), the equation of motion of a body orbiting around the Sun with tidal forces from the Galactic disk is

$$\begin{eqnarray}&&\displaystyle \frac{d{\boldsymbol{r}}}{{dt}}=-{{GM}}_{\odot }\displaystyle \frac{{\boldsymbol{r}}}{{r}^{3}}+{f}_{\mathrm{tide}},\end{eqnarray} \tag{ 61 }$$

where r is the position of the body with respect to the Sun and f_tide is the Galactic tide,

$$\begin{eqnarray}&&{f}_{\mathrm{tide}}={{\rm{\Omega }}}_{0}^{2}(x^{\prime} -y^{\prime} )-{\nu }_{0}^{2}z^{\prime} ,\end{eqnarray} \tag{ 62 }$$

where x', y', and z' give the position of the body in rotating coordinates centered on the Sun, Ω₀ is the circular frequency (i.e., the angular speed of the Sun in the Galaxy). We adopt Ω₀ = 26 km s⁻¹ (e.g., Binney & Tremaine 1987) and ρ = 0.1 M_⊙ pc⁻³ (Holmberg & Flynn 2000), which give ${\nu }_{0}^{2}/{{\rm{\Omega }}}_{0}^{2}\simeq 10$. To evaluate the analytical solutions derived in Section 2, we integrated Equation (61) for bodies initially on the ecliptic plane (i.e., i_i = 60°, Ω_i = 186°) for 4.5 Gyr with the P(EC)² Hermite scheme (Kokubo et al. 1998) and compared the orbital evolutions with the analytical solutions. The bodies are set at their perihelion at t = 0. We also performed extra numerical calculations that neglect the first term in Equation (62) to examine the effect of the radial component of the Galactic tide.

3.2. Comparison

In this section, we compare the results of numerical calculations that consider both the radial and vertical components of the Galactic tide and the analytical solutions by plotting them together in the same figures.

Figures 2 and 3 show the orbital evolutions of bodies against time for 4.5 Gyr. All bodies have (i_i, Ω_i) = (60°, 186°) but different values (color-coded) of q_i and ω_i. Bodies in the same colors on the left and right panels have the same initial orbital elements except for a_i, the semimajor axis at t = 0; a_i = 2 × 10⁴ au and 5 × 10⁴ for the left and right panels, respectively. Circles and squares are the results obtained by numerical integration of Equation (61) and the solid curves are the analytical solutions derived in Section 2. Panel (1) in Figure 2 shows the evolution of a. There are variations within ≲0.5%, but no systematic decay or increase is seen for 4.5 Gyr evolution. We find the same features in results of numerical calculations without the radial component of the Galactic tide. We conclude that the variation of a is due to the short-term effect of the Galactic tide that is neglected in the time-averaging process in the derivation of Equation (5).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Panel (2) in Figure 2 shows the evolution of e. Oscillations with various periods and amplitudes occur since the Hamiltonian of a body depends on ω_i as seen in Equation (7). The initial values are very close to 1 because e_i = 1 − q_i/a_i = 0.9985–0.99975, but can be larger if ${\sin }^{2}{\omega }_{i}\ne 1$. Panel (3) in Figure 2 shows the evolution of q. The range of time when the bodies arrive in the observable region, i.e., q ≲ 10² au, is very short compared to the period of the oscillation. For most of the time, q is outside the planetary region and planetary perturbations are negligible for those orbits. The analytical solutions and the results of numerical calculations shown in panels (2) and (3) are in good agreement. In contrast, those for the evolution of i are not in good agreement as shown in panel (4) in Figure 2. For the analytical solutions, the orbits never become retrograde as a consequence of the conservation of j. Panels (5)–(8) in Figure 2 are the same as panels (1)–(4) but for a_i = 5 × 10⁴ au. Note that only for t > 3 Gyr are results shown, except in panel (5). They have the same features as for a_i = 2 × 10⁴ au, although the agreement of the analytical solutions with the results of the numerical calculations is worse because of the shift of the oscillation phase. The shift can simply be explained by the evolution of a, which is not completely equal to a_i as shown in panel (1) in Figure 2. This makes the period of oscillation slightly longer/shorter than the period given by Equation (34). Since the variation of a is larger for large a_i and the dynamical evolution is quicker for large a_i, the shift of periods is not negligible in 4.5 Gyr for a_i = 5 × 10⁴ au. In contrast, the amplitudes of the oscillations of e and q found in the numerical calculations show rather good agreement with the analytical solutions.

The top panel in Figure 4 shows orbital evolution on the i–q plane with the same symbols and colors as in panels (1)–(4) in Figure 2. All the bodies initially have the value of j that is shown as an equi-j curve (black solid) given by Equation (8), where e = 1 − q/a is substituted. The results of numerical calculations shown with circles are scattered away from the equi-j curve for small q. Consequently, j is not conserved completely. However, for very small q (i.e., e ∼ 1), j is very small independently of i. In other words, j is approximately conserved for 4.5 Gyr. Consequently, i always reaches ≃90° when q is large. The bottom panel shows the same as the top one but for a_i = 5 × 10⁴ au. The conservation of j is worse than for a_i = 2 × 10⁴ au but still we can say it is approximately conserved even when the radial component of the Galactic tide is included.

Zoom In Zoom Out Reset image size

Panel (1) in Figure 3 shows the evolution of Ω. The evolution is characterized as a decreasing step function. The values at each step are different among the bodies although they all have the same value of Ω_i = 186°. As is clear from Equation (15), the big drop in Ω occurs when e is large, i.e., q is small. Panel (2) in Figure 3 shows the evolution of ω. Two of the six bodies, those with ω_i = 10° (black) and 210° (blue), are in the case of circulation. Their behavior of having a big change at e ≃ 1 is quite similar to that of Ω but they increase with time. The other four bodies are in the case of libration and they librate around ω = 90°/270° depending on each ω_i. Panel (3) in Figure 3 shows the evolution of L, which is expressed as a function of i, Ω, and ω (Equation (2)). Interestingly, for bodies in the case of circulation, L is almost constant beyond 4.5 Gyr. For bodies in the case of libration, their evolutions are quite similar to those of Ω. However, the phases are shifted so that the value of L is almost constant when q ≲ 10² au. Panel (4) in Figure 3 shows the evolution of $\sin B$, which is expressed as a function of i and ω (Equation (3)). For bodies in the case of circulation, $\sin \ B$ oscillates symmetrically with respect to $\sin B=0$, the Galactic plane. As given by Equation (12), the cases of circulation and libration are divided by $\sin | {B}_{i}| \,=\sqrt{1/5}\,=0.447$. In panels (1)–(4), the analytical solutions and the results of numerical calculations are in good agreement. Panels (5)–(8) in Figure 3 are the same as panels (1)–(4) but for a_i = 5 × 10⁴ au. Just as seen in Figure 2, the agreement of the analytical solutions with the results of numerical calculations is worse. In particular, in panels (5) and (7), four of six bodies show an increase in Ω and L in the numerical calculations, which is never given by the analytical solutions. From several extra numerical calculations, we confirmed that these opposite evolutions seen in Ω and L are due to the radial component of the Galactic tide that breaks the conservation of j. In panels (6) and (8), the disagreement that arises from the shifts of the periods is quite significant; however, the amplitudes of the oscillations show good agreement even after 4.5 Gyr.

From the above comparisons, we conclude that the analytical solutions are basically useful for describing the orbital evolution, except for i and Ω of comets in the observable region (i.e., q ≲ 10² au). For bodies with a = 5 × 10⁴ au, the small differences between the periods given by the analytical solutions and the periods obtained from the numerical integrations pile up and are quite large at t ∼ a few gigayears. This could be understood simply as the results of the shifts of oscillation/circulation of orbital evolution. Therefore, the time evolution normalized by the periods obtained by numerical integrations is well reproduced by the analytical solutions.

In this section, we apply the analytical solutions derived in Section 2 to fictional observable long-period comets entering the planetary region from the Oort cloud. We assume that the comets initially have very elongated orbits given by planetary scattering on the ecliptic plane. For these comets, setting e_i ≃ 1 is a good approximation and it makes the solutions simple. We call this approximation the quasi-rectilinear approximation. Since the planetary scattering does not give high inclinations (see Appendix), we assume i_i = i_⊙ = 60° and Ω_i = Ω_⊙ = 186° to be on the ecliptic plane. Using this approximation, we compare the periods of χ, Ω, and L and investigate the relation among them especially for comets in the observable region.

4.1. Preparation

In this section, we first derive the explicit expression for α₀, α₁, and α₂ (Equation (70)) and their relation, k², ${w}_{1}^{2}$, and ${w}_{2}^{2}$ (Equations (32), (37), and (51)) under the quasi-rectilinear approximation, which gives 0 < j² ≪ 1. Using the expressions, we calculate the values of Π that appear in solutions to Ω and L.

4.1.1. Solutions and Parameters

Substituting e = e_i, B = B_i, and $1-{e}_{i}^{2}={j}^{2}/{\cos }^{2}{i}_{\odot }=4{j}^{2}$ into Equation (10), we obtain

$$\begin{eqnarray}&&c=3{j}^{2}+5(1-4{j}^{2}){\sin }^{2}{B}_{i}.\end{eqnarray} \tag{ 63 }$$

Using −i_⊙ ≤ B_i ≤ i_⊙, we have the minimum and maximum values of c as

$$\begin{eqnarray}\begin{array}{lcl}{c}_{\min } & = & 3{j}^{2}\simeq 0\,{\rm{for}}\,{B}_{i}=0\\ {c}_{\max } & = & 3{j}^{2}+\displaystyle \frac{15}{4}(1-4{j}^{2})\simeq \displaystyle \frac{15}{4}\,{\rm{for}}\,{B}_{i}={i}_{\odot }.\end{array}\end{eqnarray} \tag{ 64 }$$

where j² ≪ 1 is used.

From Equation (27), the explicit expressions of the solutions are approximated as

$$\begin{eqnarray}&&{\chi }_{1}^{* }\simeq \displaystyle \frac{5{j}^{2}}{5+4{j}^{2}-c}+\displaystyle \frac{100{j}^{4}}{{\left(5+4{j}^{2}-c\right)}^{3}}+{ \mathcal O }({j}^{6})\ll 1\end{eqnarray} \tag{ 65 }$$

$$\begin{eqnarray}&&{\chi }_{2}^{* }\simeq \displaystyle \frac{5+4{j}^{2}-c}{4}-\displaystyle \frac{5{j}^{2}}{5+4{j}^{2}-c}+{ \mathcal O }({j}^{4})\end{eqnarray} \tag{ 66 }$$

where j² ≪ 1 is used but the term ${ \mathcal O }({j}^{4})$ in Equation (65) is left for the comparison with j².

To evaluate the relation between j² and χ₁*, we substitute c = c_min since the difference between j² and ${\chi }_{1}^{* }$ becomes minimum for c = c_min. It is calculated as

$$\begin{eqnarray}\begin{array}{rcl}{\chi }_{1}^{* }-{j}^{2} & \simeq & \displaystyle \frac{5{j}^{2}}{5+4{j}^{2}-{c}_{\min }}+\displaystyle \frac{100{j}^{4}}{{\left(5+4{j}^{2}-{c}_{\min }\right)}^{3}}-{j}^{2}\\ & = & \ \displaystyle \frac{{j}^{4}(75-10{j}^{2}-{j}^{4})}{{\left(5+{j}^{2}\right)}^{3}}\,\gt \,0.\end{array}\end{eqnarray} \tag{ 67 }$$

Therefore, the relation between j² and ${\chi }_{1}^{* }$ is always as ${j}^{2}\lt {\chi }_{1}^{* }$.

From Equations (26), (64), and (65), the relation between ${\chi }_{0}^{* }$ and ${\chi }_{1}^{* }$ is always as ${\chi }_{1}^{* }\lt {\chi }_{0}^{* }$. The relation between ${\chi }_{0}^{* }$ and ${\chi }_{2}^{* }$ depends on the value of c. The difference is calculated as

$$\begin{eqnarray}&&{\chi }_{0}^{* }-{\chi }_{2}^{* }\simeq c-\displaystyle \frac{5-c}{4}=\displaystyle \frac{5(c-1)}{4},\end{eqnarray} \tag{ 68 }$$

where the terms ${ \mathcal O }({j}^{2})$ are neglected. As the value of c at the separatrix (Equation (9)) is approximated as c + j² ≃ c, the relation is approximately given as

$$\begin{eqnarray}&&\begin{array}{l}{\chi }_{0}^{* }\lt {\chi }_{2}^{* }\,\mathrm{for}\,c\lt 1\,(\mathrm{circulation})\\ {\chi }_{0}^{* }\gt {\chi }_{2}^{* }\,\mathrm{for}\,c\gt 1\,(\mathrm{libration}).\end{array}\end{eqnarray} \tag{ 69 }$$

Summarizing the relations, we have

$$\begin{eqnarray}&&\begin{array}{c}\begin{array}{l}({\alpha }_{0},{\alpha }_{1},{\alpha }_{2})\\ \,=\,\left\{\begin{array}{l}({\chi }_{1}^{\ast },{\chi }_{0}^{\ast },{\chi }_{2}^{\ast })\,{\rm{f}}{\rm{o}}{\rm{r}}\,c\lt 1\,({\rm{c}}{\rm{i}}{\rm{r}}{\rm{c}}{\rm{u}}{\rm{l}}{\rm{a}}{\rm{t}}{\rm{i}}{\rm{o}}{\rm{n}})\\ ({\chi }_{1}^{\ast },{\chi }_{2}^{\ast },{\chi }_{0}^{\ast })\,{\rm{f}}{\rm{o}}{\rm{r}}\,c\gt 1\,({\rm{l}}{\rm{i}}{\rm{b}}{\rm{r}}{\rm{a}}{\rm{t}}{\rm{i}}{\rm{o}}{\rm{n}})\end{array}\right.\end{array}\end{array}\end{eqnarray} \tag{ 70 }$$

and α₁ = α₂ for the separatrix for c = 1.

Using Equations (32) and (70) and j² ≪ 1, we have

$$\begin{eqnarray}&&{k}^{2}=\left\{\begin{array}{l}\displaystyle \frac{4c}{5-c}\,\lt \,1\ \mathrm{for}\ c\lt 1\,(\mathrm{circulation})\\ \displaystyle \frac{5-c}{4c}\,\lt \,1\ \mathrm{for}\ c\gt 1\,(\mathrm{libration})\end{array}\right.\end{eqnarray} \tag{ 71 }$$

and k² = 1 for the separatrix for c = 1.

Using j² < α₀ ≪ 1 and Equations (67) and (70), the parameter ${w}_{1}^{2}$ that appears in Π for Ω given by Equation (37) is expressed as

$$\begin{eqnarray}&&{w}_{1}^{2}=\displaystyle \frac{{\alpha }_{1}-{\alpha }_{0}}{{j}^{2}-{\alpha }_{0}}\simeq -\infty ,\end{eqnarray} \tag{ 72 }$$

and another parameter ${w}_{2}^{2}$, which is for L given by Equation (51), is

$$\begin{eqnarray}&&\begin{array}{c}{w}_{2}^{2}\simeq \left\{\begin{array}{l}\displaystyle \frac{4c}{5-c}\,{\rm{f}}{\rm{o}}{\rm{r}}\,c\lt 1\,({\rm{c}}{\rm{i}}{\rm{r}}{\rm{c}}{\rm{u}}{\rm{l}}{\rm{a}}{\rm{t}}{\rm{i}}{\rm{o}}{\rm{n}})\\ 1\,{\rm{f}}{\rm{o}}{\rm{r}}\,c\gt 1\,({\rm{l}}{\rm{i}}{\rm{b}}{\rm{r}}{\rm{a}}{\rm{t}}{\rm{i}}{\rm{o}}{\rm{n}}).\end{array}\right.\end{array}\end{eqnarray} \tag{ 73 }$$

4.1.2. Complete Integrals of the Third Kind in Ω and L

Depending on the values and relation between k² and ${w}_{1}^{2}$ or ${w}_{2}^{2}$, the complete elliptic integrals of the third kind are expressed in different forms (Byrd & Friedman 1971).

For $0\lt -{w}_{1}^{2}\lt \infty$, which is called "case I" in Byrd & Friedman (1971),

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Pi }}({w}_{1}^{2},k) & = & \displaystyle \frac{{k}^{2}K}{{k}^{2}-{w}_{1}^{2}}\\ & & -\ \displaystyle \frac{\pi {w}_{1}^{2}{{\rm{\Lambda }}}_{0}(\varphi ,k)}{2\sqrt{{w}_{1}^{2}(1-{w}_{1}^{2})({w}_{1}^{2}-{k}^{2})}}\quad ({w}_{1}^{2}\lt 0),\end{array}\end{eqnarray} \tag{ 74 }$$

where

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{0}(\varphi ,k)=\displaystyle \frac{2}{\pi }[{EF}(\varphi ,k^{\prime} )+{KE}(\varphi ,k^{\prime} )-{KF}(\varphi ,k^{\prime} )],\end{eqnarray} \tag{ 75 }$$

where E = E(k) and E(φ, k') are complete and normal elliptic integrals of the second kind,

$$\begin{eqnarray}&&\sin \varphi =\sqrt{\displaystyle \frac{{w}_{1}^{2}}{{w}_{1}^{2}-{k}^{2}}},\end{eqnarray} \tag{ 76 }$$

and

$$\begin{eqnarray}&&k{{\prime} }^{2}=1-{k}^{2}.\end{eqnarray} \tag{ 77 }$$

Using Equation (72), we have $\sin \varphi =1$. Then F(φ, k') and E(φ, k') in Equation (75) become complete elliptic integrals of the first and second kinds, K' = K(k') and E' = E(k'), respectively, and

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{0}=\displaystyle \frac{2}{\pi }[{EK}^{\prime} +E^{\prime} K-{KK}^{\prime} ]=\displaystyle \frac{2}{\pi },\end{eqnarray} \tag{ 78 }$$

which is called Legendre's relation. Then ${\rm{\Pi }}({w}_{1}^{2},k)$ is approximated as

$$\begin{eqnarray}&&{\rm{\Pi }}({w}_{1}^{2},k)\simeq \displaystyle \frac{\pi }{2}\displaystyle \frac{1}{\sqrt{-{w}_{1}^{2}}}.\end{eqnarray} \tag{ 79 }$$

For ${k}^{2}\lt {w}_{2}^{2}\lt 1$, which is called "case II" in Byrd & Friedman (1971),

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{0}(\vartheta ,k)=\displaystyle \frac{2}{\pi }[{EF}(\vartheta ,k^{\prime} )+{KE}(\vartheta ,k^{\prime} )-{KF}(\vartheta ,k^{\prime} )],\end{eqnarray} \tag{ 80 }$$

$$\begin{eqnarray}&&\sin \vartheta =\sqrt{\displaystyle \frac{1-{w}_{2}^{2}}{k{{\prime} }^{2}}}\end{eqnarray} \tag{ 81 }$$

and

$$\begin{eqnarray}&&{\rm{\Pi }}({w}_{2}^{2},k)=K+\displaystyle \frac{\pi {w}_{2}[1-{{\rm{\Lambda }}}_{0}(\vartheta ,k)]}{2\sqrt{({w}_{2}^{2}-{k}^{2})(1-{w}_{2}^{2})}}.\end{eqnarray} \tag{ 82 }$$

For the special case of ${w}_{2}^{2}={k}^{2}$, which is true for the case for circulation,

$$\begin{eqnarray}&&{\rm{\Pi }}({k}^{2},k)=\displaystyle \frac{E}{1-{k}^{2}}.\end{eqnarray} \tag{ 83 }$$

For ${w}_{2}^{2}\simeq 1$ in the case of circulation, we have $\sin \,\vartheta =0$ and then

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{0}(\vartheta ,k)\simeq 0.\end{eqnarray} \tag{ 84 }$$

Therefore, ${\rm{\Pi }}({w}_{2}^{2},k)$ is approximated as

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Pi }}({w}_{2}^{2},k)\\ \ \simeq \ \left\{\begin{array}{l}\displaystyle \frac{E}{1-{k}^{2}}\ \mathrm{for}\ c\lt 1\,(\mathrm{circulation})\\ K+\displaystyle \frac{\pi {w}_{2}}{2\sqrt{({w}_{2}^{2}-{k}^{2})(1-{w}_{2}^{2})}}\ \mathrm{for}\ c\gt 1\,(\mathrm{libration}).\end{array}\right.\end{array}\end{eqnarray} \tag{ 85 }$$

4.2. Periods of χ, Ω, and L and Their Ratios

4.2.1. Mean P_χ

Assuming a uniform distribution of ω_i between 0° and 360°, the mean value of ${\sin }^{2}{\omega }_{i}$ is $\langle {\sin }^{2}{\omega }_{i}\rangle =1/2$, which corresponds to B_i ≃ 378. Using i_i = 60° and ${\sin }^{2}{\omega }_{i}=1/2$, the mean value of c given by Equation (7) is approximated as

$$\begin{eqnarray}&&\langle c\rangle \sim \displaystyle \frac{15}{8}.\end{eqnarray} \tag{ 86 }$$

Since $\langle c\rangle \gt 1$, it is in the case of libration of ω. Therefore, from Equations (70) and (71), we have

$$\begin{eqnarray}&&{\alpha }_{2}={\chi }_{0}^{* }\simeq \langle c\rangle +{j}^{2},\quad {k}^{2}=\displaystyle \frac{5-\langle c\rangle }{4\langle c\rangle }.\end{eqnarray} \tag{ 87 }$$

Substituting Equations (86) and (87) and the expansion in series of K given as

$$\begin{eqnarray}&&K=\displaystyle \frac{\pi }{2}\left[1+\displaystyle \frac{1}{4}{k}^{2}+\displaystyle \frac{9}{64}{k}^{4}\ldots \right]\quad ({k}^{2}\lt 1)\end{eqnarray} \tag{ 88 }$$

into Equation (34) and using j² ≪ 1, the mean value of P_χ for a given a and ρ₀ is estimated as

$$\begin{eqnarray}&&\langle {P}_{\chi }\rangle \simeq \displaystyle \frac{4K}{{A}_{1}\sqrt{{\alpha }_{2}}}\sim {10}^{9}\left(\displaystyle \frac{0.1}{\rho \,[{M}_{\odot }\,{\mathrm{pc}}^{-3}]}\right){\left(\displaystyle \frac{2\times {10}^{4}}{a[\mathrm{au}]}\right)}^{3/2}\,[\mathrm{yr}].\end{eqnarray} \tag{ 89 }$$

4.2.2. Ratio of ${P}_{{{\rm{\Omega }}}^{* }}$ to P_χ

Substituting Equations (36) and (79) into Equation (46), the period of Ω is approximated as

$$\begin{eqnarray}&&{P}_{{{\rm{\Omega }}}^{* }}\simeq \displaystyle \frac{16K}{{A}_{1}}\displaystyle \frac{\sqrt{({\alpha }_{1}-{\alpha }_{0})({\alpha }_{0}-{j}^{2})}}{{jc}}.\end{eqnarray} \tag{ 90 }$$

Using Equations (30), (34), and (90), the ratio of ${P}_{{{\rm{\Omega }}}^{* }}$ to P_χ is given as

$$\begin{eqnarray}&&\displaystyle \frac{{P}_{{{\rm{\Omega }}}^{* }}}{{P}_{\chi }}\simeq \displaystyle \frac{4\sqrt{{\alpha }_{0}{\alpha }_{1}{\alpha }_{2}-{\alpha }_{1}{\alpha }_{2}{j}^{2}}}{{jc}},\end{eqnarray} \tag{ 91 }$$

where j² < α₀ ≪ 1. Using Equations (26), (65), (66), and (70), the terms in the square root are calculated as

$$\begin{eqnarray}&&{\alpha }_{0}{\alpha }_{1}{\alpha }_{2}={\chi }_{0}^{* }{\chi }_{1}^{* }{\chi }_{2}^{* }\simeq \displaystyle \frac{5{{cj}}^{2}}{4}\end{eqnarray} \tag{ 92 }$$

and

$$\begin{eqnarray}&&{\alpha }_{1}{\alpha }_{2}{j}^{2}={\chi }_{1}^{* }{\chi }_{2}^{* }{j}^{2}\simeq \displaystyle \frac{{{cj}}^{2}(5-c)}{4}.\end{eqnarray} \tag{ 93 }$$

Substituting Equations (92) and (93) into Equation (91), we have

$$\begin{eqnarray}&&\displaystyle \frac{{P}_{{{\rm{\Omega }}}^{* }}}{{P}_{\chi }}\simeq 2.\end{eqnarray} \tag{ 94 }$$

Consequently, not only e, i, and ω, but also Ω displays a coupled evolution. This commensurability is established only under the quasi-rectilinear approximation. Figure 5 shows the ratios of ${P}_{{{\rm{\Omega }}}^{* }}$ to P_χ against e for B_i = 0, 378, and 60°. For any B_i, the ratio approaches 2 only when e ∼ 1. Note that Higuchi et al. (2007) already reported the commensurability but they did not show it in equations. Using the quasi-rectilinear approximation, Equation (38) is rewritten as

$$\begin{eqnarray}&&{\rm{\Omega }}\simeq {{\rm{\Omega }}}_{i}-{A}_{3}\left[{{\rm{\Pi }}}_{{\rm{\Omega }}}^{{\prime} }-{\rm{\Pi }}({u}_{0},{w}_{1}^{2},k)\right]-m\pi .\end{eqnarray} \tag{ 95 }$$

Zoom In Zoom Out Reset image size

4.2.3. Ratio of ${P}_{{L}^{* }}$ to P_χ

Substituting Equations (50) and (85) into Equation (60), the period of L is approximated as

$$\begin{eqnarray}&&\begin{array}{c}{P}_{{L}^{\ast }}\simeq \left\{\begin{array}{l}\displaystyle \frac{8\pi }{{A}_{1}j}{\left[1+\displaystyle \frac{E}{K}\right]}^{-1}\,{\rm{f}}{\rm{o}}{\rm{r}}\,c\lt 1\,({\rm{c}}{\rm{i}}{\rm{r}}{\rm{c}}{\rm{u}}{\rm{l}}{\rm{a}}{\rm{t}}{\rm{i}}{\rm{o}}{\rm{n}})\\ \displaystyle \frac{8K}{{A}_{1}\sqrt{c}}\,{\rm{f}}{\rm{o}}{\rm{r}}\,c\gt 1\,({\rm{l}}{\rm{i}}{\rm{b}}{\rm{r}}{\rm{a}}{\rm{t}}{\rm{i}}{\rm{o}}{\rm{n}}).\end{array}\right.\end{array}\end{eqnarray} \tag{ 96 }$$

Using Equations (34) and (96), we have the ratios of ${P}_{{L}^{* }}$ to P_χ as

$$\begin{eqnarray}&&\displaystyle \frac{{P}_{{L}^{* }}}{{P}_{\chi }}\simeq \left\{\begin{array}{l}\infty \ \mathrm{for}\ c\lt 1\,(\mathrm{circulation})\\ 2\ \mathrm{for}\ c\gt 1\,(\mathrm{libration})\end{array}\right.\end{eqnarray} \tag{ 97 }$$

where j² ≪ 1, 0 < E/K < 1, and α₂ ≃ c for c > 1 are used. Then Equation (54) is rewritten as

$$\begin{eqnarray}&&L\simeq \left\{\begin{array}{l}{L}_{i}\ \mathrm{for}\ c\lt 1\,(\mathrm{circulation})\\ {L}_{i}-{A}_{6}\left[{{\rm{\Pi }}}_{L}^{{\prime} }-{\rm{\Pi }}({u}_{0},{w}_{2}^{2},k)\right]-m\pi \ \mathrm{for}\ c\gt 1\,(\mathrm{libration}).\end{array}\right.\end{eqnarray} \tag{ 98 }$$

The almost fixed value of L beyond 4.5 Gyr for the case of circulation already seen in Section 2 is also explained with Equations (2), (4), and (94) as follows. Since the evolution of i is a periodic oscillation, the effect of cos i multiplied by $\tan \omega$ in Equation (4) can be assumed null when it is averaged over the period, i.e., as if $L\sim {\rm{\Omega }}+\omega$. Under the quasi-rectilinear approximation, Ω and ω for the case of circulation evolve with the same velocities but in opposite directions. Consequently, Ω and ω are canceled out and L evolves little.

For the case of libration, the second term in Equation (4) would be 0 on average over a period of libration. As a result, Ω and L evolve with the same period, which is twice P_χ. This relation is true for any e_i and does not require the quasi-rectilinear approximation.

4.3. Behavior in the Observable Region

As confirmed in Section 3, the time for q ∼ q_i is expressed as t ≃ mP_χ, where m is an integer. Using this relation and Equations (28), (95), and (98), we can use the combinations of orbital elements at t = mP_χ. as the prediction for q ≃ q_i ≲ 10² au. Table 1 summarizes the relation.

Table 1.  Orbital Elements for q ≃ q_i, i.e., at t = mP_χ under the Quasi-rectilinear Approximation

m	q	i	Ω	ω	L	sin B	Case
0	qi	ii	Ωi	ωi	Li	sin Bi	circulation/libration
odd	qi	ii	Ωi − π	ωi + π	Li	$-\sin {B}_{i}$	circulation
				ωi	Li − π	$\sin {B}_{i}$	libration
even	qi	ii	Ωi	ωi	Li	$\sin {B}_{i}$	circulation/libration

Download table as:  ASCIITypeset image

Figure 6 shows the analytical solutions to Ω, i_E, L, and $\sin B$ against q < 50 au with the results of numerical calculations for 4.5 Gyr presented in Section 3. Panel (1) in Figure 6 shows the behavior of Ω for a_i = 2 × 10⁴ au. As seen from Equation (15), the analytical solution to Ω (solid curves) drastically changes when q is near its minimum. For example, for a body with ω_i = 10° (black), Ω can have almost any value between 0 and 360° when q is small. The curves overlap every two oscillations of q since ${P}_{{{\rm{\Omega }}}^{* }}/{P}_{\chi }\simeq 2$. The disagreements between the analytical solution and the numerical calculations (circles) are not negligible in the observable region, i.e., q ≲ 10 au. Panel (2) in Figure 6 shows the behavior of i_E for a_i = 2 × 10⁴ au. Since i_E is a function of i and Ω as

$$\begin{eqnarray}&&\cos {i}_{{\rm{E}}}=\cos i\cos {i}_{\odot }+\sin i\sin {i}_{\odot }\cos {\rm{\Omega }},\end{eqnarray} \tag{ 99 }$$

the prediction of i_E in the observable region using the analytical solution is as difficult as that for i (see Figure 4) and as sensitive to q as that for Ω. Panels (5) and (6) in Figure 6 are the same as panels (1) and (2) but for a_i = 5 × 10⁴ au. The drift of the curves is seen more obviously than in panel (1) since bodies with a_i = 5 × 10⁴ au make more oscillations in 4.5 Gyr.

Zoom In Zoom Out Reset image size

Panels (3) and (4) in Figure 6 show the behavior of L and $\sin B$ for a_i = 2 × 10⁴ au, respectively, and panels (7) and (8) in Figure 6 are the same as panels (3) and (4), respectively, but for a_i = 5 × 10⁴ au. Both L and $\sin B$ are almost independent of q for q < 50 au. The agreement of the results of numerical calculations and the analytical solutions is good for a_i = 2 × 10⁴ au. For a_i = 5 × 10⁴ au, the disagreement is larger than ∼30° for some bodies beyond 4.5 Gyr but still much better than that in Ω or i_E. Therefore, we conclude that the relation among q, L, and $\sin B$ in Table 1, for any m, is safely satisfied for q in the observable region.

4.4. The Empty Ecliptic

Based on the standard scenario of the formation of the Oort cloud, the initial orbital elements of the Oort cloud comets are restricted as follows; q_i ≲ 30 au to be near a giant planet, i_i ≃ i_⊙ = 60° and Ω_i ≃ Ω_⊙ = 186° to be on the ecliptic plane. Also ω_i is uniformly distributed for 0 ≤ ω_i < 360°, L_i is given by Equation (2), and $\sin {B}_{i}$, in order to be on the ecliptic plane, is given by

$$\begin{eqnarray}&&\tan {B}_{i}=\sin ({L}_{i}-{{\rm{\Omega }}}_{\odot })\tan {i}_{\odot }.\end{eqnarray} \tag{ 100 }$$

The relation between L and $\sin B$ in Table 1 defines two planes in the Galactic coordinates. As L_i and B_i are assumed to be on the ecliptic plane the points that satisfy 0 ≤ L < 360° and $\sin B$ given by Equation (100) for L draw a curve on the $L\mbox{--}\sin B$ plane, by definition. Another set of points for an odd m that satisfy 0 ≤ L < 360° and $-\sin B$ draw a curve that defines a second plane, which is formed by a rotation of the ecliptic around the Galactic pole by 180°:

$$\begin{eqnarray}&&\tan B=-\sin (L-{{\rm{\Omega }}}_{\odot })\tan {i}_{\odot }.\end{eqnarray} \tag{ 101 }$$

We call this plane "the empty ecliptic," since it is not initially populated and this plane and the ecliptic are symmetrical about the plane perpendicular to the Galactic plane through the intersection of the ecliptic plane and the Galactic plane (just like the focus and the empty focus of an ellipse). If L and B of long-period comets in the observable region are concentrated on these two planes, it would constitute observational evidence that the comets were on the ecliptic plane at t = 0. Comets with relatively small semimajor axes such as a ∼ 10⁴ au, which satisfy P_χ = 4.5 Gyr (i.e., m = 1 for present), are predicted to be on the empty ecliptic for their first return to the planetary region.

Figure 7 shows solar system bodies with q > 1 au and a > 10³ au or e > 1 taken from the JPL Small Body Database Search Engine on 2020 June 5² on the $L\mbox{--}\sin B$ plane. 277 bodies with e ≤ 1 and 296 bodies with e > 1 are indicated with open and filled circles, respectively. The bodies are divided into three groups with L as indicated by colors. Two interstellar objects, 1I/2017 U1 ('Oumuamua) and 2I/2019 Q4 (Borisov), are additionally shown with squares for reference. We used the osculating orbital elements to calculate L and B using Equations (2) and (3). To calculate L and B for bodies with e > 1, we replace ω in Equations (2) and (3) as

$$\begin{eqnarray}&&\omega \to \omega +(\pi -{f}_{\infty }),\end{eqnarray} \tag{ 102 }$$

where

$$\begin{eqnarray}&&\cos {f}_{\infty }=-\displaystyle \frac{1}{e}.\end{eqnarray} \tag{ 103 }$$

To evaluate the concentrations on the two planes, we define the new angle ε as

$$\begin{eqnarray}&&\tan \varepsilon =\displaystyle \frac{\tan B}{\sin (L-{{\rm{\Omega }}}_{\odot })}.\end{eqnarray} \tag{ 104 }$$

The angle ε is interpreted as a longitude around the intersection of the ecliptic and the Galactic plane. For the ecliptic and empty ecliptic planes, ε = i_⊙ = 60° and ε = −i_⊙ = −60°, respectively. The solid and dashed curves in Figure 7 show the ecliptic and empty ecliptic planes, respectively. Curves for ε = 0, ±30°, and ±80° are also shown as thin dashed curves in Figure 7.

Zoom In Zoom Out Reset image size

Figure 8 shows the distribution of ε. There are two sharp peaks not exactly at the ecliptic or empty ecliptic plane but near them. An isotropic distribution would be flat in ε.

Zoom In Zoom Out Reset image size

Panel (1) in Figure 9 shows the distribution of $\sin B$. The depletions around b = 0 and b = ±90° found by Luest (1984) and Delsemme (1987) are seen (note that b = −B); however, the shape of the distribution depends on the regions of L' = L − Ω_⊙. Panels (2)–(4) in Figure 9 show the distribution of $\sin B$ for regions of L: blue (−30° < L' < 30°, 150° < L' < 210°), orange (30° < L' < 60°, 120° <L' < 150°, 210° < L' < 240°, 300° < L' < 330°), and dark orange (60° < L' < 120°, 240° < L' < 300°). The less-sharp, rather a broad peak at $| \sin B| \leqslant 0.5$ in panel (2) and the sharpest double peaks at $| \sin B| \gt 0.5$ in panel (4) are explained as a consequence of the double peaks in the distribution of ε. If the depletions are the result of the strength of the Galactic tide as Delsemme (1987) explained, the distributions are expected to be independent of L' since the strength of the Galactic tide is independent of L. Therefore, we conclude that the concentration of comets on the ecliptic and empty ecliptic planes is a better explanation than that by Delsemme (1987).

Zoom In Zoom Out Reset image size

We derived analytical solutions for L and B, the Galactic longitude and latitude of the aphelion direction of bodies orbiting around the Sun, with the perturbation from the Galactic disk in an axisymmetric approximation. We used the solutions to predict the distribution of observed long-period comets in the Galactic coordinates. To evaluate the analytical solutions, we performed numerical calculations of the orbital evolution including the effect of the radial component of the Galactic tide, which is neglected in the derivation of the analytical solutions. Our findings are summarized as follows.

• 1.  
  For bodies initially having eccentricities e ≃ 1, the analytical solutions and the results of numerical calculations show good agreement in the time evolutions normalized by the periods of the evolution. However, the Galactic inclination i and the longitude of the node in the Galactic coordinates Ω show non-negligible disagreement, especially when their perihelion distances q are small enough to be in the observable region since the vertical angular momenta of the bodies are not completely conserved.
• 2.  
  In the orbital evolution, three periods are defined: (1) P_χ, the period of oscillation of e (i.e., q), i, and B, (2) ${P}_{{{\rm{\Omega }}}^{* }}$, the mean period of circulation of Ω, and (3) ${P}_{{L}^{* }}$, the mean period of circulation of L. The period of the argument of perihelion, ω, is P_χ for the case of libration (i.e., von Zeipel–Lidov–Kozai mechanism) and 2P_χ for the case of circulation. Under the quasi-rectilinear approximation (i.e., e_i ∼ 1), the following relations are established among the analytical solutions; ${P}_{{{\rm{\Omega }}}^{* }}/{P}_{\chi }\simeq 2$ for all cases and ${P}_{{L}^{* }}/{P}_{\chi }\simeq 2$ and ${P}_{{L}^{* }}\sim \infty$ for the cases of libration and circulation of ω, respectively. Consequently, the evolutions of q and Ω are coupled in any case, the evolutions of q and L are coupled in the case of libration of ω, and L evolves very little in the case of circulation of ω.
• 3.  
  Under the quasi-rectilinear approximation, the coupled evolutions of L and B of bodies initially on the ecliptic plane with q in the planetary region draw two curves on the L–B plane when their q are small. One corresponds to the ecliptic plane and the other to the empty ecliptic defined by the longitude around the intersection of the ecliptic and the Galactic plane ε (see Equation (104)), ε = 60° and −60°, respectively. The numerical calculations showed that the coupling of L and B is quite stable at any q in the observable region and confirmed that ε would be a reliable indicator of the dynamical character of observed long-period comets. The evolution of Ω is also coupled with that of q, i, ω, and B under the rectilinear approximation; however, ε is a better indicator than Ω and others since the time variation of Ω is quite large at small q and the value of i is not nicely reproduced by the analytical solution at small q.
• 4.  
  The distribution of ε of observed solar system bodies with q > 1 au and semimajor axis a > 10³ au or e > 1 shows double peaks that might correspond to the ecliptic and empty ecliptic planes although their locations are not exactly at ε = ±60°. The concentration of the bodies on the ecliptic and empty ecliptic planes explains the depletions around B = 0 and B = ±90° (Luest 1984; Delsemme 1987) better than the explanation by Delsemme (1987).

The concentration of long-period comets from the Oort cloud on the ecliptic and empty ecliptic planes is observational evidence that the Oort cloud comets were planetesimals initially on the ecliptic plane. We expect the concentrations even when we consider the effect of passing stars. Perturbations from passing stars change the conserved quantities and may break the relation between q, B, and L more or less; however, it takes a much longer time to change the eccentricity vector (i.e., L and B) than to change i (Higuchi & Kokubo 2015). Therefore, we suggest that observers, including the space mission Comet Interceptor, focus on the ecliptic plane and/or the empty ecliptic plane to find dynamically new comets.

What we showed in Section 5 is a brief examination. An investigation of the distribution of observed small bodies has to include many factors. The bodies should be carefully chosen from the database and examined by classes defined by their original semimajor axes calculated with non-gravitational forces for active comets (e.g., Królikowska et al. 2014). The orbital elements during the last perihelion passage would be a key to the dynamical evolution if they encountered any of the planets (Kaib & Quinn 2009; Fouchard et al. 2018). Comparison with numerical calculations with all perturbations from the Galactic disk, stars, and planets is also important. The long-term behavior found in numerical calculations of comets in Fouchard et al. (2020) describes the empty ecliptic plane. Observational bias should also be taken into account. Detailed examination of the distribution of long-period comets will be our future work. The all-sky survey by the Large Synoptic Survey Telescope will provide valuable information for this study.

I am grateful to Melaine Saillenfest and Takashi Ito for their comments that greatly improved the quality of this paper and to Marc Fouchard and Eiichiro Kokubo for their comments and discussion that led to this work. I also thank David Jewitt for carefully reading the manuscript. Finally, I thank Giovanni B. Valsecchi for reviewing this paper. The numerical computations were in part carried out on PC cluster at Center for Computational Astrophysics, National Astronomical Observatory of Japan. This work was partially supported by the Programme National de Planetologie (PNP) of CNRS/INSU, co-funded by CNES.

Assume a planet in a circular orbit with semimajor axis a_planet and a body with inclination with respect to the orbital plane of the planet i. The relative velocity between the unperturbed orbits of the planet and the comet Δv is written as the two components

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{v}_{z} & = & {v}_{\mathrm{planet}}\sqrt{3-T}\sin i,\\ {\rm{\Delta }}{v}_{r} & = & {v}_{\mathrm{planet}}\sqrt{3-T}\cos i,\end{array}\end{eqnarray} \tag{ A1 }$$

where Δv_z is the component perpendicular to the orbital plane of the body, Δv_r is defined as ${\rm{\Delta }}{v}_{r}=\sqrt{{({\rm{\Delta }}v)}^{2}-{({\rm{\Delta }}{v}_{z})}^{2}}$, v_planet is the velocity of the planet, and T is the Tisserand parameter with respect to the planet defined by

$$\begin{eqnarray}&&T=\displaystyle \frac{{a}_{\mathrm{planet}}}{a}+2\sqrt{\displaystyle \frac{a}{{a}_{\mathrm{planet}}}(1-{e}^{2})}\cos i,\end{eqnarray} \tag{ A2 }$$

where a and e are the semimajor axis and eccentricity of the body. Under the two-body approximation, the velocity of the body after planetary scattering v is v = v_planet + Δv at maximum and v = v_planet − Δv at minimum for each component. Therefore, the maximum change in inclination Δi given by planetary scattering is approximated as a function of T and i,

$$\begin{eqnarray}&&\tan {\rm{\Delta }}i\simeq \displaystyle \frac{{\rm{\Delta }}{v}_{z}}{{v}_{\mathrm{planet}}-{\rm{\Delta }}{v}_{r}}\simeq \displaystyle \frac{\sqrt{3-T}\sin i}{1-\sqrt{3-T}\cos i}.\end{eqnarray} \tag{ A3 }$$

Bodies that we are interested in are those that have a' ≫ a_planet, where a' is the semimajor axis after planetary scattering. The value of T for these bodies is estimated as follows. To gain a change in velocity large enough to become nearly parabolic, a large Δv (i.e., a small T) is required. For example, for a parabolic or hyperbolic orbit, Δv needs to satisfy $| {{\boldsymbol{v}}}_{\mathrm{planet}}+{\rm{\Delta }}{\boldsymbol{v}}| \geqslant {v}_{\mathrm{esc}}$, where ${v}_{\mathrm{esc}}=\sqrt{2}{v}_{\mathrm{planet}}$ is the escape velocity at the heliocentric distance r = a_planet. This condition leads to $T\leqslant 2\sqrt{2}$. On the other hand, the chance of having an effective encounter with a planet becomes larger for smaller Δv (i.e., for larger T) since the gravitational radius of the planet is proportional to Δv⁻². Higuchi et al. (2006) numerically showed that the efficiency of having large a' is higher for smaller Δv. Therefore, we estimate that $T\simeq 2\sqrt{2}$ would be the favored value for a' ≫ a_planet.

Substituting $T\simeq 2\sqrt{2}$ and i ≪ 1 into Equation (A3), we have Δi ≪ 1. More general and detailed discussion about post-encounter inclination will be given by G. B. Valsecchi et al. (2020, in preparation.).