Gas-Giant Planet Formation

# Gas Accretion onto a Protoplanet and Formation of a Gas Giant Planet

Masahiro N. Machida11affiliation: Department of Physics, Graduate School of Science, Kyoto University, Sakyo-ku, Kyoto 606-8502, Japan; machidam@scphys.kyoto-u.ac.jp, inutsuka@tap.scphys.kyoto-u.ac.jp , Eiichiro Kokubo22affiliation: Division of Theoretical Astronomy, National Astronomical Observatory of Japan, Osawa, Mitaka, Tokyo 181-8588, Japan; kokubo@th.nao.ac.jp , Shu-ichiro Inutsuka11affiliation: Department of Physics, Graduate School of Science, Kyoto University, Sakyo-ku, Kyoto 606-8502, Japan; machidam@scphys.kyoto-u.ac.jp, inutsuka@tap.scphys.kyoto-u.ac.jp , and Tomoaki Matsumoto33affiliationmark:
###### Abstract

We investigate gas accretion onto a protoplanet, by considering the thermal effect of gas in three-dimensional hydrodynamical simulations, in which the wide region from a protoplanetary gas disk to a Jovian radius planet is resolved using the nested-grid method. We estimate the mass accretion rate and growth timescale of gas giant planets. The mass accretion rate increases with protoplanet mass for , while it becomes saturated or decreases for , where , and and are the Jovian mass and the orbital radius, respectively. This accretion rate is typically two orders of magnitude smaller than that in two-dimensional simulations. The growth timescale of a gas giant planet or the timescale of the gas accretion onto the protoplanet is about yr, that is two orders of magnitude shorter than the growth timescale of the solid core. The thermal effects barely affect the mass accretion rate because the gravitational energy dominates the thermal energy around the protoplanet. The mass accretion rate obtained in our local simulations agrees quantitatively well with those obtained in global simulations with coarser spatial resolution. The mass accretion rate is mainly determined by the protoplanet mass and the property of the protoplanetary disk. We find that the mass accretion rate is correctly calculated when the Hill or Bondi radius is sufficiently resolved. Using the oligarchic growth of protoplanets, we discuss the formation timescale of gas giant planets.

a
22affiliationtext: Faculty of Humanity and Environment, Hosei University, Fujimi, Chiyoda-ku, Tokyo 102-8160, Japan; matsu@i.hosei.ac.jp

ccretion, accretion disks — hydrodynamics — planetary systems —planets and satellites: formation— solar system: formation

## 1 Introduction

Currently, more than 400 exoplanets have been observed. Most such planets are considered to be gas giant planets. Although giant planets are preferentially observed, observations imply that gas giant planets, such as Jupiter and Saturn in our solar system, can be born around stars. Thus, it is important to understand the formation process of the gas giant planet. In general, gas planets are formed in the protoplanetary disk (or the circumstellar disk) around the protostar. However, there is a problems about the growth (or the gas accretion rate) of the gas planet, which is related to the resulting mass of a gas giant planet.

According to the core accretion scenario (Perri & Cameron, 1974; Mizuno et al., 1978; Hayashi et al., 1985), the protoplanet has a less massive hydrostatic gas envelope when a protoplanet core mass is less than where is the Earth mass, while the protoplanet captures a massive gas envelope from the protoplanetary disk, that is, a runaway gas accretion phase, to become a gas giant planet when (e.g., Mizuno et al., 1978; Mizuno, 1980; Stevenson, 1982; Bodenheimer & Pollack, 1986; Pollack et al., 1996; Ikoma et al., 2000, 2001; Hubickyj et al., 2005). Since the gas giant planet acquires almost all of its mass in the runaway accretion phase, the gas flow and mass accretion rate in this phase are important to determine the protoplanet evolution and its resulting mass. Since gas accretion onto the protoplanet may be closely related to the formation of the circumplanetary disk and acquisition process of angular momentum (Machida et al., 2008; Machida, 2009), three-dimensional calculations are required for investigating the runaway accretion phase. Machida (2009) showed that the circumplanetary disk is formed in a compact region near the protoplanet. Thus, we may have to resolve the present size of the gas giant planet to estimate the gas accretion rate onto the protoplanet. As well, to properly handle the outer boundary of the protoplanetary system (protoplanet and circumplanetary disk), the region sufficiently far from the gravitational sphere of the protoplanet, that is, the Hill radius, should also be included, since gas flows into the protoplanetary system from outside the Hill sphere. Therefore, there is a need to incorporate vastly different spatial scales, ranging from the planet current radius to the Hill radius. For example, the Hill radius of Jupiter ( cm) is about 700 times larger than Jupiter’s current radius ( cm).

So far, the gas accretion process onto the proto-gas giant planet was investigated mainly in global three-dimensional calculations (e.g., Kley et al., 2001; D’Angelo et al., 2002, 2003; Bate et al., 2003; Klahr & Kley, 2006; Dobbs-Dixon et al., 2007; Paardekooper & Mellema, 2008; Fouchet & Mayer, 2008; D’Angelo & Lubow, 2008). However, the fine structures in the proximity of the protoplanet cannot be resolved in such calculations. A compact disk is formed in the range of (Machida et al., 2008; Machida, 2009). Thus, the mass accretion rate may have to be derived by resolving the spatial resolution with at least . It should be noted that the mass accretion rate obtained from a global simulation can be correct, if the mass accretion rate is determined solely by the global structure around the Hill sphere and the circumplanetary disk contributes little to gas accretion. Even in such a case, the mass accretion rate should be investigated in calculations with higher spatial resolution to confirm the validity of the accretion rate derived in simulations with coarser spatial resolution. There are only a few studies that report the mass accretion rate onto the protoplanet with sufficiently high-spatial resolution in a local simulation (Tanigawa & Watanabe, 2002; Ayliffe & Bate, 2009). However, in these studies, the gas flow from the region outside the Hill radius to the proximity of the protoplanet was not sufficiently discussed. Furthermore, the acquisition process of the angular momentum and formation of the circumplanetary disk are not discussed in these studies.

In this paper, we focus on the mass accretion onto the protoplanet in runaway gas accretion phase, in which the gas continues to collapse onto the protoplanet without additional heating by collision of planetesimals. After the solid core (or planetary core) formation, the formation process of the gas giant planet can be divided into two phases. The quasi-static envelope slowly contracts with a timescale of yr (e.g., Ikoma et al., 2000) when the protoplanet is less massive than , while the gas rapidly collapses onto the protoplanet when . We investigate the evolution of the protoplanetary system only in the runaway gas accretion phase, based on the results of a local simulation with a sufficiently high-spatial resolution using the nested-grid method, in which the region of from the protoplanet is resolved with cells having the size of the radius of the present-day Jupiter. Since the acquisition process of the angular momentum (Machida et al., 2008) and circumplanetary disk formation (Machida, 2009) around the protoplanet were already investigated using this method, this paper will focus on gas accretion onto the protoplanet and gas flow in the proximity of the protoplanet. As a result of calculation, we found that the mass accretion rate is correctly calculated when the Hill or Bondi radius is sufficiently resolved, and it barely depends on the thermal effect around proto-gas giant planet. The structure of the paper is as follows. §2 gives the model frameworks, while §3 describes the numerical methods used. The numerical results are presented in §4 and compared with the results of previous studies in §5. §6 discusses protoplanetary growth. The conclusions of this paper are presented in §7.

## 2 Model

### 2.1 Basic Equations

A local region around a protoplanet is considered using the shearing sheet model (e.g., Goldreich & Lynden-Bell, 1965), in which the self-gravity of the protoplanetary disk is ignored. In addition, no physical viscosity is included, and the numerical viscosity can be ignored because it is sufficiently small. Thus, an inviscid gas disk model is adopted. The orbit of the protoplanet is assumed to be circular in the equatorial plane of the protoplanetary disk. Local rotating Cartesian coordinates with the origin at the protoplanet are set up, in which the -, -, and -axis are the radial, azimuthal, and vertical directions of the disk. The equations of hydrodynamics without self-gravity are solved [see, eqs. (1)-(6) of Machida et al. 2008].

For the gas, a barotropic equation of state is adopted (for details, see Machida 2009). In a local region, the protoplanetary disk has an almost constant temperature (Hayashi et al., 1985), while the gas around the protoplanet, that is, the gas envelope, has a higher temperature than the protoplanetary disk (Mizuno et al., 1978; Mizuno, 1980; Bodenheimer & Pollack, 1986; Pollack et al., 1996; Ikoma et al., 2000). Mizuno et al. (1978) studied the structure and stability of the envelope around the protoplanet, on the assumption that the envelope is spherically symmetric and in hydrostatic equilibrium. They also investigated the thermal evolution of the envelope, parameterizing the dust grain opacity, and determined the boundary between the isothermal and adiabatic regions. Using Figure 2 of Mizuno et al. (1978), the thermal evolution around the protoplanet is modeled as a function of the gas density, that is, using the barotropic equation of state, as

 P=c2s,0ρ[1−tanh(ρρcri)]+Kργtanh(ρρcri), (1)

where is the sound speed, is the adiabatic index (), and the adiabatic constant is defined as where is the critical density, wherein the gas behaves isothermally for , and adiabatically for . In this study, = (isothermal model), , 10, , and (adiabatic model) are used, where is the initial density on the equatorial plane. The hyperbolic tangent (tanh) function is used to smoothly connect the first (isothermal) and second (adiabatic) terms in Equation (1). The thermal evolution for different is plotted against the gas density in Figure 1, in which the gas temperature is constant in the isothermal model (), while it increases gradually from the initial value at in the adiabatic models ( = , 10, , and ).

In the standard disk model (Hayashi et al., 1985), at a Jovian orbit, the density is , and the temperature is  K. Thus, for example, in a model with , the gas behaves isothermally when , while it behaves adiabatically when . Comparing Figure 1 with Figure 2 of Mizuno et al. (1978), the thermal evolution of the model with (Fig. 1 broken line) corresponds to that for a gas envelope around a proto-Jovian planet (Fig.2 of Mizuno et al. 1978) with a dust opacity of  cm g. In this model setting, the critical density corresponds to the dust opacity in Mizuno et al. (1978). Mizuno et al. (1978) adopted  cm g as the most reliable parameter of a proto-Jovian planet, indicating that a more realistic gas temperature of the envelope is lower than that in the model with (the dotted line of Figure 1). In our settings, because  cm g almost corresponds to , we call the models having ‘most realistic models’. It should be noted that the critical density increases as the dust opacity decreases. Thus, in models with , the thermal energy around the protoplanet may be overestimated. On the other hand, when the isothermal equation of state is adopted, the thermal energy around the protoplanet is obviously underestimated. Therefore, it is expected that the actual thermal evolution is located between models with and .

### 2.2 Protoplanetary Disk Model and Scaling

The initial settings are similar to Miyoshi et al. (1999), Machida et al. (2006b), Machida et al. (2008), and Machida (2009). The gas flow has a constant shear in the -direction as

 \boldmathv0=(0,−3/2Ωpx,0), (2)

where is the Keplerian angular velocity of the protoplanet

 Ωp=(GMca3p)1/2, (3)

where is the gravitational constant, is the mass of the central star, and is the orbital radius of the protoplanet.

For hydrostatic equilibrium, the density is given by where () is the surface density of the unperturbed disk. The scale height is related to the sound speed by .

Using the standard solar nebular model (Hayashi, 1981; Hayashi et al., 1985), the temperature , sound speed , and gas density can be described as

 T=280(LL⊙)1/4(ap1AU)−1/2  K, (4)

where and are the protostellar and solar luminosities,

 cs=(kTμmH)1/2=1.9×104(T10K)1/2(2.34μ)1/2  cms−1, (5)

where is the mean molecular weight of the gas composed mainly of H and He, and

 ρ0=1.4×10−9(ap1AU)−11/4  gcm−3. (6)

When and are adopted, the scale height can be described as

 h=5.0×1011(ap1AU)5/4   cm. (7)

The inverse of the angular velocity is described as

 (8)

The basic equations can be normalized using unit time, , and unit length, . The density is also scalable and is normalized using . Hereafter, the normalized quantities are expressed with a tilde on top, for example, , , and . Further details can be found in Machida et al. (2008). A dimensionless description is given by Equations (14)-(19) of Machida (2009). The dimensionless quantities are converted into dimensional quantities using Equations (4)–(8). The gas flow is characterized by two parameters, the dimensionless Hill radius , and the critical density . The dimensional Hill radius is defined by

 rH=(Mp3Mc)1/3ap. (9)

In this study, a Hill radius ranging from 0.29 to 1.36 is used. As a function of the orbital radius and the mass of the central star, the parameter is related to the actual protoplanet mass in units of Jovian mass by

 MpMJup=0.12(Mc1M⊙)−1/2(ap1AU)3/4~r3H. (10)

For example, in the model with ,  AU and , the protoplanet mass is . In the parameter range of , at a Jovian orbit ( AU), protoplanets have masses of . Table 1 gives the dimensionless Hill () and Bondi () radii, the masses of protoplanets at Jovian (5.2 AU) orbit, the critical densities and sink radii for all models. Model names consist of two parts: the protoplanet mass at the Jovian orbit and the critical density. For example, model M001A3 has parameter values given as , .

## 3 Numerical Methods

### 3.1 Numerical Procedures

The purpose of this study is to investigate gas accretion onto the protoplanet in three-dimensional simulations. However, given current CPU limitations, it is impossible to calculate the complete evolution of the gas giant planet with sufficiently high-spatial resolution, that is, the evolution of the planet from a solid core () with a thin gas envelope to a protoplanet that acquires a massive atmosphere () cannot be performed. Thus, the mass accretion rate is derived using the following procedure: (i) With a fixed protoplanet mass, the evolution of the protoplanetary system until the gas flow reaches the steady state ( orbits) is calculated, (ii) Then a sink is introduced, and the mass accretion rate that is derived is averaged over a further  orbits, and (iii) Finally, steps (i) and (ii) are repeated by changing the protoplanet mass to give the growth rate of the protoplanet.

### 3.2 Nested-Grid Method

To investigate the formation of a circumplanetary disk in a protoplanetary disk, it is necessary to cover a large dynamic spatial scale range. Using the nested-grid method (for details, see Machida et al., 2005, 2006a), the regions near and remote from the protoplanet are covered with adequate resolution. Each level of a rectangular grid has the same number of cells (), but cell width depends on the grid level . The cell width is reduced by 1/2 with increasing grid level (). Eight grid levels () are used. The box size of the coarsest grid, , is , and that of the finest grid, , is . The cell width in the coarsest grid, , is , and it decreases with as the grid level increases. Thus, the finest grid has . The fixed boundary condition in the - and -direction, and the periodic boundary condition in the -direction are used.

### 3.3 Sink Cell and Smoothing Length

In the finest grid (), the cell width is . In real units, when the protoplanet is located at 5.2 AU, the cell width corresponds to  cm, or 0.8 times the Jovian radius. The evolution of the protoplanetary system is calculated using a sink. In the fiducial models, the radius of the sink is or twice the Jovian radius at Jovian orbit (for details, see §3.4). During the calculation, the gas from the region inside the sink radius is removed in each time step.

The smoothing length for the gravitational potential of the protoplanet is not explicitly used. For numerical calculations, the physical quantities are defined at the cell centre, while the origin (protoplanet’s position) is defined as the cell boundary. Thus, the region inside has a uniform gravitational potential. At a Jovian orbit, is times the Jovian radius ( or  cm).

### 3.4 Convergence Test for the Accretion Rate

To check the convergence of the mass accretion rate onto the protoplanet, the evolution of the protoplanetary system for different sink radii was calculated. The convergence of other quantities that change with the cell width or grid level with and without a sink were already investigated in Machida et al. (2008) and Machida (2009). Figure 2 shows the mass accretion rate as a function of the sink radius for models M02ISS, M02ISM, M02I, and M02ISL. The mass accretion rate is derived based on the procedure outlined in §3.1. As listed in Table 1, these models have the same Hill radius (or the same protoplanetary mass) and different sink radii . Models have at Jovian orbit ( AU), and the isothermal equation of state is used. Model M02ISL has the smallest sink radius of which corresponds to cm (0.87 ), while model M02ISL has the largest sink radius of which corresponds to cm (6 ). The protoplanetary system was calculated for  orbits, for which it was confirmed that the gas flow achieves a steady state in  orbits. Tanigawa & Watanabe (2002) also showed that steady state is reached after  orbits for a local calculation. Machida (2009) showed that the gas envelope and angular momentum around the protoplanet reach steady state after 1-10 orbits.

Figure 2 shows that the models M02ISS, M02ISM, and M02I have almost the same accretion rate of , for models, which have sink radii of (cm; 1.92 ). On the other hand, model M02ISL has a sink radius of (cm; 6 ), and has a slightly larger accretion rate of than the others. Thus, the difference in the accretion rate between M02ISS (0.2) and M02ISL (0.35) is only within a factor of 2. As a result, the accretion rate can safely be estimated when the sink radius is less than . In the following, we show results adopting (cm; ).

Since the gas inside the sink radius is removed at each time step as mentioned in §3.3, the structure, that is, the protoplanet’s atmosphere and envelope, inside the sink radius cannot be investigated. However, it is expected that the gas falling into the sink cannot escape from the protoplanet and, hence, does not influence the results obtained. In the region where , since the gravitational energy greatly dominates the thermal energy (Mizuno et al., 1978; Machida, 2009), the gas that falls into the sink cannot be pushed out by the gas pressure. In addition, the gas falling into the sink does not acquire additional angular momentum inside the sink radius (or near the protoplanet), because the angular momentum of the system is only acquired by the shearing motion of the protoplanetary disk (Machida et al., 2008). However, it is unknown how the gas trapped by the gravitational potential of the protoplanet reaches the surface of the protoplanet because the centrifugal force is expected to prevent the gas from further falling inside the sink. To understand the surface and structure of the protoplanet, the region inside the protoplanet, that is, inside needs to be resolved. Such simulation would require an enormous amount of CPU time. In the calculations, the sink radius is that is smaller than the Roche limit . Thus, it is expected that the gas inside loses its angular momentum by tidal interaction and finally falls into the gas planet. In addition, Klahr & Kley (2006) suggested that the radius of a young planet is about twice the current size of Jupiter. Therefore, when the sink radius is comparable to the size of the present planet, it is possible to safely estimate the gas accretion rate.

## 4 Gas Flow and Mass Accretion Rate

### 4.1 Spiral Pattern and a Circumplanetary Disk

Figure 3 shows the structure around the Hill sphere for each model at orbits (). The systems in Figure 3 are in steady state, because the gas flow achieves a steady state in  orbits () as shown in §3.4 (see also Miyoshi et al., 1999; Tanigawa & Watanabe, 2002; Machida et al., 2008). Models in the figure have a parameter of , for which the gas temperature increases adiabatically for . The figures show that the global structure in the models is almost the same as for the isothermal models (for details, see Machida, 2009). In calculations using the isothermal equation of state, the spiral patterns that are distributed from the upper-left to the lower-right region are observed (e.g., Lubow et al., 1999). The same patterns are seen also in Figure 3, in which the barotropic equation of state is used. This can be explained by noting that for regions far from the protoplanet, the gas density is not high enough, and the gas behaves isothermally. As well, Figure 3 shows that as the protoplanet mass increases, stronger shock waves appear and the spiral patterns become clearer. Furthermore, the density contrast between spiral and gap becomes stronger as the protoplanet mass increases. Thus, it can be concluded that the features observed in adiabatic models correspond well with those features observed in isothermal models.

To focus on the region near the protoplanet, the structures inside the Hill sphere are shown in Figure 4 in three dimensions. The grid level of Figure 4 is , while that of Figure 3 is . Thus, Figure 4 is a 16 times enlargement of the central part of Figure 3. In Figure 4, the region of is plotted using the red constant density surfaces, while that for is plotted using the orange constant density surfaces. It should be noted that the orange and red surfaces do not appear in model M001A2, because this model has a very small part that has . The density distributions on , , and plane are projected onto each wall surface. Figure 4 shows that the red and orange surfaces increase with the protoplanet mass, indicating that a massive protoplanet has a denser, or more massive, envelope. The orange surface flattens as the protoplanet mass increases, while the red surface in all models except for model M001A2 has a sufficiently flattened structure. Machida et al. (2008) and Machida (2009) showed that the specific angular momentum of the envelope increases in proportional to when the protoplanet mass is smaller than . Thus, the angular momentum also increases with the protoplanet mass. Therefore, the centrifugal radius increases with the protoplanet mass. Since the centrifugal force gradually affects the more distant region from the protoplanet as the protoplanet mass increases, the low-density region gradually flattens out owing to the centrifugal barrier. On the other hand, the red surfaces always have flattened structures, because the red surface (or higher density envelope) is located near the protoplanet that is inside the centrifugal radius even when the protoplanet is less massive .

### 4.2 Mass Accretion Rate onto a Protoplanet

Figure 5 shows the mass accretion rates onto the protoplanet for all models against the cube of the Hill radius, . With a fixed protoplanet’s orbit, can be connected to the protoplanet mass using Equation (10). As a reference, the mass accretion rate and the protoplanet mass at a Jovian orbit are given in the upper and right axes in Figure 5. The mass accretion rate and the protoplanet mass at the Jovian orbit can be converted into those at any orbit using the value of in parenthesis in the upper and right axes. In this paper, the mass accretion rate at the Jovian orbit will be used for convenience.

As shown in Figure 5, the mass accretion rates increases rapidly for protoplanetary mass in the range of . The accretion rates for isothermal and adiabatic models are almost the same in this range. The accretion rate for each model has a peak at , then it gradually decreases in the range of . As well, in this range, the accretion rates depend slightly on the equation of state used. The difference increases as the protoplanet mass increases. For , the decrease of the mass accretion rate in models with a harder equation of state, that is smaller , is less than in models with a softer equation of state, that is, larger . For example, the barotropic model of has accretion rates of  yr at , and  yr at . On the other hand, the isothermal model shows a steeper decrease of the accretion rate than the other models. The isothermal model had accretion rates of  yr at , and  yr at . Thus, the difference between isothermal and barotropic models is not so large even when the protoplanets become massive. When the protoplanet has a Jovian mass , the accretion rate in model with is only 2.4 times larger than that for the isothermal model.

As shown in §2.1, it is expected that the actual mass accretion rate is located between the isothermal and the adiabatic model with . Around the protoplanet, the thermal energy is underestimated in the isothermal model, while it is overestimated in the model with . We fitted the mass accretion rate by the solid line in Figure 5, in which, for convenience, it is fitted as a constant in the range of . It should be noted that, although the accretion rate in each model gradually decreases in this range, the rate of decrease is sufficiently small especially in adiabatic models. As well, it can be noted that, since the protoplanetary evolution in the range of was not calculated, the accretion rate may not apply in the range of . Using dimensionless quantities, the mass accretion rate is described as

 d~Mpd~t≃{[]ll0.83(~r3H)3/2 for ~r3H<0.30.14          for ~r3H>0.3. (11)

In Equation (11), the mass accretion rate for () almost corresponds to that of the average () in most realistic models having for . Since the protoplanetary evolution was calculated using a local simulation, it is possible to convert the mass accretion rate from the fixed orbit to any arbitrary orbit, arbitrary density, and temperature of the protoplanetary disk model. Using Equations (4)-(8), the dimensional mass accretion rate can be described as

 dMpdt[MJupyr]≃ ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩1.2×10−2(MpMJup)3/2(ap1AU)−13/8(LL⊙)−3/16(McM⊙)−1/4for Mp/MJup<0.036(ap1AU)3/4(McM⊙)−1/2(LL⊙)3/88.1×10−5(ap1AU)−1/2(LL⊙)3/8(McM⊙)−1for  Mp/MJup>0.036(ap1AU)3/4(McM⊙)−1/2(LL⊙)3/8. (12)

When the central protostar has a mass of , a luminosity of , and the standard model (i.e., MMSN disk model; Hayashi et al., 1985) with surface density is adopted in Equation (12), the accretion rate at any orbit can be described as

 dMpdt[MJupyr]≃ ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩1.2×10−2(MpMJup)3/2(ap1AU)−13/8(ΣΣ0)for (MpMJup)<0.036(ap1AU)3/48.1×10−5(ap1AU)−1/2(ΣΣ0)for (MpMJup)>0.036(ap1AU)3/4, (13)

where is arbitrary surface density at any orbit. Equation (13) indicates that the growth timescale of the Jovian mass protoplanet at a Jovian orbit is  yr. We will present a detailed discussion of the mass accretion rate and the growth rate of the protoplanet in §5. The following sections show the gas flows at large (§4.3) and small scales (§4.4) around the protoplanet to confirm the relationship between the accretion rate and the circumplanetary disk.

### 4.3 Gas inside the Hill Sphere

This section examines the large-scale gas flows in order to determine gas accretion from outside the Hill sphere. Each top panel in Figure 6 plots the gas streamlines falling into the protoplanet. Figure shows only models with , which is the most realistic parameter (Mizuno et al., 1978; Machida, 2009). There are no differences in the streamlines between models with different equations of state, because the gas outside the Hill sphere, that is, at a large-scale, has a low density and behaves almost isothermally. In the figure, the streamlines are inversely integrated from the region inside the Hill sphere. Thus, only the streamlines accreting onto the protoplanetary system are plotted. Colors of the streamlines indicate the velocity in the vertical direction (). For example, the gas rises upward in blue and green parts of the streamline, while it goes down in yellow and red parts. The Bondi, which is shown using a black circle,

 rB=GMpc2s,0, (14)

the doublewide Bondi, shown with a black-dotted circle, the Hill, shown with a white circle, and the doublewide Hill, shown with a white-dotted circle radius, are also plotted in the figure. Each top panel in Figure 6 shows that the gas flow falling into the protoplanetary system is controlled either by the Bondi or the Hill radius. When the Bondi radius is smaller than the Hill radius (models M001A2 and M005A2), the gas flows into the protoplanet only in the narrow band of . The bandwidth of streamlines in the region far from the Hill sphere almost corresponds to the Bondi radius (, see the lower left corner in the top panel of Fig. 6). When the protoplanet mass increases and the Bondi radius exceeds the Hill radius (M04A2, M08A2, and M1A2), the gas flows into the protoplanet in the band of , in which the streamlines in the region far from the Hill sphere have a bandwidth of the Hill radius (). It should be noted that in Figure 6, the grid-scale is different in each of the panels. Thus, it can be concluded that the gas flows into the protoplanetary system occur only in the region given as

 racc≲x≲3racc,    racc=min(rB,rH). (15)

It was also shown by D’Angelo & Lubow (2008) that the mass accretion rate is controlled by either the Bondi or the Hill radius.

The colors of the streamlines imply that the gas rises slightly in regions far from the protoplanet (see, for example, Models M04A2, M08A2, M1A2), then it gradually goes down as the streamlines approach the Bondi or Hill sphere. Near the Bondi or Hill sphere, the gas rises strongly upward near the shock front and then rapidly falls into the Bondi or Hill spheres in the vertical direction (for details, see Machida et al., 2008). The similar features of the gas flow (or streamlines) around the Hill sphere are seen in Kley et al. (2001), D’Angelo et al. (2003) and Paardekooper & Mellema (2008).

For each bottom panel in Figure 6, the plane is plotted, in which the color indicates the mass flux per unit area. In the plane, only gas passing through the gray region (or gray points) can reach the protoplanetary system (compare the gray regions in each bottom panel with the streamlines in each top panel). Most of the gray points are distributed in the region far from the equatorial plane, indicating that the gas present in the midair flows preferentially into the protoplanetary system, while the gas near the equatorial plane barely flows into it. Although the gas near the equatorial plane also flows into the Hill sphere, it flows out from the Hill sphere afterwards, as shown in Machida et al. (2008) and Paardekooper & Mellema (2008). Figure 7 shows the streamlines for model M04A2, in which the streamlines are integrated from the area of and on the plane, that is, the region denoted by the black squares. It should be noted that the streamlines in Figure 6 are inversely integrated from the region inside the Hill sphere, while, in Figure 7, they are integrated along the flow from inside the black square. The color of each streamline indicates the mass flux per unit volume at each mesh point. The gas flow through the black square in Figure 7 tracks three different paths: (i) the flow is attracted toward the protoplanet and returns to a similar orbit as the original orbit (the pass-by region), (ii) the flow is trapped by the protoplanet and falls into the protoplanetary system (the planetary atmosphere region), and (iii) the flow turns round inside the Hill sphere and goes back out (the horseshoe region). Lubow et al. (1999) showed that the flow distributed in a certain band (see, Figure 4 of Lubow et al. 1999) accretes onto the protoplanetary system in a two-dimensional simulation (see also, Tanigawa & Watanabe, 2002). However, in three-dimensional simulations, only a part of the flow in the band (or bundle) in the radial direction flows into the protoplanetary system, as shown in Figures 6 and 7. Kley et al. (2001) and Paardekooper & Mellema (2008) pointed out that the mass accretion rate derived in three-dimensional calculations is much smaller than that derived from two-dimensional calculations.

As shown in Figure 6, the bandwidth that limits the gas flowing into the protoplanetary system is proportional to min (, ). The Bondi radius is proportional to (eq. 14), while the Hill radius is proportional to (eq. 9). When the Bondi radius is smaller than the Hill radius (), the bandwidth () and accretion rate rapidly increase with the protoplanet mass because the mass dependence of the Bondi radius is strong (). On the other hand, when , the bandwidth () and accretion rate maintain an almost constant value, because the mass dependence of the Hill radius is considerably weaker (). It should be noted that the Bondi radius equals the Hill radius when . Thus, as shown in Figure 5, it is possible to qualitatively understand the protoplanet mass dependence on the accretion rate: the accretion rate rapidly increases when , while it maintains an almost constant value when . However, the power of the accretion rate ( for , see Fig. 5) derived in our calculations is slightly smaller than the Bondi solution (Bondi, 1952):

 dMpdt=4πG2M2pρ0c3s,0∝M2p. (16)

This discrepancy can be attributed to the fact that spherically symmetric flow in an isolated system is assumed in Bondi accretion. However, in the protoplanetary system, the gravitational sphere is limited to the Hill sphere, and the angular momentum that was acquired from the shearing motion in the protoplanetary disk can affect the gas accretion. Thus, it is natural that the power of the accretion rate in the protoplanetary disk does not completely match with that obtained from the Bondi solution. In fact, the accretion rate derived in this simulation is two or three orders of magnitude smaller than the Bondi accretion rate. As has been mentioned, the Bondi radius is useful in representing the gas flow when the protoplanet is less massive, indicating that the thermal pressure cannot be ignored when the thermal energy dominates the gravitational energy inside the Hill sphere. Note that, in the isothermal model, the thermal energy dominates the gravitational energy inside the Hill sphere when (for details see Machida et al., 2008).

### 4.4 Gas Flow around a Protoplanet

In this subsection, we consider the gas flow patterns near the protoplanet. Figure 8 shows the density distributions and streamlines in the neighborhood around the protoplanet for models M005A2 (left) and M1A2 (right). The streamlines are inversely integrated from the high-density region (yellow constant density surface). The left panel in Figure 8 shows that the gas spirals into the protoplanet in the vertical direction. On the other hand, in the right panel, a part of the gas accretes onto the protoplanet, while the remainder accretes onto the circumplanetary disk. The gas accreting onto the circumplanetary disk orbits in the disk and gradually falls into the protoplanet. The density contours in wall surfaces show a thick torus-like configuration of the disk in model M005A2, and a very thin disk in model M1A2.

Figure 9 shows the mass flux () for model M005A2, in which only mesh points having are plotted in color. Thus, the gas in the black region has a positive radial velocity and flows out from the protoplanetary system. Figures 9a and b indicate that a large part of the gas enters into the protoplanetary system from the vertical direction around the Hill sphere. However, the gas does not enter into the protoplanetary system on the equatorial plane. Figure 9c shows that, even when the gas has a negative radial velocity inside the Hill sphere (), it cannot enter into the region on the equatorial plane: the flow turns round and returns (see, for example, Figure 7). Therefore, although the gas enters the Hill sphere, it cannot reach the protoplanet on the equatorial plane. Paardekooper & Mellema (2008) also showed that the gas flows out from the Hill sphere on the equatorial plane (the equatorial outflow).

Machida et al. (2008) showed that the specific angular momentum of the protoplanetary system increases in proportion to . They also showed that the circumplanetary disk appears when the gravitational energy exceeds the thermal energy in the whole region of the Hill sphere, which occurs when at the Jovian orbit. Thus, when the protoplanet mass exceeds at the Jovian orbit, the circumplanetary disk appears, and then it increases its size and mass with the protoplanet mass because the massive protoplanet acquires a larger amount of angular momentum from the shearing motion in the protoplanetary disk. One can expect that the emergence of the circumplanetary disk may reduce the accretion rate. However, we confirmed that the circumplanetary disk keeps a constant mass, indicating that the circumplanetary disk does not strongly affect the gas accretion rate onto the protoplanet. In summary, although the flow pattern around a massive protoplanet is different from that around a less massive protoplanet, fine structures such as the circumstellar disk around the protoplanet does not strongly affect the accretion rate onto the protoplanet. Thus, once the gas is gravitationally captured by the protoplanet, it accretes onto the protoplanet.

### 4.5 Thermal Effects on the Mass Accretion Rate

As shown in Figure 5, when the protoplanet is less massive (), the accretion rates for models with different equations of state have almost the same value. On the other hand, when the protoplanet mass exceeds , accretion rates are slightly different for models with different equations of state. This slight difference is considered to be caused by thermal evolution around the protoplanet. The deviation from an isothermal approximation is small in a less massive system, while it becomes slight large in a massive system, as shown in Machida (2009). The density distribution around the protoplanet for models M1I (left) and M1A1 (right) is shown in Figure 10, in which the circumplanetary disks are represented by constant density surfaces. The figure indicates that the isothermal model (model M1I) has a thinner disk than the adiabatic model (model M1A1). The orange constant density surface of shows that model M1A1 has a torus-like disk that flared up outwardly, while model M1A1 has a thick disk. The red constant density surface of shows that, near the protoplanet, model M1A1 has a more compact disk than model M1I.

Figure 11 illustrates the mass of the gas envelope (or the circumplanetary disk), which is integrated from the centre. The figure indicates that a model with a harder equation of state has a more massive envelope. Since the circumstellar disk in these models is sufficiently gravitationally stable (Toomre’s Q ), the gas in the circumplanetary disk barely accretes onto the protoplanet with dynamic instability. However, as shown in Figure 10, the model with large thermal energy, that is, a harder equation of state, has a more massively flared disk. In such a disk, the effect of centrifugal and gravity forces is relatively small, because the thermal energy is relatively large. Therefore, it is likely that the accretion rate in such a disk becomes higher than that in a model with lower thermal energy in which the circumstellar disk is strongly supported by the centrifugal force, because the path of the gas streamlines accreting onto the protoplanet is changed. Finally, it follows that the massive disk has a larger accretion rate.

## 5 Comparison with Previous Simulations

The mass accretion rate onto the protoplanet (system) in the protoplanetary disk was investigated in some previous studies (Kley et al., 2001; Tanigawa & Watanabe, 2002; D’Angelo et al., 2003; Bate et al., 2003; Klahr & Kley, 2006; Dobbs-Dixon et al., 2007; Paardekooper & Mellema, 2008; Fouchet & Mayer, 2008; D’Angelo & Lubow, 2008; Ayliffe & Bate, 2009). Subsections §5.1 to 5.3 will explain the classification of the previous simulations by classifying them into three categories and discussing the salient features of each category.

### 5.1 Two-dimensional Isothermal Calculations

Tanigawa & Watanabe (2002) derived the mass accretion rate under the isothermal gas condition in a two-dimensional local simulation. The spatial resolution of their calculation is comparable to our study. However, the mass accretion rate is greatly different from the results obtained in this study. In Tanigawa & Watanabe (2002), the mass accretion rate continues to increase as a function of the protoplanet mass, while it saturates at in our three-dimensional calculation. Furthermore, although in this study, the mass accretion rate increases with the protoplanet mass in the range of , its value is much smaller than in Tanigawa & Watanabe (2002). For example, at , the accretion rate in Tanigawa & Watanabe (2002) is  yr (see, Equation [19] of Tanigawa & Watanabe 2002), while it is  yr in our isothermal model. Thus, in the range of , the accretion rate based on a two-dimensional calculation is two orders of magnitude larger than that obtained from a three-dimensional calculation. The difference becomes larger in the range of , because the accretion rate continues to increase in two-dimensional calculations, while it gradually decreases in the three-dimensional calculations. Kley et al. (2001) compared the mass accretion rate in their three-dimensional simulations with that in two-dimensional simulations (Kley, 1999) and showed that the mass accretion rate derived in three-dimensional simulations is smaller than that in two-dimensional simulations. Paardekooper & Mellema (2008) also commented on the decrease in mass accretion observed in three-dimensional simulations compared with two-dimensional simulations. The difference in mass accretion is considered to be caused by the spatial dimensions.

### 5.2 Three-dimensional Isothermal Calculations

Kley et al. (2001), D’Angelo et al. (2003), and Bate et al. (2003) estimated the mass accretion rate under the locally isothermal assumption in global simulations. The mass accretion rates derived from these global simulations corresponds well to the rate derived from our local simulations, despite the fact that the spatial resolution and size of the sink are considerably different. For example, D’Angelo et al. (2003) estimated the mass accretion rate in their global three-dimensional simulations with a spatial resolution of , which is 40 times coarser than ours, and defined the sink as the region inside one-tenth of the Hill radius (), which is 28 times larger than our sink size (see §3.3). Thus, our study can resolve the present Jovian radius, while D’Angelo et al. (2003) cannot resolve the radius. As well, in their calculation, the circumplanetary disk cannot be sufficiently resolved because the circumplanetary disk is formed in the region (Machida, 2009). In summary, the differences between D’Angelo et al. (2003) and our study are the spatial resolution and numerical setting. We calculated the evolution of the protoplanetary system in the local simulation with finer spatial resolution, while this cannot be done in D’Angelo et al. (2003)’s global simulation with coarser spatial resolution. Despite the differences, the accretion rate in our isothermal models corresponds well to their results (compare Figure 5 with Figure 7 of D’Angelo et al. 2003). The mass accretion rate has a peak value of  yr at in D’Angelo et al. (2003), while it has a peak value of  yr at in our study. Furthermore, the mass accretion rates at and are  yr and  yr in D’Angelo et al. (2003), while they are  yr and  yr in our study. It should be noted that although D’Angelo et al. (2003) investigated the evolution of the protoplanetary disk by changing the gravitational potential of the protoplanet, the mass accretion rates differed slightly in each model.

Kley et al. (2001) investigated the evolution of the protoplanetary disk only when the protoplanet has a Jovian mass () under the isothermal gas condition in their three-dimensional simulations, in which the spatial resolution is coarser than D’Angelo et al. (2003). They found that a mass accretion rate of  yr, while it is  yr at in our study (see, for example, Fig. 5). Thus, our results correspond well with their results. Bate et al. (2003) also calculated the evolution of the protoplanetary disk in a three-dimensional global simulation with a spatial resolution comparable to that of D’Angelo et al. (2003). The accretion rate derived in Bate et al. (2003) corresponds well to ours within a factor of three. In addition, it shows the same trends as in our study. That is, the accretion rate has a peak at and decreases for .

In summary, all three-dimensional calculations show almost the same mass accretion rate, though they have adopt different calculation settings (i.e., local or global calculation), spatial resolutions and sizes of the sink radius. In these calculations, the Hill radius or the region near the Hill sphere is resolved, while the region near the protoplanet, much smaller scale than the Hill radius is not always resolved sufficiently. This indicates that the mass accretion rate is regulated by the region around the Hill sphere, not by the details at smaller scale length. Thus, we conclude that the mass accretion rate is safely estimated when the Hill radius is resolved with adequate spatial resolution.

Recently, Klahr & Kley (2006), Paardekooper & Mellema (2008), Fouchet & Mayer (2008) and Ayliffe & Bate (2009) investigated the evolution of the protoplanetary disk using three-dimensional, radiation-hydrodynamical simulations. Klahr & Kley (2006) estimated the mass accretion rate of  yr when the protoplanet has a Jovian mass. They concluded that the mass accretion rate in a radiative model is larger than that in an isothermal model. Paardekooper & Mellema (2008) calculated the evolution of the protoplanetary disk with several models with different protoplanetary mass under the locally isothermal approximation, while they calculated it with the radiation-hydrodynamical simulations only with protoplanet mass of and . The accretion rates in their isothermal calculation are identical to our results. As well, the accretion rate in their radiation-hydrodynamical simulations is also comparable to our results. It should be noted that, since the range of the protoplanetary mass in their radiation-hydrodynamical simulations is different from our study, a direct comparison of the results cannot be performed. They concluded that although the accretion rate in the radiation-hydrodynamical simulation is smaller than that in the isothermal simulation, the difference is not dramatic. It should be noted that, for radiation-hydrodynamical simulations, Paardekooper & Mellema (2008) only calculated the disk evolution with very small protoplanetary masses of and . These masses are considerably different from that used by Klahr & Kley (2006).

Fouchet & Mayer (2008) investigated the evolution of the protoplanet disk including both the self-gravity and radiation physics for a protoplanet with a Jovian mass, . They also calculated the protoplanetary disk under the isothermal approximation and compared them with the radiation-hydrodynamical simulation. They obtained an accretion rate of  yr in the isothermal simulation and  yr in the radiation-hydrodynamical simulation. These rates are well identical to Klahr & Kley (2006) and our results of  yr at . Ayliffe & Bate (2009) also estimated the mass accretion rate in their radiation-hydrodynamical simulation and found an accretion rate similar to those derived in previous studies.

However, the radiative effect on the mass accretion rate is controversial. Paardekooper & Mellema (2008) and Fouchet & Mayer (2008) showed that the accretion rate in a nonisothermal disk is reduced compared to the isothermal disk (see also, Ayliffe & Bate, 2009), while Klahr & Kley (2006) showed the mass accretion rate in radiative model is larger than in the isothermal model. In our study, the accretion rate for the isothermal model is smaller than for the adiabatic models. However, the difference in the accretion rate between the isothermal and radiative models is not so large. This can be attributed to the gravitational energy of the protoplanet dominating the thermal energy (Machida, 2009). Thus, it is expected that the accretion rate barely depends on the thermal evolution in the protoplanetary disk. For example, when the Jovian mass planet is adopted at Jovian orbit, the accretion rate is in the range  yr (Kley et al., 2001; D’Angelo et al., 2003; Bate et al., 2003; Klahr & Kley, 2006; Fouchet & Mayer, 2008; D’Angelo & Lubow, 2008). This indicates that the growth timescale of the Jovian planet at Jovian orbit is   years.

### 5.4 Disk Viscosity and Effect of the Gap on the Mass Accretion Rate

In this subsection, we discuss the relation between gap formation and the mass accretion rate onto the protoplanet. Some previous studies pointed out that the reduction of the mass accretion rate for massive protoplanets is related to the gap formation (e.g., D’Angelo et al., 2003; Bate et al., 2003). The density gap in the protoplanetary disk is the result of the competition between torques exerted on the disk by the planet and by the disk itself. The planet gives angular momentum to the outer part of the disk, and it takes angular momentum from the inner part of the disk (Goldreich & Tremaine, 1980). Thus, a protoplanet tends to open a gap. On the other hand, the disk kinematic viscosity makes the gas back to the gap. As a result, the gap formation depends on both the disk kinematic viscosity and mass of the protoplanet (Lin & Papaloizou, 1986; Bryden et al., 1999). A massive protoplanet in the disk with a smaller viscosity tends to show a clear gap. However, the formation condition, size, and width of the gap are not clearly understood.

In the inviscid () gas disk model as adopted in this study, it is considered that the gap formation occurs when the Hill radius exceeds the scale height of the protoplanetary disk, i.e., (Lin & Papaloizou, 1986; Bryden et al., 1999; Crida et al., 2006). Using equations (7) and (9), this condition can be described as

 MpMJup≳0.12(ap1AU)3/4, (17)

which indicates that, at Jovian orbit, the density gap begins to appear when the mass of the protoplanet exceeds . As shown in §4.3, we found that the gas flows into the planetary system in the region when the protoplanet is more massive than at Jovian orbit. Thus, when the gap width becomes larger than , the mass accretion rate onto the protoplanet is expected to decline. The decline of the mass accretion rate at in global simulations (e.g., D’Angelo et al., 2003) is caused by a wide gap formation. On the other hand, our local simulation shows no clear gap, though the low-density region is shallower than that in the global simulations. Local simulations are not appropriate to treat gap formation because of the radial boundary condition (Miyoshi et al., 1999; Tanigawa & Watanabe, 2002), since the gap properties, such as the gap depth and width, depend on the size of the simulation box. Therefore, when the mass of protoplanets exceeds at Jovian orbit, it is expected that the mass accretion rate derived in the local simulation is smaller than that in the global simulation owing to the gap. On the other hand, in the range of , the mass accretion rate derived in the local simulation is always applicable, because no clear gap forms in this mass range. In addition, the mass accretion rate derived in this study would be valid even for when the disk viscosity is sufficiently large, because no clear gap appears in such a viscous disk. Moreover, the mass accretion rate in the local simulation well agrees with that in the global simulation in the range of as denoted in §5.2. This may be because the density gap barely affects the gas accretion in this mass range. However, global simulations show a rapid decline of the mass accretion rate for , in which the density gap seems to significantly affect the mass accretion. In such a case, the mass accretion rate derived in the local simulation cannot be applicable.

In this study, the mass accretion rate is derived in the local simulation. This accretion rate may be overestimated when the protoplanet is sufficiently massive and the protoplanetary disk has a sufficiently small viscosity, because the gap formation decreases the mass accretion rate. Thus, the mass accretion rate obtained here corresponds to the possible maximum value that is attainable without forming a density gap in the disk. However, even when the density gap is formed, the formula of the mass accretion rate is applicable with the reduced disk surface density owing to the gap, because we investigated the gas flow in a dimensionless form, and all physical quantities are scalable, as denoted in Equation (13).

## 6 Growth Timescale of Gas Giant Planets

As mentioned in §1, to form a gas giant planet, the solid core with a mass of needs to acquire the gas component from the protoplanetary disk. At the end of the “oligarchic growth” stage (Kokubo & Ida, 1998), several protoplanets composed of the solid component, which will be referred to as the planetary core, appear with an orbital separation of about 10 Hill radii (Kokubo & Ida, 2002). Only when a planetary core acquires sufficient gas from the protoplanetary disk before the dissipation of the protoplanetary gas disk, the gas giant planet is possible to form. The mass and growth timescale of the planetary core were investigated in detail using -body simulations. In the remainder of this section, using the mass accretion rate onto the gas giant planet obtained in this study, as well as the mass and growth timescale of the planet core obtained from -body simulations (Kokubo & Ida, 2002), an investigation of the final mass and growth timescale of the gas giant planet is performed.

Kokubo & Ida (2002) derived an (isolation) mass and growth timescale of the planetary core as a function of disk model. Adopting (the orbital separation normalized by the mutual Hill radius of the cores), (the power index of the radial surface density distribution), and (the mass of the central star) in Equations (17) and (26) of Kokubo & Ida (2002), the isolation mass and growth timescale of the planetary core can be described as

 Miso≃0.16(ficeΣ110)3/2(ap1AU)2M⊕, (18)
 τcore≃3.2×105f−1/2ice(Σ110)−9/10(ap1AU)59/20  yr, (19)

where (ice factor) and (reference solid surface density at 1 AU) can be related to the solid surface mass density as

 Σsolid=ficeΣ1(ap1AU)−3/2   g cm−3. (20)

The ice factor expresses the increase of solid by ice condensation over the snow-line and (=2.7 AU) and 4.2 (). As the reference solid surface density, we adopt the minimum-mass disk model (Hayashi et al., 1985) with . In the following, the solid surface density is parameterized as = 1/2, 1 and 2. We set the mass of the planetary core as the isolation mass .

In a conventional scenario of a gas planet formation (e.g., Mizuno, 1980; Stevenson, 1982; Bodenheimer & Pollack, 1986), the planetary core is heated by the continuous accretion of planetesimals, and has a hydrostatic gas envelope when a planetary core is less massive than . Then, the planetary core cannot sustain the gas envelope and the rapid gas accretion is triggered (i.e., the runaway gas accretion phase) after it grows up to . However, the current planetary accretion theory (e.g., Kokubo & Ida, 1998, 2002) indicates that the accretion of planetesimals almost halts after the isolation core formation (or after the oligarchic growth stage). In this case, the planetary core cannot have a hydrostatic envelope owing to the absence of the heat source even when the planetary core is . Therefore, the runaway gas accretion is triggered just after the isolation core formation.

Figure 12a shows the mass of the planetary core for different based on Equation (18). In the figure, the sudden rise of the core mass at 2.7 AU is caused by the increase of solid component at the snow-line. In the standard model (, the red line), the planetary core has a mass of at the Jovian orbit ( AU) and at the Saturnian orbit ( AU).

Figure 12b shows the growth timescales of the planetary core (thin solid lines, Equation [19]) and Jovian-mass gas planet (thick solid lines) in each orbit. Using Equation (12), the mass accretion rate is integrated as a function of time until the planetary core reaches the Jovian mass in each orbit using

 MJup=Mcore+∫τgas0dMpdtdt,  Mp(t=0)=Mcore, (21)

where , and are used in Equation (12). Thus, represents the time required to form a Jovian-mass gas planet after the formation of the isolation core. The figure indicates that in any orbit. For example, in the standard model, that is, , the growth timescale of the planetary core is  yr at the Jovian orbit and  yr at the Saturnian orbit, while that of a Jovian gas planet is  yr at the Jovian orbit and  yr at the Saturnian orbit. Therefore, the gas-accretion timescale is two or three orders of magnitude shorter than the core growth timescale.

As shown in Kokubo & Ida (2002), gas giant planets can form in the limited range of the protoplanetary disk where the planetary core is large enough to capture the gas component before the dissipation of the gas disk. The solid curves of Figures 12c and d show the planet mass for the disk lifetime of   yr (c) and  yr (d). In these figures, using Equation (12), the planet mass at each orbit is derived by

 Mp=Mcore+∫tgasdMpdtdt,  Mp(t=0)=Mcore. (22)

The dotted lines of Figures 12c and d indicate the maximum mass in the ring of the protoplanetary disk around the planet (Hayashi et al., 1985), which is described as

 Mmax=2πapδrΣg, (23)

where is the ring width and is expressed by , where is the parameter representing the ring width. Adopting and in equation (9), can be described as Thus, Equations (23) gives

 Mmax=(4πa2pfHΣg)3/2(3M⊙)−1/2. (24)

In equation (24), is adopted as the gas surface density of protoplanetary disk, where is related to the solid surface density as because the constant ratio of the gas to solid surface density is assumed. The fiducial gas surface density is given by the standard disk model (Hayashi et al., 1985) as

 ΣHg=1.7×103(ap1AU)−3/2  gcm−2. (25)

Equation (24) is the mass of the ring in the protoplanetary disk in the range of , where the origin of is the position of the planet. As shown in §4.3, our calculation indicates that the gas in the range of flows into the protoplanet system. Thus, is adopted in equation (24). Since the protoplanet cannot acquire the gas component exceeding the mass in this ring, the mass means the maximum mass of the gas component acquired by the protoplanet at each orbit (for details see, Hayashi et al., 1985).

The thick gray curves of Figures 12c and d are the planet mass attainable at each orbit. Figure 12c indicates that the planet mass is determined by the dotted line (i.e., eq. [24]) for  AU, while by the solid curve (i.e., eq. [22]) for  AU. For AU, since the growth timescale of the gas planet is shorter than the lifetime of the gas disk, the planetary core accumulates all gas in the ring with in  yr. On the other hand, for  AU, since the gas disk dissipates before the protoplanet accumulates all gas in the ring, the planet mass is determined by equation (22). As shown in equation (12), since the mass accretion rate is a decreasing function of the orbital radius , the planet mass decreases with for  AU. The thick gray curve in Figure 12c shows that the planet, respectively, has almost the Jovian and Saturnian mass at Jovian and Saturnian orbit. For  AU, the mass of the gas component is less massive and comparable to the solid core mass (). In addition, the gas planet formed at  AU has a core mass of . Note that Saumon & Guillot (2004) studied the internal structure of Jupiter and Saturn, and estimated the core mass of the Jupiter in the range of , and core mass of the Saturn in the range of . These features (mass of gas and solid components at each orbit) well correspond to those of giant planets in our solar system (Jupiter, Saturn, Uranus, and Neptune).

Although we assumed that the gas disk dissipates in  yr after the isolation core formation in Figure 12c, observations indicate that the protoplanetary disk has a lifetime of Myr (Haisch et al., 2001; Hartmann, 2005; Silverstone et al., 2006; Flaherty & Muzerolle, 2008). Figure 12d shows the planet mass when the gas disk exists for  yr after the isolation core formation. In this figure, the planet mass is determined only by the dotted line (i.e., the maximum mass in the ring, eq. [24]). As shown in Figure 12b, since the growth timescale of a gas planet is shorter than  yr at any orbit, the protoplanet can accumulate all gas component around it when the gas disk exists for  yr. The thick gray line in Figure 12d shows that although the Jovian mass planet is formed at Jovian orbit, more massive planets than Jovian mass are formed at Saturnian, Uranian, and Neptunian orbits, which contradicts our solar system.

Figures 12c and d indicate that the gas giant planet with the core mass of can be formed in the orbital range of  AU. They show that the mass of planets exceeds Jovian mass in the range of  AU when the gas disk exists for  yr after the isolation core formation. Thus, to realize our solar system in the present model, it is expected that the gas component of the protoplanetary disk needs to dissipate just after the isolation core formation.

We need to pay attention to viscous evolution of the disk for estimating the planet mass. In Figures 12c and d, the planet mass in the range of  AU for disk with  yr and in whole orbital range for disk with  yr is determined by the disk surface density at the location of the planetary core formation when the gas accretion timescale is shorter than the disk lifetime. In other words, the planet mass is determined by in situ surface density (or in situ disk mass). However, in these figures, the disk viscous evolution is ignored. If the gap is refilled by the viscous evolution, the planet can acquire the gas component exceeding . Thus, the planet mass becomes more massive than gray line in these figure, when the disk has a significantly large viscosity.

At last, we comment on the hydrostatic envelope. As shown in the above, to estimate the formation time of a gas giant planet, we assumed no planetesimal accretion onto the planetary core after the isolation core formation. Although this assumption is supported by recent numerical simulations, the formation timescale of a gas giant planet may be considerably different when the planetesimal accretion does not completely halt even after the isolation core formation. With the planetesimal accretion, the planetary core can sustain a hydrostatic envelope for  yr at the maximum (Ikoma et al., 2000). In this situation, the formation time (or gas accretion time) for the gas giant planet is prolonged. Thus, in Figure 12b gives the minimum time necessary for the formation of gas giant planets.

## 7 Summary and Discussion

In this study, a calculation of gas accretion onto a planetary core including the thermal effect was performed, in which the mass accretion rate and growth timescale of a gas giant planet were estimated. l,

The mass accretion rate derived in this study quantitatively agrees with those derived in previous simulations. It is surprising that the same value for the accretion rate is derived from three-dimensional simulations with different parameter values, including parameters such as spatial resolution, the size of the sink, and the treatment of the thermal effect, including locally isothermal approximation, barotropic equation of state, and radiative hydrodynamics. This suggests that the mass accretion rate is only determined by ‘the protoplanet mass’ and ‘the properties of the protoplanetary disk.’ In our calculation, we resolved the Jovian radius of the cell width. On the other hand, for example, Kley et al. (2001) derived the mass accretion rate by using the sink radius of . In §4.4, we showed that almost all gas falling into the half Hill radius () can accrete onto the surface of the protoplanet. Although the radiative effect can affect the mass accretion rate, the difference is less than a factor of . As a result, this study, like previous studies, shows that the mass accretion rate can be accurately estimated for resolution higher than half Hill radius. In other words, the high spatial resolution that resolve Jovian radius in not necessary to investigate the mass accretion rate and growth times scale of gas giant plants.

A comparison of our results with the -body simulations for the solid core aggregation shows that the gas giant planet can form in the orbital range of in a short timescale of  yr after the planetary core inside the gas giant planet is formed in a longer timescale of  yr. In each orbit, the growth timescale of the gas giant planet is two orders of magnitude shorter than the core aggregation timescale. This may indicate that, in our solar system, the planetary core with began to capture the gas component to form gas giant planets such as Jupiter and Saturn just before the dissipation of the protoplanetary disk was completed. Otherwise, the protoplanet accumulates all gas around it before the gas dissipation when the gas disk exists for  yr. In the latter case, the mass of the gas giant planet may be determined only by disk properties (i.e., the surface density and density profile of the protoplanetary disk). Thus, to investigate the formation of the gas giant planet, we need to precisely estimate the properties of the protoplanetary disk, and the gas dissipation process (or timescale).

In this study, it was not possible to investigate gap formation in detail, since the simulations were limited in a local region around the protoplanet. However, it was possible to resolve the present Jovian radius and the circumplanetary disk. To understand gap formation and the mass accretion rate, it is necessary to simulate a protoplanetary system using a global simulation that can resolve the protoplanet. However, such a calculation requires a large amount of CPU time. As well, issues concerning disk viscosity would need to be better understood before further investigation of both the ‘protoplanetary’ and the ‘circumplanetary disk’ can be performed. In general, it is considered that disk viscosity is caused by the magneto-rotational-instability (MRI). Thus, to properly calculate the viscous disk (circumplanetary and circumstellar disks), it is necessary to include magnetic effects. Although it is currently possible to calculate the evolution of the magnetized disk with a lower spatial resolution (Machida et al., 2006b), a simulation with a higher spatial resolution (or more CPU power) is necessary to investigate the disk viscosity correctly. For such a simulation, the next generation of supercomputers is required. However, the mass accretion rate derived in the present local simulation is consistent with those obtained in global simulations. Thus, we could link the mass accretion rate in a global simulation to that in a local simulation. Therefore, we can safely estimate the growth of the gaseous protoplanet.

## Acknowledgments

We have greatly benefited from discussion with S. Watanabe, H. Tanaka, T. Tanigawa, and T. Muto. We are very grateful to an anonymous referee for a number of very useful suggestions and comments. Numerical calculations were carried out using a Fujitsu VPP5000 at the Center for Computational Astrophysics, the National Astronomical Observatory of Japan. This work is supported by the Grants-in-Aid from MEXT (18740104, 20540238, 21740136).

## References

• Ayliffe & Bate (2009) Ayliffe, B. A., & Bate, M. R. 2009, MNRAS, 393, 49
• Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214
• Bate et al. (2003) Bate, M. R., Lubow, S. H., Ogilvie, G. I., & Miller, K. A. 2003, MNRAS, 341, 213
• Bodenheimer & Pollack (1986) Bodenheimer, P., & Pollack, J. B. 1986, Icarus, 67, 391
• Bondi (1952) Bondi, H. 1952, MNRAS, 112, 195
• Bryden et al. (1999) Bryden, G., Chen, X., Lin, D. N. C., Nelson, R. P., & Papaloizou, J. C. B. 1999, ApJ, 514, 344
• Crida et al. (2006) Crida, A., Morbidelli, A., & Masset, F. 2006, Icarus, 181, 587
• D’Angelo et al. (2002) D’Angelo, G., Henning, T., & Kley, W. 2002, A&A, 385, 647
• D’Angelo et al. (2003b) D’Angelo, G., Henning, T., & Kley, W. 2003b, ApJ, 599, 548
• D’Angelo et al. (2003) D’Angelo, G., Kley, W., & Henning, T. 2003, ApJ, 586, 540
• D’Angelo & Lubow (2008) D’Angelo, G., & Lubow, S. H. 2008, ApJ, 685, 560
• Dobbs-Dixon et al. (2007) Dobbs-Dixon, I., Li, S. L., & Lin, D. N. C. 2007, ApJ, 660, 791
• Flaherty & Muzerolle (2008) Flaherty, K. M., & Muzerolle, J. 2008, AJ, 135, 966
• Fouchet & Mayer (2008) Fouchet, L., & Mayer, L. 2008, arXiv:0806.3975
• Goldreich & Tremaine (1980) Goldreich, P., & Tremaine, S. 1980, ApJ, 241, 425
• Goldreich & Lynden-Bell (1965) Goldreich, P., & Lynden-Bell, D. 1965, MNRAS, 130, 125
• Haisch et al. (2001) Haisch, K. E., Jr., Lada, E. A., & Lada, C. J. 2001, ApJl, 553, L153
• Hartmann (2005) Hartmann, L. 2005, The Nature and Evolution of Disks Around Hot Stars, 337, 3
• Hayashi (1981) Hayashi, C. 1981, Prog. Theor. Phys. Suppl., 70, 35
• Hayashi et al. (1985) Hayashi, C., Nakazawa, K., & Nakagawa, Y. 1985, in Protostars and Planets II, ed. D. C. Black & M. S. Matthews (Tucson: Univ. Arizona Press), 1100
• Hubickyj et al. (2005) Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2005, Icarus, 179, 415
• Ikoma et al. (2000) Ikoma, M., Nakazawa, K., & Emori, H. 2000, ApJ, 537, 1013
• Ikoma et al. (2001) Ikoma, M., Emori, H., & Nakazawa, K. 2001, ApJ, 553, 999
• Klahr & Kley (2006) Klahr, H., & Kley, W. 2006, A&A, 445, 747
• Kley (1999) Kley, W. 1999, MNRAS, 303, 696
• Kley et al. (2001) Kley, W., D’Angelo, G., & Henning, T. 2001, ApJ, 547, 457
• Kokubo & Ida (1998) Kokubo, E., & Ida, S. 1998, Icarus, 131, 171
• Kokubo & Ida (2002) Kokubo, E., & Ida, S. 2002, ApJ, 581, 666
• Lin & Papaloizou (1993) Lin, D. N. C., & Papaloizou, J. C. B. 1993, Protostars and Planets III, 749
• Lin & Papaloizou (1986) Lin, D. N. C., & Papaloizou, J. 1986, ApJ, 309, 846
• Lubow et al. (1999) Lubow, S. H., Seibert, M., & Artymowicz, P. 1999, ApJ, 526, 1001
• Machida et al. (2005) Machida, M. N., Matsumoto, T., Tomisaka, K., & Hanawa, T. 2005, MNRAS, 362, 369
• Machida et al. (2006a) Machida, M. N., Matsumoto, T., Hanawa, T., & Tomisaka, K. 2006a, ApJ, 645, 1227
• Machida et al. (2006b) Machida, M. N., Inutsuka, S., & Matsumoto, T., 2006b, ApJl, 649, L129
• Machida et al. (2008) Machida, M. N., Kokubo, E., Inutsuka, S.-i., & Matsumoto, T. 2008, ApJ, 685, 1220
• Machida (2009) Machida, M. N. 2009, MNRAS, 392, 514
• Miyoshi et al. (1999) Miyoshi, K., Takeuchi, T., Tanaka, H., & Ida, S. 1999, ApJ, 516, 451
• Mizuno et al. (1978) Mizuno, H., Nakazawa, K., & Hayashi, C. 1978, Progress of Theoretical Physics, 60, 699
• Mizuno (1980) Mizuno, H. 1980, Prog. Theor. Phys., 64, 544
• Papaloizou et al. (2007) Papaloizou, J. C. B., Nelson, R. P., Kley, W., Masset, F. S., & Artymowicz, P. 2007, Protostars and Planets V, 655
• Paardekooper & Mellema (2006) Paardekooper, S.-J., & Mellema, G. 2006, A&A, 459, L17
• Paardekooper & Mellema (2008) Paardekooper, S.-J., & Mellema, G. 2008, A&A, 478, 245
• Perri & Cameron (1974) Perri, F., & Cameron, A. G. W. 1974, Icarus, 22, 416
• Pollack et al. (1996) Pollack, J. B., Hubickyj, O., Bodenheimer, P., Lissauer, J. J., Podolak, M., & Greenzweig, Y. 1996, Icarus, 124, 62
• Saumon & Guillot (2004) Saumon, D., & Guillot, T. 2004, ApJ, 609, 1170
• Silverstone et al. (2006) Silverstone, M. D., et al. 2006, ApJ, 639, 1138
• Stevenson (1982) Stevenson, D. J. 1982, Planet. Space Sci., 30, 755
• Tanigawa & Watanabe (2002) Tanigawa, T., & Watanabe, S. 2002, ApJ, 580, 506