Anisotropy of Long-period Comets Explained by Their Formation Process

Long-period comets coming from the Oort cloud are thought to be planetesimals formed in the planetary region on the ecliptic plane. We have investigated the orbital evolution of these bodies due to the Galactic tide. We extended Higuchi et al. (2007) and derived the analytical solutions to the Galactic longitude and latitude of the direction of aphelion, $L$ and $B$. Using the analytical solutions, we show that the ratio of the periods of the evolution of $L$ and $B$ is very close to either 2 or $\infty$ for initial eccentricities $e_i\simeq 1$, as is true for the Oort cloud comets. From the relation between $L$ and $B$, we predict that Oort cloud comets returning to the planetary region are concentrated on the ecliptic plane and a second plane, which we call the"empty ecliptic."This consists in a rotation of the ecliptic around the Galactic pole by 180$^\circ$. Our numerical integrations confirm that the radial component of the Galactic tide, which is neglected in the derivation of the analytical solutions, is not strong enough to break the relation between $L$ and $B$ derived analytically. Brief examination of observational data shows that there are concentrations near both the ecliptic and the empty ecliptic. We also show that the anomalies of the distribution of $B$ of long-period comets mentioned by several authors are explained by the concentrations on the two planes more consistently than by previous explanations.


INTRODUCTION
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, 1992Brasser 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.
Using the orbital elements, the unit vector of the direction of aphelion in the Galactic coordinates r Q is written as − cos ω cos Ω + sin ω sin Ω cos i − cos ω sin Ω − sin ω cos Ω cos i − sin ω sin i    (1) Then L and sin B are written as where θ = atan sin ω cos i cos ω for cos ω < 0 π + atan sin ω cos i cos ω for cos ω > 0. (4)

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 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 R = − ν 2 0 4 a 2 sin 2 i 1 − e 2 + 5e 2 sin 2 ω , where ν 0 = √ 4πGρ 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: c = sin 2 i 1 − e 2 + 5e 2 sin 2 ω .
The simplified z-component of the angular momentum, which is a conserved quantity under the axisymmetric approximation of the potential, is written as j = 1 − e 2 cos i.
For e 2 > 0, the sufficient condition on B for libration is given as Matese & Whitman (1989) defined the value B = asin 1/5 ≃ ±26.6 • as a barrier that the latitude of perihelion cannot migrate across.

Time variations
Substituting Equation (6) into Lagrange's planetary equations (Murray & Dermott 1999), we obtain the time variations of e, i, Ω, and ̟ as 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 ̟ = ω + Ω, From Equation (2), we have From Equation (4), we have Substituting Equations (4), (3), (14), and (17) into Equation (19), we have Substituting Equation (20) into Equation (18), we have Introducing χ = 1 − e 2 , Equation (13) can be rewritten as Using Equations (7) and (8), ω can be rewritten as Substituting sin ω and cos ω from Equation (23) into Equation (22), where The solution to Equation (24) is expressed using a Jacobian elliptic function, sn (Byrd & Friedman 1971), where we define where F is a normal elliptic integral of the first kind with modulus k and the amplitude φ 0 , The sign of F (φ 0 , k) in Equation (31) depends on ω i , the value of ω at t = 0; it is negative for sin(2ω i ) > 0 and positive for sin(2ω i ) < 0. From Equation (28), we learn that χ oscillates between the maximum (α 1 ) and minimum (α 0 ) values according to the parameter u, with the period of u = 2K, where K = K(k) is a complete elliptic integral of the first kind with modulus k. The period given in time is (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 4 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 (eq.(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 χ . 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 2 ω and sin B are obtained from equations (23) and (3), respectively. To calculate the real ω from sin 2 ω, 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|.

Solution to Ω
Substituting Equations (7) and (8) into Equation (15) and replacing 1 − e 2 with χ, we have where The integration of (35) with t is rewritten as where Ω i is the value of Ω at t = 0 and Since sn 2 (u, k) oscillates with a period of u = 2K, we split u ′ as 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 where Π(K, w 2 1 , k) is a complete elliptic integral of the third kind and The period of Ω can be estimated by the linear approximation as where The period of Ω * is obtained as .
(46) Figure 1 shows P Ω * in orange against B i for a = 2 × 10 4 au, q i = 10 au, and i i = 60 • . The behavior of P Ω * is quite similar to that of P χ . The approximate relation between P Ω * and P χ is shown in Section 4.

Solution to L
From Equations (3), (8), and (23), sin 2 B is expressed with c, j, and χ as Then, tan 2 B is written as Substituting Equation (28) into Equation (48), tan 2 B is expressed as where Substituting Equation (48) into Equation (18), we have where As well as Ω, the integration of (52) with t is expressed using an elliptic integral of the third kind as where L i is the value of L at t = 0 and The period of L is also estimated in the same manner as Ω using the linear approximation, where The period of L * is obtained as The black dashed curve in Figure 1 shows P L * plotted against B i for a = 2 × 10 4 au, q i = 10 au, and i i = 60 • . The behavior of P L * looks identical to that of P Ω * for B i > 26 • .6. In contrast, for B i < 26 • .6, P L * suddenly becomes 10 2−5 times larger than that for B i > 26 • .6. For this example, P L * for c < 1 is much longer than the age of the solar system.

EVALUATION OF ANALYTICAL SOLUTIONS WITH NUMERICAL CALCULATIONS
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.

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 where r is the position of the body with respect to the Sun and f tide is the Galactic tide, where x ′ , y ′ , and z ′ give the position of the body in rotating coordinates centered on the Sun, Ω 0 is the circular frequency (i.e., the angular speed of the Sun in the Galaxy), ν 0 = √ 4πGρ is the vertical frequency, and ρ is the total density in the solar neighborhood. We adopt Ω 0 = 26 km s −1 (e.g., Binney & Tremaine 1987) and ρ = 0.1 M ⊙ pc −3 (Holmberg & Flynn 2000), which give ν 2 0 /Ω 2 0 ≃ 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) 2 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.

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 4 au and 5 × 10 4 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).
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 Figure 2 shows the evolution of q. The range of the time when the bodies arrive in the observable region, i.e., q 10 2 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 4 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 4 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 4 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 4 au. The conservation of j is worse than for a i = 2 × 10 4 au but still we can say it is approximately conserved even when the radial component of the Galactic tide is included.
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 ω (eq. (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 2 au. Panel (4) in Figure 3 shows the evolution of sin B, which is expressed as a function of i and ω (eq.(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 | = 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 4 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 2 au). For bodies with a = 5 × 10 4 au, the small differences between the periods given by the analytical solutions and the ones obtained from the numerical integrations pile up and are quite large at t ∼ a few Gyr. 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.

QUASI-RECTILINEAR APPROXIMATION AND APPLICATION TO LONG-PERIOD COMETS
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 as 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.

solutions and parameters
Substituting e = e i , B = B i , and 1 − e 2 i = j 2 / cos 2 i ⊙ = 4j 2 into Equation (10), we obtain Using −i ⊙ ≤ B i ≤ i ⊙ , we have the minimum and maximum values of c as where j 2 ≪ 1 is used. From Equation (27), the explicit expressions of the solutions are approximated as where j 2 ≪ 1 is used but the term O(j 4 ) in Equation (65) is left for the comparison with j 2 .
To evaluate the relation between j 2 and χ * 1 , we substitute c = c min since the difference between j 2 and χ * 1 becomes minimum for c = c min . It is calculated as Therefore, the relation between j 2 and χ * 1 is always as j 2 < χ * 1 .   From Equations (26), (64), and (65), the relation between χ * 0 and χ * 1 is always as χ * 1 < χ * 0 . The relation between χ * 0 and χ * 2 depends on the value of c. The difference is calculated as where the terms O(j 2 ) are neglected. As the value of c at the separatrix (eq. (9)) is approximated as c + j 2 ≃ c, the relation is approximately given as Summarizing the relations, we have and α 1 = α 2 for the separatrix for c = 1. Using Equations (32) and (70) and j 2 ≪ 1, we have and k 2 = 1 for the separatrix for c = 1. Using j 2 < α 0 ≪ 1 and Equations (67) and (70), the parameter w 2 1 that appears in Π for Ω given by Equation (37) is expressed as and another parameter w 2 2 , which is for L given by Equation (51), is w 2 2 ≃ 4c 5−c for c < 1 (circulation) 1 for c > 1 (libration) .

Complete integrals of the third kind in Ω and L
Depending on the values and relation between k 2 and w 2 1 or w 2 2 , the complete elliptic integrals of the third kind are expressed in different forms (Byrd & Friedman 1971).
For 0 < −w 2 1 < ∞, which is called "case I" in Byrd & Friedman (1971), where where E = E(k) and E(ϕ, k ′ ) are a complete and normal elliptic integrals of the second kind, and Using Equation (72), we have sin ϕ = 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 which is called Legendre's relation. Then Π(w 2 1 , k) is approximated as For k 2 < w 2 2 < 1, which is called "case II" in Byrd & Friedman (1971), and For the special case of w 2 2 = k 2 , which is true for the case for circulation, For w 2 2 ≃ 1 in the case of circulation, we have sin ϑ = 0 and then Therefore, Π(w 2 2 , k) is approximated as for c > 1 (libration). Assuming a uniform distribution of ω i between 0 • and 360 • , the mean value of sin 2 ω i is sin 2 ω i = 1/2, which corresponds to B i ≃ 37 • .8. Using i i = 60 • and sin 2 ω i = 1/2, the mean value of c given by Equation (7) is approximated as Since c > 1, it is in the case of libration of ω. Therefore, from Equations (70) and (71), we have Substituting Equations (86) and (87) and the expansion in series of K given as into Equation (34) and using j 2 ≪ 1, the mean value of P χ for a given a and ρ 0 is estimated as

Ratio of PΩ * to Pχ
Substituting Equations (36) and (79) into equation (46), the period of Ω is approximated as Using Equations (30), (34), and (90), the ratio of P Ω * to P χ is given as where j 2 < α 0 ≪ 1. Using Equations (26), (65), (66), and (70), the terms in the square root are calculated as and Substituting Equations (92) and (93) into Equation (91), we have 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 Ω * to P χ against e for B i = 0, 37 • .8, 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

Ratio of PL * to Pχ
Substituting Equations (50) and (85) into equation (60), the period of L is approximated as Using Equations (34) and (96), we have the ratios of P L * to P χ as where j 2 ≪ 1, 0 < E/K < 1, and α 2 ≃ c for c > 1 are used. Then Equation (54) is rewritten as 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 ω in Equation (4) can be assumed null when it is averaged over the period, i.e., as if L ∼ Ω + ω. 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.

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 2 au.  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 4 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 Ω * /P χ ≃ 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 4 au. Since i E is a function of i and Ω as cos i E = cos i cos i ⊙ + sin i sin i ⊙ cos Ω, the prediction of i E in the observable region using the analytical solution is as difficult as much as that for i (see Figure  4) and as sensitive to q as much as that for Ω. Panels (5) and (6) in Figure 6 are the same as panels (1) and (2) but for a i = 5 × 10 4 au. The drift of the curves is seen more obviously than in panel (1) since bodies with a i = 5 × 10 4 au make more oscillations in 4.5 Gyr. Panels (3) and (4) in Figure 6 show the behavior of L and sin B for a i = 2 × 10 4 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 4 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 4 au. For a i = 5 × 10 4 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.

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 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 − 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 • : We call this plane as "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 values such as a ∼ 10 4 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.
5. OBSERVATIONAL DATA Figure 7 shows solar system bodies with q > 1 au and a > 10 3 au or e > 1 taken from the JPL Small Body Database Search Engine on 2020 June 5 1 on the L − 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 where To evaluate the concentrations on the two planes, we define the new angle ε as The angle ε is interpreted as a longitude around the intersection of the ecliptic and the Galactic plane.  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. 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 ε.
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| ≤ 0.5 in panel (2) and the sharpest double peaks at | sin B| > 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).

SUMMARY AND DISCUSSION
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 ascending 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 Ω * , 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 Ω * /P χ ≃ 2 for all cases and P L * /P χ ≃ 2 and P L * ∼ ∞ 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 the semimajor axis a > 10 3 au or e > 1 shows the 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 an 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) is the one that 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.

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

A. THE INITIAL INCLINATION GIVEN BY PLANETARY SCATTERING
Assume a planet in a circular orbit with the semimajor axis a planet and a body with the 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 is written as the two components where ∆v z is the component perpendicular to the orbital plane of the body, ∆v r is defined as ∆v r = (∆v) 2 − (∆v z ) 2 , v planet is the velocity of the planet, and T is the Tisserand parameter with respect to the planet defined by T = a planet a + 2 a a planet (1 − e 2 ) cos i,