A&A 410, 711-723 (2003)
DOI: 10.1051/0004-6361:20031248

Enhanced collisional growth of a protoplanet that has an atmosphere

S. Inaba1 - M. Ikoma2


1 - Department of Earth and Planetary Sciences, Tokyo Institute of Technology, 2-12-1 Oookayama, Meguro-ku, Tokyo 152-8551, Japan
2 - Interactive Research Center of Science, Tokyo Institute of Technology, 2-12-1 Oookayama, Meguro-ku, Tokyo 152-8551, Japan

Received 14 February 2003 / Accepted 12 May 2003

Abstract
Once a protoplanet becomes larger than about lunar size, it accumulates a significant atmosphere that surrounds the solid core. When a planetesimal approaches the protoplanet, it interacts with the atmosphere. If enough energy of the planetesimal is lost by gas drag of the atmosphere, it is captured in the atmosphere even if its original trajectory would not lead to a direct collision with the solid core of the protoplanet. This increases the collision rate, resulting in faster growth of the protoplanet. We have derived the analytical calculations for the collision rate, and calculated the structure of the atmosphere and the trajectories of the planetesimals in the atmosphere. As a result of their large gas drag, small planetesimals are easily captured, resulting in a large rate of collision with the protoplanet. A collision rate of a protoplanet of Earth size with a planetesimal of 100 m radius is, for example, enhanced by a factor of ${\sim} 10$. These effects play an essential role in the study of formation of solid cores of gas giant planets by the core accretion model.

Key words: solar system: formation

1 Introduction

Planets form in a circumstellar disk composed of gas and dust. A gas giant planet in particular has to be born in such a disk to acquire a large amount of gas. We are assuming the standard core accretion model in which formation of gas giant planets proceeds as follows. A number of planetesimals of $\sim$km radius form by collisions and sticking of small dust particles (Weidenschilling 1997a), or gravitational instability of a dust disk (Goldreich & Ward 1973; Goodamn & Pindor 2000). Runaway growth of planetesimals occurs in the course of accumulation of gravitationally interacting planetesimals, resulting in a bimodal mass distribution of planetesimals, where a small number of protoplanets and a large number of small planetesimals coexist (Wetherill & Stewart 1989). Protoplanets larger than about the size of the moon can attract gas gravitationally and can have atmospheres. The solid part of the protoplanet is often called the core. The atmospheric mass increases with the mass of the core. When the atmospheric mass is comparable to the mass of the core, runaway gas accretion occurs, and the protoplanet acquires a large amount of gas from a circumstellar disk to become a gas giant planet (Mizuno 1980; Bodenheimer & Pollack 1986).

One of the challenging problems in the core accretion model is the timescale-problem for the formation of Jupiter and Saturn. In a pioneering contribution, Pollack et al. (1996) simulated the formation of gas giant planets calculating the gas and planetesimal accretion rates in a self-consistent manner and constructed one plausible pathway by which Jupiter and Saturn could have formed in their current positions before the disappearance of the disk gas. They adopted a simple model of planetesimal accretion, where a single core grows by accumulation of planetesimals whose velocities are constant, and the effects of fragmentation of planetesimals was not considered. Weidenschilling (1997b) reexamined the accretion process of the solid core under the same initial condition except that velocities of planetesimals are not constant. The core could not grow to sufficiently large mass to acquire a large amount of gas within the probable lifetime of the disk gas.

In the early stage of the planetesimal accumulation, the velocities of the planetesimals are determined by the gravitational interactions with other planetesimals and are so small that fragmentation does not occur when they collide with each other. On the other hand, in the late stage of the planetesimal accumulation, protoplanets enhance the velocities of the planetesimals (Ida & Makino 1993). Collisions of the planetesimals in the high velocities lead to their fragmentation. The created fragments are vulnerable to radial migration to the central star by gas drag of the disk gas (Wetherill & Stewart 1993). A portion of the mass of solid is removed as fragments. As a result, the growth rates of the cores decrease. Inaba & Wetherill (2001) calculated planetesimal accumulation with fragmentation and realized that even the largest core could not reach the so-called critical core mass, beyond which the runaway gas accretion occurs, and thus failed to become a gas giant planet .

Inaba & Wetherill (2001) assumed that the planetesimals were captured by the protoplanet only when they collided with the core of the protoplanet. The rate of these collisions has been examined in detail by many researchers (e.g., Ida & Nakazawa 1989; Greenzweig & Lissauer 1992). However, this collision rate is not useful if the protoplanet has an atmosphere. Because the atmosphere surrounds the core, planetesimals interact with the atmosphere first. Planetesimal energies are reduced by gas drag of the atmosphere, and if enough energy is lost, they can be captured in the atmosphere. The gas drag depends not only on the gas density of the atmosphere but also the size of the planetesimal. As mentioned above, many small fragments are produced by collisions of planetesimals. The small fragments have large gas drags and are easily captured in the atmosphere. As a result, fragments that would not be captured by base cores can be captured in the atmosphere. Because the radius of the atmosphere can become much larger than that of the core, a large increase in the collision rate is expected. In this paper we will calculate this collision rate.

Previously, Podolak et al. (1988) studied the interaction of planetesimals with atmospheres, including detailed physics of gas drag and ablation of the planetesimals. They calculated the trajectories of the planetesimals and obtained the range of the impact parameters with which the planetesimals are captured in the atmosphere. The range of the impact parameter is strongly dependent on the radius of the planetesimal. Their choice of parameters is so limited that it is difficult to use their result in the simulation of planetesimal accumulation. In this study we will derive an analytical expression of the collision rate, which is useful for a wide range of the parameters.

The effects of the atmospheric gas drag on the accumulation of planetesimals are also important to explain the present compositions of the envelopes of Jupiter and Saturn. They are known to be enriched in heavy elements (at least by a factor of 5 for Saturn) compared with the solar value (Wuchterl et al. 2000). Some captured planetesimals would be ablated and produce small particles that are vaporized in an inner part of the atmosphere (i.e., the proto-envelope). This would increase the amount of heavy elements in Jupiter's and Saturn's envelopes.

In the present paper, we will investigate what size of planetesimals are captured in what location of the atmosphere of the protoplanet. If the planetesimals are captured far away from the surface of the core, it helps the growth of the protoplanet and enhances the total mass of heavy elements in the atmosphere. In Sect. 2, we first describe the calculation methods of the atmospheric structure and motion of planetesimals in the atmosphere. Next, we find that the effective radius of a planet capturing a planetesimal is much larger than the radius of the solid core. After that, we derive an approximate expression of the effective radius, which gives us insight into its dependence on several parameters concerning the structure of the atmosphere as well as the dynamics of planetesimals. In Sect. 3, we obtain the enhanced collision rate of a protoplanet that has an atmosphere.

2 Enhanced radius of a protoplanet that has an atmosphere

2.1 Definition of the enhanced radius

The gas density of the atmosphere is much larger close to the center of the protoplanet (see Fig. 1). A planetesimal is more easily captured in the close region because of the large gas drag. If the distance between the planetesimal and the protoplanet is smaller than a critical distance, the planetesimal is captured in the atmosphere. Otherwise, the planetesimal passes through the atmosphere with a loss of energy. In this subsection, we will define the critical distance.


  \begin{figure} \par\resizebox{11cm}{!}{\includegraphics[clip]{MS3608f1.eps}}\end{figure} Figure 1: Gas density of the atmosphere is plotted as a function of distance from the center of the planet. The solid and dashed lines are the numerical and analytical solutions, respectively. The detail of the former is described in the Sect. 2.2 and that of the latter in Appendix. a) The dependence on the core mass, $M_{\rm core}$, and b) the dependences on the luminosity, L and the grain depletion factor, f, are shown.
Open with DEXTER

Before the planetesimal enters the atmosphere, the energy of the planetesimal is constant and determined by $v_{\infty}$ where $v_{\infty}$ is the relative velocity between the planetesimal and the protoplanet when their distance is infinite. We assume that the mass of the planetesimal, $m_{\mbox{\scriptsize p}}$, is much smaller than that of the protoplanet, $M_{\mbox{\scriptsize p}}$. Then the impact energy is given by $E_0 = 1/2 m_{\mbox{\scriptsize p}} v_{\infty}^2$. Once the planetesimal gets into the atmosphere, the gas drag starts removing the energy, $E = 1/2 m_{\mbox{\scriptsize p}} v^2 - G M_r m_{\mbox{\scriptsize p}}/r$, where v is the relative velocity, r is the relative distance between the planetesimal and the protoplanet, and Mr is the total mass contained within r. In the two body system composed of the planetesimal and the protoplanet, the planetesimal is captured by the protoplanet if E<0.

In order to include the effect of a central star approximately in this two body system, we introduce the Hill radius of the protoplanet, $r_{\mbox{\scriptsize H}}$, within which the gravity of the protoplanet is larger than that of the central star, defined by (Ida & Nakazawa 1989)

 \begin{displaymath}r_{\mbox{\scriptsize H}} =\left( \frac{M_{\mbox{\scriptsize p... ...box{\scriptsize c}}} \right)^{1/3} A_{\mbox{\scriptsize p}}, \ \end{displaymath} (1)

where $M_{\mbox{\scriptsize c}}$ is the mass of the central star and $A_{\mbox{\scriptsize p}}$ is the semimajor axis of the protoplanet. In this system, the planetesimal is captured by the protoplanet if the apocentre of the planetesimal that orbits the protoplanet is within the Hill radius. Otherwise, the planetesimal is not captured by the protoplanet but orbits the central star. The largest energy of the captured planetesimal is given by $E_{\mbox{\scriptsize c}} = - G M_r m_{\mbox{\scriptsize p}}/(r_{\mbox{\scriptsize H}} + r)$. In general, r is much smaller than $r_{\mbox{\scriptsize H}}$, $E_{\mbox{\scriptsize c}} = - G M_r m_{\mbox{\scriptsize p}}/r_{\mbox{\scriptsize H}}$. If E becomes smaller than $E_{\mbox{\scriptsize c}}$ by the gas drag, the planetesimal is captured by the protoplanet. The critical distance, $R_{\mbox{\scriptsize c}}$, is defined by the value of r when $E = E_{\mbox{\scriptsize c}}$, that is,

 \begin{displaymath}\frac12 v^2 - \frac{GM_{\scriptsize r}}{r} = - \frac{GM_{\sc... ...ace{3ex} \mbox{at} \hspace{1ex} r = R_{\mbox{\scriptsize c}}. \end{displaymath} (2)

Because the planetesimal is captured by the protoplanet if r is smaller than $R_{\mbox{\scriptsize c}}$, $R_{\mbox{\scriptsize c}}$ reduces to the radius of the core, $R_{\mbox{\scriptsize core}}$, if the protoplanet does not have the atmosphere. As shown later, $R_{\mbox{\scriptsize c}}$ is much larger than $R_{\mbox{\scriptsize core}}$. It is adequate to call $R_{\mbox{\scriptsize c}}$ the enhanced radius of the protoplanet. In Sect. 3, we will show that the discussion for the two body system approximately hold true for the three body system that includes the central star.

2.2 Structure of the atmosphere

As mentioned above, the density of the atmospheric gas is one of the key factors for the enhanced radius, because the gas drag is dependent on it. Thus, we should know the density profile of the atmosphere as accurately as reasonable. In this study, we calculate the structure of the atmosphere.

We consider a spherically-symmetric protoplanet composed of a rigid core and a chemically-uniform atmosphere. The mass fractions of hydrogen, helium, and the other heavy elements are 0.740, 0.243, and 0.017, respectively, that are the same to those used by Pollack et al. (1996). The atmosphere is assumed to be in purely hydrostatic equilibrium. The outer boundary of the atmosphere (i.e., the radius of the protoplanet, R) is assumed to be the smaller of the Hill radius, $r_{\rm H}$, defined by Eq. (1) and the Bondi radius, $r_{\rm B}$, defined by

 \begin{displaymath}r_{\rm B} = \frac{GM_{\rm p}}{c_{\rm S}^2}, \end{displaymath} (3)

where $c_{\rm S}$ is the sound velocity of the disk gas. The atmosphere is assumed to connect smoothly with the disk at radius R. These assumptions are basically the same as those of Mizuno (1980).

The basic equations are the usual set of equations for study of stellar structure (e.g., Kippenhahn & Weigert 1990):

 \begin{displaymath}\frac{{\rm d}P}{{\rm d}r} = - \frac{GM_r\rho}{r^2}, \end{displaymath} (4)


 \begin{displaymath}\frac{{\rm d}M_r}{{\rm d}r} = 4 \pi r^2 \rho, \end{displaymath} (5)


 \begin{displaymath}\frac{{\rm d}T}{{\rm d}r} = \left\{ \begin{array}{ll} \dis... ..._r\rho}{P r^2 T} } & \mbox{otherwise}. \end{array} \right. \end{displaymath} (6)

Here P is the pressure, r is the distance from the center of the core, Mr is the mass inside a sphere of radius r, $\rho$ is the density, T is the temperature, $\kappa$ is the opacity, $a_{\rm r}$ is the radiation density constant, c is the light velocity, $\Gamma_2$ is the second adiabatic exponent, and Lr is the energy per unit time passing through a sphere of radius r. In this study, we assume

 \begin{displaymath}L_r = L = \mbox{constant}, \end{displaymath} (7)

where L is called the luminosity. This assumption might be oversimplification, because it is likely that energy sources are distributed by planetesimals throughout the atmosphere. As for the situation under consideration, the energy deposition can occur at various positions, depending on planetesimal size. However, the dependence of the enhanced radius on the luminosity is small, as will be found later.

The opacity is given by

 \begin{displaymath}\kappa (\rho, T) = \kappa_{\rm gas} (\rho, T)+f\kappa_{\rm gr} (\rho, T), \end{displaymath} (8)

where $\kappa_{\rm gas}$ is the gas opacity given by Alexander & Ferguson (1994), $\kappa_{\rm gr}$ is the grain opacity with an interstellar size distribution given by Pollack et al. (1985), and f is the so-called grain depletion factor introduced by Mizuno (1980). We treat L and f as parameters. The equation of state used here is one including non-ideal effects given by Saumon et al. (1995).

The inner boundary condition is

 \begin{displaymath}\begin{array}{ll} M_r = M_{\rm core} & \mbox{at} \hspace{1e... ...rm core}}{4 \pi \rho_{\rm core}} \right)^{1/3} }, \end{array}\end{displaymath} (9)

where $\rho_{\rm core}$ is the mean density of the core. The outer boundary conditions are

\begin{displaymath}\begin{array}{llll} T = T_{\rm neb} & \mbox{and} & \rho = \rho_{\rm neb} & \mbox{at} \hspace{1ex} r = R, \end{array}\end{displaymath} (10)

where $T_{\rm neb}$ and $\rho_{\rm neb}$ are the temperature and the density of the disk gas, respectively. The values of $\rho_{\rm core}$, $T_{\rm neb}$, and $\rho_{\rm neb}$ adopted in this study are $3.2~ \rm g~cm^{-3}$, $150~ \rm K$, and $5.0 \times 10^{-11}~ \rm g~cm^{-3}$, respectively, that are the same used by Pollack et al. (1996).

In Fig. 1, we show examples of the density profiles for some different values of $M_{\rm core}$, L, and f. In this figure, we also compare the results of our numerical calculations with the analytical solutions derived in Appendix. The density increases with decreasing distance from the planet's center: its dependence being expressed as $\rho \propto 1/r^{3-4}$. The density at a given radius increases rapidly with the mass of the core (see Fig. 1a), because the gravity becomes strong. On the other hand, the larger the luminosity Land/or the grain depletion factor f (i.e., the opacity) becomes, the smaller the density becomes especially in deep regions of the atmosphere (see Fig. 1b). All of this behavior can be understood by the help of the analytical solutions obtained in Appendix.

2.3 Numerical simulation of planetesimal's motion

In the numerical simulation, we consider the gravity of the core and the atmosphere, the gas drag of the atmosphere, and ablation. We assume that the planetesimal has no atmosphere. The equations that describe the time evolution of the position and mass of the planetesimal are given by the following two equations (Borovicka & Spurný 1996):

 \begin{displaymath}\displaystyle \left\{ \begin{array}{@{~}ll} \displaystyle m... ... D \pi r_{\mbox{\scriptsize p}}^2 \rho v^3, \end{array}\right. \end{displaymath} (11)

where ${\vec v}$ and ${\vec x}$ are the relative velocity and position vectors, respectively, Mr is the total mass of the core and the atmosphere within the orbit of the planetesimal, D is the gas drag coefficient, and $\sigma$ is the ablation coefficient. We adopt the gas drag coefficient given by Podolak et al. (1988). Their gas coefficient is similar to that obtained by Adachi et al. (1976). We calculate Eq. (11) using a 4th order Runge-Kutta scheme. The gas density of the atmosphere is obtained numerically as described in the Sect. 2.2.

We study the two cases with and without ablation. When we consider the effect of ablation, we use the ablation coefficient, $\sigma = 0.1 ~\mbox{s}^2 ~\mbox{km}^{-2}$even though the theoretical value of $0.001~ \mbox{s}^2~ \mbox{km}^{-2}$(Field & Ferrara 1995) is much smaller (Borovicka & Spurný 1996). In the next subsection, we will show that the effect of ablation is small to obtain the enhanced radius of the protoplanet.

The energy and mass of the planetesimal are constant before entry into the atmosphere. The outer boundary of the atmosphere is determined by Eq. (1) or (3). Because we have the analytical solution to the motion of the planetesimal before the planetesimal gets into the atmosphere, we carry out the simulation after the distance between the planetesimal and the protoplanet is equal to the radius of the atmosphere, R. The parameters used in the simulation are shown in Table 1.

Table 1: Parameters used in the simulations.

2.4 Numerical results: Enhanced radius


  \begin{figure} \par\resizebox{11cm}{!}{\includegraphics[clip]{MS3608f2.eps}}\end{figure} Figure 2: a) The enhanced radii of the protoplanet with $M_{\mbox{\scriptsize core}} = 1~ M_{\oplus}$, $L = 1 \times 10^{-8} ~L_{\odot}$, and f=0.01 when the velocities of planetesimals are $1 r_{\mbox{\scriptsize H}} \Omega_{\mbox{\scriptsize K}}$and b) the dynamical pressure the planetesimals take. The cases without and with ablation are shown by the solid lines and the filled circles. The strength of the planetesimal given by Benz & Asphaug (1999) is shown by the dashed curve as well. Planetesimals with the radii larger than $5\times 10^8$ cm cannot be captured in the atmosphere.
Open with DEXTER

In Fig. 2a, we show the enhanced radii normalized by the radius of the core as a function of the radius of the planetesimal for the case of $M_{\mbox{\scriptsize core}} = 1~ M_{\oplus}$, $L = 1 \times 10^{-8} ~L_{\odot}$, f=0.01, and $v_{\infty} = 1 r_{\mbox{\scriptsize H}} \Omega_{\mbox{\scriptsize K}}$, where $M_{\oplus}$ is Earth mass, and $\Omega_{\mbox{\scriptsize K}}$ is Keplerian angular velocity and given by $ \sqrt{G M_{\odot}/A_{\mbox{\scriptsize p}}^3}$. With an increase in the radius of the planetesimal, the enhanced radius decreases. The large planetesimals have to get into a deep region of the atmosphere to be captured. The planetesimal with the radius larger than $5\times 10^8$ cm cannot be captured in the atmosphere even if it takes the largest gas drag produced in the deepest region of the atmosphere. A planetesimal of this size can only be captured by the protoplanet by colliding with the core. The enhanced radii for the small planetesimals are larger than the radius of the core. A planetesimal with radius of $1\times 10^4$ cm is captured at $R_{\mbox{\scriptsize c}} \sim 10 R_{\mbox{\scriptsize core}}$. This increases the collisional cross section between the protoplanet and the planetesimal and, as a result, the growth rate of the protoplanet increases. The mass loss of planetesimals by ablation is small. We have nearly the same results for the cases with and without ablation.

Planetesimals might lose mass by fragmentation. If dynamical pressure ($\rho v^2$) exceeds the strength of the planetesimal, fragmentation of the planetesimal occurs. The largest dynamical pressure that the planetesimal experiences is shown in Fig. 2b. An increase in the dynamical pressure is caused by a sharp increase in the gas density of a deep region of the atmosphere. It is still smaller than the strength of the planetesimal that was obtained by Benz & Asphaug (1999) with use of a smooth particle hydrodynamics method. When the planetesimals are captured, they are not fragmented by the dynamical pressure.


  \begin{figure} \par\resizebox{11cm}{!}{\includegraphics[clip]{MS3608f3.eps}}\end{figure} Figure 3: The enhanced radii of protoplanets with various parameters. a) The enhanced radii of the protoplanets for the case with $L = 1 \times 10^{-8} ~L_{\odot}$, $v = 1 r_{\mbox{\scriptsize H}} \Omega_{\mbox{\scriptsize K}}$, and $M_{\mbox{\scriptsize core}} = 4, 2, 1, 0.1 ~M_{\oplus}$. b) Those for the cases with $M_{\mbox{\scriptsize core}} = 1~ M_{\oplus}$, $L = 1 \times 10^{-8} ~L_{\odot}$, and $v = 0.1, 1, 10 r_{\mbox{\scriptsize H}} \Omega_{\mbox{\scriptsize K}}$, c) those for the cases with $M_{\mbox{\scriptsize core}} = 1~ M_{\oplus}$, $v = 1 r_{\mbox{\scriptsize H}} \Omega_{\mbox{\scriptsize K}}$, and $L = 1\times 10^{-7}, 1\times 10^{-8}, 1\times 10^{-9} ~L_{\odot}$. In every cases, the grain depletion factor is f=0.01. Comparison between the enhanced radii obtained from the numerical simulation (solid lines) and the analytical formula (dashed lines) is shown as well.
Open with DEXTER

We show the enhanced radii of the protoplanets with the different parameters (the mass of the core, the luminosity, the initial velocity) in Fig. 3 by the solid lines. The atmosphere of the protoplanet with a massive core is dense because of the large gravity of the core. The large gas drag due to the dense atmosphere increases the enhanced radius (Fig. 3a). Planetesimals with the large initial velocities have more energy to be captured. By getting into a deep region of the atmosphere, they take large gas drag and are captured (Fig. 3b). The dependence of the initial velocity on the enhanced radius is small specially for $v < 1 r_{\mbox{\scriptsize H}} \Omega_{\mbox{\scriptsize K}}$ because the velocity is determined by $\sqrt{GM_{\mbox{\scriptsize p}}/R}$ just before the planetesimal enters the atmosphere. With a decrease in the luminosity of the atmosphere, the gas density of the atmosphere increases as shown in the Sect. 2.2. The gas drag becomes more effective and the enhanced radii increase (Fig. 3c).

2.5 Approximate solution for the enhanced radius

In this subsection, we derive an approximate solution for the enhanced radius, $R_{\mbox{\scriptsize c}}$. Because the effect of the ablation is shown to be small in the last subsection, we consider only the gas drag. The gas drag of the atmosphere decrease the relative velocity from v to $v-\Delta v$, which results in the decrease in the energy as $\Delta E = m_{\mbox{\scriptsize p}} v \Delta v$. The decrease in the velocity is given approximately by

 \begin{displaymath}\Delta v = \pi r_{\mbox{\scriptsize p}}^2 \rho(r) v^2 \Delta t /m_{\mbox{\scriptsize p}}, \end{displaymath} (12)

where $r_{\mbox{\scriptsize p}}$ is the radius of the planetesimal, $\rho(r)$ is the gas density of the atmosphere at r, and $\Delta t$ is the period while the planetesimal passes the atmosphere. We assume that the velocity of the planetesimal decreases only at $r=R_{\mbox{\scriptsize c}}$. This assumption is reasonable because the gas drag is effective in a deep region by a sharp increase in the gas density of the atmosphere. Then the velocity is given by the sum of the relative velocity at infinity and the escape velocity:

 \begin{displaymath}v^2 = v_{\infty}^2 + \frac{2 G M_{\mbox{\scriptsize core}}}{R_{\mbox{\scriptsize c}}}\cdot \end{displaymath} (13)

Because the mass of the atmosphere is much smaller than that of the core, we replace the total mass of the core and the atmosphere with the mass of the core. The gas density is given approximately by that at $R_{\mbox{\scriptsize c}}$ as

 \begin{displaymath}\rho(r) = \rho(R_{\mbox{\scriptsize c}}). \end{displaymath} (14)

Because the planetesimal orbits the protoplanet at $R_{\mbox{\scriptsize c}}$while the gas drag is effective, we assume that the period, $\Delta t$, is given by

 \begin{displaymath}\Delta t = R_{\mbox{\scriptsize c}}/ v. \end{displaymath} (15)

It will be shown that the approximate solutions can reproduce the numerical results well. With the use of the above equations, we have

 \begin{displaymath}\Delta E = \frac34 m_{\mbox{\scriptsize p}} \left( v_{\infty... ...}} \frac{R_{\mbox{\scriptsize c}}}{r_{\mbox{\scriptsize p}}}, \end{displaymath} (16)

where $\rho_{\mbox{\scriptsize p}}$ is the material density of the planetesimal.

For the planetesimal to be captured by the protoplanet, the energy of the planetesimal have to be smaller than $E_{\mbox{\scriptsize c}}$ (see Sect. 2.1). If $\Delta E$ is equal to $E_0 - E_{\mbox{\scriptsize c}}$, the planetesimal is captured by the protoplanet. Using Eq. (16), we have the radius of the captured planetesimal as a function of the gas density of the atmosphere and $R_{\mbox{\scriptsize c}}$:

 \begin{displaymath}r_{\mbox{\scriptsize p}} = \frac32 \frac{\displaystyle v_{\i... ...e c}})}{\rho_{\mbox{\scriptsize p}}} R_{\mbox{\scriptsize c}}. \end{displaymath} (17)

If $v_{\infty}$ is much smaller than the escape velocity of the core, $v_{\mbox{\scriptsize esc}} =\sqrt{2 G M_{\mbox{\scriptsize core}}/R_{\mbox{\scriptsize core}}}$, $r_{\mbox{\scriptsize p}}$ is independent of $v_{\infty}$ and given by

 \begin{displaymath}r_{\mbox{\scriptsize p}} = \frac32 \frac{\rho(R_{\mbox{\scri... ...e c}})}{\rho_{\mbox{\scriptsize p}}} r_{\mbox{\scriptsize H}}. \end{displaymath} (18)

As seen in Sect. 2.2, the gas density increases with proximity to the center of the protoplanet. Equation (18) shows that large planetesimals can be captured in a deep region of the atmosphere, while small ones in an outer region.

We show the comparison between the numerical results (solid lines) and the approximate solution given by Eq. (17) with the approximate solution of the atmospheric structure given in the Appendix (dashed lines) in Fig. 3 as well. The enhanced radii for the planetesimal with radius larger than $1\times 10^9$ cm are different from each other (a factor of less than two). This is caused by the difference of the gas density obtained analytically and numerically. However, for the wide range of the size of the planetesimal, both numerical and analytical results are consistent with each other.

2.6 Comparison with the previous work

Pollack et al. (1996) simulated formation of Jovian planets, including the effects of the enhanced radius. In their simulation, the core of the protoplanet grew by accumulation of a specific size of planetesimals (100 km radius). In Fig. 4, we show the enhanced radius given by Pollack et al. (1996) as a function of the core mass (dashed curve). With an increase in the core mass, the enhanced radius increases as well. The enhanced radius becomes larger than the radius of the core after the core mass increases to ${\sim} 2 ~M_{\oplus}$. That is, the planetesimal of 100 km radius was captured in the atmosphere of the protoplanet by the gas drag. For the protoplanet with the core mass of $15 ~M_{\oplus}$, the enhanced radius became about 10 times larger than the radius of the core.

It is valuable to compare their result with that obtained in this study. By the use of the derived analytical solution (Eq. (17)), we have the enhanced radius of a protoplanet, given the core mass, the luminosity, and the radius of a planetesimal. For the comparison, we adopt the same core mass and luminosity to those of Pollack et al. (1996). The enhanced radius given by this study is shown by the solid curve in Fig. 4. The difference of the enhanced radii for the core mass larger than $12 ~M_{\oplus}$ comes from the fact that the contribution of the mass of the atmosphere to the total mass becomes large. In Pollack et al. (1996), for example, the total mass of the protoplanet is $26 ~M_{\oplus}$when the core mass is $15 ~M_{\oplus}$. On the other hand, we assume the protoplanet has a negligible atmosphere. However, we obtain the similar enhanced radius to that of Pollack et al. (1996) for the wide range of the core mass.


  \begin{figure} \par\resizebox{8cm}{!}{\includegraphics*[clip]{MS3608f5.eps}}\end{figure} Figure 4: The enhanced radii of the protoplanet normalized by the core radius given by case J1 of Pollack et al. (1996) (dashed curve) and this study (solid curve). The radius of the planetesimal is 100 km.
Open with DEXTER

3 Enhanced collision rate of the protoplanet

In this section, we will show the relationship of the enhanced radius to the collision rate. A collision rate of a protoplanet without an atmosphere (a core) and a planetesimal under the gravity of the sun was studied (Ida & Nakazawa 1989; Greenzweig & Lissauer 1992). If the distance between the protoplanet and the planetesimal is smaller than the radius of the core of the protoplanet, they will collide. Using the analytical and numerical methods, the collision rate as a function of the relative velocity and the radius of the core of the protoplanet was obtained. Inaba et al. (2001) compiled their data and used the collision rate for the simulation of planetesimal accumulation. For simplicity, we consider the case that the eccentricities and inclinations of a protoplanet and a planetesimal are zero, i.e., circular and non-inclined orbits. The semimajor axis is the only parameter. The initial relative position is determined by the difference between the two semimajor axes, $b= a_{\mbox{\scriptsize p}} - A_{\mbox{\scriptsize p}}$, where $a_{\mbox{\scriptsize p}}$ and $A_{\mbox{\scriptsize p}}$ are the semimajor axes of the planetesimal and the protoplanet, respectively. Because the orbit with $b < 1.9 r_{\mbox{\scriptsize H}}$ corresponds to the horseshoe orbit, the planetesimal cannot approach the protoplanet (Ida & Nakazawa 1989). A planetesimal with $b > 2.6 r_{\mbox{\scriptsize H}}$, is not influenced by the gravity of the protoplanet. The planetesimal with $1.9r_{\mbox{\scriptsize H}} < b < 2.6r_{\mbox{\scriptsize H}}$can get to the close place to the center of the protoplanet. Ida & Nakazawa (1989) summed the orbits, of which the closest distance is smaller than the radius of the core, and obtained the collision rate of the core and the planetesimal.


  \begin{figure} \par\resizebox{11cm}{!}{\includegraphics[clip]{MS3608f4.eps}}\end{figure} Figure 5: The closest distances of the planetesimals with radius of 100 m from the protoplanet with $M_{\mbox{\scriptsize core}} = 1~ M_{\oplus}$, $L = 1 \times 10^{-8} ~L_{\odot}$, and f=0.01for the cases without the atmosphere (the dashed line) and with the atmosphere (the solid line). The enhanced radius, $R_{\mbox{\scriptsize c}}$, obtained from Fig. 2 and the radius of the core, $R_{\mbox{\scriptsize core}}$, are shown by the horizontal dotted lines.
Open with DEXTER

If the protoplanet has an atmosphere, the situation is different. Even if the closest distance is larger than the radius of the core, it can be smaller than the enhanced radius. Then we expect the planetesimal is captured in the atmosphere even if it does not collide with the core. We carry out two simulations of protoplanet-planetesimal interaction, one for a protoplanet with an atmosphere and the other for a no-atmosphere condition. The latter case is similar to the previous studies (Ida & Nakazawa 1989). Gravity of the sun was included in both cases. The mass of the core of the protoplanet is $M_{\mbox{\scriptsize core}} = 1~ M_{\oplus}$and the radius of the planetesimal is 100 m. When we consider the atmosphere of the protoplanet, the luminosity of the atmosphere is $L = 1 \times 10^{-8} ~L_{\odot}$and the grain depletion factor is f=0.01. In Fig. 5, the closest distances between the protoplanet and the planetesimal are shown. The solid horizontal line corresponds to the initial distance, $b/r_{\mbox{\scriptsize H}}$, of the planetesimal from the protoplanet that is at the origin. The planetesimal with the small (<1.9) and large (>2.6) initial distances are not considered because they cannot interact with the atmosphere. The dashed line corresponds to the closest distances between the planetesimal and the protoplanet that does not have the atmosphere, while the solid line gives the closest distances of the planetesimal from the protoplanet that has the atmosphere. If the closest distance is equal to the radius of the core that is given by the dotted line with $R_{\mbox{\scriptsize core}}$, the planetesimal collides with the core. Because of the gas drag of the atmosphere, the planetesimal lose the energy. The closest distances of some planetesimals from the center of the protoplanet become small for the case that the protoplanet has an atmosphere. When the protoplanet has an atmosphere, the planetesimals that have initial distances in the ranges of $2.04 {-} 2.17 r_{\mbox{\scriptsize H}}$ and $2.34 {-}2.43 r_{\mbox{\scriptsize H}}$lose energy and collide with the core. The boundaries of the two ranges are approximately determined by the closest distance for the case that the protoplanet does not have the atmosphere and the enhanced radius. Once the planetesimal approaches the protoplanet within the enhanced radius, it loses energy and is captured by the effective gas drag. Therefore, it is a good approximation to replace the radius of the core, $R_{\mbox{\scriptsize core}}$, with the enhanced radius, $R_{\mbox{\scriptsize c}}$, in order to obtain a collision rate between a protoplanet that has the atmosphere and a planetesimal.

The average number of collisions between a protoplanet and planetesimals per unit time, $\langle N_{\mbox{\scriptsize col}} \rangle_{Mm}$, is given by Nakazawa et al. (1989) as

 \begin{displaymath}\langle N_{\mbox{\scriptsize col}} \rangle_{Mm}= r_{\mbox{\sc... ...ptsize s}}(m) \langle P_{\mbox{\scriptsize col}} \rangle_{Mm}, \end{displaymath} (19)

where $n_{\mbox{\scriptsize s}}(m)$ is the surface number density of planetesimals with mass of m, and $\langle P_{\mbox{\scriptsize col}} \rangle_{Mm}$ is the non-dimensional mean collision rate. Inaba et al. (2001) obtained an approximate formula of $\langle P_{\mbox{\scriptsize col}} \rangle_{Mm}$ for the case without an atmosphere. In Fig. 5, it is shown that the replacement of $R_{\mbox{\scriptsize core}}$ with $R_{\mbox{\scriptsize c}}$ is a good approximation. According to Inaba et al. (2001), we propose a new formula of the mean collision rate when we consider the atmosphere of the protoplanet. For the high velocity regime ( $e^{*}A_{\mbox{\scriptsize p}}/r_{\mbox{\scriptsize H}} > 2$), where e* is the dispersion of the eccentricity of the planetesimals, the mean collision rate is expressed as

 \begin{displaymath}\langle P_{\mbox{\scriptsize col}} \rangle_{\mbox{\scriptsize... ...l G}(i^*/e^{*})} {(e^{*} A_{\mbox{\scriptsize p}})^2} \right), \end{displaymath} (20)

where i* is the dispersion of the inclination of the planetesimals. The functions ${\cal F}(x)$ and ${\cal G}(x)$ are given, respectively, by

  \begin{displaymath}{\cal F}(x) \equiv 8 \int^{1}_{0} {\rm d}\lambda \frac{x^{2}E... .../2\right]} {\left[x^{2}+\left(1-x^{2}\right)\lambda^2\right]}, \end{displaymath} (21)

where K and E are the complete elliptic integrals of the first and second kinds, respectively. While for the low velocity regime, ( $e^{*}A_{\mbox{\scriptsize p}}/r_{\mbox{\scriptsize H}} < 0.2$),

 \begin{displaymath}\langle P_{\mbox{\scriptsize col}} \rangle_{\mbox{\scriptsize... ...11.3 \sqrt{R_{\mbox{\scriptsize c}}/r_{\mbox{\scriptsize H}}}, \end{displaymath} (22)

and for the medium regime ( $0.2 < e^{*}A_{\mbox{\scriptsize p}}/r_{\mbox{\scriptsize H}} < 2$),

 \begin{displaymath}\langle P_{\mbox{\scriptsize col}} \rangle_{\mbox{\scriptsize... ...{\mbox{\scriptsize H}}}{R_{\mbox{\scriptsize c}}} \right)\cdot \end{displaymath} (23)

With the use of the above equations, the mean collision rate is given by

 \begin{displaymath}\langle P_{\mbox{\scriptsize col}} \rangle_{Mm} =\min\left[\l... ...}} \rangle_{\mbox{\scriptsize low}}^{-2}\right)^{-1/2}\right]. \end{displaymath} (24)

4 Conclusions and discussion

The presence of an atmosphere surrounding a protoplanet increases the radial distance over which planetesimals are captured by the protoplanet. Planetesimals approaching the protoplanet interact with the atmosphere. Their velocities are reduced by the gas drag. If they lose enough energy, they are captured in the atmosphere. The enhanced radius of the protoplanet is much larger than the radius of the core of the protoplanet. This distance reduces to the radius of the core if the protoplanet does not have the atmosphere. We call this distance the enhanced radius of the protoplanet.

The enhanced radius depends on the gas density of the atmosphere. The gas density increases with an increase in the mass of the core because of the large gravity. With increases in the luminosity of the atmosphere and/or the grain opacity of the atmosphere, the gas density decreases. Approximate solutions to the gas density of the atmosphere have been derived as a function of mass of the core, luminosity, and grain opacity of the atmosphere. The solutions reproduce the gas density obtained by the numerical simulations well. The enhanced radius also depends on the size of a planetesimal and increases with a decrease in the size of the planetesimal because gas drag is more effective for small planetesimals. We have derived the approximate solutions to the enhanced radius of a protoplanet as a function of the size of the planetesimal and the gas density. With the use of the approximate solutions, we have derived the collision rate between the protoplanet that has an atmosphere and the planetesimal.

We study the effect of ablation on the enhanced radius. By ablation, planetesimals lose their mass in the atmosphere. If the planetesimal lose mass, gas drag becomes more effective so that the protoplanet can capture the planetesimal. Even if we adopt the relatively large value of the ablation constant ( $0.1~ \mbox{s}^2~ \mbox{km}^{-2}$), only small amount of mass is lost and, as a result, the enhanced radius does not change. Furthermore, we examine if fragmentation by the dynamical pressure occurs while planetesimals are captured in the atmosphere. The dynamical pressure that the planetesimal experiences is smaller than the strength of planetesimals and the planetesimals are not fragmented.


  \begin{figure} \par\resizebox{8cm}{!}{\includegraphics[clip]{MS3608f6.eps}}\end{figure} Figure 6: The growth time, $T_{\mbox{\scriptsize g}}$, of the protoplanet with $10 ~M_{\oplus}$by accumulation of uniform planetesimals with radius of $r_{\mbox{\scriptsize p}}$without (dashed curve) and with (solid curve) the effect of enhanced radius.
Open with DEXTER

Based on the results obtained in this study, we simply estimate the growth time of the protoplanet with $10 ~M_{\oplus}$. With the use of Eq. (19), we have the growth rate of the protoplanet as

 \begin{displaymath}\frac{{\rm d}M}{{\rm d}t} = \Omega_{\mbox{\scriptsize K}} \Si... ...iptsize H}}^2 \langle P_{\mbox{\scriptsize col}} \rangle_{Mm}, \end{displaymath} (25)

where $\Sigma$ is the surface density of solid. As shown in Eq. (20), the collision rate, $\langle P_{\mbox{\scriptsize col}} \rangle_{Mm}$, is dependent on the random velocities (i.e., eccentricities) of planetesimals. When the protoplanet is large, its gravitational interactions with planetesimals determine the eccentricities of the planetesimals (Ida & Makino 1993):

 \begin{displaymath}e^{*} = 11 \left( \frac{A_{\mbox{\scriptsize p}}}{5.2~\mbox{A... ...ight)^{1/6} r_{\mbox{\scriptsize H}}/A_{\mbox{\scriptsize p}}. \end{displaymath} (26)

Using Eqs. (20), (25), and (26), the growth time of the protoplanet ( $T_{\mbox{\scriptsize g}}= M/({\rm d}M/{\rm d}t)$) is approximately given by

 \begin{displaymath}T_{\mbox{\scriptsize g}} = 1.5 \times 10^7 \left( \frac{A_{\m... ...{R_{\mbox{\scriptsize c}}} \right) \hspace{3ex} \mbox{years}. \end{displaymath} (27)

The growth time is dependent on the enhanced radius. The protoplanet with a large enhanced radius has a short growth time. As seen in Fig. 3c, the enhanced radius decreases with an increase in the luminosity. Assuming that the luminosity is given by the gravitational energy released by the incoming planetesimals, that is $L= G M/R_{\mbox{\scriptsize core}} {\rm d}M/{\rm d}t$, we have

 \begin{displaymath}L = 8.1 \times 10^{-8} \left( \frac{A_{\mbox{\scriptsize p}}}... ...{R_{\mbox{\scriptsize core}}} \right) \hspace{3ex} L_{\odot}. \end{displaymath} (28)

Therefore, the luminosity and the enhanced radius are dependent on each other. After several trials, we find the enhanced radius that does not contradict the luminosity. Using the obtained enhanced radius and Eq. (27), we have the growth time of the protoplanet with $10 ~M_{\oplus}$, as shown in Fig. 6. By the inclusion of the enhanced radius, the growth time reduces by factors of 3 and 10 when planetesimals with radii of $1 \times 10^8$ cm and $1\times 10^4$ cm, respectively, are accumulated by the protoplanet. If we consider that the protoplanet grows by accumulating planetesimals of $1\times 10^4$ cm radius, the growth time of the protoplanet becomes less than 105 years. This growth time is much smaller than the probable lifetime of a disk and might explain the formation of Jovian planets. However, we have to be very careful to have the conclusion because this estimate is too simple. This estimate does not include the size distribution of planetesimals or the decrease of the surface density of solid due to migration of planetesimals toward the central star by the gas drag of a circumstellar disk. The detailed calculations of formation of Jovian planets including the above effects will be studied in Inaba et al. (2003).

Acknowledgements
We would like to express our gratitude to G. W. Wetherill for extremely helpful discussions. We also thank M. Podolak for motivating this study, T. Fagan for a comment on this paper and H. Tanaka for discussion. S. Inaba and M. Ikoma received financial support from Research Fellowships of the Japan Society for Promotion of Science for Young Scientists (Nos. 06423 and 06405).

Appendix A: Approximate solutions for the atmospheric structure

We make additional assumptions to derive an approximate solution. First, we assume

\begin{displaymath}M_r = M_{\rm core} = \mbox{constant}, \end{displaymath} (A.1)

since the atmospheric mass is much smaller than the mass of the core before the critical core mass. From this assumption, Eq. (5) can be eliminated. Next, we use the equation of state for ideal gas given by

 \begin{displaymath}P = \frac{k}{\mu m_{\rm H}} \rho T \end{displaymath} (A.2)

with assumption of constant mean molecular weight, $\mu$ (k: Boltzmann constant; $m_{\rm H}$: mass of a hydrogen atom). Finally, we assume that the opacity is constant in radiative zones and the second adiabatic exponent is constant in convective zones; $\Gamma_2$ is equal to the usual adiabatic exponent $\gamma$ in the ideal approximation.

For convenience, we introduce the following dimensionless quantities:

\begin{displaymath}p \equiv \frac{P}{P_0}, \hspace{1em} \theta \equiv \frac{T}{... ...space{1em} \mbox{and} \hspace{1em} q \equiv \frac{r}{r_0}, \end{displaymath} (A.3)

where the quantities with suffices 0 are some reference values. Furthermore, we define two non-dimensional numbers such as

\begin{displaymath}V_0 \equiv \frac{GM_{\rm core} \rho_0}{r_0 P_0} \hspace{1em}... ...c{3 \kappa L P_0}{4 \pi a_{\rm r} c G M_{\rm core} T_0^4}\cdot \end{displaymath} (A.4)

Using these, we can rewrite Eq. (4) as

 \begin{displaymath}\frac{{\rm d}p}{{\rm d}q} = - V_0 \frac{\sigma}{q^2}, \end{displaymath} (A.5)

Eq. (6) as

 \begin{displaymath}\left\{ \begin{array}{ll} \displaystyle{ \frac{{\rm d} \th... ...rac{\sigma}{q^2} } & \mbox{otherwise}, \end{array} \right. \end{displaymath} (A.6)

and Eq. (A.2) as

 \begin{displaymath}p = \sigma \theta. \end{displaymath} (A.7)

First, we derive an analytical solution for a radiative atmosphere. Dividing the first equation of Eq. (A.6) by Eq. (A.5) and integrating the obtained equation, we have

 \begin{displaymath}\theta^4 = 1 + W_0 (p - 1). \end{displaymath} (A.8)

Inserting Eq. (A.7) into Eq. (A.5) and using Eq. (A.8) leads to

\begin{displaymath}-\frac{V_0}{q^2} {\rm d}q = \left\{4 + \frac{4 \left(1-W_0\right)}{\theta^4+W_0-1} \right\} {\rm d}\theta. \end{displaymath} (A.9)

This equation can be integrated analytically to be

 \begin{displaymath}\frac{1}{q} = 1+\frac{1}{V_0} \left\{ 4(\theta-1) + f(\theta; w_0) \right\}, \end{displaymath} (A.10)

where

\begin{displaymath}w_0 = \left\vert 1 - W_0 \right\vert^{1/4} \end{displaymath} (A.11)

and

\begin{displaymath}f(\theta; w_0) = \left\{ \begin{array}{ll} \displaystyle{ ... ...}w_0}{w_0^2-1}} \right) \right\} } & % \end{array} \right. \end{displaymath} (A.12)

(the upper and the lower equations for W0 < 1 and W0 > 1, respectively). It should be noted that this solution reduces to that derived by Mizuno (1980) and Stevenson (1982) for $\theta^4 \gg 1$ and $p \gg 1$, that is,
 
                  $\displaystyle \theta$ = $\displaystyle \frac{V_0}{4} \frac{1}{q},$ (A.13)
p = $\displaystyle \frac{1}{W_0} \left(\frac{V_0}{4}\right)^4 \frac{1}{q^4},$ (A.14)
$\displaystyle \sigma$ = $\displaystyle \frac{1}{W_0} \left(\frac{V_0}{4}\right)^3 \frac{1}{q^3}\cdot$ (A.15)

Next, we consider convective cases. From Eq. (A.5) and the second equation of Eq. (A.6), we obtain the familiar adiabatic relation such as

 \begin{displaymath}p = \sigma^\gamma. \end{displaymath} (A.16)

Substituting this in Eq. (A.5) and integrating the obtained differential equation, we have

 \begin{displaymath}\frac{1}{q} = 1 + \frac{\gamma}{(\gamma-1)V_0} \left( \sigma^{\gamma-1} - 1\right). \end{displaymath} (A.17)

In summary, we use Eqs. (A.7), (A.8), and (A.10) for radiative zones; we use Eqs. (A.7), (A.16), and (A.17) for convective zones.

We divide the atmosphere into three layers in order to reproduce the numerical solutions as accurately as possible using our approximate solutions. This requirement arises from the fact that the grain opacity drops abruptly because of evaporation (typically of ice and silicate), so that we cannot replace $\kappa$ with a constant value throughout the atmosphere. From our experience, we decide to calculate in the following way.

We choose the normalized temperature $\theta$ as the independent variable. The two boundaries are the points at which $T = 170 ~ \rm K$ and $T = 1700 ~ \rm K$. These temperatures are typical evaporation temperatures of ice and silicate. We calculate our approximate solutions inward from the outer boundary and change the reference values at the boundaries. That is,

1.
for $T_{\rm neb} < T < 170 ~ \rm K$, we choose

\begin{displaymath}T_0=T_{\rm neb}, \hspace{1em} P_0=P_{\rm neb}, \hspace{1em} ... ...} \mbox{and} \hspace{1em} \kappa = 4 f ~{\rm cm^2 ~ g^{-1}}, \end{displaymath}

2.
for $170 ~ {\rm K} < T < 1700 ~ \rm K$,

\begin{displaymath}T_0=170~{\rm K}, \hspace{1em} P_0=P_{(T~=~170~{\rm K})}, \hs... ...} \mbox{and} \hspace{1em} \kappa = 2 f ~{\rm cm^2 ~ g^{-1}}, \end{displaymath}

3.
for $1700 ~ {\rm K} < T$,

\begin{displaymath}T_0=1700~{\rm K}, \hspace{1em} P_0=P_{(T~=~1700~{\rm K})}, \... ... \mbox{and} \hspace{1em} \kappa = 0.01 ~{\rm cm^2 ~ g^{-1}}. \end{displaymath}

Only at each boundary, we judge whether radiative or convective solutions should be used. If $W_0 < 4(1-1/\gamma)$ ( $\gamma = 7/5$ here), then we use the radiative solution; otherwise, the convective solution.

In Fig. 1, we compare the above-obtained solutions with our numerical solutions described in the Sect. 2.2. As seen from these figures, our analytical solutions reproduce the numerical solutions well, except for the model with $M_{\rm core} = 10 ~M_\oplus$ in Fig. 1a. The difference between the analytical and numerical solutions for the model of $M_{\rm core} = 10 ~M_\oplus$ originates mainly from neglecting the contribution of the atmospheric mass to the planetary radius R (see Eq. (3)); the atmospheric mass is about 1/3 times as large as the core mass at the critical core mass (see Stevenson 1982). However, this deviation is not important in this study, because the information required here is not the density at a given radius, but the radius for a given density. In this sense, our analytical solutions coincide with the numerical ones within the accuracy of 50%.

References


Copyright ESO 2003