{'fraction_non_alphanumeric': 0.03852101785463236, 'fraction_numerical': 0.018854366808816407, 'mean_word_length': 4.362989025718354, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Accretion Rates of Planetesimals by Protoplanets Embedded in Nebular Gas

When protoplanets growing by accretion of planetesimals have atmospheres, small planetesimals approaching the protoplanets lose their energy by gas drag from the atmospheres, which leads them to be captured within the Hill sphere of the protoplanets. As a result, growth rates of the protoplanets are enhanced. In order to study the effect of an atmosphere on planetary growth rates, we performed numerical integration of orbits of planetesimals for a wide range of orbital elements and obtained the effective accretion rates of planetesimals onto planets that have atmospheres. Numerical results are obtained as a function of planetesimals' eccentricity, inclination, planet's radius, and non-dimensional gas-drag parameters which can be expressed by several physical quantities such as the radius of planetesimals and the mass of the protoplanet. Assuming that the radial distribution of the gas density near the surface can be approximated by a power-law, we performed analytic calculation for the loss of planetesimals' kinetic energy due to gas drag, and confirmed agreement with numerical results. We confirmed that the above approximation of the power-law density distribution is reasonable for accretion rate of protoplanets with one to ten Earth-masses, unless the size of planetesimals is too small. We also calculated the accretion rates of planetesimals averaged over a Rayleigh distribution of eccentricities and inclinations, and derived a semi-analytical formula of accretion rates, which reproduces the numerical results very well. Using the obtained expression of the accretion rate, we examined the growth of protoplanets in nebular gas. We found that the effect of atmospheric gas drag can enhance the growth rate significantly, depending on the size of planetesimals.


Introduction
Planets are thought to be formed in circum-stellar gas disks by accumulating a large number of planetesimals. When a protoplanet reaches about the Moon mass, it starts to attract surrounding gas of the disk by its gravity, which would lead to the formation of a thick atmosphere, if the growth of the protoplanet proceeds in the nebular gas. Thus, in the process of planetesimal accretion onto planets, interactions between planetesimals and atmospheres should be considered. Interactions between planetary atmospheres and planetesimals have been studied mainly in the context of giant planet formation. Podolak et al. (1988) studied the interaction of planetesimals with atmospheres in order to examine the ablation and dissolution of incoming planetesimals due to gas drag by calculating trajectories and mass-loss rates of the planetesimals. Pollack et al. (1996) and Hubickyj et al. (2005) followed the method of Podolak et al. (1988) to calculate energy and mass deposition on the atmosphere when planetesimals are passing through it as a part of their formation model of giant planets. Also, in relation to the origin of natural satellites, the interaction between primordial atmospheres and incoming objects was considered, in order to study the possibility of capture of satellites by the gas drag, as well as the orbital evolution of captured satellites (Pollack et al. 1979;Ćuk and Burns 2004).
Interactions between planetesimals and atmospheres also affect the growth rate of protoplanets. In protoplanetary disks, solid protoplanets form first, and then solid protoplanets with a sufficiently large mass (about 10 Earth masses) acquire a large amount of gas to become gas-giant protoplanets (Mizuno 1980;Bodenheimer and Pollack 1986;Pollack et al. 1996;Ikoma et al. 2000;Hubickyj et al. 2005). This scenario is widely accepted, but the growth timescale of solid planets predicted by results of theoretical studies of planetary accretion (e.g., Tanaka and Ida 1999;Kokubo and Ida 2000) are too long in two senses.
First, the predicted growth timescale is longer than the typical timescale of inward migration of protoplanets predicted by current theories including linear analyses (e.g., Ward 1986Ward , 1997Korycansky and Pollack 1993;Tanaka et al. 2002) and numerical simulations (e.g., D'Angelo et al. 2002;Masset et al. 2006), which results in the loss of forming planets into the central star before the completion of their growth. Second, planets cannot reach a critical mass that is necessary for the onset of runaway gas accretion before dissipation of disk gas, which means that gas giant planets are difficult to be formed. Current theories predict that solid planets grow in proto-planetary disks via runaway growth phase (Wetherill and Stewart 1989) followed by oligarchic growth phase (Kokubo and Ida 1996), and in both phases, protoplanets grow mainly by capturing objects that are much smaller than the protoplanets. Thus, in order to estimate the growth timescale of solid planets, one has to investigate the stage where a large number of small objects are accreted by a small number of large objects.
Collision rates of small objects onto protoplanets have been studied in detail by numerical orbital integration and analytical calculations (Nishida 1983;Wetherill and Cox 1985;Ida and Nakazawa 1989;Lissauer 1990, 1992;Dones and Tremaine 1993). The early stage of planetary accretion is likely to have proceeded in the presence of the protoplanetary gas disk, but most of the previous studies did not consider the effect of atmosphere around protoplanets, which likely enhances effective accretion rates through gas drag. Kary et al. (1993) and Kary and Lissauer (1995) studied accretion rates when planetesimals are migrating radially inward by gas drag from the nebular gas, but they did not consider the effect of planets' atmospheres.  first examined the effect of atmospheric gas drag on planetary accretion rates in detail. They showed that the presence of atmospheres around planets can enhance the accretion rate greatly depending on the size of objects that are passing through the atmospheres. The results were used in their simulation of planetary accretion in the Jovian region . Chambers (2006a,b) also used the results of  in his semi-analytic model of oligarchic growth of planets, and showed that the enhancement of the accretion rate can shorten the growth timescale so that gas giant planets can be formed before the dispersal of the gas disk, if incoming planetesimals are small ( 100m). Although  demonstrated the importance of the atmosphere for planetary accretion rates, most of their results were based on calculations neglecting the solar gravity, and three-body orbital integration was performed only in the case of initially circular orbits. In addition, they used a realistic model for the atmospheric structure rather than a simplified one, which makes the study of parameter-dependence or comparison with analytic calculations rather difficult.
In the present work, we obtain accretion rates of planetesimals onto protoplanets that have atmospheres, by using three-body orbital integration for a wide range of parameters. We use a simple power-law density distribution model for the atmosphere, which allows a systematic study of the accretion rate and its parameter-dependence. From the obtained accretion rates, we derive a semi-analytical expression of the accretion rates for planets with atmospheres and examine growth of protoplanets in the nebular gas. In §2, we first show the basic formulation and some analytic estimations which are the basis of our numerical simulation and analysis to obtain the accretion rates. In §3, we show the results of our numerical simulations to obtain the accretion rates, together with some demonstrative calculations, which are useful in understanding the basic behavior of planetesimals' motion under the gas drag of planetary atmospheres. In §4, we derive a semi-analytic formula for accretion rates under the assumption of power-law density profile for atmospheres. In §5, we examine the validity of the power-law approximation for density profile of atmospheres. Using the obtained semi-analytic formula for accretion rates, we examine the growth of protoplanets with atmospheres in §6. Our conclusions and discussion are presented in §7.

Basic equations
When the masses of planetesimals and the planet are much smaller than the solar mass and their orbital eccentricities and inclinations are sufficiently small, their motions in a rotating coordinate system can be described by linearized equations, called Hill's equations. We use a coordinate system centered on the planet, where the x-axis points radially outwards, the y-axis points in the direction of the planet's orbital motion, and the z-axis is normal to the x-y plane. We scale the distance by the mutual Hill radius r H ≡ ah with h ≡ {(M + m)/3M * } 1/3 (M and m are the masses of the planet and a planetesimal, respectively; a is the semi-major axis of the planet; and M * is the mass of the central star) and the time by Ω −1 K (Ω K is the planet's orbital angular frequency). The non-dimensional equation for the relative motion between the planet and a planetesimal can then be written as whereṽ is velocity,t is time, andΦ is the Hill potential given bỹ wherer = x 2 +ỹ 2 +z 2 . Variables with tildes denote non-dimensional quantities. The acceleration due to gas dragã drag is described bỹ where F drag = −(C D /2)πr 2 s ρ g ∆u∆u is the drag force for a planetesimal with radius r s , m is the mass of the planetesimal, C D is the gas drag coefficient of order unity, ρ g is gas density, ρ s is the internal density of planetesimals,r s is the normalized physical radius of the planetesimals, and ∆u is the velocity of the objects relative to the gas.

Atmospheric Structure
In the present work, we assume that the density structure of the atmosphere is spherically symmetric. For sufficiently massive solid protoplanets (> 10 Earth masses), runaway gas accretion occurs and the density distribution of the inflowing gas significantly deviates from spherical symmetry. However, spherically symmetric atmosphere is a good approximation for atmospheres of low-mass protoplanets because hydrostatic (or quasi-hydrostatic) equilibrium is realized in such cases (e.g., Mizuno 1980;Bodenheimer and Pollack 1986;Pollack et al. 1996;Ikoma et al. 2000). In addition, in order to study the gas drag effect systematically, we assume a power-law function for the atmospheric structure as where ρ b is the density at the bottom of the atmosphere, or equivalently, the density at the surface of the solid core,r is the distance from the center of the planet,r p is the radius of the solid part of the planet, and the exponent α is assumed to be constant.
For a purely radiative atmosphere with a constant opacity, the density distribution is given by the above power-law with α = 3 (e.g., Stevenson (1982); see Appendix A). In more realistic atmospheric structures, on the other hand, the slope of the density distribution changes near the sublimation points of ice and silicate dust, especially when opacity (or dust/gas ratio) is large (e.g., Rafikov 2006), and outside of Bondi radius where gas is no longer bounded by the planet gravity. However, the overall structure in most radial regions of atmospheres can still be well approximated by a power-law distribution. In the present work, we adopt the above simple model for atmospheric structures in order to examine the effect of atmosphere systematically for a wide range of parameters and to better understand physics behind it. Effects of a more realistic atmospheric structure based on the analytic solution described in  will be discussed in Section 5.
Once a power-law function is assumed for gas density, effects of gas drag can be described by two parameters: the exponent α, and the non-dimensional gas drag parameter ξ, which will be defined in the next subsection.

Non-dimensional gas drag parameter ξ
Assuming that density profile is expressed by Eq. (4), Eq. (3) can be rewritten in the formã where ξ is a non-dimensional parameter representing the strength of gas drag, and is defined by In the above, M ⊕ is the Earth mass, ρ core and ρ ⊕ are the internal density of planets and the Earth, respectively, and we assumed α = 3 as mentioned above. For a purely radiative atmosphere with large optical thickness, ρ b ∝ M 3 L −1 f −1 κ is a good approximation (see Appendix A or Sasaki and Nakazawa (1990)), where L is the released energy at the solid surface of the planet per unit time due to accretion of incoming planetesimals, and f κ is the opacity depletion factor, that is, f κ = 1 when dust to gas ratio is the same as the ratio in interstellar clouds assuming dust optical property is the same as that of interstellar dust; f κ decreases when dust itself is depleted or when dust opacity decreases (due to growth of dust, for example). Thus ρ b , which is a function of the planet mass, can be eliminated from the above equation. Using the analytic formula for atmospheric structures given by , we found that ρ b can be approximated by where we used accretion energy L = GMṀ /r p andṀ is the accretion rate of planetesimals onto the planet.

Normalized accretion rate
Following the definition of collision rate given by Nakazawa et al. (1989), we define the accretion rate of planetesimals per unit surface number density as whereẽ andĩ are planetesimals' orbital eccentricity and inclination divided by h, andb is the initial value of the normalized semimajor axis of an incoming planetesimal relative to the planet atỹ = ∞. p acc (ẽ,ĩ,b) is the accretion probability for a given set of the orbital elements (ẽ,ĩ,b) defined by where τ and ̟ are the initial phase angles for epicyclic motion in horizontal and vertical directions, respectively, and ϕ acc (ẽ,ĩ,b, τ, ̟) is a judgment function whether an object is accreted by the planet: unity if the object is accreted, and zero otherwise. In general, accretion can happen in two ways: (i) a planetesimal hits the planet regardless of whether it loses enough energy to become bound, or (ii) a planetesimal loses enough energy through gas drag to become gravitationally bound. We thus define the following two quantities in a similar way: P col is the normalized accretion rate due to direct collision with the solid surface of a planet without an atmosphere, which is the definition identical to that of Nakazawa et al. (1989); and P cap is the normalized rate of capture of planetesimals within the Hill sphere of the planet due to gas drag, with the assumption that the planet is a point mass (i.e., direct collision onto the planet is neglected). We also define p col and p cap in a similar way.
When a planet does not have an atmosphere and the random velocity of planetesimals is large enough, the rate of direct collision onto the solid surface of the planet with radius r p is analytically given by (Greenzweig and Lissauer 1990) where K and E are the complete elliptic integrals of the first and second kinds, respectively, with k ≡ 3ẽ 2 /4(ẽ 2 +ĩ 2 ). Collision probability p col can be written as 2.5. Analytic estimate of accretion rates due to atmospheric drag Without energy dissipation, objects coming into the Hill sphere always escape from the sphere in the end unless they collide with the planet, but they can be captured within the Hill sphere if there is energy dissipation such as atmospheric gas drag.  introduced "enhanced radius" of a planet due to atmospheric gas drag; planetesimals which come closer than this radius lose a significant amount of kinetic energy and cannot escape out of the Hill radius, thus it can be regarded as an effective radius of the planet for capturing planetesimals. In the present study, we assume a power-law density distribution for atmosphere, which enables us to obtain the expression for the enhanced radius analytically.
We first estimate energy dissipation of objects that pass through the atmosphere in a way similar to . We assume that most of energy dissipation in one encounter occurs near the point of closest approach to the planet, because the gas density increases steeply with decreasing distance from the planet. We will confirm the validity of this assumption later, using orbital integration ( §3.2). Energy dissipation due to gas drag ∆Ẽ is then estimated by the work acting on a unit mass of the object: whereã drag = |ã drag |; ∆l is the path length of the orbit near the point of closest approach and is substituted by f crmin (f c is a correction factor of order unity, andr min is the minimum distance from the origin); andṽ ∞ is the relative velocity between the planet and an approaching planetesimal before acceleration due to the planet's gravity. When the relative velocity is dominated by the random velocity rather than Kepler shear,ṽ ∞ can be approximated by the relative velocity at the origin in the unperturbed, non-gravitating solution to the Hill's equation  given as -9 - The total specific energy of a particle on the Hill coordinate system is given bỹ whereṽ is the velocity of the particle, andẼ corresponds to the Jacobi integral . Note that the orbital elementsẽ,ĩ, andb in the last expression of Eq.(14) are the initial unperturbed values, thus the term for mutual gravity (−3/r) is neglected in this expression. Since theΦ = 0 contour surface is closed, a planetesimal cannot escape from the planet's Hill sphere, if the total energyẼ becomes negative due to gas drag (Ohtsuki 1993). Therefore, the condition for capture by gas drag is given by 1 The critical value of the minimum distancer min obtained by equating both sides of Eq. (15) is denoted byR c , which can be regarded as the enhanced radius of a planet inside which incoming objects get effectively captured by the gas drag . Thus, using Eqs. (12) and (15), the enhanced radius is obtained as the solution of and is shown in Fig. 1 as a function ofṽ ∞ . Depending onṽ ∞ ,R c can be divided into three regimes and approximated bỹ In the low velocity regime, the effect of solar gravity is important. In the intermediate velocity regime, the two-body approximation is valid (e.g., Greenzweig and Lissauer 1990) andṽ ∞ is smaller than the escape velocity atr =R c as well. In the high velocity regime, v ∞ is larger than the escape velocity atr =R c .
Once we obtain the enhanced radius, the capture probability p cap can be obtained by replacingr p in Eq. (11) withR c , namely Then, the capture rate P cap can be obtained by integrating Eq.(18) with respect tob, as in the case of P col (Eq. (10)). Note that the integration cannot be done analytically becausẽ R c is not a simple function but a solution of Eq. (16), and is a function ofb throughṽ ∞ (Eq. (13)).

Numerical Method
We integrate Eq. (1) for particles with various initial orbital elements, using the eighthorder Runge-Kutta integrator. For a given (ẽ,ĩ), we calculate a phase volume for accretion in the (b, τ, ̟) phase space by finding orbits that lead to accretion. For high velocity cases, the phase volume is much smaller than the total volume, and we need to calculate a large number of orbits to obtain P acc with sufficient accuracy. We thus adopt a kind of adaptive mesh refinement method for τ and ̟ mesh; phase volume for accretion is narrowed down by several steps of calculations from a coarser mesh with a larger target to a finer mesh with a smaller target (Ohtsuki 1993). In our numerical code, the initial azimuthal distance is set toỹ 0 = 50; we confirmed that this is large enough to obtain accretion rates with sufficient accuracy. We stop our orbital integration when one of the following conditions is met: (i) r >ỹ 0 +b 0 + 2ẽ 0 +ĩ 0 (the subscript 0 denotes initial values), (ii)r <r p , (iii)Ẽ < 0. 2 Hereafter, we set α = 3 unless otherwise stated.
Gas velocity is assumed to be zero on the rotating coordinate system, i.e., ∆ũ =ṽ wherẽ v is the velocity vector of particles, since the gas velocity at a distancer from the planet is much smaller than the Keplerian velocity around the planet (= 3/r) at the location, which is the typical velocity of a planetesimal moving around the planet under the strong effect of the planet's gravity. As we consider atmospheres in hydrostatic equilibrium before runaway gas accretion, gas velocity is much smaller than the sound speed. Also, the sound speed in an atmosphere well inside of the Bondi radius is always smaller than the Keplerian velocity. Thus, the neglect of gas velocity is consistent with the assumption of atmospheres in hydrostatic equilibrium. We do not consider ablation of incoming particles when they are passing through an atmosphere. Although ablation may be important for atmospheres of massive protoplanets (≥ 10M ⊕ ; Benvenuto and Brunini 2008), its effect can be neglected for atmospheres before runaway gas accretion considered in the present work .

Effects of gas drag
Before discussing numerical results of accretion rates in detail, first we show some examples of orbital calculations, which help us understand the basic dynamical behavior of planetesimals affected by gas drag in atmospheres. Figure 2 shows examples of coplanar orbits under gas drag of atmospheres, with three different values of the gas drag parameter ξ. In the case with ξ = 1.66 × 10 −5 (left panel), the orbit is slightly affected by gas drag, but the energy dissipation was not significant and the object escapes from the Hill sphere. When ξ = 1.67 × 10 −5 (middle panel), the particle loses enough energy by gas drag to be captured within the Hill sphere. When ξ = 2 × 10 −4 , the particle loses a large amount of energy at the first encounter, and spirals toward the planet quickly.
Next, we examine the change of energy of a particle. Figure 3 shows the change ofẼ for the captured orbit described in the middle panel of Fig. 2. The left and right panels show the plots ofẼ as a function of time and the distance from the planet, respectively. We find from the right panel that the energy reduces greatly when the particle passes through the dense part of the atmosphere near the planet's surface (i.e.,r ≤ 0.1 in the case shown in Fig.  3), whereas it is almost constant when it is far from the planet. Also, the left panel shows that significant energy decrease occurs instantaneously, because most of energy dissipation occurs near pericenters, where particles move with the fastest velocity in an orbit. Note that, for the purpose of demonstration, the condition (iii) for stopping integration described above was not applied to the orbital integrations shown in Figs. 2 and 3. Figure 4 shows the total change of energyẼ due to gas drag through an encounter with the planet until one of the truncation conditions for orbital integration (i)-(iii) is met, as a function of the minimum approach distance from the planet's center. Two groups of points with ∆Ẽ 1 × 10 −9 (red pluses and green crosses) show numerical results for the low-and high-velocity cases, respectively. We can see that the energy dissipation is well described by the analytic estimates given by Eq. (12), which are shown by the two lines (we set f c = 2 in Eq. (12)). For the low velocity case, ∆Ẽ ∝r −α min , which corresponds to the second term of Eq. (12) and means that the particles are well accelerated by the planet gravity atr ∼r min . For the high velocity case, ∆Ẽ ∝r 1−α min , which corresponds to the first term of Eq. (12) and means that the particle velocity aroundr ∼r min is not much larger thanṽ ∞ . We also plot the changes ofẼ for calculations without gas drag, which corresponds to the numerical error in our numerical integration. We find that the numerical error is much smaller than the calculated energy changes due to gas drag, even in the case when gas drag is very weak (ξ = 1 × 10 −9 ).
In Fig. 5, we show the total change ofẼ in the same way as Fig. 4 in the case with r p = 0.005 to examine the dependence on ξ andẽ. Note that there are cases when the second or a later close approach reduces energy by an amount comparable to or larger than that due to the first one even when planetesimals are not captured within the Hill sphere yet, but such events are rare unless planetesimals are captured by the first close approach. Thus, the values of ∆Ẽ shown in Fig. 5 practically represent the energy change due to the first close approach to the planet, which can be used as an indicator whether planetesimals are captured or not from the point of view of energetics. The left panel shows the dependence on ξ for a low-velocity case. In this cases, we find ∆Ẽ ∝r −α min for most points, as explained above. In the weak gas drag case, planetesimals reduce their energy due to single close encounter and escape from the Hill sphere. In the strong gas drag case, points shift upward in accordance with stronger energy dissipation. The points in this case fall on a horizontal band with ∆Ẽ ∼ 2-3 atr min 0.015; in this regime, particles are captured by gas drag rather than direct collision onto the planet, i.e., orbital integration is stopped by the conditionẼ < 0. Thus, values of ∆Ẽ in this horizontal band in the plot show the amount of energy that needs to be dissipated for capture, and the upper end of this horizontal band (r min ∼ 0.015) corresponds toR c (see Eq. (17)). When gas drag is weaker, energy dissipation becomes smaller andR c becomes smaller accordingly, and eventually the trend ∆Ẽ ∝r −α min does not reach the energy needed to be captured (∆Ẽ ∼ 2-3 in this case) before it reaches the surface of the planet (r =r p ); in this case, particles are accreted by direct collision. Note that the points withr min 0.005 represent cases where orbital integrations were stopped by the condition of direct collision to the planet. The right panel shows the dependence on random velocity (i.e.,ẽ andĩ). The low velocity case is identical to the strong drag case in the left panel. With increasing random velocity, the energy needed to be captured increases because initial kinetic energy is larger, andR c decreases accordingly. In the high velocity case, the numerical results fall on a line with a nearly constant slope which extends all the way to r min =r p , which means that gas drag is not strong for particles with such a high random velocity and they are accreted only by direct collision with the planet's solid surface.
3.3. Accretion rates for the case withẽ = 2ĩ Next, we examine accretion rates in the case ofẽ = 2ĩ. Figure 6 shows the plots of p col (ẽ,ĩ,b) and p cap (ẽ,ĩ,b) as a function ofb, for two different sets of (ẽ,ĩ). Crosses and pluses represent numerical results for p col (Eq.(11)) and p cap (Eq. (18)), respectively, and the dotted and solid lines show the analytic results. Note that we empirically set f c = 1.83 for the analytic results of p cap ; changing f c only shifts the lines in the vertical direction and does not change the functional form. In the high (ẽ,ĩ) case, we find that the numerical results are in good agreement with the analytic solutions over a wide range ofb, except for the small-b region where some small deviations are observed. In the low (ẽ,ĩ) case, numerical results deviate from the analytic solution significantly. This is becauseṽ ∞ is not large enough to neglect the three-body effect in this case. Figure 7 shows numerical results of P acc for several values of ξ, along with P col (bottom line). We also plot numerical results for P col,Hill , which is the collision rate onto the Hill sphere (top line). All these curves can be roughly divided into two parts: a horizontal part in the low-velocity regime and a part with a negative slope in the high-velocity regime. In the lower eccentricity regime, previous studies showed that P col is independent ofṽ ∞ and depends onr p as Inaba et al. 2001) This corresponds to the horizontal part of the curve for P col in Fig. 7, even though the horizontal part can be seen only in the lowestẽ regime displayed in the left figure. In the high eccentricity regime, P col can be well described by Eq. (10) and is thus proportional tõ e −2 ; this corresponds to the negative-slope part of the above curve.
In the case for planets with atmospheres, this general trend still holds, but changes in some parts. In the lower eccentricity regime, P cap is expected to be proportional to ξ 1/2α , (17)). On the other hand, Fig. 7 shows a rather weak dependence on ξ in the low-velocity regime, especially when ξ ≥ 10 −5 . In order to examine the ξ-dependence in the low-velocity regime in detail, we calculated capture rates as a function of ξ forẽ =ĩ = 0 (Fig. 8). In the low ξ regime (ξ 10 −5 when α = 3), P cap is indeed ∝ ξ 1/2α . In the high ξ regime (ξ > 1), on the other hand, P cap is constant and is equal to the collision rate onto the Hill sphere, because particles entering the Hill sphere immediately get captured due to the strong gas drag. In between the above two extreme regimes of ξ in Fig. 8, there is a part where P cap is well approximated by ξ 1/10α at 10 −5 ξ 1 when α = 3. We thus obtain the following empirical formula for capture rates in the low velocity regime P cap0 = min(P cap0,lowξ , P cap0,medξ , P col0,Hill ), where P cap0,lowξ is the capture rate in the low ξ regime, P cap0,medξ is the one for the intermediate ξ regime, and P col0,Hill is the collision rate onto the Hill sphere, respectively, and these variables are defined in the case withẽ =ĩ = 0, given by P cap0,lowξ = 10.7 ξ 0.082 P cap0,medξ = P col,Hill ξ 1/10α , P col0,Hill = 4.392.
As we can see in Fig. 8, the expressions (20)-(23) reproduce our numerical results very well, even when α = 2 and 4.
In the high-velocity regime, we can see from Fig. 7 that, whenẽ ∼ 10 and ξ is small, theẽ-dependence of P cap is stronger than that of P col (∝ẽ −2 ), becauseR c decreases with increasing random velocity in this regime (see Eq. (17)). Thus, with increasingẽ, P acc for small ξ approaches P col at a certainẽ (ẽ ∼ 5 when ξ = 10 −7 andr p = 0.005, for example). In §4, we will derive an analytic expression of the capture rate in this regime. The expression shows that the capture rate is expected to be proportional toẽ −8/3 when α = 3 andẽ = 2ĩ, which agrees well with the above numerical results.

Accretion rates averaged over the distribution of eccentricities and inclinations
Planetesimals in protoplanetary disks have a distribution of eccentricities and inclinations that is well described by the Rayleigh distribution function (Ida and Makino 1992), thus we next show P acc for incoming planetesimals that have the Rayleigh distribution in eccentricity and inclination.
Normalized accretion rate for planetesimals that have a Rayleigh distribution inẽ and i with the root mean square valuesẽ * andĩ * is given by where f R is the Rayleigh distribution function given by In order to obtain P acc from Eq. (24) for certain regimes ofẽ * andĩ * (10 −1 ≤ẽ * ≤ 10 5/4 and 10 −1 ≤ĩ * ≤ 10 5/4 in this paper), we need values of P acc for a wide range ofẽ andĩ including the region outside of the regime ofẽ * andĩ * (Ohtsuki 1999). To do that, we follow the method described in Ohtsuki (1999) with a slight modification.
We first obtain P acc (ẽ,ĩ) by direct orbital integration described in the previous sections at the 10 × 10 mesh points in the square region 10 −1 ≤ẽ ≤ 10 5/4 and 10 −1 ≤ĩ ≤ 10 5/4 with even interval in the logarithmic space. We also perform orbital integration at supplementary coarse grid points in the low inclination regime, i.e., 3×2 mesh points withẽ = 10 −1 , 10 0 , 10 1 andĩ = 10 −3 , 10 −2 . Values outside of the region are obtained by extrapolation from the calculated region and by the analytic estimation for the high-velocity regime (see §4). We confirmed that P acc for 10 −1 ≤ẽ * ≤ 10 5/4 does not depend sensitively on the P acc values for the outside of the square region 10 −1 ≤ẽ ≤ 10 5/4 and 10 −1 ≤ĩ ≤ 10 5/4 . Fig. 9 shows numerical results for P acc as a function ofẽ * whenẽ * = 2ĩ * . While this figure looks similar to Fig. 7, the results shown in Fig. 9 change more smoothly with increasing random velocity than those of Fig. 7, and the absolute values in Fig. 9 tend to be slightly higher than those of Fig. 7 in the high-velocity regime. These two differences come from the average over the Rayleigh distribution, and similar results were obtained for the Rayleigh distribution average of collision rates (Greenzweig and Lissauer 1992).

A semi-analytic expression for accretion rates
Next, we derive a semi-analytical expression for P acc based on the numerical and analytical calculations for accretion rates described above. As shown in Figs. 7 and 9, with decreasing strength of gas drag (i.e., ξ), accretion rates decrease and approach the collision rates onto the solid surface of the planet without an atmosphere. Thus P acc can be approximated by the larger of P col and P cap : 3 P acc (r p , ξ, α,ẽ * ,ĩ * ) = max( P col (r p ,ẽ * ,ĩ * ), P cap (ξ, α,ẽ * ,ĩ * )).
As for the capture rate averaged over the Rayleigh distribution, we found that the following expression reproduces our numerical results very well: with η = 2/3. The accretion rate for the low-velocity regime can be simply written in terms of P cap0 (Eq.(20)) as P cap low = P cap0 , because P cap0 is independent ofẽ andĩ. We further divide P cap high into two functions as We define that P cap high1 is for a relatively low-velocity regime in the sense that the relative velocity is low enough for the collision cross-section to be enhanced by gravitational focusing, whereas it is still high enough to hold the two-body approximation. On the other hand, P cap high2 is defined for a higher velocity regime where gravitational focusing is negligible.
In the former regime, P cap,high1 can be obtained by the integration of Eq. (18) as in Eq. (8): where γ ≡ (α + 2)/2α, and R(k, γ) is an integral similar to the complete elliptic integral, and is defined by The Rayleigh distribution average of P cap,high1 can be calculated as where Γ is the gamma function, 2 F 1 is the hypergeometric function, and we assumeẽ * >ĩ * and α > 1. This gives P cap high1 ∝ (ẽ * ) −(2γ+1) = (ẽ * ) −2(α+1)/α whenẽ * = 2ĩ * , whose dependence onẽ * (i.e., ∝ (ẽ * ) −2(α+1)/α ) agrees with that of P cap,high1 onẽ whenẽ = 2ĩ. For the typical case of α = 3 andẽ * = 2ĩ * , we have P cap high1 ≃ 1.52 × 10 2 (f c ξ) 1/3 (ẽ * ) −8/3 . In deriving Eq. (38), we treated R(k, γ) as a constant 4 , because R(k, γ) is a slowly varying function of k as long as k is not close to unity (note that 0 ≤ k ≤ 3/4 from its definition). For the still higher velocity regime, thus The curves in Fig. 9 show P acc (r p , ξ, α,ẽ * ,ĩ * ) obtained using the above semi-analytic expressions. We can see that P acc is well approximated by the above semi-analytical expressions within a factor of two for the entire regions of random velocity and ξ. Whenr p is large, P acc tends to approach P col because planetesimals tend to collide with the planet more easily before losing a sufficient amount of energy by gas drag to get captured. For example, whenr p = 0.005 (left panel) and ξ = 10 −9 , P acc = P col , which can be deduced by Fig. 1. On the other hand, P acc can be approximated by P cap when ξ is large. In summary, the effect of atmosphere is statistically negligible in the entire velocity regime when ξ is so small that P col high > P cap high1 is satisfied, where P col high and P cap high1 are given by Eqs.(30) and (38), respectively. This condition can be written in terms of ξ as ξ < 4 × 10 −9 r p 0.005 3 . ( To derive the condition, we setẽ * = 2 as the lowest limit of the range for the functions. It should be noted that when the mass of planetesimals are small ( 10 14 g), the assumption ofẽ * = 2ĩ * breaks down if protoplanetary gas disk still exists and runaway bodies significantly perturb the orbit of such small bodies (Ohtsuki et al. 2002). This is because, once inclination is damped toĩ ≤ 1, the stirring rate of inclination by the large bodies is much smaller than that of eccentricity in such a velocity regime, and only eccentricity tends to be enhanced (Ida 1990).
Using analytic calculation of planetesimals' energy dissipation due to gas drag under the two-body approximation neglecting the solar gravity,  derived a relation between the size of planetesimals and the enhanced radius of protoplanets, for a given atmospheric structure (their Eq. (17)). Substituting this enhanced radius for the planetary radius in the expressions for collision rates derived by Inaba et al. (2001), semianalytic expressions of accretion rates for protoplanets with atmospheres can be obtained (Inaba and Ikoma 2003, their Eqs.(20)- (24)). We find that P acc with the effect of atmospheric gas drag obtained from these expressions derived by  ( P col in their notation) is smaller by a factor of about two as compared to ours in the high-velocity regime, even when the same atmospheric structure is assumed (Fig.10). This discrepancy is caused by the more simplified procedures in deriving the expression in  as compared to ours; the cause for the discrepancy is further discussed in detail in Appendix B. On the other hand,  obtained the enhanced radius using more realistic atmospheric structures, while we assumed a power-law density distribution. Effects of this assumption will be discussed in Section 5.

Comparison with the case of more realistic atmospheric structures
Here, we examine the validity of our assumption of the simplified atmospheric structure, i.e., power-law function (Eq.(4)). As shown above, the assumption of power-law approximation for gas density in atmospheres allowed a systematic study and understanding of the dependence of the accretion rate on various physical parameters. In order to check the validity of this assumption, we compared the results shown previously in this paper with those obtained using a more realistic atmospheric structure based on the analytic solution shown by . Although the index α should be properly chosen depending on the parameter values for a given atmospheric structure, we found that the use of power law functions is a good approximation in most cases to describe realistic atmospheres when opacity depletion factor f κ is much smaller than unity (f κ 10 −2 ), because optically thin atmospheres tend to be radiative rather than convective.
The assumption of small f κ would be reasonable, because a large fraction of dust in protoplanetary disks in this stage would be depleted after the accretion of large bodies such as planetesimals or protoplanets. In addition, small f κ also corresponds to the situation when dust size is much larger than submicron size, and this would be the case for large protoplanets (M M ⊕ ), because dust in the atmospheres grows quickly to reduce opacity by a factor of 100 from that of sub-micron size dust (Movshovitz and Podolak 2008).
However, the power-law distribution would not be a good approximation at outer region of the atmospheres in some cases. The density profile fitted by a power-law function at the bottom of an atmosphere tends to underestimate the gas density at the outer region when planet mass is small (M ∼ 0.1M ⊕ ) and overestimate it when planet mass is large (M ∼ 10M ⊕ ) for typical disk parameters. In the case of small planets, however, the atmosphere is not very effective in capturing planetesimals because it is not massive enough. On the other hand, the effect of atmospheric gas drag on accretion rates is most important for such large bodies as M ∼ 10M ⊕ , because they have thick and massive atmospheres . Also, the formation timescale of protoplanets is determined by the later stage of accretion (i.e., when M ∼ 1-10M ⊕ for ∼ 5AU) where their growth timescale is long. Thus here we concentrate on the cases of large protoplanets, and compare accretion rates for atmospheres with power-law density profiles and those for more realistic atmospheric structures. As examples, we consider cases of 1M ⊕ and 10M ⊕ planets at 5AU when disk gas temperature is 125K and density is 1.5 × 10 −8 kg m −3 . We set the density of solid part of the planets to be 3.4 × 10 3 kg m −3 (r p = 0.001), which represents the density of 10M ⊕ planets that consist of ice and rock for roughly even mass (e.g., Fortney et al. 2007). Figure 11 shows the density profiles derived by the method given by  together with fitted power-law profiles. As we can see, the fitted power-law functions overestimate the gas density in the outer region (r 0.01), by a factor of up to ∼ 10 2 for the 10M ⊕ case. As for the 1M ⊕ case, the line for the power-law distribution crosses the realistic one and underestimates the gas density at the outermost region (r > 0.2) because the region is already outside of the outer boundary of the analytic solution given by  and the temperature and density are constant. The density profile at the outer region also depends on the outer boundary condition, which will be briefly discussed later in this section. Figure 12 shows ∆Ẽ obtained by orbital integration as a function ofr min for the two density profiles for the 10M ⊕ planet for various planetesimal sizes (10 -10 7 m) withẽ = 0.1 and 10 in order to examine the dependence ofR c on r s (or ξ) andẽ whenẽ = 2ĩ (the values ofr min for the upper end of the horizontal band for each curve in this plot corresponds toR c , as we discussed in §3.2). In the case ofẽ = 0.1 (left panel), we can see from the shapes of the curves that the difference of ∆Ẽ directly reflects the assumed density profiles (Fig. 11).
When planetesimals are large (r s = 10 4 km),R c of the two atmospheres are almost the same because the two curves eventually agree where the curves reach ∆Ẽ ∼ 2 − 3, which is necessary to be captured within the Hill sphere whenẽ ∼ 0.1. But as planetesimals become smaller, the curves of ∆Ẽ, which is dissipation energy per unit mass, shift upward, which yields two differentR c , the values for the realistic atmospheric model being smaller. When r s =1km, the difference in the value ofR c for the two atmospheric profiles is about a factor of three. In the case ofẽ = 10 (right panel), energy dissipation required for capture is larger. Thus, in order to be captured, planetesimals have to pass through the inner dense region of the atmosphere, where the difference of the two density profiles are smaller, yielding better agreement inR c forr s ≥ 1km. Figure 13 shows the plots of P acc in the case ofẽ = 2ĩ for the two atmospheres for various planetesimal sizes. In the case of M = 10M ⊕ (left panel), P acc for the fitted powerlaw density profile overestimates the accretion rate for small planetesimal sizes, as easily deduced from the results of Figs. 11 or 12, whereas the two values of P acc for the two atmospheric models are almost the same when planetesimals are large. Smaller planetesimal cases tend to give a larger difference in P acc because planetesimals are captured at outer atmospheres, where the difference of the two density profiles is larger. The values of P acc differ by a factor of 2-3 for 1km planetesimals, and the difference is about a factor of 15 for 10m planetesimals. As for the dependence on eccentricity, there are two effects that shift P acc in opposite directions. One is theẽ-dependence of the amount of energy dissipation required for capture, which we already discussed above and the difference in P acc for the two atmospheric profiles becomes smaller for planetesimals with higher eccentricities (Fig. 12). This tendency can be seen in the case of r s = 10 2 km, where the difference between the results for the two atmospheric profiles becomes smaller with increasingẽ. The other is the effect of gravitational focusing by the planet. When the initial random velocity of incoming planetesimals is small (ẽ 1), velocities of planetesimals atR c (≪ 1) are determined by the planet's gravity, and thus the effective collision cross section is proportional toR c (cf. Eq. (18)). Therefore 5 , the difference in P acc is roughly proportional to the difference of R c . However, when the initial random velocity is large, planetesimals are not significantly accelerated by the planet's gravity at locations withr 6/ẽ 2 . Thus, whenR c 6/ẽ 2 , P acc tends to be proportional to the geometrical cross section, i.e., ∝R 2 c , which amplifies the difference in P acc . 6 This tendency can be seen in the case of r s = 10m, where the difference between the results for the two atmospheric profiles becomes larger with increasingẽ. In the case of M = 1M ⊕ , the difference between the values of P acc for the two atmospheric profiles can be explained basically in a similar way. In both 1M ⊕ and 10M ⊕ cases, the difference between the results for the atmospheric profiles is rather small in the low-velocity regime withẽ ≤ 1. One difference from the M = 10M ⊕ case is that P acc with the power-law model for the smallest planetesimal case (r s = 0.1m) underestimates the accretion rate wheñ e 10, which simply reflects the underestimation of the density at the outer region (see Fig. 11). However, it should be noted that the assumption of our numerical method (i.e., planetesimals move on a Keplerian orbit while weakly perturbed by gas drag) breaks down for such small planetesimals. Their motion is expected to be strongly coupled with the gas flow, and the accretion rate of such small planetesimals is determined by the flow pattern of the nebular gas in the vicinity of a planet, rather than the structure of the atmosphere near the planet's surface. Therefore, the difference in the accretion rates for the two atmospheric profiles is practically negligible in the M = 1M ⊕ case, if we take into account the fact that sufficiently small planetesimals are coupled with the gas drag flow.
In summary, the power-law atmospheric density profile is a reasonable approximation for accretion rates of large (∼ 10M ⊕ ) protoplanets when planetesimals' size is not too small. In the case of 1M ⊕ , the difference of the two atmospheric profiles practically does not cause a significant difference in P acc , if we take account of the fact that the motion of planetesimals with very small sizes ( 0.1m) are strongly coupled with the nebular gas.
We further note that the outer boundary condition of atmospheres is important for the structure of the outer region of atmospheres, where the difference between realistic density profiles and fitted power-law profiles can be significant, although the inner dense region, which account for most of the atmospheric mass, is insensitive to the boundary condition. For example, when the temperature of the outer boundary is decreased to 50K, which is a typical disk temperature at 5AU for a disk that is optically thin in the vertical direction and optically thick in the horizontal direction, the agreement between the realistic density profile and the fitted power-law function is significantly improved (see Fig. 11); in this case, the difference in accretion rates for the two atmospheric profiles is expected to be reduced greatly. Conversely, when the density at the outer boundary decreases due to the dispersal of the nebular gas, for example, the difference between the two profile increases (Fig. 11).
roughly ∝R 2 c , planetesimals need a close encounter that minimum distance to the protoplanet have to be closer to the enhanced radius, which means that planetesimals would be captured by the encounter.

Growth rates of protoplanets
In order to examine the consequence of the obtained formula of accretion rates, we calculate the time evolution of planet mass using a simple semi-analytic model with the accretion rate obtained above. We basically follow the calculation model of Chambers (2006a), and also compare our results with his results.
We consider a protoplanet embedded in a disk of single-size planetesimals, whose random velocity is determined by the equilibrium between damping by gas drag and excitation by the protoplanet. The growth rate of a protoplanet is given by where Σ p is the surface density of planetesimals and T is the orbital period of the protoplanet.
To integrate this equation, we need formulas for the surface density and P acc . The surface density is given by where Σ tot is the initial surface density of planetesimals. This expression was derived under the assumption that the protoplanet does not migrate radially. In this case, the final mass of the planet becomes the isolation mass, (2πa 2b f Σ tot ) 3/2 /(3M * ) 1/2 , whereb f is the width of the feeding zone of the planet scaled by the Hill radius. As for the normalized accretion rate, we use our semi-analytic expression for the high-velocity regime (Eq.(38)) in accordance with Chambers (2006a), who used the expression based on  for the same velocity regime. Equation (38) includes the gas drag parameter ξ, which is a function of the accretion rate (see Eq. (7)). Thus we solve simultaneous equations with respect to dM/dt (Eq.(42) with Eq.(38)) and ξ self-consistently. We set Σ tot = 100 kg m −2 and a = 5AU, and also assume that mean molecular weight is 2.8, material density is 2.0 × 10 3 kg m −3 , gas density of the disk is 3.9 × 10 −8 kg m −3 , opacity is 1.0 × 10 −3 m 2 kg −1 , and α = 3 for the atmospheric gas density distribution. Figure 14(a) shows the growth of the protoplanet. We first calculated the growth with the formula used in Chambers (2006a) to check our calculation, and confirmed that their results in the case of planetesimals with 10km diameter can be reproduced (thin dashed line in Fig. 14(a); see also Fig.8 of Chambers (2006a)). As demonstrated by the semi-analytic calculation of Chambers (2006a), acceleration of the growth of protoplanets by the effect of atmospheric gas drag becomes notable when M > M ⊕ , because at this stage the atmospheres of the protoplanets become dense enough to capture 10-km-diameter planetesimals by gas drag.
On the other hand, the three thick lines in Fig. 14(a) show results obtained using our new formula for accretion rates with the effect of atmospheric gas drag. We also plot the change of values of ξ in the course of planetary growth in these three cases in Fig. 14(b). We can see that the growth based on our formula of the accretion rate for planetesimals with 10km is faster than that obtained by Chambers (2006a) by a factor of about two, which reflects the difference between our formula and that of , as discussed in §4 and Appendix B. In addition to the case of 10-km-diameter planetesimals, we also performed calculations with smaller planetesimal sizes (1km and 100m). The acceleration of the protoplanet's growth is much more significant in the case of small planetesimals for two reasons. First, in the early stage where M 10 −2 M ⊕ , the growth rate is enhanced for the small planetesimal cases because their equilibrium random velocity is smaller due to strong drag by the protoplanetary disk gas than that in the case of the larger planetesimal case. Second, in the late stage where M 0.1M ⊕ , the enhancement of accretion rates by atmospheric gas drag is more significant in the small planetesimal case, because corresponding values of ξ are larger for small planetesimals as shown in Fig. 14(b). As a result, the growth is accelerated by more than a factor of 10 in the case of 100-m-diameter planetesimals as compared to the case of 10-km-diameter. We do not intend to draw many conclusions from our growth calculation presented in Fig. 14, because our model is rather simple and does not include other important effects such as fragmentation or migration. However, our calculation shows the importance of atmospheric gas drag in planetary accretion, as first demonstrated and emphasized by  and , and it also demonstrates the usefulness of the new formula for accretion rates derived in the present work. A more detailed growth model, which includes fragmentation of planetesimals, non-equilibrium eccentricity, and radial migration as well as planetesimal capture by atmospheric gas drag can be found in Chambers (2006aChambers ( , 2008.

Conclusions and Discussion
In the present work, we examined accretion rates of planetesimals onto protoplanets that have atmospheres by analytic calculation and numerical orbital integration with gas drag. Assuming that the radial distribution of the atmospheric gas density can be approximated by a power-law and that most of energy dissipation of a planetesimal passing through a protoplanet's atmosphere occurs near the point of closest approach to the protoplanet, we analytically estimated dissipation of kinetic energy of the planetesimal, and confirmed agreement with orbital integration. We performed three-body calculation for a large number of orbits and obtained accretion rates for a wide range of parameters. We also examined the case when planetesimals have the Rayleigh distribution in orbital eccentricities and in-clinations, and derived a semi-analytic expression of the accretion rate with the effect of atmospheric gas drag.
We found that our semi-analytic formula given by Eq. (26) can reproduce our numerical results for the accretion rate of planetesimals with the Rayleigh distribution in eccentricities and inclinations by planets with atmospheres very well. The formula is basically described as a function of five non-dimensional parameters: the normalized radius of solid surface of planetsr p , the non-dimensional gas-drag coefficient ξ, the exponent of power-law function for the gas density profile α, and r.m.s. eccentricityẽ * and inclinationĩ * of planetesimals. Many other physical parameters with real dimension such as semi-major axis, planet mass, radius of incoming particles, etc. are all reduced to ξ as in Eq. (7), thus our results can be applied to a wide range of situations through the above non-dimensional parameters.
We also performed orbital calculation with gas drag from an atmosphere with a realistic density distribution, and compared with the above results for an atmosphere with a powerlaw density distribution. We found that the results using these two different atmospheric profiles agree well with each other, for example, when the protoplanet's mass is 10M ⊕ and planetesimals are not too small (r s 1km), while the results using the power-law density distribution tends to overestimate the accretion rate when planetesimals' size is smaller and their random velocity is large (ẽ ≫ 1). We also found that the degree of deviation depends on the outer boundary condition of the atmosphere because the structure of the outer atmosphere is sensitive to the boundary.
Using a simple semi-analytic model and the above results of accretion rate with atmospheric gas drag, we performed calculation for the growth of protoplanets. We confirmed the results of previous studies that atmospheric gas drag can significantly enhance the growth rate of protoplanets Chambers 2006a). We also found that the acceleration of the growth is significantly enhanced for small planetesimals. Since our semi-analytic formula is expressed in a general form in terms of non-dimensional parameters, it can be used in the study of planetary accretion in various situations.
(MEXT-16077202), and CPS running under the auspices of the MEXT Global COE Program entitled "Foundation of International Center for Planetary Science". Numerical simulation was performed at the Global Scientific Information and Computing Center of Tokyo Institute of Technology and at the Center for Planning and Information Systems of Japan Aerospace Exploration Agency.
Appendix A: Gas density distribution of a wholly radiative atmosphere near the planetary surface Equations of hydrostatic equilibrium and diffusion approximation of radiative transfer are given by 16σT 3 3κρ where P is pressure, G is the gravitational constant, M is planet's mass, r is distance from the planet center, ρ is gas density, σ is the Stefan-Boltzmann constant, T is gas temperature, L is luminosity, and κ is opacity. Eliminating r from equations (44) and (45), we have Integrating this equation, we have where temperature and pressure at the outer boundary (photosphere) is assumed to be much smaller than those of atmospheres near the planetary surface. Using the equation of state P = ρR ′ T (R ′ is the gas constant), Substitution of this into Eq. (45) yields Integrating this, we have Substituting this into Eq. (48), we have Since M = (4/3)πr 3 p ρ core , where r p and ρ core are radius and mean density of the planet, respectively, the density at the surface of the planet can be written in the form: Appendix B: Comparison with Inaba and Ikoma (2003) As shown in Fig. 10, accretion rates derived in the present work are systematically larger than those obtained by using the formulas derived by . Here, we compare our results given by Eq. (26) with their results, which are obtained mainly by a simpler way (i.e., two-body orbital integration and rough estimation of dissipation energy of gas drag), and we examine the cause of the difference.
Calculating the amount of planetesimals' kinetic energy dissipated by gas drag under the two-body approximation neglecting the solar gravity,  derived a relation between the enhanced radius of a protoplanet and the size of planetesimals, which also depends on the atmospheric density distribution (their Eq. (17)). Substituting the obtained enhanced radius into the semi-analytic expressions of collision rates derived by Inaba et al. (2001), accretion rates for protoplanets with atmospheres are obtained (Inaba and Ikoma 2003, their Eqs.(20)-(24)).
Since P acc is a quantity that is integrated with respect toẽ,ĩ, andb, P acc is not an explicit function ofẽ,ĩ, andb. However, P acc of  ( P col in their notation) is a function ofR c , which explicitly depends onẽ,ĩ, andb. Here we assume thatṽ ∞ = (ẽ * ) 2 + (ĩ * ) 2 forR c in calculating their P acc . In the following comparison, we assume an atmospheric structure described by a power-law function (i.e., Eq.(4)). Figure 10 shows normalized accretion rates obtained by our semi-analytic expression and that of  when ξ = 10 −7 andr p = 10 −3 . In this case, all accreting objects are practically captured by the gas drag of the atmosphere, rather than direct collision. We can see that the overall features of the two results are similar to each other; the accretion rates in the low velocity regime roughly agree with each other and those in the high velocity regime have the same slope. However P acc for ξ = 10 −7 of our result is about a factor of 2.3 higher than theirs. We find that the difference of the factor of 2.3 between the two formulas comes from the following three effects: (i) The amount of energy dissipation for incoming particles in an atmosphere is estimated by their Eq.(16) in . This estimate is smaller by a factor of about two compared to the more accurate estimate in the present work. The correction factor is denoted by f c in the present work and was estimated empirically to be 1.83. This factor appears in the expression of the accretion rate as P acc ∝ f 1/α c (Eq.(38)), thus it enhances P acc by ∼ 1.83 1/3 = 1.22. (ii)  does not seem to fully take account of theb-dependence ofR c in obtaining the accretion rate averaged over the Rayleigh distribution. If it is correctly taken into account, this enhances the accretion rate by 1.17. (iii) Instead of fully taking account of theẽ-and i-dependence ofR c and the averaging over the Rayleigh distribution,ṽ ∞ = (ẽ * ) 2 + (ĩ * ) 2 seems to have been used in their calculation of P col . If we fully take account of this, the accretion rate is enhanced by about 1.6. Although all of their calculation procedures are not explicitly written in , the product of the above three factors is about 2.3, which is consistent with the difference observed in Fig. 10, thus we believe that the difference comes from the three factors described above. Based on the consideration, we can deduce the reason why the two accretion rates in the low velocity regime roughly agree with each other. As described above, two causes out of the three are related to the dependence ofR c onẽ,ĩ andb, whileR c does not depend on them in the low velocity regime. Fig. 1.-Enhanced radiusR c as a function ofṽ ∞ when α = 3 and f c = 2. Four lines correspond to the cases with ξ = 10 −3 , 10 −5 , 10 −7 , 10 −9 from top to bottom. Core radius in this normalization is typically about 0.005 for 1AU and 0.001 for 5AU.   -Total change of the energyẼ due to gas drag through an encounter with the planet for each orbit, as a function of the minimum distance to the planet center in the case with α = 3.113. Red pluses and green crosses show the low-velocity (ẽ = 3) and high-velocity (ẽ = 30) cases, respectively. For each value ofẽ, numerical results for ξ = 1 × 10 −9 with a wide range of initial phase angle τ and impact parameterb in the two-dimensional case (ĩ = 0) are shown. Blue asterisks show the case without gas drag, that is, the error of the numerical integration. Long-dashed line shows the approximate analytic results (Eq.(12)) for the high-velocity regime (i.e., the term proportional tor 1−α min ) and short-dashed line shows that for the low-velocity regime (the term proportional tor −α min ). We set f c = 2. Left panel shows the dependence on ξ when (ẽ,ĩ) = (0.1, 0.05): red pluses, green crosses, and blue asterisks are the cases with ξ = 2.56 × 10 −9 (weak), 2.56 × 10 −8 (intermediate), and 2.56 × 10 −7 (strong), respectively. Right panel shows the dependence on initial orbital elements when ξ = 2.56 × 10 −7 : red pluses, green crosses, and blue asterisks are the case with (ẽ,ĩ) = (10, 5), (1, 0.5), and (0.1, 0.05), respectively. The vertical dashed line shows the planet's radiusr p = 5 × 10 −3 .    21) and (22) are also shown. Fig. 9.-Normalized accretion rates averaged over the Rayleigh distribution forr p = 0.005 (left panel) andr p = 0.001 (right panel), as a function ofẽ * (= 2ĩ * ). Marks show simulation results for ξ = 10 −1 , 10 −3 , 10 −5 , 10 −7 , and 10 −9 , and lines show the semi-analytical formula given by Eq. (26). The Rayleigh distribution average of collision rate (Eq. (27)) is also shown by the long dashed line for comparison in the case ofr p = 0.001 (left panel). Note that, in the case of r p = 0.005 and ξ = 10 −9 , there is no effect of the atmosphere on P acc , i.e., P acc = P col Fig. 10.-Normalized accretion rate averaged over the Rayleigh velocity distribution as a function ofẽ * (= 2ĩ * ) for ξ = 10 −7 andr p = 0.001. Solid line shows the accretion rate obtained by our semi-analytic expression (Eq.(26)), and dashed line is obtained by using the formulas derived by . In both cases, a power-law density distribution is assumed for the atmospheric structure (α = 3). Collision rate for the case without an atmosphere is also shown by the dotted line for comparison. Fig. 11.-Density profiles of a realistic atmosphere derived by the analytic method given in  for 1M ⊕ and 10M ⊕ planets at 5AU (thick-solid lines) in the case when ρ core = 3.4 × 10 3 kg m −3 , L = 10 −7 L ⊙ , and f κ = 10 −2 , and power-law functions fitted to the realistic one at the bottom (thick-dashed lines). The exponents for the fit lines are α = 3.094 for 1M ⊕ and α = 3.027 for 10M ⊕ . Thick lines show the case when the temperature and gas density at the outer boundary are 125K and 1.5 × 10 −8 kg m −3 , respectively. In order to see the dependence on the boundary condition, we also show density profiles when the temperature at the boundary is lower (50K, thin dot-dashed line) and the density is lower (1.5 × 10 −9 kg m −3 , thin dotted line) in the 10M ⊕ case. Fig. 12.-Change of the energyẼ due to gas drag as a function of minimum approach distance (the same as Fig. 5) for a realistic atmosphere model (red pluses) and a power-law model (green crosses) when M = 10M ⊕ at 5AU. Left and right panels show the case when (ẽ,ĩ) = (0.1, 0.05) and (10, 5), respectively. Four different planetesimal sizes are assumed for each of the two atmosphere models; r s = 10 4 km, 10 2 km, 1km, and 10m, from lower-left to upper-right.