Orbital Dynamics with the Gravitational Perturbation due to a Disk

The secular behavior of an orbit under the gravitational perturbation due to a two-dimensional uniform disk is studied in this paper, through analytical and numerical approaches. We develop the secular approximation of this problem and obtain the averaged Hamiltonian for this system first. We find that, when the ratio of the semimajor axes of the inner orbit and the disk radius takes very small values ($\ll1$), and if the inclination between the inner orbit and the disk is greater than the critical value of $30^\circ$, the inner orbit will undergo the (classical) Lidov-Kozai resonance in which variations of eccentricity and inclination are usually very large and the system has two equilibrium points at $\omega=\pi/2,3\pi/2$ ($\omega$ is the argument of perihelion). The critical value will slightly drop to about $27^\circ$ as the ratio increases to 0.4. However, the secular resonances will not occur for the outer orbit and the variations of the eccentricity and inclination are small. When the ratio of the orbit and the disk radius is nearly $1$, there are many more complicated Lidov-Kozai resonance types which lead to the orbital behaviors that are different from the classical Lidov-Kozai case. In these resonances, the system has more equilibrium points which could appear at $\omega=0,\pi/2,\pi,3\pi/2$, and even other values of $\omega$. The variations of eccentricity and inclination become relatively moderate, moreover, in some cases the orbit can be maintained at a highly inclined state. In addition, a analysis shows that a Kuzmin disk can also lead to the (classical) Lidov-Kozai resonance and the critical inclination is also $30^\circ$.


INTRODUCTION
Disk-like structures are ubiquitous in the universe, from the galactic disks in spiral galaxies to the accretion disks in active galactic nuclei (AGN) and X-ray binaries, to protoplanetary disks, debris disks, and planetary rings in the formation and evolution process of planetary systems (Sellwood 1989;Latter et al. 2018). In our solar system, all outer planets have planetary rings and the disk-shaped Kuiper belt has been observed. In addition, there is a theoretical Oort cloud considered as comet reservoir in the outer solar system, including a disk-shaped inner Oort cloud and a spherical outer Oort cloud (Oort et al. 1950;Hills 1981;Levison et al. 2001), which has not yet been confirmed by observation. Recently, Sefilian & Touma (2019) proposed that a debris disc of icy material exists outside the orbit of Neptune, with a combined mass around 10 times that of Earth, whose self-gravity could be responsible for the strange orbital architecture of Trans-Neptunian Objects (TNOs) in the outer solar system (Trujillo & Sheppard 2014;Batygin & Brown 2016).
The gravitational effect of a astronomical disk plays an important role in the formation of the dynamical architecture of various systems. Nagasawa et al. (2003) showed that the self-gravity of a dissipative protoplanetary disk could have a significant impact on the planetary eccentricity in the extrasolar multiple planetary system. In binary systems, many exoplanets have higher eccentricities than that of the planets in the solar system, some of which could be induced via the Lidov-Kozai effect (or resonance; mechanism) (Lidov 1962;Kozai 1962;Holman et al. 1997;Innanen et al. 1997;Wu & Murray 2003;Takeda & Rasio 2005). However, when the inclination angle between the protoplanetary disk plane, in which planetesimals are embedded, and the binary plane is too large, the growth of the kilometer-sized planetesimals could be inhibited due to the Lidov-Kozai effect (Marzari et al. 2009). This is a challenge for the current planetary formation theories. In order to understand planet formation in stellar binaries, many authors have investigated the influence of the protoplanetary disk gravity on planetesimal dynamics. Some found that the fast apsidal precession on planetesimal orbit induced by the gravitational effect of an axisymmetric protoplanetary disk can effectively suppress the Lidov-Kozai effect or the excitation of planetesimal eccentricity, which is conducive to the planetesimal growth resulting in the formation of planetary embryo (Batygin et al. 2011;Rafikov 2013a,b). However, it was shown that the gravity of an eccentric disk will instead excite planetesimal eccentricity to high values, leading to high impact velocities, and therefore prevent their growth (Marzari et al. 2012(Marzari et al. , 2013Lines et al. 2016). Zhao et al. (2012) studied the Lidov-Kozai effect on planetesimal dynamics with the perturbations from both the companion star and the circumprimary disk in the inclined binary system. They showed that the Lidov-Kozai effect will be similarly suppressed if the gravitational effect of disk is included, but the Lidov-Kozai effect can work at arbitrarily low inclinations in the Kozai-on region in which planetesimal eccentricities can be excited to extremely high values (∼ 1). Hence the planetesimal with very high orbital eccentricity may become a "hot planetesimal" as the shrink of the planetesimal orbit due to the gas drag damping of the gaseous disk.
Terquem & Ajmia (2010) found that the secular perturbation of an annulus disk will lead to the Lidov-Kozai effect for a planet if the planetary orbit is well inside the disk inner cavity. Furthermore, if the planetary orbit crosses the disk but most of the disk mass is beyond the orbit, they found that the oscillations of both the planetary eccentricity and inclination were not observed when the initial inclination of the orbit is below a critical value, which is significantly smaller than 39.2 • . The authors pointed out that the critical value could be 30 • in some case (case A in the paper). The analogous discussion for a three-dimensional disk and a warped disk can be found in the literature Teyssandier et al. (2012) and Terquem (2013), respectively.
In galactic dynamics, a series of papers (Vokrouhlickỳ & Karas 1998;Šubr et al. 2004;Šubr & Karas 2005;Karas &Šubr 2007) investigated the secular evolution of the stellar orbits in a galactic center surrounded by a massive accretion disk. They indicated that the stellar orbits near the black hole will undergo the Lidov-Kozai resonance under the perturbation of the accretion disk. This is helpful for the formation of highly eccentric stellar orbit in the vicinity of the black hole and increasing the star-capture rate of the black hole. Hass andŠubr also studied the Lidov-Kozai resonance in an eccentric stellar disk around a supermassive black hole (Haas &Šubr 2016;Šubr & Haas 2016).
In this paper, we aim to understand the orbital dynamics under the secular perturbation of a disk through both analytical and numerical methods. We analytically demonstrate for the first time that the gravitational perturbation from a uniform disk will induce the Lidov-Kozai effect (or resonance). In Section 2, we describe the dynamical model of a massless test particle under the disk perturbation, focusing on the multipole expansion of the disturbing function and its averaging. In Section 3, we provide a analytic study on the secular problem. Successive numerical studies are presented in Section 4. And results are summarize and discussed in Section 5.

THE DYNAMICAL MODEL
We consider a massless test particle orbiting around a central body of mass M which is surrounded by a disk of mass m d . We assume m d M , hence the motion of the particle is dominated by the central body, and the particle's orbit is a Keplerian ellipse but slightly perturbed by the gravitational potential of the disk. The Hamiltonian for this system is written as follows: where µ = GM (G is the gravitational constant), a is the semimajor axis of the particle's orbit; V is the gravitational potential of the disk. Note that the Hamiltonian has the opposite sign relative to the standard form. For simplicity, in our case the disk is considered to be a two-dimensional uniform disk. The gravitational potential exerted by the uniform disk on the particle is given by (Alberti & Vidal 2007) where R and σ are the radius and constant mass density of the disk, respectively. r is the distance between the particle and the disk's center (or the central body), ϕ is the angle between the position vector r of the particle and the disk plane.
Since we are interested in the secular behaviour of the particle's orbit under the perturbation from the disk, we would like to average the perturbing potential V (or the Hamiltonian) over the mean anomaly M of the particle's orbit, and this results in the elimination of the short-period terms in the perturbing potential. This process is know as secular approximation. Unfortunately, it is almost impossible to average Equation (2) directly and then obtain an analytical expression even for a constant mass density σ. However, when r/R < 1 (or r/R > 1), Equation (2) can be expanded in r/R (or R/r) by means of Legendre's polynomials P n . This result in r/R < 1 : where m d = σπR 2 , and r/R > 1 : The detailed derivation of Equations (3),(4) can be seen in Appendix A.
We take the disk plane to be the equatorial plane of the central body, and the inclination of the particle's orbit is measured with respect to this plane. Thus, we have where i, ω, f are the inclination, argument of perihelion, and true anomaly of the particle's orbit respectively (moreover, we use the most common variables e, Ω to denote the eccentricity and longitude of ascending node of the orbit in this paper). In fact, the key step in the process of averaging Equation (3) is to obtain the average value of the formula r| sin(f + ω)|, the average value and details are presented in Appendix B. Finally, averaging the potential in Equations (3),(4) over the mean anomaly M , we get r/R < 1 : 4 sin i π 1 + 1 2 e 2 − e 2 cos 2ω + a R 2 1 2 1 + 3 2 e 2 1 − 3 2 sin 2 i + 15 8 e 2 sin 2 i cos 2ω We drop the constant term independent of the orbital elements in the above expansion. The above expansion is at the quadrupole level of approximation. The first (second) term in the expansion is called the dipole (quadrupole), which is proportional to a/R ((a/R) 2 ). Note that Equation (6) only applies to the inner orbits whose apocenter distances are smaller than the disk radius R (geometrically, the inner orbit is located inside the sphere of radius R). Likewise r/R > 1 : Equation (7) only applies to the outer orbits whose pericenter distances are greater than the disk radius R. The outer orbit is outside the sphere of radius R. A closed form of the potential of uniform disk had been derived in Lass & Blitzer (1983), involving complete elliptic integrals of three kinds (Byrd & Friedman 1971). The closed form is numerically equivalent to the integral form in Equation (2), and Equations (3)(4) can also been obtained by expanding the closed form in the appropriate limits, but the derivations are very complicated and tedious (particularly for the case of r/R > 1). On the other hand, complete elliptic integrals as well as the potential in closed form can be computed precisely and fast, in comparison with the integral form, hence we will adopt the closed form in our full model where the potential of uniform disk is neither approximated nor averaged (see details in Section 4).

QUALITATIVE ANALYSIS OF SECULAR DYNAMICAL BEHAVIOR
In this section we star the qualitative study of the dynamics of a particle's orbit under the secular perturbation from the uniform disk. We consider two secular problems: The first is about the inner orbit and the second is about the outer orbit.

Dynamics for the inner orbit
It is convenient to understand the dynamical behavior of the orbit using the canonical Delaunay variables (Brouwer & Clemence 1961): In the inner orbit problem, the averaged Hamiltonian for the system considered here is where k = (4Gm d /πR). V di and V quad are the dipole term and the quadrupole term of the potential V in Equation (6), respectively. is a dimensionless control parameter which takes a value of 0 (1) in the dipole (quadrupole) approximation of the potential/Hamiltonian.
The averaged Hamiltonian does not depend on the mean anomaly M , nor the longitude of ascending node Ω, thuṡ L,H are constant of motion, which implies that the semimajor axis a, the z component of the angular momentum of the orbit are conserved in the secular problem. Apparently, H/L also remains constant, it follows that where J z is a constant (Kozai integral). This indicates oscillations of e and i are coupled and in antiphase. Since L,H and the Hamiltonian itself are all constant, the degree of freedom for this system is reduced to one, related to the couple (G, g). Thus the system is analytically integrable in principle.
The equations of motion about the canonical variables G, g are given bẏ When a is much smaller than R (i.e., a/R 1), the quadrupole term effect is negligible and hence the dipole approximation can describe the dynamical behaviour of the system well. Next, we first consider the above equations of motion in the dipole approximation.

Dipole approximation ( = 0)
The equilibrium point of the system satisfies the following equations In the dipole approximation ( = 0), solvingĠ = 0, then we get: ω = 0, π/2, π, 3π/2. Substitution of these values of ω into Equation (14) yieldsġ > 0, at ω = 0, π (cos 2ω = 1) Thus Equations (15) has no solution at ω = 0, π. At ω = π/2, 3π/2 (cos 2ω = −1), Equation (14) becomeṡ The above equation is expressed in terms of the Delaunay variables. As mentioned above, where L as well as H are constant. Solving the equationġ = 0, one obtains Therefore, when the above inequality is satisfied,ġ = 0 as well as Equations (15) have solutions, and the system has two equilibrium points at ω = π/2, 3π/2. This implies that the Lidov-Kozai resonance (or effect) will occur for the inner orbits of |J z | ≤ √ 3/2. When |J z | > √ 3/2, Equations (15) do not have solution and hence there is no equilibrium point for the system. The inner orbits of |J z | > √ 3/2 do not undergo any secular resonances. Apparently, the critical value of J c for the occurrence of the Lidov-Kozai resonance is √ 3/2 (for the prograde orbits), i.e. J c = √ 3/2. Figure 1 shows the trajectories in the (e, ω) phase space for the inner orbits with different values of J z . In our numerical calculations, we take such a set of dimensionless parameters: G = 1, M = 1, m d = 0.01, R = 100 (hereinafter the same). In Figure 1, all orbits have the semimajor axis a = 10. Figure 1(a)-(c) correspond to J z = 0.2, 0.5, 0.8 respectively. In Figure1(a)-(c), the Lidov-Kozai resonance occurs, and there is a stable equilibrium point at ω = π/2 surrounded by the libration island (closed trajectories). Variations of e are usually very large when the Lidov-Kozai resonance is active, and small eccentricity can even be excited to near 1 for J z = 0.2 (see Figure 1(a)). Moreover, it is predictable that variations of i are also usually very large in the Lidov-Kozai resonance, because √ 1 − e 2 cos i remains constant and e, i oscillate in anti-phase. In Figure 1(d), J z = 0.9 (> √ 3/2), the Lidov-Kozai resonance does not occur for the orbits, there is no equilibrium point and variations of e are small.
The condition making Inequality (20) to hold for any e is: In other words, as long as the inclination of the inner orbit is larger than 30 • , the orbit will undergo the Lidov-Kozai effect in which both e and i oscillate dramatically. When the inclination is below 30 • , the Lidov-Kozai effect does not work and the oscillations of both e and i are very small. Figure  The phase space trajectories in ω ∈ (π, 2π) are identical to those in (0, π), and they are not be presented here. The blue and magenta lines represent the circulating trajectories in (0, 2π) and librating trajectories around ω = π/2, respectively. Lidov-Kozai effect becomes very significant, e dramatically oscillates between 0.01 and 0.83, whereas i between 60 • and 25 • (approximately). The oscillations of both e and i are very large. The Lidov-Kozai effect in the disk problem is similar to that in the restricted three-body problem in nature. The difference is that the critical inclination in the disk problem is 30 • (in the dipole approximation), and in the restricted three-body problem is 39.2 • (Kozai 1962;Innanen et al. 1997;Naoz et al. 2013).
Following the analysis of Innanen et al. (1997) for the Lidov-Kozai effect in the restricted three-body problem, one can obtain the maximum value reached by the eccentricity and the evolution time T evol to reach the maximum value starting from a small initial eccentricity e 0 for the disk problem. We briefly present here the derivation. According to Equations (12),(13), we have where α = (k a/µ)/R. For the orbit with very small initial eccentricity e 0 and large initial inclination i 0 , the i remains almost constant before e is excited to a large value because di/dt has the factor e 2 . And when the Lidov-Kozai effect works, ω will quickly move to a value which makesω = 0, thus according to Equation (14) we have 2 sin 2 i 0 (1 − cos 2ω) = 1. Taking only the first order of e, de/dt becomes Solving Equation (23), one gets the time T evol it takes to reach e max starting from e 0 by where T is the orbital period. Note that we must have 4 sin 2 i 0 > 1, namely i 0 > 30 • , for the increase of e, which is consistent with the previous analysis (Equation (21)). If the initial inclination is smaller than 30 • , the actual growth of eccentricity is very small. For very small initial eccentricity e 0 and large initial inclination i 0 (> 30 • ), since √ 1 − e 2 cos i remains constant, the eccentricity grows from e 0 to e max simultaneously as the inclination drops from i 0 to i min , that is 1 − e 2 max cos i min = 1 − e 2 0 cos i 0 According to Equation (23), the minimum value i min the inclination can drop to is 30 • . Thus, ignoring the small quantity of e 2 0 , one obtains Two examples illustrating the values of T evol and e max predicted by Equations (24),(27) are shown in Figure 2 (see middle and right panel). Generally, when a takes small values, Equation (27) can provide rather good values for e max , but the expected time T evol given by Equation (24) is a little less than the actual time required to reach the maximum eccentricity. When a takes large values, the quadrupole term effect becomes significant, and hence Equations (24)
At ω = π/2, 3π/2, Equation (14) becomeṡ For a certain value of k 0 (or a/R), we can solve Equation (31) numerically and then obtain the values of J z which makė g = 0 has solutions at ω = π/2, 3π/2. Consequently, corresponding to these values of J z , the system has equilibrium points only at ω = π/2, 3π/2 when a/R < 0.42, and the orbits undergo the classical Lidov-Kozai resonance as shown in Figure 1. Figure 3 provides the critical value J c and the corresponding critical angle i c for occurrence of the classical Lidov-Kozai resonance as a function of a/R ranging from 0 to 0.4. If J z is smaller than the critical value J c , the Lidov-Kozai resonance occurs. If J z is larger than the critical value J c , there are no any secular resonances.
Comparison of the red line and the points in Figure 3 shows that the quadrupole approximation agrees very well with the full model for a/R below 0.4 (in the full model, the potential of the uniform disk is neither approximated nor averaged). In the dipole approximation, J c and i c do not depend on the value of a/R. But in the quadrupole approximation, J c slightly increases from √ 3/2 (≈ 0.866) to 0.896 (approximately) as a/R increases from 0 to 0.4, meanwhile, i c drops from 30 • to 26.4 • (in the full model closer to 27 • ).
When a/R > 0.42,ġ = 0 still has solutions at ω = π/2, 3π/2 for some values of J z . However, as mentioned above,ġ = 0 may also have solutions at ω = 0, π. Consequently, the equilibrium points of the system may appear at ω = 0, π/2, π, 3π/2. In this case, the phase space structure as well as the dynamical behaviors of the orbits are different from that of the classical Lidov-Kozai case of a/R < 0.42. Figure 4(a) shows a new (e, ω) phase space structure with the equilibrium points at ω = 0, π/2, π, 3π/2. One observes that the small eccentricities cannot be pumped to large values even at very high inclinations, and the corresponding inclination variations are also very small (see Figure  4(b)). Hence, in this case the small eccentricity orbits can be maintained at a highly inclined configuration. This is a significant difference from the classical Lidov-Kozai case, in which the small eccentricities with high inclinations will be excited to large values and the inclined orbits are unstable.
In fact, when a/R > 0.42, there are many other types of the Lidov-Kozai resonance. One of them has been shown in Figure 4, and more types can be seen in Section 4.

Dynamics for the outer orbit
In the outer orbit problem, the averaged Hamiltonian is Similarly, L,H and J z remain constant. AndĠ where k 1 = 15Gm d R 2 /16. Apparently, G is constant, and hence e as well as i are also constant. This implies that variations of e and i are small in the full model.ġ remains as a non-zero constant (except for H 2 /G 2 = 1/5 or i = 63.4 • ), which means the precession of ω from 0 to 2π is linear. Thus, the trajectories of the outer orbits are circulating in the phase space, and the outer orbits do not undergo secular resonances. Figure 5 shows the phase space trajectories for a = 300, J z = 0.6 in the full model, and in Figure 5 only the trajectories of the outer orbits are presented. For these outer orbits, the variations of eccentricity and inclination are small. We have run many cases of the outer orbit in the full model, and the results show that the variations of e and i are small for these outer orbits (even if a/R is only a little greater than 1). Hence the outer orbits have strong stability.
In this section we perform our numerical study based on the full model. We introduce our full model first and show the validity of secular approximation within limits by comparisons with the full model. Then we focus more on the dynamics of the orbit of a/R ∼ 1 that is difficult to investigate by analytical methods.

Full model
Consider a massless test particle moving under the gravitational field of a central body and a uniform disk, the equation of motion for the massless particle is given bÿ where r is the position vector of the particle. In the full model, the potential V is neither approximated nor averaged. In order to easily compute the acceleration vector ∇V , in the above equation of motion we adopt the closed form of the potential of uniform disk derived in Lass & Blitzer (1983), instead of the integral form in Equation (2). Accordingly, ∇V can also be written in closed form in terms of complete elliptic integrals, which can be computed precisely and easily. The detailed expressions of ∇V and its computational approaches can be found in Krogh et al. (1982) and Fukushima (2010). It needs to remark that the acceleration of the particle becomes infinite at the boundary of the uniform disk, hence we will terminate the calculations once the particle passes through the boundary. Equation (34) is integrated using the Runge-Kutta-Fehlberg 7(8) integrator. In most cases, the integrator conserves the z-component of angular momentum and the energy of the orbit within the relative error of 10 −4 .
A comparison between full model and dipole approximation on (e, ω) phase portraits for a/R = 0.1 is illustrated in Figure 6. In these phase portraits, the equilibrium points and the trajectories given in the dipole approximation are almost identical with that of the full model. This indicates that the dipole approximation can sufficiently describe the dynamical behaviour of the orbit when a/R takes a small value.  Figure 7. One observes: The quadrupole approximation agrees very well with the full model even a/R = 0.5, which suggests that the quadrupole approximation is still valid for large values of a/R. The dipole approximation and the full model are in good overall agreement up to a/R = 0.4, except for the oscillation period in the case a/R = 0.4. In the case a/R = 0.5, the dipole approximation is significantly inconsistent with the full model. The eccentricity and inclination variations in the dipole approximation are very large, and the initial small eccentricity is still excited to a large value due to the dipole-level Lidov-Kozai effect which does not depend on the value of a/R. However, as illustrated in Section 3.1.2, when a/R > 0.42 the excitation of the small eccentricity could be "suppressed" induced by the quadrupole effect. As a result, in the full model/quadrupole approximation the eccentricity and inclination variations are small for a/R = 0.5.

Dynamics for the orbit of a/R ∼ 1
We have analytically surveyed the dynamics of the inner orbit and outer orbit under the secular perturbation of a uniform disk. However, in the previous analysis, we only considered the limiting case of a/R taking small (or large) values. For the case of a/R ∼ 1, octupole and higher order terms become nonnegligible and would need to be taken into account. On the other hand, once the orbits cross the sphere surface r = R, our Equations (6)(7) break down. Thus, it is rather difficult to study the dynamical behaviour of the orbit of a/R ∼ 1 through analytical methods, we have to resort to numerical methods.
We have carried out massive numerical calculations for the case a/R ∼ 1 using the full model. And the results show that the orbits could undergo many more complicated Lidov-Kozai resonances which are different from the classical type (depicted in Figure 1). Figure 8 illustrates seven resonant phase space structures, each of which corresponds to a special Lidov-Kozai resonance type. These resonances shown in Figure 8 (a),(b),...,(g) will be called Type a,b,...,g, respectively. One observes that these resonances all have equilibrium points at ω = 0, π. In Type d and Type e there are two equilibrium points at ω = π/2, and in Type g the equilibrium point disappears at ω = π/2. In particular, in Type b,d the equilibrium points appear at other values of ω (besides ω = 0, π/2, π) and these values are not fixed. We have demonstrated that the quadrupole term effect can lead to the equilibrium points at ω = 0, π in Type a. Consequently, it is plausible to suppose that the appearances of the equilibrium points in the other types are also attributed to the high-order term effects.
In fact, which type of resonance the orbit will undergo under the uniform disk perturbation is determined by the parameters a/R and J z . We obtain the distributions of the all resonance types in the parameter space a/R × J z = (g) a = 120, Jz = 0.8 (the disk radius R is 100). The black lines correspond to librating trajectories, and the blue lines to the circulating trajectories. The equilibrium points appear at ω = 0, π/2, π, and other positions, for example, at ω = π/6, 5π/6 approximately (see panel (b)).
[0, 3] × [0, 1] through global numerical calculations (see Figure 9). In our runs, the parameter space in a/R is covered in steps of 0.01 and J z in steps of 0.01. In Figure 9, the red region gives the distribution of the classical type (as shown in Figure 1(a)), the yellow region gives the distribution of Type a, and green: Type b; purple: Type c; white: Type d; gray: Type g; brown: Type e,f. The blue region is non-resonance region where the Lidov-Kozai resonance does not occur (see Figure 1(d)). As Type e and Type f both have the same number of the stable and unstable equilibrium points, we merge them into the same distribution for simplicity. The differences between (e, ω) phase space portraits caused by varying a/R or J z in the same distribution region are only the size of libration island, the oscillation amplitudes of e and ω, and the positions of the equilibrium points. Figure 9. Distribution of all resonance types in the parameter space a/R × Jz, for the secular problem of the uniform disk perturbation. The red region represents the distribution of the classical Lidov-Kozai resonance as shown in Figure 1(a). The blue region represents the distribution of the non-resonance type shown in Figure 1(d), and the other regions, yellow: Type a; green: Type b; purple: Type c; white: Type d; brown: Type e and Type f; gray: Type g. The points corresponding to the phase space portraits shown in Figure 8 are highlighted in black stars.
As mentioned in the quadrupole approximation, when a/R > 0.42 there are such Lidov-Kozai resonances in which the system has equilibrium points at ω = 0, π. The numerical results show the critical value of a/R in the full model is about 0.43, which corresponds to the value of a/R of the boundary point between the red region and yellow region at J z = 0 in Figure 9. The analytical value of 0.42 is rather close to the numerical value of 0.43.
In the case of a/R ∼ 1, the most main resonance is Type b corresponding to the green region in Figure 9. Liking in Type a, the small eccentricities in Type b cannot grow to large values even at high inclinations, and hence the corresponding inclinations also cannot drop to very low values (as √ 1 − e 2 cos i remains constant). As a result, in Type b the small eccentricity orbits can be maintained at a highly inclined state. Overall, in the case of a/R ∼ 1, the libration islands are a little small in comparison with the case of a/R 1, the orbital resonances are not very dramatic, and the variations of the eccentricity as well as the inclination are relatively moderate.
In the region of a/R > 1, the resonances shown in Figure 8 gradually fade away as a/R increases, and the blue nonresonance region becomes larger and larger (see Figure 9), more and more orbits no longer undergo secular resonances.
Although when a/R is large, such as a/R = 3, the orbital resonances still exist for some values of J z , the libration islands are rather small in these resonances and appear at high eccentricities (see Figure 10). The most low eccentricity (e 0.6) orbits (actually the outer orbits) do not undergo secular resonance, their trajectories are circulating from 0 to 2π and the eccentricity variations are small.

CONCLUSION AND DISCUSSION
In this paper, we have studied the secular behaviour of a particle's orbit under the gravitational perturbation from a uniform disk. By averaging the multipole expansion of the disturbing potential over the orbit, we develop the secular approximation for the secular problem. We can analytically derive some properties of the system based on the secular approximation.
For the inner orbit problem, we first consider the dipole level of the secular approximation, i.e., the dipole approximation. We demonstrated that when the Kozai integral J z ≤ √ 3/2, the Lidov-Kozai resonance occurs for the inner orbits and the system has two equilibrium points at ω = π/2, 3π/2. The critical value √ 3/2 corresponds to the critical inclination 30 • above which the eccentricity and inclination variations of the orbits are usually large. The maximum eccentricity e max reached by the eccentricity depends only on the initial inclination i 0 , and e max increases as i 0 increases (for the prograde orbits). For the very large value of i 0 , the eccentricity can be excited to near 1 due to the Lidov-Kozai effect. The oscillation period or evolution time T evol ∝ 1/m d . When a/R 1, the dipole approximation agrees well with the full model. When a/R takes larger values, the dipole approximation is inadequate to describe the behaviour of the system, and hence we need to take into account the quadrupole approximation.
In the quadrupole approximation, the critical value of J c for the occurrence of the Lidov-Kozai resonance slightly increases from √ 3/2(≈ 0.866) to 0.896 as a/R increases from 0 to 0.4 due to the quadrupole effect, and the corresponding critical inclination i c drops from 30 • to about 26.4 • (the value in the full model is closer to 27 • ). When a/R > 0.42, besides ω = π/2, 3π/2, the equilibrium points of the system could also appear at ω = 0, π, which leads to the behaviours of the orbits different from that in the classical Lidov-Kozai resonance of a/R < 0.42.
For the outer orbit problem, we find that the outer orbits do not undergo the Lidov-Kozai resonance under the secular perturbation of the uniform disk. The variations of the eccentricity and inclination for the outer orbits are small, and the outer orbits have strong stability.
We investigate the secular dynamics of the orbits with a/R ∼ 1 through the numerical methods. We find that there are many more complicated Lidov-Kozai resonance types in which the equilibrium points of the system appear at ω = 0, π/2, π, 3π/2, even other values of ω. The eccentricity (as well as inclination) oscillations in these types are relatively moderate on the whole. In particular, in some resonance types the highly inclined orbits are stable.
We also find that the multipole expansion of the potential due to a Kuzmin disk (Kuzmin 1956) is similar to that of the uniform disk in form (see Appendix C). And the difference between them is only in the numerical factor which does not affect the results of qualitative analysis. This implies that the orbit under the secular perturbation of the Kuzmin disk has the similar dynamical behavior with that under the uniform disk. As a result, for the orbit located at the central region of the Kuzmin disk, the Lidov-Kozai effect kicks in for the orbit if its inclination is larger than the critical value 30 • (in the dipole approximation). In addition, Terquem & Ajmia (2010) mentioned that under the perturbation of the disk with the decreasing surface density σ(ρ) ∝ ρ −1/2 , inner radius R i = 1AU and outer radius R o = 100AU (case A in the paper), the critical inclination i c for the eccentricity growth is also about 30 • .
For the problem of the annulus disk with a non-zero inner radius R i , Terquem & Ajmia (2010) found that the Lidov-Kozai effect of i c = 39.2 • will occur for the orbit that is well inside the inner cavity of the annulus (i.e. a R i ). We have checked the behaviour of the orbit under the perturbation of the unform annulus with a small radius ratio R i /R o (below 0.5). When a R i , the orbit is indeed subject to the Lidov-Kozai effect of i c = 39.2 • . For the orbit crossing the annulus, we find that the behaviour of the orbit is similar to that in the uniform disk case. In fact, since the radius ratio R i /R o is small, to a certain extent the cumulative effect of the perturbation of the annulus is equivalent to that of the disk with R i = 0. Thus, the annulus perturbation will not significantly change the previous dynamical behavior of the orbit.

ACKNOWLEDGMENTS
We thank the anonymous referee for helpful comments and suggestions improving the paper. This work is supported by the Chinese Academy of Sciences, the National Natural Science Foundation of China (NSFC) (Nos. 11673053,11673049), and The Youth Innovation Promotion Association of Chiness Academy of Sciences (2019265).

APPENDIX
A. MULTIPOLE EXPANSION OF THE GRAVITATIONAL POTENTIAL OF UNIFORM DISK 1. When r/R < 1, it is difficult to expand the potential function in Equation (2) directly into a power series of r/R. For this we need to perform some operations on Equation (2) first. By integrating with respect to ρ, the potential function can be decomposed into: where I 1 = π 0 R 2 + r 2 − 2rR cos ϕ cos θ dθ (A2) cos θ ln R − r cos ϕ cos θ + R 2 + r 2 − 2rR cos ϕ cos θ dθ (A3) I 3 = π 0 cos θ ln r(1 − cos ϕ cos θ) dθ (A4) Now, I 1 can be expanded in term of the Legendre polynomials as follows r R n P n (cos ϕ cos θ)dθ where P n are Legendre polynomials. We take P n up to n = 2, which is I 1 = π 0 R 1 + r R 2 − 2 r R cos ϕ cos θ · 1 + r R P 1 (cos ϕ cos θ) + r R 2 P 2 (cos ϕ cos θ) + O r R 3 = Rπ 1 + 1 4 r R

C. KUZMIN DISK
Kuzmin disk is a classical model for the razor-thin disk galaxy, and it is of infinite radial extent but has finite mass as the surface density decreases fast with radius. The potential-density pairs of Kuzmin disk is given by (Binney & Tremaine 2008) where M k is the total mass of Kuzmin disk, R 2 k = x 2 + y 2 . C(> 0) is the radial scale length of the disk. Equation (C18a) can also be written as where r 2 = R 2 k + z 2 = x 2 + y 2 + z 2 and sin ϕ = |z|/r. When r/C < 1, we can expand Equation (C19) in r/C using Legendre polynomials. Here we take Legendre polynomial up to P 2 , and then we get The average of V K over the mean anomaly M is V K = GM K C a C 2 sin i π 1 + 1 2 e 2 − e 2 cos 2ω + a C 2 1 2 1 + 3 2 e 2 1 − 3 2 sin 2 i + 15 8 e 2 sin 2 i cos 2ω (C21) The above equation omits the constant term not involving the orbital elements. By comparing Equation (C21) and Equation (6), it is easy to find that they are identical in form and only different in the numerical factor.