2025 Volume 101 Issue 3 Pages 143-176
In 1962, Yoshihide Kozai reported his findings on the secular dynamics of asteroids moving in orbits with high inclination and eccentricity. In contrast to the classic understanding of the stability of planetary motion in the solar system, Kozai showed that asteroids can significantly change their orbital shape over a long timescale in an oscillatory manner between nearly circular orbits and highly elliptic orbits. An anti-correlated variation between orbital inclination and eccentricity characterizes this oscillation. The importance of Kozai’s work in understanding the dynamical evolution of various systems was recognized decades later, including the fields of irregular satellites of planets, Oort Cloud, extrasolar planets, binary star systems, type Ia supernovae, planet climate, merging black hole systems, and so on.
In this paper, we review Yoshihide Kozai’s work on three-body dynamics and its increasing impact on various fields of astronomy. We describe the perturbation theory of the three-body problem in Section 1 and summarize the essence of Kozai’s work so that readers can trace the mathematical procedures in a comprehensive and self-contained way in Section 2. Section 3 presents examples where Kozai’s work has recently been highlighted. Section 4 describes his academic life. We have added a discussion in Section 5.
For this article, we consulted two standard textbooks1),2) and a monograph3) among many past studies in this field.1
We begin by considering the case where two objects exert gravity on each other. This circumstance is called the gravitational two-body problem or simply the two-body problem. In this paper, we focus on the motion of objects under Newtonian dynamics and do not consider relativistic effects.
Suppose the two objects’ masses are M and m, respectively. The equations of motion for each object can be described by two sets of ordinary differential equations. However, we can reduce the two-body problem to the relative motion of a single body by describing the system on the coordinates concerning one of the objects. In this case, we can consider only one equation of motion that describes the relative motion.
Assume r is the position vector from mass point M to mass point m, and t is time. Then, the equation for the relative motion of mass m with respect to mass M is as follows:
| \[\frac{d^2\boldsymbol{r}}{dt^2} = - \mu \frac{\boldsymbol{r}}{r^3},\] | [1] |
where
| \[\mu = {\cal G} (M + m),\] | [2] |
and \({\cal G}\) as the gravitational constant. The right-hand side of Eq. [1] denotes the mutual gravity working between the mass points M and m.
The solution of the equations of motion of the two-body system [1] is the well-known Keplerian motion (or the Keplerian orbit). The Keplerian motion follows the three laws of Kepler, empirically formulated in the early 17th century by analyzing the positions of the solar system planets observed by Tycho Brahe over decades. There are three types of shapes in the orbits of Keplerian motion: ellipses (including the circle), parabolas, and hyperbolas. The object’s orbital motion is around the other object located at a focus of these quadratic curves.
Six variables, called orbital elements, define the Keplerian motion. Here is a brief description of the orbital elements, particularly for the elliptic orbits. Consult Fig. 1 for the geometric relationship of each element (except for the sixth time-dependent one).
(Color online) Keplerian orbital elements that describe an elliptic orbit. The panel a explains the elements in the orbital plane. a is the semimajor axis and b is the semiminor axis defined by \(b = a \sqrt{1-e^2}\) using eccentricity e. f denotes an angle referred to as the true anomaly. r is the position vector from mass M to mass m. The point P denotes the location of the pericenter. In this panel, the mass point M is located at one of the ellipse foci. The distance between the ellipse’s center and the foci is ae. The panel b explains the elements concerning the reference plane. Ω (or h) denotes the longitude of ascending node, and ω (or g) denotes the argument of pericenter. In this panel, the mass point M is still located in one of the foci and not in the ellipse’s center. The point A denotes the location of the ascending node and B denotes the location of the descending node. We assume that mass m orbits counterclockwise in this figure. Υ denotes a reference direction fixed in the inertial space. The reference plane is also fixed in the inertial space.
・Semimajor axis a
This element describes the size of the orbit. It represents the longer radius of the ellipse, i.e., half of its longest diameter. This quantity is directly related to the orbital period of Keplerian motion. If we write the orbital period as P and define the mean angular velocity (the mean motion) n as follows:
| \[n = \frac{2\pi}{P},\] | [3] |
Then, Kepler’s third law states the following:
| \[n^2 a^3 = \mu .\] | [4] |
The semimajor axis is also directly related to the orbital energy of the two-body system. The total orbital energy (the sum of kinetic energy and potential energy) is given as follows:
| \[E = -\frac{\mu}{2a} .\] | [5] |
・Eccentricity e
This element describes the orbit’s shape and indicates the degree of deviation from a complete circle. The closer e is to 0, the more circular the orbit is. As e approaches 1, the ellipse becomes more elongated. When e ＝ 1, the orbit is a parabola. When e ＞ 1, the orbit is a hyperbola. When e ≥ 1, the orbit is not closed, and the motion is not periodic. The eccentricity e is connected to the angular momentum of the two-body system and is proportional to the quantity \(\sqrt{\mu a ( 1-e^2)}\) .
・Inclination I
This element is the angle between the Keplerian orbit plane and the reference plane. For planets in the solar system, the reference plane is typically the ecliptic plane (Earth’s orbital plane around the Sun) at a given epoch. Inclination I is measured in radians, ranging from 0 to π. When \(I < \frac{\pi}{2}\) , the orbital motion is called prograde. When \(I > \frac{\pi}{2}\) , the orbital motion is called retrograde. When \(I = \frac{\pi}{2}\) , it is called the polar orbit. The inclination I is related to the vertical component of the angular momentum vector of the two-body system and is proportional to the quantity \(\sqrt{\mu a ( 1-e^2)} \cos I\) .
・Longitude of ascending node Ω (or h)
This element describes the orientation of the Keplerian orbit in three-dimensional space. The range of this element is from 0 to 2π. This defines the angle between a reference direction (usually the direction of the vernal equinox Υ fixed in inertial space) and the orbit’s ascending node. The ascending node is where the orbiting body crosses the reference plane from south to north.
・Argument of pericenter ω (or g)
This element describes the orientation of the orbit within its orbital plane. The range of this element is from 0 to 2π. Specifically, it is the angle between the ascending node of the orbit and the point closest to the central mass (the pericenter or periastron). When the central mass is the Sun or Earth, this is called the argument of perihelion or the argument of perigee, respectively.
The above five variables are time-independent constants for Keplerian orbits. The remaining sixth orbital element is time-varying.
・A variable directly connected to time
The true anomaly, the angle f in Fig. 1, is the angular distance from the pericenter to the object. Its rate of change is not proportional to time unless eccentricity e is zero. The most commonly used element directly related to time t is the mean anomaly l, which is defined by:
| \[l = n (t - t_0),\] | [6] |
where n is the mean motion defined by Eq. [3] and it is a constant. t_0 is the time when l ＝ 0 is realized, i.e., the time the object passes its pericenter. The mean anomaly l is connected to an angle called the eccentric anomaly u as follows:
| \[l = u - e \sin u .\] | [7] |
The eccentric anomaly u is connected to the true anomaly f through the geometric relationship shown in Fig. 2. For example, the following relationships are often useful (e.g., Ref. 1, their p. 24-25):
| \[\tan \frac{f}{2} = \sqrt{\frac{1+e}{1-e}} \tan \frac{u}{2} ,\] | [8] |
and
| \[\cos f = \frac{\cos u - e}{1 - e \cos u}, \quad \sin f = \frac{\sqrt{1-e^2}\sin u}{1 - e \cos u} .\] | [9] |
(Color online) Geometric relationship between eccentric anomaly u and true anomaly f. Here, O is the ellipse’s center (the black curve). The outer blue circle is the circumscribed circle of the orbital ellipse, within which the eccentric anomaly is defined.
Equation [7] is called Kepler’s equation, and it links the position and time of a celestial object in its orbit. However, Kepler’s equation is transcendental and cannot be solved algebraically for u as u ＝ u(l). Instead, numerical methods (iterative techniques) or polynomial approximations are typically used to find the solution of u ＝ u(l) for a given l. Danby (1992) (Ref. 6, his Section 6.6) explains in detail how to employ numerical methods to solve Kepler’s equation in an accessible manner for beginning learners.
The distance between the two bodies in the Keplerian orbit r is described using true anomaly f as follows:
| \[r = \frac{a \left( 1-e^2 \right)}{1 + e \cos{f}} .\] | [10] |
It can also be expressed as follows using eccentric anomaly u:
| \[r = a (1 - e\cos u) .\] | [11] |
For more detailed discussions on the general principles of the two-body problem as described in this section, readers should consult standard, classical textbooks such as Smart (1953),7) Brouwer and Clemence (1961),1) Danby (1992),6) Boccaletti and Pucacco (1996, 1998),8),9) Murray and Dermott (1999),2) Merritt (2013),10) or Tremaine (2023).11)
1.2. Three-body problem.The two-body problem is integrable. It has an exact analytic solution, the Keplerian motion. However, by just adding one object to the system, the system exhibits quite a different behavior: The three-body problem is not integrable (e.g., Ref. 12, their Section 3), and we have no exact analytic solution except in very few exceptional cases.13),14)
Although the three-body problem is generally non-integrable, the nearly integrable hierarchical three-body system described below serves as a valuable platform for studying the motion of objects in the solar system (e.g., Ref. 15). A hierarchical three-body system comprises a massive central primary mass (M) accompanied by a smaller secondary mass (m' ＜ M), as well as an even smaller tertiary mass (m ＜ m'). Figure 3a schematically shows the relative locations of the three masses. The tertiary mass can orbit inside or outside the (M,m') binary as shown in Fig. 3b. Unless the orbits of the secondary and tertiary around the primary get too close or intersect each other, the two binaries, (M,m') and (M,m), usually both behave almost like the Keplerian motion. In this case, we can obtain their orbital solution primarily through perturbation methods (e.g., Ref. 16). The binary (M,m') would make a pure Keplerian motion if m ≪ m', not being disturbed by m at all, while the motion of m is affected by the binary (M,m'). This is called the restricted three-body problem (denoted hereafter R3BP). In particular, when the orbit of m' in the binary (M,m') is circular, the system results in the circular restricted three-body problem (CR3BP). This is the problem Kozai (1962)17) intended to study.
(Color online) Schematic illustrations of the relative motion of the hierarchical three-body system. a: Relative geometric configuration of a three-body system (M, m', m) with respect to the primary mass M. The vectors r and r' are position vectors from the mass M to the masses m and m', respectively. We presume that the masses have the relation of M ＞ m' ＞ m here. b: Two typical patterns of the restricted three-body problem. (i) The inner problem where the orbit of the perturbed body m (the green ellipse) lies inside the orbit of the perturbing body m' (the black ellipse). (ii) The outer problem where the perturbed body’s orbit (the green ellipse) lies outside the orbit of the perturbing body (the black ellipse). This figure is modified and adopted from Ito and Ohtsuka (Ref. 3, their Fig. 1).
Following the derivation described in standard textbooks (e.g., Refs. 1, 2 and 6), we have the equation of motion of the mass point m relative to the central mass M as follows:
| \[\frac{d^2\boldsymbol{r}}{dt^2} + \mu\frac{\boldsymbol{r}}{r^3} = \nabla R .\] | [12] |
On the other hand, the equation of motion of the mass point m' relative to the central mass M is expressed as follows:
| \[\frac{d^2\boldsymbol{r}'}{dt^2} + \widetilde{\mu} \frac{\boldsymbol{r}'}{{r'}^3} = \nabla R' ,\] | [13] |
where
| \[\widetilde{\mu} = {\cal G} (M + m') .\] | [14] |
R in Eq. [12] as well as R' in Eq. [13] are called disturbing functions. The function forms are as follows:
| \[R = \frac{{\cal G} m'}{\Delta}-{\cal G} m' \frac{\boldsymbol{r} \cdot \boldsymbol{r}'}{{r'}^3} ,\] | [15] |
and
| \[R' = \frac{{\cal G} m }{\Delta}-{\cal G} m \frac{\boldsymbol{r}' \cdot \boldsymbol{r}}{r^3} ,\] | [16] |
with the conventional notation for mutual distance
| \[\Delta \equiv \left| \boldsymbol{r}' - \boldsymbol{r} \right|=\left| \boldsymbol{r} - \boldsymbol{r}' \right| .\] | [17] |
The disturbing functions represent the gravitational interaction between mass m and mass m'. The first term of the disturbing function [15] or [16] is called the direct part and represents the major component of the mutual perturbation between the two masses. The second term is called the indirect part. The indirect part vanishes or becomes constant when we average the system over the orbit (e.g., Ref. 2).
We expand the disturbing function, R or R', in an infinite series of orbital elements using Legendre polynomials. Applying the cosine formula to the triangle m-M- m' with angle S in Fig. 3a, we obtain
| \[\left| \boldsymbol{r}' - \boldsymbol{r} \right|^2 = r^2 + {r'}^2 - 2 r r' \cos S .\] | [18] |
Using Eq. [18], we can expand the inverse of Δ in Eq. [17] through the Legendre polynomials Pj (cos S). When r ＜ r', the expansion becomes
| \begin{aligned} \frac{1}{\Delta} &= \frac{1}{r'} \left( 1 - 2 \frac{r}{r'} \cos S + \left( \frac{r}{r'} \right)^2 \right)^{-\frac{1}{2}} \\ & = \frac{1}{r'} \sum_{j=0}^{\infty} \left( \frac{r}{r'} \right)^j P_j (\cos S) . \end{aligned} | [19] |
On the other hand, when r ＞ r', it becomes
| \begin{aligned} \frac{1}{\Delta} &= \frac{1}{r} \left( 1 - 2 \frac{r'}{r} \cos S + \left( \frac{r'}{r} \right)^2 \right)^{-\frac{1}{2}} \\ & = \frac{1}{r} \sum_{j=0}^{\infty} \left( \frac{r'}{r} \right)^j P_j (\cos S) . \end{aligned} | [20] |
In the remainder of this section, we discuss the restricted three-body problem (R3BP) and assume m ≪ m'. Let us remember here that m is the mass of the perturbed object, and m' is the mass of the perturbing object. The orbital condition r ＜ r' required for the expression of 1/Δ in Eq. [19] indicates that we are dealing with the inner problem (see Fig. 3b(i)). In this case, the term with j ＝ 0 in Eq. [19] does not depend on r and disappears through its differentiation (e.g., Ref. 2, their Section 6.3, p. 229). Therefore, we can ignore the j ＝ 0 term in Eq. [19]. In addition, we have the following relationship,
| \[\boldsymbol{r} \cdot \boldsymbol{r'} = r r' \cos S = r r' P_1 (\cos S),\] | [21] |
and the indirect part of the disturbing function cancels out the term of j ＝ 1 in Eq. [19]. Therefore, to construct the disturbing function of the inner problem, we consider only the j ≥ 2 terms in the expansions of 1/Δ in Eq. [19] as
| \[\frac{1}{\Delta} = \frac{1}{r'} \sum_{j=2}^{\infty} \left( \frac{r}{r'} \right)^j P_j (\cos S) .\] | [22] |
For the outer problem where condition r ＞ r' holds (see Fig. 3b(ii)), the expression of 1/Δ in Eq. [20] should be employed. The treatment of the outer problem is more complicated than the inner problem. However, the principle is the same. Readers can consult previous works (e.g., Refs. 18 and 19) or textbooks2) for more details on the treatment of outer problems. Yoshihide Kozai also provided a concise review of the treatment of the outer problem and its application to Pluto’s motion, publishing it in the Proceedings of the Japan Academy, Ser. B.20)
In the circular restricted three-body problem (CR3BP), the length of the position vector of the perturber (r') with respect to the primary mass has a constant value equivalent to its semimajor axis a' because the perturber’s eccentricity is zero (e' ＝ 0). When r' ＝ a', we do not need to consider odd terms (j ＝ 3, 5, 7, …) in the expansion of Eqs. [19] or [20] because they vanish after the averaging procedure (see Appendix A for more details). Therefore, the disturbing function of the inner CR3BP is simplified from Eqs. [15] and [19] to:
| \[ R = \frac{{\cal G} m'}{a'} \sum_{n=1}^{\infty} \left( \frac{r}{a'} \right)^{2n} P_{2n} (\cos S) . \] | [23] |
This is the disturbing function that Kozai (1962)17) considered. Section 2 summarizes his discussion based on this function.
Note that in the hierarchical N-body problem, including the three-body problem addressed in this paper, a set of canonical coordinates known as the Jacobi coordinates are commonly used (e.g., Refs. 21 and 22). In the Jacobi coordinates, the position and velocity of object k are given relative to the barycenter of objects 0 through k － 1 (e.g., see Fig. 9.20 of Ref. 2). Usually, index 0 is assigned to the central mass (e.g., the Sun), and other indices are assigned to the remaining objects in the order of increasing semimajor axis. For the Sun-Jupiter-asteroid system, Jupiter’s coordinate origin is the Sun, consistent with the relative coordinates employed in this paper. However, the asteroid’s coordinate origin was the barycenter of the Sun-Jupiter binary rather than the Sun itself. Because the Jacobi coordinates are intrinsically canonical, the Hamiltonian can be described in a more compact form than relative or barycentric coordinates (e.g., Ref. 11). In addition, the Jacobi coordinates facilitate the development of a highly efficient algorithm for numerical orbit integration in a weakly perturbed N-body system.22) It is straightforward to show that 1/Δ included in the direct part of the disturbing function for the inner problem as seen in Eq. [19] or in Eq. [22] can be derived from the general Hamiltonian described using the Jacobi coordinates (Ref. 3, their p. 7). However, note again that we do not use the Jacobi coordinates in this paper and instead stay on the relative coordinates, as Kozai employed the relative coordinates in his study.
The leading publication we introduce here is Kozai (1962),17) published in The Astronomical Journal. This paper was more widely and quickly recognized than other studies on the same subject, as mentioned later. Kozai (1962)17) was selected as one of the 53 “Selected Fundamental Papers Published this Century in the Astronomical Journal and the Astrophysical Journal”,23) and was republished in 1999 on the centennial anniversary of the American Astronomical Society (AAS). A short editorial by Brian Geoffrey Mardsen, the AAS president at that time, is available.24)
In Kozai (1962),17) the disturbing function of CR3BP in the Hamiltonian form was first described. Next, the system’s degrees of freedom were reduced using the averaging technique, and the system was brought into an integrable form. In this section, we first describe his method and then summarize his conclusions. Additionally, we compare Kozai’s calculations with the results of the numerical integration from a modern perspective.
2.1. Averaging of disturbing function.Having the disturbing function R in the form of Eq. [23] in hand, Kozai carried out the double averaging of R by the mean anomalies of the perturbed and perturbing bodies. This procedure is an important part of reducing the degrees of freedom of the entire system (see Section 2.2). When the motion of an object around the central mass is close to the Keplerian motion, the variation rate of its mean anomaly is orders of magnitude larger than those of the other orbital elements. Hence, it is justified to eliminate its mean anomaly by averaging, assuming that the other orbital elements of the object do not significantly change over a period of mean anomaly, namely one orbital motion around the central body. Eliminating the mean anomaly using the averaging procedure can be regarded as part of the canonical transformation that divides the system’s Hamiltonian into periodic and secular parts (e.g., Ref. 9). Historically, this procedure was devised by Delaunay (1860, 1867),25),26) and substantially developed by von Zeipel (1916, 1917).27)-30)
To perform the averaging procedure, we assume no significant resonant relationship exists between the mean motions of the perturbed and perturbing bodies. Bearing this assumption in mind, we select the n-th term of the disturbing function R for the inner problem in Eq. [23], and refer to it as R2n. Defining a mass factor μ' as
| \[\mu' = {\cal G} m' ,\] | [24] |
and use it for simplicity, we obtain
| \[R_{2n} = \frac{\mu'}{a'} \left( \frac{r}{a'} \right)^{2n} P_{2n} (\cos S) .\] | [25] |
First, we average R2n by the mean anomaly l' of the perturbing body. Using the symbols ＜ and ＞ for averaging, it is
| \[\left< R_{2n} \right>_{l'} = \frac{\mu'}{a'} \left( \frac{r}{a'} \right)^{2n} \left< P_{2n} \right>_{l'},\] | [26] |
where
| \[\left< P_{2n} \right>_{l'} = \frac{1}{2\pi} \int^{2\pi}_0 P_{2n} (\cos S) dl'.\] | [27] |
The angle S is expressed by the orbital angles through a spherical trigonometric relationship (e.g., Ref. 17, Eq. (7) on p. 592)
| \[\cos S = \cos (f+g) \cos (f'+g') + \cos I \sin (f+g) \sin (f'+g'),\] | [28] |
where f, f' are true anomalies, and g, g' are the arguments of pericenter of the perturbed and perturbing bodies, respectively, and I is the mutual inclination angle of the two orbital planes measured at the node of the two orbits. As for how the angle S is defined in space, we made Fig. 4 and showed its geometry there. The orbital plane of the perturbing body was selected as the reference plane to measure g from the mutual node. The argument of pericenter of the perturbing body (g') is not defined in CR3BP. Therefore, in Eq. [28], we regard f' ＋ g' as a single fast-oscillating variable represented by f'. In practice, we can replace the integral over the mean anomaly of the perturbing body ∫dl' for an integral over its true anomaly ∫df' in the discussion here because f' ＝ l' in CR3BP.
(Color online) The spatial configuration of the angle S between perturbed and perturbing bodies in CR3BP.
To obtain 〈P2n〉l' of Eq. [27], we calculate the time average of cos2n S by l' as follows:
| \[\left< \cos^{2n} S \right>_{l'} = \frac{1}{2\pi} \int^{2\pi}_0 \cos^{2n} S dl' .\] | [29] |
Then, we average 〈R2n〉l' of Eq. [26] by the mean anomaly l, of the perturbed body as
| \[\left< \left< R_{2n} \right>_{l'} \right>_{l} = \frac{\mu'}{a'} \left( \frac{a}{a'} \right)^{2n} \frac{1}{2\pi} \int^{2\pi}_0 \left( \frac{r}{a} \right)^{2n} \left< P_{2n} \right>_{l'} dl .\] | [30] |
Equation [30] is the general form of the n-th term of the doubly averaged disturbing function for the inner CR3BP. If we define the ratio of the semimajor axes as α ＝ \(\frac{a}{a'}\) ＜ 1, this term has the magnitude of O (α2n).
2.2. Reduction of degrees of freedom.The averaging procedure of the disturbing function described in Section 2.1 is a part of the procedure to reduce the degrees of freedom of the CR3BP system and transform it into an integrable one. Thus, we can obtain a comprehensive understanding of the motion of a perturbed object in CR3BP without directly solving the equations of motion. The following is a summary of Kozai’s explanation of how the degrees of the disturbing function can be reduced.
In his 1962 paper,17) Kozai denoted m as the mass of the asteroid and m' as the mass of Jupiter. The solar mass is set to unity. Instead of Kepler orbital elements, Kozai used the so-called Delaunay elements. The Delaunay elements are a set of canonical variables that describe the Keplerian motion, and their relationship to the Kepler orbital elements are as follows (e.g., Ref. 8):
| \begin{aligned} L &= \sqrt{\mu a}, \\ &\hskip3em l = \; \mbox{mean anomaly}, \\ G &= L \sqrt{1-e^2}, \\ &\hskip3em g= \; \mbox{argument of perihelion}\; (\omega), \\ H &= G \cos I, \\ &\hskip3em h = \; \mbox{longitude of ascending node}\; (\Omega), \\ \end{aligned} | [31] |
where l, g, h are the canonical coordinates, and L, G, H are the corresponding conjugate momenta. The element G is proportional to the magnitude of the angular momentum vector of the perturbed object (see the description about eccentricity e in Section 1.1), and H is proportional to its vertical component (see the description regarding inclination I in Section 1.1).
Let us symbolically denote the Hamiltonian F for R3BP using the Delaunay elements as F(L, G, H, l, g, h; L', G', H', l', g', h'). Here, L, G, H, l, g, h are the canonical variables of the perturbed object, and L', G', H', l, g', h' are those of the perturbing object. The Hamiltonian F is the sum of the disturbing part (practically equivalent to the disturbing function R) and the orbital energy part (derived from the object’s Keplerian motion), and it originally had six degrees of freedom. However, the Hamiltonian contains h and h' (longitude of ascending node) only in the combination h － h' if we choose the invariable plane as the reference plane. In other words, the Hamiltonian is symmetric regarding the rotation around the system’s total angular momentum vector. This characteristic is called Jacobi’s elimination of the nodes (e.g., Refs. 31-34) and a typical outcome is h － h' ＝ π (e.g., Refs. 35 and 36). This relationship allows us to remove h and h' from the Hamiltonian, and the conjugate momenta H and H' become constants of motion. Consequently, the original Hamiltonian with six degrees of freedom is converted into F(L, G, l, g; L', G,' l,' g') with four degrees of freedom.
Then, Kozai averages the Hamiltonian (in other words, disturbing function) using the method described in Section 2.1 and further reduces its degrees of freedom. This procedure eliminates mean anomalies of the perturbed body (l) and the perturbing body (l') from the Hamiltonian. Consequently, their conjugate momenta L and L' become constants of motion, which yields another result through the definition of L and L' in Eq. [31]: the semimajor axes of both objects ( a and a') are constant. The Hamiltonian with four degrees of freedom is now transformed into a new Hamiltonian F(G, g; G', g')x-math> with two degrees of freedom. Here, we recall that Kozai (1962) addressed the CR3BP where the pericenter of the perturbing body is not defined (or g' can be considered constant). \(G' = L'\sqrt{1-{e'}^2}\) is also a constant parameter. Therefore, we do not need to consider g' or G' in the Hamiltonian, and the function form of the Hamiltonian is F(G, g) with just one degree of freedom.
Because we do not consider energy dissipation, the Hamiltonian F is a constant of motion. This indicates that any combination of (G, g) uniquely determines the dynamical state of the motion of a perturbed body. In other words, the doubly averaged CR3BP is integrable.
Recall that the Hamiltonian F can be divided into a Keplerian part (denoted \({\cal K}\) ) and a perturbation part (denoted W). In general, it can be written as follows:
| \[F = {\cal K} + m' W .\] | [32] |
The factor m' multiplied by W on the right-hand side of Eq. [32] denotes that the magnitude of the perturbation is proportional to the mass of the perturber. In general, \({\cal K}\) is a function of one of the Delaunay elements \(L = \sqrt{\mu a}\) , and is expressed as \({\cal K} \propto -L^{-2}\) . Since the semimajor axis of the perturbed body a is constant after the averaging operation, \({\cal K}\) is also constant. Therefore, in the following, we discuss only the perturbation part of the Hamiltonian, W. Now that we know W depends only on G and g, the canonical equations of motion of this system are as follows:
| \[\frac{dG}{dt} = -m'\frac{\partial W}{\partial g}, \quad \frac{dg}{dt} = m' \frac{\partial W}{\partial G} .\] | [33] |
Due to the traditional sign difference between celestial mechanics and analytic mechanics, W is equivalent to disturbing function R with the sign reversed (i.e., W ＝ －R).
Let us note that for the expression of K \((\propto -L^{-2})\) in Eq. [32], readers can refer to the expression of F0 on p. 553 of Brouwer and Clemence (1961)1) except for the sign change from ＋ to －. Additionally, see the expression of H0 in Eq. (10.90) on p. 285 of Boccaletti and Pucacco (1998),9) with the sign change from ＋ to －. Note also that the canonical equations of motion [33] can be directly and numerically integrated, and a time-dependent solution can be obtained in the form of G(t) and g(t). This line of studies using the truncated disturbing function is described as being semi-analytic, and there is a great deal of literature on it (e.g., Ref. 37). However, we do not mention it in this paper because Kozai did not do it either.
The conservation of one of the Delaunay elements of the perturbed body, H, defined in Eq. [31], plays a crucial role in Kozai’s study. This quantity is proportional to the vertical component of the angular momentum of the perturbed body. Kozai introduced a new, constant parameter Θ which is proportional to the square of H as follows:
| \[\Theta = \left( 1-e^2 \right) \cos^2 I .\] | [34] |
The fact that Θ is a constant means that eccentricity e and inclination I oscillate in an anti-correlated manner.
2.3. Stationary point.Kozai (1962)17) then searched for conditions under which the perturbed object’s argument of pericenter librates in this Hamiltonian system. More specifically, he searched for the location of the stationary points that the perturbation part of the doubly averaged disturbing Hamiltonian W can have. Note that the libration represents an oscillatory motion of an angle around a particular fixed value. The angle variable in libration does not complete a full circle (2π). On the other hand, circulation is a continuous motion of an angle in one direction that completes a full circle. At these stationary points, g and e (so G) of the perturbed body are supposed to be constant. To search for stationary points, we must know the specific function form of the disturbing function. Let us remember that α has been defined as the ratio of the semimajor axis of the perturbed object a and the perturbing object a' as α ＝ a/a'. In the inner CR3BP where α ＜ 1, Kozai selected the terms of the disturbing function in leading order, O(α2). This is equivalent to taking terms up to 2n ＝ 2 from the general form of the disturbing function described in Eq. [30]. Subsequently, the disturbing part of the Hamiltonian W (= －R) becomes as follows:
| \begin{aligned} W_{O\left(\alpha^2\right)} &= -\frac{1}{16} \frac{{\cal G}}{a'} \left( \frac{a}{a'} \right)^2 \big[ 15 e^2 \sin^2 I \cos 2 g\notag\\ &\hskip2em - \big(3e^2+2 \big) \big( 3\sin^2 I -2 \big) \big]. \end{aligned} | [35] |
Readers can refer to Appendix B for details on how \(W_{O\left(\alpha^2\right)}\) in Eq. [35] was derived from Eq. [30]. Note that the truncation of the disturbing function on the order of α2, as in Eq. [35], is commonly referred to as the quadrupole-level (or quadrupole-order) approximation. In contrast, truncation at the order of α3 is known as the octupole-level (or octupole-order) approximation (e.g., Ref. 38, see her p. 443).
Applying \(W_{O\left(\alpha^2\right)}\) in Eq. [35] to W in Eq. [33], we can rewrite the canonical equations of motion for G and g as follows:
| \begin{aligned} \frac{dG}{dt}&=-m' \frac{\partial W_{O\left(\alpha^2\right)}}{\partial g}\\ &=\frac{\mu'}{a'} \left( \frac{a}{a'}\right)^2\frac{15}{8} e^2 \sin^2 I \sin 2g, \end{aligned} | [36] |
| \begin{aligned} \frac{dg}{dt} &= m' \frac{\partial W_{O\left(\alpha^2\right)}}{\partial G}\\ &=-\frac{\mu'}{a'} \left( \frac{a}{a'}\right)^2\frac{3}{8G} \big[ \big( 5\cos^2 I - \big(1-e^2\big) \big)\\ &\hskip1em -5\big( \cos^2 I - \big(1-e^2\big) \big) \cos 2g \big] . \end{aligned} | [37] |
Let us determine where the stationary points of the disturbing Hamiltonian W are located at the quadrupole-level approximation using Eqs. [36] and [37]. At the stationary points, both of the following equations must be satisfied:
| \[\frac{dG}{dt} = 0, \quad \frac{dg}{dt} = 0 .\] | [38] |
From Eq. [36] and the first equation of Eq. [38], we know that the stationary points can be realized when sin 2g = 0. This means cos 2g ＝ ＋1 or cos 2g ＝ －1. When cos 2g ＝ ＋1, from Eq. [37] and the second equation of Eq. [38] we have
| \[4 \left( 1 - e^2 \right) = 0,\] | [39] |
which yields e ＝ 1 . However, we do not consider such a particular case e ＝ 1.
When cos 2g ＝ －1, from Eq. [37] and the second equation of Eq. [38] we have
| \[10 \cos^2 I - 6 \left( 1-e^2 \right) = 0 .\] | [40] |
Using Eq. [34], Eq. [40] yields the following relation:
| \[\Theta = \frac{3}{5} \left( 1-e^2 \right)^2 .\] | [41] |
As long as we deal with elliptic orbits where 0 ≤ e ＜ 1, we have
| \[0 ＜ \left( 1-e^2 \right)^2 \leq 1 ,\] | [42] |
and thus Eq. [41] gives
| \[\Theta \leq \frac{3}{5} .\] | [43] |
Equation [43] means that the doubly averaged disturbing Hamiltonian \(W_{O\left(\alpha^2\right)}\) in Eq. [35] can have stationary points only when the condition [43] is satisfied ― in other words, the parameter Θ has an upper limit beyond which the perturbed body’s argument of pericenter cannot librate. Equation [43], together with the definition of Θ in Eq. [34], indicates that there is a lower limit for the perturbed object’s inclination for the libration of its argument of pericenter to occur. We obtain this value by setting e ＝ 0, and Eq. [43] becomes
| \[\cos^2 I \leq \frac{3}{5} .\] | [44] |
By assuming I ＜ 90° (the prograde motion of the object), Eq. [44] can be translated as follows:
| \[I \geq \cos^{-1} \sqrt{\frac{3}{5}} = 39^\circ . 231 \cdots.\] | [45] |
Equation [45] is a representation of the necessary condition for the perturbed object’s argument of pericenter to librate at the quadrupole-level approximation ― when the object’s inclination I satisfies this condition at e ＝ 0, the object’s argument of pericenter g can librate.
Note that the range of I in Eq. [45] or the range of Θ in Eq. [43], which describes the necessary condition for the occurrence of the libration of the argument of the pericenter of the perturbed object, actually depends on the ratio of the semimajor axis of the perturbed and the perturbing objects, α ＝ a/a'. When deriving Eqs. [44] and [45], we assumed a ≪ a', therefore α → 0. This is equivalent to considering only the lowest order term (2n ＝ 2) of the averaged disturbing function [30]. In other words, we limited ourselves to use \(W_{O\left(\alpha^2\right)}\) in Eq. [35] for the disturbing Hamiltonian. When α is not negligibly small, we either use higher-order expansion of the disturbing function, or give up the analytical expansion and resort to numerical quadrature (see Section 2.5). The numerical quadrature is to perform a double-averaging integral of the direct term of the disturbing function R in Eq. [15] over the mean anomaly of the perturbed object (l) and that of the perturbing object (l') as follows:
| \[\frac{\mu'}{4\pi^2} \int_0^{2\pi} \int_0^{2\pi} \frac{1}{\Delta} dl dl' .\] | [46] |
By numerically evaluating the definite double integral in Eq. [46], Kozai (1962)17) discovered that as α increases, the upper limit of Θ, below which the perturbed body’s argument of pericenter is allowed to librate, also increases. In other words, when α is large, objects with an even lower orbital inclination than I ～ 39°.231 can exhibit the libration of the argument of pericenter. This circumstance is later illustrated as Fig. 7a in Appendix D.
Note that if we assume I ＞ 90° (the retrograde motion of the object), we can derive the following conclusion from Eq. [44], instead of Eq. [45], as the range of inclination that allows the libration of the argument of the pericenter to happen:
| \[I \leq \cos^{-1} \left( -\sqrt{\frac{3}{5}} \right) = 140^\circ . 768 \cdots.\] | [47] |
Later in Section 3.1, the implication that Eq. [47] turns out to be quite significant.
2.4. Numerical integration.So far, the general framework of the theory of motion of small bodies that Kozai (1962)17) presented has been described. We present a numerical example of the orbital evolution of three typical objects to examine the validity of this theory. The results are presented in the five upper rows of Fig. 5. Here, we consider CR3BP, in which the perturbing object (a proxy of Jupiter) has the same mass and semimajor axis as Jupiter in our solar system (a' ＝ 5.2042 au). We consider the inner problem and place three asteroid proxies as perturbed bodies orbiting inside the orbit of the Jupiter proxy. Each of the considered asteroid proxies is intended to represent the following actual asteroids: (4690) Strasbourg, (1373) Cincinnati, and (3040) Kozai. Then, we numerically integrate their orbital evolution over 100 kyr (105 years) in the future direction by directly integrating the equations of motion [12]. The nominal step size in the numerical integration is 1 day, and the data output interval is 100 years. For the numerical integration scheme, we employed the so-called Wisdom-Holman symplectic map (Refs. 22 and 39) implemented as the SWIFT package.40) We have added a slight improvement to this code (e.g., Refs. 41-43).
(Color online) (Upper five rows) Time evolutions of the orbital elements of the three asteroid proxies. Left: (4690) Strasbourg. Middle: (1373) Cincinnati. Right: (3040) Kozai. (The bottom row) The equipotential contours of the three asteroids are plotted in black on the (e cos g, e sin g) plane of the CR3BP framework. They are obtained by using the numerical quadrature defined in Eq. [46]. The (Θ, α) values of (4690) Strasbourg are (Θ, α) ＝ (0.9045, 0.3726). For (1373) Cincinnati (Θ, α) ＝ (0.5325, 0.6569), and for (3040) Kozai (Θ, α) ＝ (0.452, 0.354). The green curves denote the results of the direct numerical integration described in the upper five rows. In other words, the numerical results plotted in the second top row for e and in the fourth top row for g (black curves) are converted into green curves on the (e cos g, e sin g) plane in the bottom panels. In the numerical integration, the initial values of the orbital elements of the asteroids’ were taken from the JPL Horizons web interface as of June 7, 2017. The initial locations of the asteroids are as follows: (e, g) ＝ (0.1089756, 105°.515364) for Strasbourg, (e, g) ＝ (0.3151321, 99°.948105) for Cincinnati, and (e, g) ＝ (0.2005303, 288°.967682) for Kozai. See the main text for details of the Cincinnati panel’s red and blue partial circles. This figure was adopted from Ito and Ohtsuka (Ref. 3, their Figs. 2, 8 and 9) with modification.
Among the three columns of the panels in Fig. 5, the motion of (4690) Strasbourg shown in the left panels exhibits the most common behavior of the objects in the inner CR3BP. We found several noticeable characteristics in this case:
・Semimajor axis (a) remains almost constant, although it exhibits a short-term oscillation with a small amplitude.
・Eccentricity (e) and inclination (I) exhibit regular, anti-correlated oscillations. When e becomes large (or small), I becomes small (or large).
・Argument of pericenter (g) circulates in the prograde direction. The circulation period of this system correlates with the e-I oscillation.
・Longitude of ascending node (h) circulates in the retrograde direction. Its circulation period does not appear to have any particular correlation with e, I, or g.
The first property (a being almost constant) originates from the general fact that the semimajor axis of the perturbed body remains constant in the doubly averaged CR3BP. The second and third characteristics (the regular and correlated oscillations of e, I, and g) reflect the constant parameter Θ.
On the other hand, the motions of the other two asteroids shown in Fig. 5 (the middle and right columns) are qualitatively different from that of (4690) Strasbourg. Their argument of pericenter (g) librates around \(\frac{\pi}{2}\) or \(\frac{3\pi}{2}\) , instead of circulating from 0 to 2π. The oscillation still correlates with the e-I coupling. Regarding (3040) Kozai, the argument of pericenter g seems to librate around \(\frac{3\pi}{2}\) with a similar correlation.
The key to understanding the difference lies in their Θ values: Θ ～ 0.91 for (4690) Strasbourg, Θ ～ 0.55 for (1373) Cincinnati, and Θ ～ 0.45 for (3040) Kozai. In other words, as we showed in Section 2.3, the condition for Θ (Eq. [43]) determines whether the argument of the pericenter g of each object can librate or not.
Let us note a relevant but slightly different subject. In the circular restricted three-body problem at the quadrupole-level approximation, we can employ the following constant parameter44):
| \[c_2 = e^2 \left( 1 - \frac{5}{2} \sin^2 I \sin^2 g \right) .\] | [48] |
The argument of pericenter g of the perturbed object librates when c2 ＜ 0, and it circulates when c2 ＞ 0. When c2 ＝ 0, it corresponds to a special case where the perturbed body moves along a separatrix that divides the motion of the argument of pericenter into libration and circulation in phase space. Actually, parameters equivalent to c2 in Eq. [48] have been independently devised and used under different notations in later studies, such as Cse in Kinoshita and Nakai (2007),45) CK in Katz et al. (2011),46) or CKL in Antognini (2015).47) Furthermore, Vashkov’yak (2021)48) recently discussed the terminology for parameters such as Θ or c2 in this line of study. He argued that the constant parameter c2 should be referred to as the Lidov integral, as it was Lidov (1961)44) who first derived this parameter and explicitly used it in the analysis of the doubly averaged CR3BP.
2.5. Numerical quadrature.To facilitate the reader’s understanding, we briefly recap the difference between the numerical quadrature that Kozai (1962)17) dealt with and the direct numerical integration described in Section 2.4. Numerical quadrature is generally not directly related to differential equations. In general, it is an operation to find the value of a definite integral \(\int_a^b f(x) dx\) , where x can be time t or other independent variable. When the disturbing function is given by Eq. [15], its numerical quadrature corresponds to the double averaging operation given by Eq. [46]. On the other hand, the direct numerical integration of the equations of motion used in this paper corresponds to a procedure of numerically solving the ordinary differential equation \(\frac{dy}{dt}=f(y;t)\) into the form of a time-dependent solution y ＝ y(t), where t is time. Direct numerical integration is an authentic method for solving differential equations; however, it often incurs high computational costs. In this regard, the numerical quadrature is a useful and inexpensive method.
Unlike the direct numerical integration of the equations of motion, we cannot obtain information about the time variation of motion from the numerical quadrature. However, this operation allows us to draw motion trajectories in the phase space, such as in the bottom panels of Fig. 5, particularly when the system is integrable. This allows us to obtain a global picture of the motion without directly solving the equations of motion. Another important aspect is that the calculation cost of the numerical quadrature in this type of study is usually much lower than that of the numerical integration. This is because the numerical quadrature does not pursue a time-dependent solution and it is applied to systems with reduced degrees of freedom.
We now compare the results obtained through the numerical quadrature (definite integral) of the disturbing function adopted by Kozai (1962)17) with those obtained by the direct numerical integration of the equations of motion. We perform the numerical quadrature of the disturbing function along Eq. [46] and draw the equipotential contours in the bottom panels of Fig. 5. Then, we draw the result of direct numerical integration over the equipotential contours. For simplicity, one can draw equipotential contours using the analytically obtained disturbing function truncated at a finite order. However, numerical quadrature is expected to yield more accurate results than the truncated disturbing function, mainly when α is large. Kozai (1962)17) carried out the numerical quadrature under the name of “numerical harmonic analysis” (see his Fig. 1 in p. 593). To draw the equipotential contours, we use the coordinates (e cos g, e sin g) which is a variant of the Poincaré coordinates (e.g., Refs. 49 and 50).
The trajectories obtained by direct numerical integration of the equations of motion Eq. [12] are shown as green curves in each of the bottom panels of Fig. 5. The equipotential contours derived by the quadrature (black) are similar to the trajectories obtained from the direct numerical integration of the equations of motion (green). The equipotential contours of Cincinnati (bottom middle panel) possess stationary points. This is because its Θ values are smaller than the critical value \(\left(\Theta = \frac{3}{5}\right)\) in the quadrupole-level approximation. This is also the case for the asteroid Kozai (see bottom right panel). In contrast, Θ in Strasbourg is as large as 0.9. With this value, the equipotential contours (bottom left panel) do not possess any local extremums except for a minimum at the origin (0, 0). The argument of pericenter of this asteroid g circulates regularly from 0 to 2π. Kozai’s theoretical prediction is well demonstrated here.
We may find slight differences between the black and green trajectories. This discrepancy can be attributed to the approximations used in the double-averaging operation of the disturbing function, which entirely ignores mean motion resonances. By carrying out the orbital averaging procedure, the disturbing function becomes equivalent to an axisymmetric ring potential corresponding to an approximation of the perturbing planet as a ring along its orbit. As a result, any mean motion resonances (commensurability) between the perturbed and perturbing bodies will lose their effect. Although none of the three-body systems shown in Fig. 5 is trapped in a strong mean motion resonance, they may be involved with weak and high-order resonances, which the double-averaging operation neglects.
The red and blue partial circles in the bottom panel for Cincinnati in Fig. 5 represent the conditions under which the orbits of the perturbed and perturbing bodies intersect each other at the ascending node (red) and descending node (blue) of a perturbed body (e.g., Refs. 51-53). In these curves, the disturbing function becomes non-analytic and forms a set of borders on the (e cos g, e sin g) plane (note that a non-analytic function is a function that is not locally given by a convergent power series). The curves of orbit intersection for the other two asteroids are outside the panel ranges and are not shown. The black dashed circles represent the theoretically largest eccentricity (emax) of the perturbed body in each CR3BP system: \(e_\mathrm{max} = \sqrt{1-\Theta}\) which is realized when cos I ＝ 1.
2.6. About the asteroid (3040) Kozai.We note that Yoshihide Kozai’s prediction that the asteroid Cincinnati’s g librates was later confirmed by a direct numerical integration of equations of motion, including perturbations from four giant planets and Pluto (Ref. 54, his p. 210, although no figure or table was given). We also note that when Kozai published his paper in 1962, no asteroid was known to have the predicted librating behavior of orbits. Years later, James G. Williams found that the argument of perihelion of a stony asteroid 1979 BA, which is on a tilted orbit with high eccentricity that exceeds Mars’ orbit and the innermost region of the asteroid belt, is in libration. He proposed it to be named (3040) Kozai because it falls in the category that Yoshihide Kozai predicted to have the librating behavior of the orbit. See MPC (Minor Planet Circulars) 9770 (1985 July 2) or Milani et al.55) for more details.
Although the Kozai oscillation was discovered in Kozai (1962),17) it has remained of interest among a somewhat limited community of celestial mechanics for nearly three decades. However, with the rapid progress of observational astronomy, it gradually became evident that the Kozai oscillation is essential in understanding the dynamical evolution of various astronomical systems. In the following subsections, we discuss some examples. Regarding the statistics of the various terms used to represent this phenomenon and our own opinion on selecting the relevant term, we prepared the Appendix C.
3.1. Irregular satellites of giant planets.Giant planets in the solar system have many satellites. The Hill radius rH of a giant planet with mass mp orbiting at a distance ap from the Sun Msun, where the gravity of the Sun and the planet rivals, is defined as \(r_\mathrm{H} = a_\mathrm{p} \left( \frac{m_\mathrm{p}}{3M_\mathrm{sun}} \right)^\frac{1}{3}\) . Suppose r is the planetocentric distance of a satellite from its mother planet. The satellites near the giant planet \((r \lesssim 0.05 \: r_\mathrm{H})\) are generally orbiting in circular orbits in the equatorial plane of the giant planet in the same direction (prograde) as the planet’s spin. These regular satellites are considered to have formed in the protoplanetary disk, which has been fragmented and absorbed into planetesimals and planets. There are many additional irregular satellites in the outer region of the giant planet, \(r \gtrsim 0.05 \: r_\mathrm{H}\) , revolving in orbits with high eccentricity and inclinations, including some revolving in retrograde orbits. A typical theoretical work is Nesvorný et al. (2003),56), which studied the evolution of 60,000 test satellite orbits by considering the Kozai oscillation into account. They concluded that the satellite orbits highly inclined to the ecliptic are unstable due to the Kozai oscillation. This is because it radially stretches the satellite orbits and makes them either escape from the planetary Hill sphere, collide with massive inner moons, or impact the parent planet.
As of August 5, 2024, Jupiter is known to have 64 irregular satellites among the 72 confirmed satellites.2 These are probably minor bodies captured by giant planets. When the orbits of irregular satellites affected by the Kozai oscillation become highly eccentric, their pericenter may approach the planet, leading to potential collisions. Alternatively, their apocenter can be far enough to experience significant gravitational perturbations from other planets, causing them to deviate from their CR3BP orbits. There are 23 irregular satellites with an inclination angle of 40° ＜ I ＜ 140° that are probably unstable and likely to disappear in the long term due to the Kozai oscillation.
3.2. Galactic tidal force and the Oort Cloud.There is supposed to be a reservoir of comets as a spherical shell, extending to the solar system’s farthest reaches at distances between approximately 1,000 and 100,000 au, known as the Oort Cloud.57) Because the galactic plane is inclined at about ～ 60° to the ecliptic plane, the galactic tidal force influences the comet’s motion similarly to Jupiter’s gravitational influence that affects the motion of asteroids. It can be treated as a perturbation of the Keplerian motion of the cometary system (e.g., Refs. 58 and 59). The formula of the vertical component of the galactic tidal force (relative to the galactic plane) acting on cometary motion is fundamentally equivalent to what is derived from Kozai’s theory17) about the perturbation that asteroids receive from Jupiter (e.g., Ref. 60). The galactic tidal force has played a crucial role in making the outer part of the Oort Cloud dynamically isotropic, aided by random close encounters with nearby stars. The tidal influence effectively regulates the flux and the orbital distribution of the incoming comets on a time scale of 109 years (e.g., Refs. 61 and 62).
3.3. Shrinking binary orbits.The size of the main sequence stars shrank upon their birth from protostars that were much larger. The presence of main sequence stars with companion stars orbiting only a few stellar radii away, and recent discoveries of hot Jupiters imply that their binary orbits have shrunk by at least an order of magnitude after their formation. Fabrycky and Tremaine (2007)63) studied stellar binary orbit shrinkage by the Kozai oscillation due to a distant third companion on binary orbits. They demonstrated by direct numerical integration that binaries with periods as short as 0.1-10 days can be produced on a time scale of 109 years from binaries with much more extended periods of 10 to 105 days by the combined effect of the Kozai oscillation and tidal friction. This result is consistent with most short-period binaries having distant third companions. They also argued that this process explains the presence of a significant number of hot Jupiter population.
3.4. Type Ia supernovae.Type Ia supernovae are explosions caused by merging events of two white dwarfs in close binaries that result in merged mass exceeding the Chandrasekhar limit for a stable white dwarf. However, shrinking distant binary orbits close enough to merge within Hubble time has been an enigma. This is because such shrinking is feasible only for a few binaries with P ＜ 0.3 days. Thompson (2011)64) found that the Kozai oscillation in triple stars helps shrink white dwarf binaries with an orbital period P ＜ 300 days within Hubble time if the third star of the solar mass is in a high-inclination prograde orbit. This effect is even more efficient for systems with a retrograde tertiary star.
3.5. Blue stragglers.Blue stragglers in globular clusters are a group of enigmatic young stars populating above the turnoff point in the color-magnitude diagram. The theory of stellar evolution shows that massive stars leave the main sequence and move onto the giant branch in color-magnitude diagrams. Among all stars born simultaneously in a globular cluster, this process occurs from the most massive stars to the least massive stars as time passes. Therefore, the turnoff point represents the most massive stars remaining in the main sequence of a globular cluster. The blue stragglers should have different origins from single stars. They might be stars that recently gained mass from stellar collisions or a merging process/mass accretion of binary stars. Perets and Fabrycky (2009)65) suggested that the Kozai oscillation in triple stars could be the dominant process for forming close binaries and explaining many of the characteristics of blue stragglers, including binary fraction, period-eccentricity distribution, and their location above the turnoff point of their host clusters.
3.6. Orbit obliquity distribution of exoplanets.Spectroscopic time-series observations of exoplanet transits enable the identification of the transiting path, whether occulting from the approaching side or the receding side of the rotating host star surface. Albrecht et al. (2012)66) argued that the orbital rotation vectors of 14 hot Jupiters and the spin vectors of their host stars are categorized in two groups: well-aligned and misaligned from the measurements of the Rossiter-McLaughlin effect (e.g., Ref. 67). The well-aligned (low obliquity) system could result from the star-planet tidal interactions induced by the Kozai oscillation, for which the expected tidal timescale is short. On the other hand, the high-obliquity (misaligned or retrograde) systems have long tidal timescales and retain their initial random directions.
3.7. Long-term climate change on planets.As we have discussed, the Kozai oscillation produces significant variations in the eccentricity and orbital inclination of perturbed object. Its obliquity (the angle between the object’s equatorial plane and its orbital plane) also varies. If the perturbed object is a planet, particularly if it has climatic components such as atmosphere or ice sheet, this variation can change the distribution of incoming energy flux from its host star at different seasons and latitudes, causing climatic changes on the planet’s surface. This phenomenon is well known as the so-called Milankovitch Cycles for the Earth (Refs. 68-70, see Watanabe et al. (2023)71) for the recent research status). However, the influence of the Kozai oscillation on the Milankovitch Cycles on the Earth is weak because its orbit is almost circular, and its orbital inclination remains low. This phenomenon would be particularly remarkable on planets with large orbital eccentricities, whether they are inside or outside our solar system. For example, Pluto is a typical manifestation of the Kozai oscillation within our solar system with significant eccentricity and obliquity variations on timescales of several million years (e.g., Refs. 72-74). As a result, variations in the incoming solar radiation flux and its geographic distribution on Pluto induce considerable changes on its surface, together with the local accumulation or erosion of volatile ice layers.75) This kind of dynamical process will also have a significant impact on the climate of extrasolar planets with large eccentricities (e.g., Refs. 76 and 77). In recent years, studies have shown that significant variations in eccentricity and inclination that are caused by the Kozai oscillation could lead to the occurrence of a drastic bifurcation and associated hysteresis of the planetary surface environment. This type of phenomenon could have caused the snowball status on Earth (e.g., Ref. 78).
3.8. Planet around white dwarf.White dwarfs are the degenerate cores of low-mass stars that remain after the death of these stars. It was previously thought that the survival of nearby planets was unlikely. Stephan et al. (2021)79) pointed out that observations of accreting material from an ice giant planet onto a white dwarf (WD) and a Jupiter-sized planet transiting another WD both indicate that these planets are in close orbits well inside the stellar envelope radius of the host stars during their red giant phase. This implies that the planets must have migrated to their current orbits after their host stars had become WDs. Furthermore, the former WD should be hot and have a short cooling time, which indicates a fast migration mechanism. They demonstrated that the eccentric Kozai oscillation (see Section 5.1), combined with the stellar evolution onto WD and tidal effects, can produce the observed orbital configurations, assuming the presence of distant third stellar companions.
3.9. Black hole binaries Gravitational wave.Several galaxies are expected to have a central massive black hole (MBH), numerous stellar black holes (BH), and BH-BH binaries in their centers. These binaries form a hierarchical triple system with a central MBH, and gravitational perturbations from the MBH can cause high-eccentricity excitation in the BH-BH binary orbit. The pericenter distance may become sufficiently small during this process, so gravitational wave emission drives the BH-BH binary to merge. Hoang et al. (2018)80) made simulations considering the Kozai oscillation in BH-BH binaries and demonstrated that BH-BH merging can occur. Antonini et al. (2017)81) considered the formation of BH mergers through the evolution of triple stars in densely populated fields. Simulations that follow self-consistently the evolution of massive triple stars with the Kozai oscillations produced close pericenter passage binaries that could merge with gravitational wave emission.
The gravitational waves from the merging black hole binaries observed by LIGO/VIRGO showed highly asymmetric masses and significant spin that can be explained if one or both binary components are remnants of previous mergers.82) The binaries that survived the supernovae explosions generally have orbital separations too wide to merge by themselves but could merge with the aid of an external third companion, which gives rise to the Kozai oscillation.
Yoshihide Kozai (1928-2018) is a Japanese astronomer famous for his various works on the dynamics of small solar system bodies, planetary satellites, and artificial satellites around the Earth. Several oral history publications are available for Kozai’s research and life, such as DeVorkin (Ref. 83, an online publication by American Institute of Physics). More recently, Kozai summarized his academic career and personal life84) including his work on the present subject.
Kozai studied celestial mechanics under the guidance of Yusuke Hagihara (1897-1979), a member of the Japan Academy, who published nine-volume textbooks on celestial mechanics (Refs. 85-93). Kozai educated himself by reading classic textbooks on celestial mechanics (e.g., Refs. 94-98).
Kozai (1954)99) discussed the secular perturbation of asteroids by analyzing the effect of higher-order terms in the disturbing function, which marked his first paper on the subject. He wrote several papers on the orbits of Saturn’s inner satellites and obtained a doctorate from the University of Tokyo in 1958. Figure 6 is a photo of Yoshihide Kozai using a hand computer (Tiger calculator) that he extensively used for his thesis work. He took pride in his exceptional skill to turn the handle rapidly.
(Color online) Yoshihide Kozai demonstrates how to use a hand computer (Tiger calculator). The Japanese characters at the back of the instrument state that this is a property of the Tokyo Astronomical Observatory. This photo was taken in around 1955 and later given to Fumi Yoshida when he published Kozai (2016).84) Yoshida then worked at the National Astronomical Observatory of Japan. This photo was transcribed from Ref. 84 (his Fig. 1) under the permission of the publisher (Credit: Research in Astronomy and Astrophysics).
The International Geophysical Organization has proposed a campaign called “International Geophysical Year (IGY)” to organize international collaboration to measure the Earth from July 1957 through December 1958. It was also when the USSR and the USA competed to develop rocket-launching technology for defense. The USSR launched Sputnik’s first artificial satellite to orbit Earth in October 1957, followed by Explorer 1 from the USA in January 1958. Both satellites carried detectors that were useful for the IGY campaign. The rapid development of the IGY campaign created an urgent demand for theoretical researchers qualified to conduct orbital analyses of artificial satellites. In the United States, the director of the Smithsonian Astrophysical Observatory (SAO), Fred Whipple, led a project to track the orbits of artificial satellites by deploying 12 wide-field Baker-Nunn satellite tracking cameras worldwide. One of them was installed at the Tokyo Astronomical Observatory, and Kozai worked on this campaign. Whipple recognized the need to recruit specialists in satellite orbit studies and invited Yoshihide Kozai, who was relatively young then, to the SAO. Kozai joined the SAO satellite tracking group in October 1958.
Kozai wrote several influential articles during his four-year stay in the United States. In Kozai (1959),100) he proposed formulas for the perturbed satellite orbits due to non-spherical components of the Earth’s gravitational potential. He analyzed the measured motion of satellite 1958 β2, the Vanguard satellite. He found that the third-order spherical harmonics term of the Earth’s gravitational potential, corresponding to a north-south asymmetric “pear-shaped” geopotential term, can reproduce well the observed variations of orbital elements with the 80-days period.100),101)
In Kozai (1959),100) he proposed formulas to calculate the gravitational perturbations of the Moon and Sun to the artificial satellites orbiting the Earth. He used these formulas to predict that the lifetime of a satellite with a large eccentricity, 1959 Delta 2, for which NASA expected a lifetime as long as 20 years, would be as short as two years as its perigee height is quickly reduced due to solar perturbation.102) This satellite re-entered the Earth’s atmosphere after a gradual decrease in its orbital altitude about 1.96 years after launch, almost exactly as Kozai and Whitney (1959)102) predicted.3
The technology for tracking satellite orbits evolved from optical Baker-Nunn cameras in the 1950s, radio Doppler observations in the 1960s, and laser ranging in the 1970s. The satellite orbit calculation started with celestial mechanics formulas developed by Kozai and others, and was refined along with the development of computer technology during the following decades. Nevertheless, the basic principles of understanding satellite orbits remain similar to those of the classic works of celestial mechanists including Kozai.
4.2. As a community leaderIn addition to his scientific achievements, Kozai showed his eminent leadership in the astronomical community. He served as the director of the Tokyo Astronomical Observatory (TAO) of the University of Tokyo (1981-1988) and modernized its operation to be open to the academic community. In addition, he led the reorganization of TAO into an inter-university national institute. He founded the National Astronomical Observatory (NAO) and served as its first director general (1988-1994), enabling the construction of the 8.2 m Subaru Telescope on Mauna Kea. The wide-field cameras of the Subaru Telescope, Suprime-Cam104) and Hyper Suprime-Cam,105) allowed us to find many transnepturian objects that were relevant samples to verify the dynamics described in this paper (e.g., Refs. 106-108). These cameras were also successful in finding the most distant galaxies to study the early history of the universe (e.g., Ref. 109). The Subaru Telescope revolutionally enhanced the presence of Japan’s astronomical community.110) Kozai also helped establish the gravitational wave detection group in Japan.
He served as the president of the International Astronomical Union (IAU) from 1988 to 1991. He was awarded the Imperial Prize and the Japan Academy Prize in 1979 for his work on celestial mechanics. He has been a lifelong member of the Japan Academy since 1980, and served as its vice president from 2016 to 2018. The American Astronomical Society awarded him the Brouwer Award in 1989. Other honors include the Order of the Sacred Treasure, Gold and Silver Star (2002), and the Person of Cultural Merit (2009).
The study of celestial mechanics has a long history. In particular, many prominent scientists have worked on the three-body problem for hundreds of years, producing various achievements. Among them, Kozai’s work occupies an epoch-making status. However, other scientists have previously worked on similar subjects. The historical background of Kozai’s work is summarized in Appendix D. Although the three-body problem is a fundamental and classic area of research, many issues remain unsolved. Moreover, because it is fundamental, its development could have implications for mechanics and astronomy in a broader sense. As for the anti-correlated e-I oscillation discussed in this paper, there are at least two remaining theoretical challenges to be addressed. One is the effect of eccentricity of the perturber (i.e., when e' ＞ 0). The other is the influence of the mean motion resonance between the perturbed and perturbing bodies.
5.1. Effect of perturbing body's eccentricity.Unlike CR3BP where the degrees of freedom can be reduced to unity through the averaging procedure, it is not easy to reduce the degrees of freedom of the system when the eccentricity of the perturbing object is non-zero (e' ＞ 0). This is because the non-zero eccentricity of the perturber makes the disturbing potential nonaxisymmetric, and the vertical component of the perturbed body’s angular momentum \(\left( \propto \sqrt{1-e^2} \cos I \right)\) would not remain constant anymore even in the doubly averaged system. More specifically, the perturber’s argument of pericenter g' can be regarded as g' ＝ π - h in the disturbing function of the restricted three-body problem when e' ＞ 0 (Ref. 111, their Appendix A on their p. 7). Therefore, the perturbed body’s longitude of ascending node h remains in the disturbing function even after the quantity h － h' is eliminated through Jacobi’s elimination of the nodes. Hence the conjugate momentum of h, \(H = \sqrt{\mu a \left( 1-e^2 \right)} \cos I\) , could not become constant even after the double averaging operation. As a result, the perturbed body’s dynamical behavior can be significantly different when e' ＞ 0. The phenomenon is now called the eccentric Kozai oscillation. It provides answers to many questions in the solar system dynamics that the classical framework which assumes e' ＝ 0 could not answer, such as the origin of objects on retrograde orbits (e.g., Refs. 36, 38 and 46). The eccentric Kozai oscillation in the outer problem is also formulated (e.g., Refs. 112-114).
Here, we highlight another issue. Even when the orbit of the perturbing body is eccentric, its argument of pericenter g' does not appear in the disturbing function if we truncate the disturbing function to the quadrupole order (α2). In this case, the system remains integrable. For more detailed discussions, refer to Naoz et al. (2013, their Section 3)115) or Shevchenko (2017, his Section 4.4.5).116) Interestingly, Lidov and Ziglin (1976)117) referred to this situation as “a happy coincidence”. However, this coincidence no longer holds if we include the octupole-order (α3) or higher terms in the disturbing function (e.g., Refs. 111 and 118). It is worth noting the statement made by Lithwick and Naoz (2011, their footnote 2 on p. 2)111) regarding this happy coincidence: “However, it is perhaps more of an unhappy coincidence because it has misled some researchers into neglecting the role of the [perturbing] planet’s eccentricity.”
5.2. Interaction with mean motion resonance.Mean motion resonance in astronomy refers to a circumstance where two objects’ mean motion ( n and n' in the notation in this paper) is in an integer ratio, such as 1:2 or 5:3. These two objects are described as commensurable. This condition makes the configuration of their orbital motions repeat relative to each other after a certain number of orbits. Mean motion resonances, particularly lower-order ones, profoundly influence the dynamics of objects in the solar and extrasolar systems (e.g., Ref. 119).
In the present paper, we assume that no mean motion resonance is working between the orbital motions of the perturbed and perturbing bodies. However, mean motion resonances, ubiquitous in the solar system and many other dynamical systems, are a phenomenon that continues to intrigue and engage astronomers and astrophysicists. It is known that the Kozai oscillation can coexist with mean motion resonance, and they exhibit synergy to create more interesting dynamical phenomena (e.g., Refs. 120-122). A typical example is Pluto’s motion, which is in the 3:2 mean motion resonance with Neptune. There is a large body of literature on the long-term orbital motion of Pluto (e.g., Ref. 123-126), but we do not fully understand its characteristics yet. The Kozai oscillation in the outer problem is at work in this system. Still, Pluto’s inclination (～ 17°) is not high enough to trigger the oscillation in the standard theory we introduced in the present paper.127) The motion of Plutinos (a group of trans-Neptunian objects, including Pluto, in the 3:2 mean motion resonance with Neptune) exhibits even greater variety. The Kozai oscillation is also involved in many other mean motion resonances working at the outer part of the solar system and in various exoplanetary systems. Here, we uncover a wealth of potential for theoretical research in this area, inspiring and motivating future studies.
We thank the two anonymous reviewers for their detailed and critical feedback. Their suggestions significantly improved the presentation and quality of the article. This research was supported by a JSPS KAKENHI Grant (JP22K03679). The numerical calculations used to create most figures were performed at the Center for Computational Astrophysics (CfCA), National Astronomical Observatory of Japan. We employ the GNU Scientific Library (GSL) for the numerical quadrature presented in this paper, constituting a collection of routines for high-precision numerical computing. The authors acknowledge the FreeBSD Project for providing a reliable, open-source operating system that facilitated the computational tasks and data analyses required for this study. The authors used Overleaf to provide a collaborative and efficient online LaTeX environment, facilitating the manuscript’s preparation and formatting. This study used NASA’s Astrophysics Data System (ADS) Bibliographic Services.
In Section 1.3, the expansion of the disturbing function in Eq. [23] contains only even terms, 2j (j ＝ 1, 2, …). We omitted the odd terms, (2j － 1), in the expansion of Eq. [19]. This omission is justified because, in CR3BP where e' ＝ 0 (or r' ＝ a'), all odd terms of the direct part of the disturbing function R (i.e., 1/Δ in Eq. [22]) vanish after the averaging procedure. In this section, we demonstrate how these terms vanish and why our method is justified.
We start from a generalized expression of R in Eq. [19] where the disturbing function retains both even and odd terms as follows. We do not assume e' ＝ 0 yet:
| \[ R = \frac{\mu'}{r'} \sum_{j=1}^{\infty} \left( \frac{r}{r'} \right)^{j} P_{j} (\cos S) . \] | [49] |
As before, we omitted the term j ＝ 0 from the summation in Eq. [49] because it does not depend on the perturbed body’s radial distance (. We average R in Eq. [49] over the perturbing body’s mean anomaly l' (or its mean longitude λ'). This can be expressed as follows:
| \begin{aligned} \left< R \right>_{l'} &= \frac{1}{2\pi} \int^{2\pi}_0 R dl' \\ &= \frac{1}{2\pi} \int^{2\pi}_0 \frac{\mu'}{r'} \sum_{j=1}^{\infty} \left( \frac{r}{r'} \right)^{j} P_{j} (\cos S) dl' \end{aligned} | [50] |
Now, let us consider CR3BP and assume e' ＝ 0. This means that r' ＝ a' (a constant), and the mean anomaly l' of the perturbing body is equivalent to its true anomaly f' and its eccentric anomaly u'. We also assume that the perturbing body’s orbit is in the reference plane such that its inclination is I' ＝ 0.
Here, let us note the following facts:
・The perturbed body’s radial distance (r) from the central mass is independent of the perturbing body’s mean anomaly l' or its mean longitude λ'.
・The Legendre polynomials of the odd order P2j－1(x) are odd functions of x. More specifically, they are polynomials that consist of terms containing x, x3, x5, …, x2j－1. Consequently, P2j－1(cos S) turn out to be polynomials composed of terms containing cos S, cos3 S, cos5 S, …, cos2j－1 S. See Jackson (1975, his Section 3, p. 85)128) for the detailed properties of the Legendre equations and the Legendre polynomials.
・In general, cos2j－1 θ can be expressed as the sum of the terms that consist of cos (2j － 1) θ (j ＝ 1, 2, …) through the power-reduction formulae of spherical trigonometry, such as
| \begin{aligned} \cos^3 \theta &{=} \frac{ 3 \cos \theta {+} \cos 3 \theta}{ 4}, \\ \cos^5 \theta &{=} \frac{10 \cos \theta {+} 5 \cos 3 \theta {+} \cos 5 \theta}{16}, \\ \cos^7 \theta &{=} \frac{35 \cos \theta {+} 21 \cos 3 \theta {+} 7 \cos 5 \theta {+} \cos 7 \theta}{64} . \end{aligned} | [51] |
As a result, P2j－1(cos S) turn out to be functions of cos S, cos 3S, cos 5S, …, cos(2j － 1) S.
・In CR3BP, cos S in Eq. [49] can be expressed as follows (Eq. (7) of Ref. 17, note that the symbol s is equivalent to our cos S):
| \begin{aligned} \cos S &= \cos (f+g) \cos (f'+g') \\ &\hskip2em + \cos I \sin (f+g) \sin (f'+g') , \end{aligned} | [52] |
where I is the orbital inclination of the perturbed body’s orbit relative to that of the perturbing body. In CR3BP, we can regard the quantity f' + g' in Eq. [52] as being equivalent to the mean longitude of the perturbing body, λ'. Therefore, we can rewrite Eq. [52] as follows:
| \[\cos S = \cos (f+g) \cos \lambda' + \cos I \sin (f+g) \sin \lambda' .\] | [53] |
At this point, let us pick the odd terms of the averaged R from Eq. [50] when e' ＝ 0. Denoting them as \(\left< R_\mathrm{odd} \right>_{\lambda'}\) , we have:
| \begin{aligned} & \left< R_\mathrm{odd} \right>_{\lambda'}\\ &= \frac{1}{2\pi} \int^{2\pi}_0 \frac{\mu'}{a'} \sum_{j=1}^{\infty} \left( \frac{r}{a'} \right)^{2j-1} P_{2j-1} (\cos S) d\lambda' \\ &= \frac{\mu'}{2\pi a'} \int^{2\pi}_0 \left[ \frac{r}{a'} P_{1} (\cos S) + \left( \frac{r}{a'} \right)^{3} P_{3} (\cos S)\right.\\ &\qquad\left.+ \left( \frac{r}{a'} \right)^{5} P_{5} (\cos S) + \cdots \right] d\lambda' . \end{aligned} | [54] |
The disturbing function R depends on several variables. However, based on the definition of cos S in Eq. [53], we can regard that f, g, I are constant in the integral in the right-hand side of Eq. [54]. This is because all of f, g, I pertain to the perturbed body’s motion and not to that of the perturbing body. Consequently, we can regard that cos S in Eq. [53] depends solely on λ' within the integral in Eq. [54]. In this case, cos S has the periodicity of 2π. Evidently, cos 3S, cos 5S, …, cos (2j － 1)S are functions only of λ' with period of \(\frac{2\pi}{3}, \frac{2\pi}{5}, \ldots, \frac{2\pi}{2j-1}\) in this integral. This means that P2j－1(cos S) has a periodicity of \(\frac{2\pi}{2j-1}\) (j ＝ 1, 2, …) in the definite integral in Eq. [54] because P2j－1(cos S) is a function of cos S, cos 3S, cos 5S, …, cos(2j－1)S. As a result, integrating all terms on the right-hand side of Eq. [54] from 0 to 2π with respect to λ' will yield zero for any j. This gives us the conclusion that
| \[\left< R_\mathrm{odd} \right>_{\lambda'} = 0 .\] | [55] |
Equation [55] validates our earlier statement that all odd terms of R vanish after the averaging procedure in the CR3BP where e' ＝ 0.
In this appendix, we derive the functional form of the doubly averaged disturbing function at the quadrupole level ( \(W_{O\left(\alpha^2\right)}\) in Eq. [35]) used in Section 2.3.
We first carry out the averaging operation of the disturbing function R with respect to the perturbing body’s mean anomaly l', followed by averaging with respect to the perturbed body’s mean anomaly l. In CR3BP where the orbit of the perturbing body is circular (e' ＝ 0) and planar (I' ＝ 0), we cannot define its ascending node or pericenter. Therefore, the averaging operation over l' is equivalent to that over the mean longitude of the perturbing body, λ'. Since we presume that there is no mean motion resonance acting on the considered system, the distance between the central mass and the perturbed body ( and that between the central mass and the perturbing body (r') are independent. In addition, r' is assumed to be a constant; r' ＝ a'. This explains why the average of R2n in Eq. [25] over l' has the form of Eq. [26], where the operands 〈 and 〉l' do not work on the coefficients \(\frac{\mu'}{a'} \left( \frac{r}{a'} \right)^{2n}\) on the right-hand side of Eq. [26].
Now we consider the specific function form of the average of 〈P2n(cos S)〉l' that appears on the right-hand side of Eq. [26] and is defined in Eq. [27]. To proceed, we must establish two key points: the specific dependence of cos S on l' (or λ') and the specific form of each P2n (n ＝ 1, 2, …). For the function form of cos S, we consult Kozai (1962, his Eq. (7)).17) We already showed its form in Appendix A as Eq. [52] using f' + g' and also as Eq. [53] using λ'.
Many textbooks describe the characteristics of the Legendre polynomials P2n(x) in detail (e.g., Ref. 128, his Section 3.2, p. 85). When 2n ＝ 2, the function takes the following form:
| \[P_2 (x) = \frac{1}{2} \left( 3 x^2 - 1\right) ,\] | [56] |
which gives
| \[P_2 (\cos S) = \frac{1}{2} \left( 3 \cos^2 S - 1\right) .\] | [57] |
By averaging the both sides of Eq. [57] with respect to l', we obtain
| \begin{aligned} \left< P_2 (\cos S) \right>_{l'} & = \frac{1}{2} \left( 3 \left< \cos^2 S \right>_{l'} - 1 \right) \\ & = \frac{3}{2} \int_0^{2\pi} \cos^2 S d l' - \frac{1}{2} . \end{aligned} | [58] |
As mentioned above, we can replace dl' on the right-hand side of Eq. [58] with dλ' in the considered system. Thus, we can rewrite Eq. [58] as follows:
| \begin{aligned} \left< P_2 (\cos S) \right>_{\lambda'} & = \frac{1}{2} \left( 3 \left< \cos^2 S \right>_{\lambda'} - 1 \right) \\ & = \frac{3}{2} \int_0^{2\pi} \cos^2 S d \lambda' - \frac{1}{2} . \end{aligned} | [59] |
We now calculate the specific form of cos2 S using Eq. [53] and prepare for its averaging as follows
| \begin{aligned} \cos^2 S &= \left( \cos (f+g) \cos \lambda' + \cos I \sin (f+g) \sin \lambda' \right)^2 \\ &= \cos^2 (f+g) \cos^2 \lambda' \\ & \quad\quad + 2 \cos (f+g) \sin (f+g) \cos I \cos \lambda' \sin \lambda' \\ & \quad\quad\quad + \sin^2 (f+g) \cos^2 I\sin^2 \lambda' . \end{aligned} | [60] |
Using the expression of cos2 S in Eq. [60] and employing the usual trigonometric formulas of λ' such as
| \begin{aligned} &\cos^2 \lambda' = \frac{1 + \cos 2 \lambda'}{2}, \quad \sin^2 \lambda' = \frac{1 - \cos 2 \lambda'}{2}, \quad\notag\\ &\cos \lambda' \sin \lambda' = \frac{\sin 2 \lambda'}{2} , \end{aligned} | [61] |
we obtain the following expression for the average of cos2 S with respect to λ':
| \begin{aligned} & \left< \cos^2 S \right>_{\lambda'}\\ &\quad = \frac{1}{2\pi} \int_0^{2\pi} \cos^2 S d \lambda \\ &\quad = \frac{1}{2\pi} \int_0^{2\pi} \left[ \cos^2 (f+g) \frac{1+\cos 2\lambda'}{2} \right.\\ &\quad \quad \quad + \cos (f+g) \sin (f+g) \cos I \sin 2 \lambda' \\ &\quad \quad \quad \left. + \sin^2 (f+g) \cos^2 I \frac{1-\cos 2\lambda'}{2} \right] d\lambda' \\ &\quad = \frac{1}{2\pi} \int_0^{2\pi} \frac{1}{2} \left[ \cos^2 (f+g) (1+\cos 2\lambda')\right.\\ &\quad \quad \quad \left.+ \cos I \sin 2(f+g) \sin 2 \lambda' \right. \\ &\quad \quad \quad \left. + \sin^2 (f+g) \cos^2 I (1-\cos 2\lambda') \right] d\lambda' \\ &\quad = \frac{1}{2} \left[ \cos^2 (f+g) + \sin^2 (f+g) \cos^2 I \right] \\ &\quad = \frac{1}{4} \sin^2 I \cos 2 (f+g) - \frac{1}{4} \sin^2 I + \frac{1}{2} . \end{aligned} | [62] |
By substituting Eq. [62] into Eq. [59], we obtain
| \begin{aligned} & \left< P_2 (\cos S) \right>_{\lambda'}\\ &\quad = \frac{3}{2} \left< \cos^2 S \right>_{\lambda'} - \frac{1}{2} \\ &\quad = \frac{3}{8} \sin ^2 I \cos 2 (f+g) - \frac{3}{8} \sin^2 I + \frac{1}{4} . \end{aligned} | [63] |
We now perform the averaging operation of the disturbing function with respect to the perturbed body’s mean anomaly l (or by its mean longitude λ). By setting 2n ＝ 2 and substituting the result of Eq. [63] into Eq. [30], we have the following relationship:
| \begin{aligned} & \left< \left< R_{2} \right>_{\lambda'} \right>_{l}\notag\\ & = \frac{\mu'}{a'} \left( \frac{a}{a'} \right)^{2} \frac{1}{2\pi} \int^{2\pi}_0 \left( \frac{r}{a} \right)^{2} \left< P_{2} (\cos S) \right>_{\lambda'} dl\notag \\ & = \frac{\mu'}{a'} \left( \frac{a}{a'} \right)^{2} \frac{1}{2\pi} \int^{2\pi}_0 \left[ \left( \frac{r}{a} \right)^{2} \right. \notag\\ & \qquad \left. \cdot \left( \frac{3}{8} \sin ^2 I \cos 2 (f+g) - \frac{3}{8} \sin^2 I + \frac{1}{4} \right) \right] dl \notag\\ & = \frac{\mu'}{a'} \left( \frac{a}{a'} \right)^{2} \frac{1}{2\pi} \int^{2\pi}_0 \left( \frac{r}{a} \right)^{2} \frac{3}{8} \sin ^2 I \cos 2 f \cos 2 g dl \notag\\ & \quad - \frac{\mu'}{a'} \left( \frac{a}{a'} \right)^{2} \frac{1}{2\pi} \int^{2\pi}_0 \left( \frac{r}{a} \right)^{2} \frac{3}{8} \sin ^2 I \sin 2 f \sin 2 g dl \notag\\ & \quad \quad + \frac{\mu'}{a'} \left( \frac{a}{a'} \right)^{2} \frac{1}{2\pi} \int^{2\pi}_0 \left( \frac{r}{a} \right)^{2} \left( - \frac{3}{8} \sin^2 I + \frac{1}{4} \right) dl . \end{aligned} | [64] |
In what follows, we evaluate each term on the right-hand side of Eq. [64]. For this purpose, we convert the integral variable from the mean anomaly l to the eccentric anomaly u by employing the differential relationship between l and u (e.g., Ref. 1, their p. 22)
| \[\frac{dl}{du} = \frac{r}{a} .\] | [65] |
We also use the geometric relationship between the eccentric anomaly and the true anomaly given in Eq. [9] at Section 1.1, as well as the relationship between r, a, and u, which was also given in Eq. [11].
Among the three terms on the right-hand side of Eq. [64], first let us pick the term that includes cos 2f as a factor. Employing the usual trigonometric formulas such as
| \[\cos^3 u = \frac{3 \cos u + \cos 3 u}{4}, \quad \cos^2 u = \frac{1 + \cos 2 u}{2} ,\] | [66] |
the essential part of the cos 2f term can be expressed as follows:
| \begin{aligned} & \frac{1}{2\pi} \int^{2\pi}_0 \left( \frac{r}{a} \right)^{2} \cos 2 f dl \\ &= \frac{1}{2\pi} \int^{2\pi}_0 \left( \frac{r}{a} \right)^{2} \left( 2 \cos^2 f -1 \right) \frac{dl}{du} du \\ &= \frac{1}{2\pi} \int^{2\pi}_0 \left( \frac{r}{a} \right)^{2} \left[ 2 \left( \frac{\cos u - e}{1 - e \cos u} \right)^2 - 1 \right] \left( \frac{r}{a} \right) du \\ &= \frac{1}{2\pi} \int^{2\pi}_0 \left( 1 - e \cos u \right)^3 \left[ 2 \left( \frac{\cos u - e}{1 - e \cos u} \right)^2 - 1 \right] du \\ &= \frac{1}{2\pi} \int^{2\pi}_0 \left[ 2 e^2 - 1 - \left(2 e^3 + e \right) \cos u \right. \\ & \qquad\qquad\quad \left. + \left( e^2 + 2 \right) \frac{1+\cos 2u}{2} \right. \\ & \qquad\qquad\quad\quad \left. + \left( e^3 - 2e \right) \frac{3 \cos u + \cos 3u}{4} \right] du \\ &= 2 e^2 - 1 + \frac{e^2 + 2}{2} \\ &= \frac{5 e^2}{2} . \end{aligned} | [67] |
Next, we deal with the second term on the right-hand side of Eq. [64] which includes sin 2f as a factor. The essential part of this term results in the following expression:
| \begin{aligned} & \frac{1}{2\pi} \int^{2\pi}_0 \left( \frac{r}{a} \right)^{2} \sin 2 f dl \notag\\ &= \frac{1}{2\pi} \int^{2\pi}_0 \left( \frac{r}{a} \right)^{2} 2 \sin f \cos f \left( \frac{r}{a} \right) du \notag\\ &= \frac{1}{2\pi} \int^{2\pi}_0 \left[ 2 \left( \frac{r}{a} \right)^{2} \left( \frac{\sqrt{1-e^2}\sin u}{1 - e \cos u} \right) \right. \notag \\ & \left. \quad\quad\quad\quad\quad\quad \times \left( \frac{\cos u - e}{1 - e \cos u} \right) \left( \frac{r}{a} \right) \right] du \notag\\ &= \frac{1}{2\pi} \int^{2\pi}_0 \left[ 2 \left( \frac{r}{a} \right)^{2} \left( \frac{a}{r} \right)^{2} \right. \notag \\ & \left. \quad\quad \times \sqrt{1-e^2} \sin u \left( \cos u - e \right) \left( \frac{r}{a} \right) \right] du \notag \\ &= \frac{2\sqrt{1-e^2}}{2\pi} \int^{2\pi}_0 \left[ \left( 1 - e \cos u \right) \right. \notag\\ & \left. \quad\quad\quad\quad\quad \times \left( \sin u \cos u - e \sin u \right) \right] du \notag \\ &= \frac{2\sqrt{1-e^2}}{2\pi} \int^{2\pi}_0 \left( \sin u \cos u - e \sin u \right.\notag\\ &\left. \quad\quad - e \sin u \cos^2 u + e^2 \sin u \cos u \right) du \notag\\ &= \frac{2\sqrt{1-e^2}}{2\pi} \int^{2\pi}_0 \left( -\frac{5e}{4} \sin u \right.\notag\\ &\quad\quad \left. +\frac{1+e^2}{2} \sin 2u -\frac{e}{4} \sin 3u \right) du \notag\\ &= 0 . \end{aligned} | [68] |
The last term on the right-hand side of Eq. [64] does not contain f or u. The essential part of this term results in the following expression:
| \begin{aligned} & \frac{1}{2\pi} \int^{2\pi}_0 \left( \frac{r}{a} \right)^{2} dl \\ &= \frac{1}{2\pi} \int^{2\pi}_0 \left( \frac{r}{a} \right)^{2} \left( \frac{r}{a} \right) du \\ &= \frac{1}{2\pi} \int^{2\pi}_0 \left( \frac{r}{a} \right)^{3} du \\ &= \frac{1}{2\pi} \int^{2\pi}_0 \left( 1 + e \cos u \right)^{3} du \\ &= \frac{1}{2\pi} \!\! \int^{2\pi}_0 \!\!\!\! \left( 1 + 3e \cos u + 3e^2 \cos^2 u + e^3 \cos^3 u \right) du \\ &= \frac{1}{2\pi} \int^{2\pi}_0 \left( 1 + 3e \cos u + 3e^2 \frac{1+\cos 2u}{2} \right. \\ & \qquad\qquad\qquad\qquad\qquad \left. + e^3 \frac{3\cos u + \cos 3u}{4} \right) du \\ &= \frac{1}{2\pi} \int^{2\pi}_0 \left[ 1 + \frac{3e^2}{2} + \left( 3e + \frac{3e^2}{4} \right) \cos u \right. \\ & \qquad\qquad\qquad\quad\quad \left. + \frac{3e^2}{2} \cos 2u + \frac{ e^2}{4} \cos 3u \right] du \\ &= 1 + \frac{3 e^2}{2} . \end{aligned} | [69] |
Substituting the results given in Eqs. [67], [68], and [69] into the right-hand side of Eq. [64], we obtain the final expression of 〈〈R2〉λ'〉l as follows:
| \begin{aligned} & \left< \left< R_{2} \right>_{\lambda'} \right>_{l} \notag \\ &= \frac{\mu'}{a'} \left( \frac{a}{a'} \right)^{2} \! \frac{1}{2\pi} \left( \frac{3}{8} \sin ^2 I \cos 2 g \right) \int^{2\pi}_0 \!\! \left( \frac{r}{a} \right)^{2} \! \cos 2 f dl \notag\\ & -\frac{\mu'}{a'} \left( \frac{a}{a'} \right)^{2} \! \frac{1}{2\pi} \left( \frac{3}{8} \sin ^2 I \sin 2 g \right) \int^{2\pi}_0 \!\! \left( \frac{r}{a} \right)^{2} \! \sin 2 f dl\notag \\ & +\frac{\mu'}{a'} \left( \frac{a}{a'} \right)^{2} \! \frac{1}{2\pi} \left( - \frac{3}{8} \sin^2 I + \frac{1}{4} \right) \int^{2\pi}_0 \!\! \left( \frac{r}{a} \right)^{2} \! dl \notag\\ &= \frac{\mu'}{a'} \left( \frac{a}{a'} \right)^{2} \! \left[ \frac{3}{8} \sin^2 I \cos 2g \cdot \frac{5e^2}{2} \right.\notag\\ &\left. \qquad\qquad\qquad\qquad - \left( \frac{3}{8} \sin^2 I - \frac{1}{4} \right) \left( 1 + \frac{3e^2}{2} \right) \right] \notag\\ &= \frac{1}{16} \frac{\mu'}{a'} \left( \frac{a}{a'} \right)^2 \! \left[ 15 e^2 \sin^2 I \cos 2 g \right. \notag \\ & \left. \qquad\qquad\qquad\qquad -\left(3e^2+2 \right) \left( 3\sin^2 I -2 \right) \right] . \end{aligned} | [70] |
The result expressed as 〈〈R2〉λ'〉l in Eq. [70] is equivalent to \(W_{O\left(\alpha^2\right)}\) in Eq. [35], differing only in the leading sign and the mass factor \((\mu' = {\cal G} m')\) .
The higher-order terms of the disturbing function R can be derived in a similar manner for both inner and outer problems. Ito (2016)19) obtained the function forms of the doubly averaged disturbing function up to 2n ＝ 12 for the inner CR3BP and up to 2n ＝ 14 for the outer CR3BP.
This section briefly considers the most relevant term to describe the phenomenon Kozai investigated. Ito and Ohtsuka (2019, their subsection 6.2.5)3) conducted an extensive study on the choice of words describing this phenomenon in the literature. Their survey was among the abstracts of refereed and non-refereed publications in astronomy and physics registered in NASA’s Astrophysics Data System (ADS, https://ui.adsabs.harvard.edu/). They searched for the keywords, “Kozai ―” and “Lidov ―” where “―” represents one of the following terms: mechanism, resonance, cycle, oscillation, effect, dynamics, perturbation, or libration. Although it was a crude statistical method, they found a general tendency as to which term the authors in the field prefer to describe these secular dynamics. Among the words searched, the most popular terms were mechanism, resonance, cycle, oscillation, and effect listed in order of frequency. In conducting the survey, Ito and Ohtsuka3) recommended using oscillation, despite mechanism being most frequently used to express the subject in past literature. They reasoned that Kozai’s work primarily focuses on the apparent phenomenon itself (oscillation) rather than on a detailed investigation of its underlining (mechanism). The terms oscillation and cycle express a phenomenon, whereas mechanism and effect are generic words lacking a concrete image. We followed their opinion and have used oscillation in the title of this paper.
On the other hand, although the term resonance has also been commonly used (e.g., Refs. 24, 64 and 114), Ito and Ohtsuka (2019)3) highlighted that this term is not entirely appropriate for denoting this subject. This is because a resonance is supposed to work between a frequency associated with one object and a frequency associated with another object. However, the primary focus is the motion of the argument of pericenter of the perturbed body, whereas no such motion can be defined for the perturbing body’s circular orbit in CR3BP. Namely, only the frequency of a single object is involved here; therefore, it is not appropriate to call this oscillation a resonance. However, some claim that it is indeed a resonance. A reason for this is as follows. As described, the typical outcome of this oscillation is that the argument of the pericenter of the perturbed object, g, undergoes libration. This can be expressed as follows:
| \[\frac{dg}{dt} \sim 0 .\] | [71] |
This equation can be rewritten using the longitude of the ascending node h and the longitude of pericenter \(\varpi = g + h\) as follows:
| \[\frac{d\varpi}{dt} \sim \frac{dh}{dt} .\] | [72] |
Equation [72] indicates that the variation rates of \(\varpi\) and h are in the 1:1 commensurability, and one may regard that they are in a resonance.
As is well known now, an equivalent work to what Kozai (1962)17) did was published slightly earlier by a Soviet Union scientist, Michail L’vovich Lidov, in 1961. Both Lidov (1961)44) and Kozai (1962)17) worked on the inner problem of CR3BP, with neither addressing the outer problem in which the perturbed body orbits outside the perturber. Lidov (1961)44) dealt only with the quadrupole-order approximation, but he provided a detailed analysis of the system behavior at the border between the circulation and libration of g. His analysis result is summarized as a plot which is now called the Lidov diagram (e.g., Ref. 3, their Fig. 10 in their p. 35). This diagram is quite useful to understand the dynamical nature of the doubly averaged CR3BP.
Yoshihide Kozai and Michail Lidov met once in 1961 at a conference in Moscow (e.g., Ref. 129), and discussed this subject (Ref. 3, their p. 90). Kozai (1962)17) cites the presentation given by Lidov at the conference, so it would be fair to say that Kozai’s work and Lidov’s work are not independent. For this reason, the e-I oscillation that we discussed in this paper was called the Kozai-Lidov or Lidov-Kozai oscillation.
However, a recent extensive survey of literature from the late 19th to early 20th century3) revealed that pioneering work on this subject had already been established during that era. A Swedish astronomer, Edvard Hugo von Zeipel, accomplished this. von Zeipel’s achievement at that time encompasses most of the fundamental ingredients Kozai and Lidov published ～ 60 years later. Let us briefly summarize the work by von Zeipel.
・von Zeipel brought the Lindstedt series (Refs. 98, 130-137) into canonical perturbation theory. The Lindstedt series is also known as the Lindstedt-Poincaré method. In general, direct application of the classical perturbation theory to a weakly perturbed system produces artificial secular terms that tend to increase indefinitely. The Lindstedt series provides a theoretical framework for systematically suppressing the occurrence of the secular terms (e.g., Ref. 138). von Zeipel applied the Lindstedt series to the doubly averaged CR3BP to locate the possible local extremums of the disturbing function in addition to the origin (e cos g, e sin g) ＝ (0, 0). It is probably a work product of his later famous works, currently known as the von Zeipel method in the canonical perturbation theory.27)-30)
・von Zeipel dealt both with inner and outer problems because his subject was the orbital motion of long-period comets, which can take either of the two forms. Subsequent to von Zeipel, no one accomplished the detailed treatment of the outer problem until the 1990s when trans-Neptunian objects and exoplanets were discovered.
・von Zeipel considered the case of the orbit intersection between perturbed and perturbing objects. This subject was analytically considered again only at the end of the 20th century.
・von Zeipel gave general discussions on the approximation higher than the quadrupole-order (α2) and yielded solutions that are as accurate as modern studies.
・von Zeipel showed equipotential trajectories in his 1910 article. Although they appear to be handwritten, they are quite accurate.
・von Zeipel predicted the existence of stable orbits of small bodies around Neptune’s orbit. We can interpret this as a prediction of trans-Neptunian objects or Centaurs. We should remember that even Pluto was not discovered in his lifetime.
Let us now quantitatively compare von Zeipel’s theory with Kozai’s and more recent works to examine its accuracy. An example subject chosen is the dependence of the maximum (threshold) value of Θ that allows the libration of the argument of pericenter g on the ratio of the semimajor axis, α ＝ a/a'. Kozai (1962)17) estimated this dependence in the inner problem by numerical quadrature (see the third column of Kozai’s Table I in his p. 592, as well as Kozai’s Fig. 1 in his p. 593). von Zeipel (1910)139) used higher-order approximation and carried out similar calculations. von Zeipel (1910)139) used a notation k2 instead of Θ. The comparison results are presented in the panel a of Fig. 7. This tells us that von Zeipel’s k2 almost closely aligns with Kozai’s Θ up to α ＝ 0.9.
We should remark here that the dependence of Θ’s threshold value on α in the inner problem was theoretically established earlier by von Zeipel (1901).140) For the threshold value Θ or k2 in the limit of α → 0 (that is, \(\Theta = \frac{3}{5}\) ), von Zeipel (1898)141) had already found this value at the end of the nineteenth century.
(Color online) a: Dependence of Kozai’s Θ and von Zeipel’s k2 on α in the inner problem. These values denote the threshold of the libration of the perturbed object’s g. At each α, the perturbed object’s g can librate if the object’s Θ (or k2) is lower than the plotted values. The black curves with ＋ denote the numerically calculated values of Θ in Kozai (1962).17) The red open circles denote the values of k2 tabulated in von Zeipel (1910).139) b: The threshold value of the libration of the perturbed object’s g in the outer problem represented as k'22.0 (the black solid curve) and k'22.0 (the magenta solid curve) along with their dependence on α' that Ito and Ohtsuka (2019)3) calculated. The tabulated results in von Zeipel (1910)139) are converted into k'22.0 (the filled red circle) and k'20.2 (the open blue square with a central dot), and plotted them in this panel. The panel b is taken from Ito and Ohtsuka (2019, their Fig. 22a).3)
A significant feature of von Zeipel’s work is that he dealt with the outer problem with the same precision as the inner problem. Not only did Kozai and Lidov shook off to construct a theory along these lines for the outer problem, but it is even fair to say that, after von Zeipel’s work, no proper theory was published until the 1990s. From the theoretical viewpoint, the outer problem is more complex and challenging than the inner problem. For example, in the inner problem, as we have already mentioned, the equilibrium point of the argument of the pericenter g occurs only when cos g ＝ －1 (i.e., \(g=\pm \frac{\pi}{2}\) ). On the other hand, in the outer problem, the equilibrium point of g also occurs when cos g ＝ ＋1 (i.e., g ＝ 0 or g ＝ π). Therefore, the thresholds that produce the libration of g are also divided into two branches. von Zeipel (1910)139) denoted these threshold values as k'22.0 (for cos g ＝ ＋1) and k'20.2 (for cos g ＝ －1). Figure 7b shows the dependence of k'22.0 and k'20.2 on \(\alpha' (\equiv \alpha^{-1} = a'/a)\) . For comparison, we also plotted the result of a recent numerical quadrature3) for the same quantities. The panel b indicates that the result of von Zeipel agrees with the modern numerical results excellently, with only a slight deviation in the k'20.2 branch for α' ＞ 0.80.
D.2. Citation analysis.All of the works by Kozai, Lidov, and von Zeipel deal with the same subject. However, there is a significant difference in the paths these three works take toward recognition by the academic community. As a prominent example highlighting the difference, in Fig. 8 we present the annual citation counts of their works since 1960 based on data stored on NASA’s Astrophysics Data System (ADS). Note that the vertical axis of Fig. 8 has a logarithmic scale.
Time series of the citation frequencies of Kozai’s publications (cyan bars), Lidov’s publications (red bars), and von Zeipel’s publications (blue bars), between 1961 and 2024 (and partly 2025) based on the ADS as of November 25, 2024. Note that the statistics in 2025 at the bottom right part of the plot are incomplete because of the period we retrieved the citation records. For the citation data for “Kozai,” we bundle the citation frequencies of the two publications: The main paper17) and a meeting abstract by the same author with the same title published in the same issue of the same journal (abstracts of papers presented at the 111th Meeting of the American Astronomical Society at Yale University, New Haven, Connecticut, August 26-29, 1962, on p. 579 of The Astronomical Journal, 67, 1962). Similarly, for the citation data for “Lidov,” we grouped the citation frequencies of closely relevant publications: The original publication,44) its English translations,142),143) and papers on the same topic.144),145)
As we mentioned, Lidov’s achievement on doubly averaged CR3BP is equivalent to Kozai’s work. However, Fig. 8 indicates that Lidov’s work was not cited as frequently as Kozai’s. We presume this is mainly due to the disparity in the popularity of the journals that published their papers, although Lidov’s original paper written in Russian was quickly translated into English and published in two different journals.142),143) To our knowledge, relying on ADS, the first publication citing both Lidov (1962)142) and Kozai (1962)17) is Lowrey (1971).146) But for the next 33 years since Lowrey (1971),146) it seems that not only the practical equivalence of these two works but Lidov’s work has been largely forgotten. Subsequently, after the publication of a paper on the secular dynamics of irregular satellites of giant planets147) that cites Lidov (1962)142) and Kozai (1962)17) at the same time for the first time in the twenty-first century, Lidov’s work began rapidly gaining attention. Today, more and more people have come to know the equivalence of Lidov’s work and Kozai’s work.
In contrast to Kozai’s or Lidov’s work, von Zeipel’s work has been inconspicuous and perhaps even ostracized for a long time. The only citation of von Zeipel (1910)139) in the modern literature was made by Bailey and Emel’yanenko (1996)148) until 2019. However, Bailey and Emel’yanenko’s (1996)148) viewpoint was not on the phenomenon of the e-I correlated oscillation or the conditions under which argument of pericenter g of the perturbed body librates in the averaged CR3BP. Since Bailey and Emel’yanenko (1996)148) did not mention the whole picture of von Zeipel’s (1910) work either, von Zeipel’s work continued to be buried even after 1996.
Another reason may be that all of von Zeipel’s relevant publications (Refs. 139-141) were written in French. It is no exaggeration to say that English is the only language accepted in the modern academic world. Literature written in languages other than English is not considered, even if it may contain important findings. Ito and Ohtsuka (2019)3) happened to be aware of the existence of the work by von Zeipel in the process of reading Bailey and Emel’yanenko (1996).148) Then, Ito and Ohtsuka scrutinized the contents of von Zeipel (1910)139) with great interest. As a result, Ito and Ohtsuka found that the accuracy and completeness of von Zeipel’s series of works have the same merit as the later works by Lidov and Kozai, and are even better in some aspects. In other words, von Zeipel’s work is fully upward compatible with Lidov’s and Kozai’s work. Based on these facts, we assert that the prefix “von Zeipel-Lidov-Kozai” should be used to designate the theoretical framework of this line of studies.
Among the three independent works, we summarize the unique aspects of Kozai’s work as follows: Kozai used electronic computers to perform a precise numerical quadrature of the disturbing function of several actual asteroids and discovered that (1373) Cincinnati’s argument of pericenter g will librate. He continued his work along this line and later identified the locations of the stationary values of (g, e) by drawing equipotential plots of the asteroid’s disturbing Hamiltonian that are in typical mean motion resonance with Jupiter (e.g., Ref. 149, he employed numerical quadrature here again).
At the end of this section, let us note a new fact that we found while preparing this paper. First, Hagihara (1972)86) gave a detailed summary of Kozai (1962)17) in his Section 9.25 “Orbit with High Inclination and Eccentricity”. Then, we found that Hagihara (1974)89) had pointed out the identity of Kozai’s work on the doubly averaged inner CR3BP with what von Zeipel published at the end of the 19th century. Hagihara did not cite von Zeipel’s main work139) but cited the previous ones140),141) in his book. Hagihara begins his Section “14.6 von Zeipel’s theory” with the following statement:
“Consider with von Zeipel (1898; 1901) the three-body problem in which one of the masses is large compared with the other two. […]” (his p. 538)89)
Hagihara further writes:
“The discussion, and hence the theory of Lindstedt, fails […] if the orbital mutual inclination I is greater than 39°14' for Jupiter and Saturn. If the inclination is large enough, secular terms may appear that endanger the system’s stability by the presence of linear terms in t. The condition has been further studied by von Zeipel (1901), and recently by Kozai (Section 9.25).” (his p. 550)89)
The above statements indicate that Hagihara was aware of von Zeipel’s achievement in this line of study not long after the publication of Kozai (1962),17) including the threshold inclination value in the quadrupole level approximation: I ＝ 39°14'. This value is obviously derived from the theoretical framework described in Section 2.3 (see Eq. [45]). Although it seems Hagihara was not aware of the additional advanced work by von Zeipel (1910),139) the fact that Hagihara had recognized von Zeipel’s early work already in the 1970s is quite impressive.
As far as we are aware, Hagihara’s (1974)89) recognition of von Zeipel’s early work did not receive subsequent attention from the broader academic community. We browsed through Hagihara (1974)89) in the course of preparing our descriptions in Appendices A and B, when we happened to come across the above description. Hagihara’s series of treatises, including Hagihara (1974), is not available in electronic form and is currently not easily accessible to the general public. However, his works are so qualitatively detailed and quantitatively voluminous that they could be considered an encyclopedia of celestial mechanics. By reviewing Hagihara’s works and other historical works published in the early 20th century or the 19th century such as Tisserand94)-97) or Poincaré,98),150)-155) we suspect that there is still a possibility of discovering new insights into the research history of the dynamical mechanisms corresponding to the von Zeipel-Lidov-Kozai oscillation.
Contributed by Masanori IYE, M.J.A.; Edited by Takaaki KAJITA, M.J.A.
Correspondence should be addressed to: M. Iye, National Astronomical Observatory of Japan, Mitaka, Tokyo, Japan (e-mail: masanori.iye2@gmail.com).
1 For readers who understand Japanese, additional relevant articles are available (e.g., Refs. 4 and 5).
2 https://www.nao.ac.jp/new-info/satellite.html
3 This satellite is better known as Explore VI (Explore 6) or the “Paddle wheel” as Kozai and Whitney (1959) 102) mentioned in the Appendix. However, the official designation of this satellite is 1959 Delta 1, not Delta 2 (e.g., https://nssdc.gsfc.nasa.gov/nmc/spacecraft/display.action?id=1959-004A). Kozai also discussed the solar radiation pressure that works to drag down artificial satellites. 103)
Masanori Iye, born in Sapporo in 1949, graduated from the University of Tokyo in 1972 and completed his Ph.~D. in 1977 on a theoretical work to explain the origin of the spiral structures of galaxies. He became an assistant professor in the Department of Astronomy at the University of Tokyo in 1978. He extended his study at the Institute of Astronomy, University of Cambridge, UK, from 1982 to 1983. He then spent another year at the European Southern Observatory in Munich conducting observational research of galaxies using the 3.6 m telescope in Chile. Upon his return from Europe, he chaired a study group to design an 8 m class telescope and proposed the active optics system enabling a lightweight telescope. Simultaneously, he developed the first cryogenic CCD camera for astronomy and was appointed associate professor at the Tokyo Astronomical Observatory in 1986. He devoted himself to constructing the Subaru Telescope and its scientific instruments in Maunakea, Hawaii, until 2002. He then developed the laser-guide-star adaptive optics system and enhanced the telescope vision 10 times in 2011. He discovered the most distant galaxy, 13 billion light years away, using Subaru Telescope and elucidated the epoch of cosmic reionization in 2006. Since 2005, he has also promoted the Thirty Meter Telescope (TMT) project and served as the first vice-chair of the TMT International Observatory Governing Board. He received the Nishina Memorial Prize, the Toray Science and Technology Prize, the Medal with Purple Ribbon, the Order of the Sacred Treasure, the Gold and Silver Star, the Japan Academy Prize, and other honors. Iye is a Member of the Japan Academy, a professor emeritus of NAOJ, a member of the International Astronomical Union, and a Fellow of the International Society for Optics and Photonics.
Takashi Ito was born in the Yamagata prefecture, Japan, in November 1967, and graduated from the department of geophysics, faculty of science, the University of Tokyo, in March 1991. He then went on to study planetary science at graduate school, but withdrew in order to take up a position as an assistant professor at the National Astronomical Observatory of Japan (NAOJ) in April 1995. He subsequently received his Ph.D. from the University of Tokyo in March 2000 for his contributions to theoretical and numerical studies of the stability of planetary systems. Ever since taking up his position at NAOJ, he has been working on the operation and management of the open-use computing facilities of NAOJ for astrophysics. His activity encompasses research in the solar system dynamics. This includes his term as a visiting scholar at Lunar and Planetary Laboratory (LPL) of the University of Arizona, USA, from 2002 through 2004. His research in the solar system dynamics is mainly focused on analytical theory and numerical simulation. During his theoretical study, in 2019 he became aware of the equivalence between Hugo von Zeipel's 19th-century work on cometary dynamics and the ideas put forward by Yoshihide Kozai in the 1960s, as outlined in this edition of Proceedings of the Japan Academy, Series B. Ito has recently begun participating in ground-based observations of small solar system bodies. Since the spring of 2020, he has had the privilege of serving as a selected member of the observation team for NASA's New Horizons space mission. Currently, he is an associate professor (senior lecturer equivalent) at Center for Computational Astrophysics (CfCA), NAOJ. He is also a member of several academic organizations such as International Astronomical Union (IAU), Asia Oceania Geosciences Society (AOGS), Japan Geoscience Union (JpGU), Astronomical Society of Japan (ASJ), and Japanese Society for Planetary Sciences (JSPS). In recognition of his contributions to the field, asteroid 6527 Takashiito was named in his honor.
