Orbital evolution of planetesimals in gaseous disks

Planets are formed from collisional growth of small bodies in a protoplanetary disk. Bodies much larger than approximately $1$\,m are mainly controlled by the gravity of the host star and experience weak gas drag; their orbits are mainly expressed by orbital elements: semimajor axes $a$, eccentricities $e$, and inclinations $i$, which are modulated by gas drag. In a previous study, $\dot a$, $\dot e$, and $\dot i$ were analytically derived for $e \ll 1$ and $i \ll H/a$, where $H$ is the scale height of the disk. Their formulae are valid in the early stage of planet formation. However, once massive planets are formed, $e$ and $i$ increase greatly. Indeed, some small bodies in the solar system have very large $e$ and $i$. Therefore, in this paper, I analytically derive formulae for $\dot a$, $\dot e$, and $\dot i$ for $1-e^2 \ll 1$ and $i \ll H/a$ and for $i \gg H/a$. The formulae combined from these limited equations will represent the results of orbital integration unless $e \geq 1$ or $i>\pi - H/a$. Since the derived formulae are applicable for bodies not only in a protoplanetary disk but also in a circumplanetary disk, I discuss the possibility of the capture of satellites in a circumplanetary disk using the formulae.


Introduction
Planets are formed in a circumstellar disk composed of gas and solid materials (solids are of the order of 1% in mass). The solid material is initially sub-micron grains, which are controlled by an aerodynamical frictional force that is much stronger than the gravity of the central star (Adachi et al. , 1976, hereafter AHN). As solid bodies grow, gas drag becomes relatively less important. Once bodies get much larger than 1 m, they have Keplerian orbits around the central star that are slightly altered by gas drag; then, their orbits are characterized by orbital elements: semimajor axes a, eccentricities e, and inclinations i. These bodies grow via collisions, and the collisional rates are given by relative velocities determined by e and i (e.g., Inaba et al. , 2001). Damping due to gas drag and stirring by the largest body in each annulus of the disk mainly control e and i, which evolve in the protoplanetary disk during planet formation. In addition, radial drift due to gas drag, which is expressed byȧ, reduces small bodies, which stalls the growth of bodies (e.g., Kobayashi et al. , 2011Kobayashi et al. , , 2010. Therefore, the time derivative of a, e, and i (ȧ,ė, andi) caused by gas drag are very important for planet formation.
Protoplanets are formed out of collisions with kilometer-sized or larger bodies called planetesimals. While protoplanets grow, e and i of planetesimals are determined by the Hill radius of the protoplanets, and their e and i are smaller than 0.3 unless the protoplanets are greater than ten Earth masses (see equation 15 of Kobayashi et al. , 2010). Therefore, AHN derived formulae ofȧ,ė, andi due to gas drag for a body with low e < ∼ 0.3 and i ≪ 0.1. However, e and i may possibly increase when more massive planets are formed. Indeed, in the solar system, some comets, asteroids, and Kuiper belt objects have very high e and i (e.g., Kobayashi et al. , 2005). In addition, if inclined and eccentric orbits of irregular satellites around Jovian planets are originated from captures due to interaction with circumplanetary disks (e.g., Fujita et al. , 2013), these captured bodies with high e and i evolve their orbits in the disks. Therefore, analytic formulae forȧ,ė, andi for bodies with high e and i are helpful for the analysis of small bodies in the late stage of planet formation.
In this paper, I first introduce a model for gaseous disks such as protoplanetary and circumplanetary disks, and then, I revisit the derivation of Adachi et al. (1976) for the analytic formulae ofȧ,ė, andi. Next, I deriveȧ,ė, andi for bodies with high e and/or high i. By combining these limited solutions, I construct approximate formulae forȧ,ė, andi, which are applicable for all e and i unless e ≥ 1 or i > π − H/a. Lastly, I discuss the orbital evolution of satellites captured by circumplanetary disks using the derived analytic formulae forȧ,ė, andi.

Nebula disk model and gas drag law
In order to evaluate the drag force due to nebula gas, the disk model is set as follows. A gaseous disk rotates around a central object with mass M * , which is axisymmetric and in a steady state. In a cylindrical coordinate system (r, θ, z), the gas density ρ is defined from the force equilibrium in the z direction in a vertical isothermal disk as where σ(r)(= ∞ −∞ ρdz) is the surface density of the nebula disk, H(r) = √ 2c/Ω K is the disk scale height, Ω K = (GM * /r 3 ) 1/2 is the Keplerian angular velocity, and G is the gravitational constant. For simplicity, the r-dependences of σ and c are assumed as σ ∝ r −α , c ∝ r −β , respectively. This relations give ρ ∝ r −α ′ , where α ′ = α − β + 3/2. In the minimum-mass solar nebula model (Hayashi et al. , 1985), for example, α = 3/2 and β = 1/4. The angular gas velocity Ω gas is obtained from the force equilibrium in the r direction as (Tanaka et al. , 2002) In Equation (2), the terms of O(z 4 /r 4 ) and higher are ignored. This treatment is valid even for investigation of the gas drag effect on highly inclined orbits because the gas drag (and the nebula gas) is negligible at a high altitude (z ≫ H).
At the midplane of the disk, the relative velocity difference between the gas motion and the Keplerian rotation is given by For a body with mass m and radius d, gas drag force per unit mass can be written as AHN where C D is the dimensionless gas drag coefficient, u is the relative velocity vector between the body and the gas, u =| u |, and A = C D πd 2 /2m. Although C D depends on Mach number M and Reynolds number Re, C D is a constant for high Re (d > ∼ 1 km in the minimum-mass solar nebula) or for M ≫ 1 (e or i is much larger than H/a) (AHN).
General expressions for the change in a, e, and i In this paper, I investigate the time variations of semimajor axis a, eccentricity e, and inclination i of a body due to gas drag for constant C D (and then constant A). The time derivatives of a, e, and i are given by AHN as where f and ω are the true anomaly and the argument of pericenter, respectively, κ = Ω gas ρ 0 is the midplane density at r = a, and v k = (GM * /a) 1/2 is the Keplerian velocity. If the variation timescales of a, e, and i are much longer than the orbital time, the evolution of a, e, and i follows the averaged rate. The orbital averaging is taken as where T K = 2π(a 3 /GM * ) 1/2 is the Keplerian period. The same averaging is taken for e and i.
The derived formulae are in good agreement with the results of orbital integrations for e ≪ 1 and i ≪ H(a) 2 /a 2 . While Adachi et al. (1976) provided the term of i 3 inȧ, they did not take into account the vertical dependence of ρ, which includes other i 3 terms. Since the sum of these i 3 terms is negligible, I thus exclude the i 3 term derived by AHN. Inaba et al. (2001) found that the mean root squares of these limited solutions are in better agreement with the results of orbital evolution than the simple summation by Adachi et al. (1976). The averaged variation rates of a, e, and i are therefore given by where K = 2.157 and E = 1.211 are the first and second complete elliptic integrals of argument 3/4, respectively, and τ 0 = (Aρ 0 v K ) −1 is the stopping time due to gas drag for u = v K . Note that I corrected an error in the factor of the η 2 term forė in Adachi et al. (1976), which was pointed out by Kary et al. (1993).
For i = 0.01, Equations (11) to (13) are compared with the results of orbital integrations in Figure 1. These formulae are valid unless e > 0.2. Moreover, the i dependence in these formulae are valid for i < H(a)/a (see Figure 2).

Case of high eccentricity and low inclination
Here, let us consider the case where e is almost equal to unity and i is much smaller than H(a)/a. Expanding Equations (5) to (7) with respect to (1 − e 2 ) under the assumption of i ≪ H(a)/a, keeping only the lowest-order terms of (1 − e 2 ), and applying the orbital averaging such as Equation (10) to these equations, where Ψ = 1 2π 2π 0 (1 + cos f ) α−β+1 (2 − 1 + cos f ) 3 − 2 1 + cos f df, The dependences ofȧ andė on f are seen in the integral Ψ, while a term proportional to sin 2f sin 2ω ini vanishes by the orbital averaging because of an odd function of f . The integrals of Ψ, Φ 1 , and Φ 2 are functions of α − β. In the minimum-mass solar nebular model, α − β is 5/4, and then, Ψ = 0.79, Φ 1 = 0.71, and Φ 2 = −0.16.
The e dependences in these formulae are applicable for e > 0.9 as shown in Figure 1. Although the effective range of these formulae is limited, the e dependences improve the high e parts in Equations (11) to (13) as shown below.
For this derivation, Ω gas = Ω K , since the relative velocity is mainly determined by inclination.
In order to apply averaging over half an orbit around the ascending node,ȧ,ė, andi are integrated from f = −ω − π/2 to f = −ω + π/2. Sinceȧ,ė, andi are Gaussian functions as shown in Equations (20) to (22), they are negligible for large (f + ω) 2 and the integral is thus approximated to be that over interval [−∞, ∞] as follows: where H = H(a). Using this, Equations (20) to (22) are integrated around the ascending node, which results in the averaged variation rates of a, e, and i in half an orbit.
The variation rates due to the penetration near the descending node (f ≈ −ω − π) are obtained in the same way as above. Summing up the changes at two penetrations, the averaged changes are given by The validity of Equations (30) to (32) is shown in Figures 2 and 3. These formulae are applicable for i > 2H/a.

Combined equations
The variation rates of a, e, and i in two limited cases for i ≪ H/a are derived above. The formulae for low e do not well reproduce the variation rate in e ∼ 1, while high-e formulae overestimate the values for low e. Combination of low-eccentricity formulae of Equations (11) to (13) with the 1 − e 2 dependence derived in Equations (14) to (16) gives These formulae are given in a very simple way, but they are in good agreement with the results of orbital integration if i < H/2a (see Figures 1 to 3).
If H/2a < i < H/a, the variation rates of a, e, and i are obtained from combination of the low-i where MIN(D, E) is the smaller of D and E.
In conclusion, the variation rates for a, e, and i are approximately given by • Equations (30) to (32) for i > 2H/a.
In the intermediate i, the formulae tend to deviate from the right values but the accuracies are within a factor of 1.5 (see Figures 1 to 3). It should be noted that these formulae are not applicable to the case of i > π − H/a where a body experiences gas drag with relative velocity ∼ v K not only around the nodes but also for a whole orbit.

Application to captured satellites
Jovian planets have many satellites, which may be formed in circumplanetary disks. Satellites close to planets mainly have circular and coplanar orbits and may be formed in the disks. However, distant satellites tend to have inclined orbits. Here, I discuss the possibility of the capture of satellites in the disks because the formulae forȧ,ė, andi that I derive in this paper are applicable to bodies with high e and i.
Orbital evolution of bodies with high e is predicted from these analytic formulae. When a body is captured by gas drag in a circumplanetary disk, e of the captured body is approximately 1. For e > 0.9, |ė|/e and |ȧ|/a are very large. Variation rate of the pericenter distance q is much smaller than those of a and e. Indeed,q = (1 − e)ȧ − aė is estimated to be zero in Equations (14) and (15). The result is caused by the neglect of the higher terms of (1 − e 2 ), and these higher (1 − e 2 ) terms giveq/q a positive value butq/q is much smaller than |ȧ|/a and |ė|/e. Therefore, the orbital evolution occurs along with almost constant q. With decreasing e, the orbital evolution changes. Since |ȧ|/a becomes smaller than |ė|/e for e < 0.5 to 0.6, e decreases with almost constant a. Once e ≪ η,ȧ becomes dominant for the orbital evolution; the body drifts to the host planet in the timescale of τ 0 /2η 2 .
The bodies that will be satellites are temporally captured by a planet at first Suetsugu et al. , 2011), and the apocenter distances of the bodies decrease to less than the Hill radius of the host planet during the temporal capture of bodies (e.g., Fujita et al. , 2013). The change of orbital eccentricity in an orbit around the host planet is given by ∆e ≈ ė T K . The body is fully captured by gas drag if f cap ∆e ∼ 1 during the temporal capture, where f cap is the number of close encounters with the planet during the temporal capture. Using the combined formulae (Equations 30 to 38) at e = 1, ∆e is given by C 1 (i)T K (q)/τ 0 (q), where T K (q) and τ 0 (q) are T K and τ 0 at the pericenter distance q, respectively. Therefore, the necessary condition for capture is given by where the interior density of bodies, ρ d , is assumed to be 1 g cm −3 , the Hill radius of Jupiter is applied to q, and f cap is possibly approximately 100 Suetsugu et al. , 2011). As shown in Figure 4, C 1 (i) is mainly 0.1 to 10. This density is comparable to or less than the 'minimum mass subnebula' disk that contains a mass in solids equal to the mass of current Jovian satellites and gas according to the solar composition (Canup and Ward , 2002). It should be noted that the temporally captured bodies are significantly affected by the central star. However, the temporally captured bodies rotate around the host planet, which means that the perturbation by the central star is roughly canceled out in a temporally captured orbit. Therefore, the energy loss due to gas drag estimated above may lead to bound orbits.
Inclination decreases during the full capture by gas drag, which is estimated as C 2 (i) = [( i /i)(e/ ė )] e=1 in Figure 4. The initial inclination is damped during capture for 20 • < i < 30 • , while inclinations remain high after capture for other i.
However, inclinations keep decreasing due to gas drag after capture. A dissipation time of the disk, T disk , that is shorter than the damping time of inclination is thus necessary for the formation of high-inclination satellites: where M p is the host planet mass. Since the dissipation processes of circumplanetary disks are not clear yet (Fujii et al. , 2014), it is difficult to discuss the dissipation timescale. However, the dissipation timescale needed to form high-inclination satellites seems too short. Therefore, the capture of highinclination satellites might have occurred in the timescale estimated in Equation (40) before the disk dissipation and the resulting satellites tend to have retrograde orbits (see Figure 4).

Summary
I have investigated the time derivatives of orbital semimajor axis a, eccentricity e, and inclination i of a body orbiting in a gaseous disk.
• I have derivedȧ,ė, andi for e > 0.9 and i < H/2a (Equations 14 to 16) and for i > 2H/a (Equations 30 to 32). In addition, I have modified the formulae derived by AHN; Equations (11) to (13) are valid for e < 0.2 and i < H/2a, where H is the disk scale height.
• I have combined the formulae in the limited cases and have constructed approximate formulae foṙ a,ė, andi (Equations 30 to 38), which are applicable unless e ≥ 1 or i > π − H/a.
• Using these formulae, I have discussed the orbital evolution of satellites captured by a circumplanetary disk. High-inclination satellites are formed if the bodies are captured in approximately 10 4 years before the disk dissipation.