Accretion onto Circumplanetary Disks

Distribution of Accreting Gas and Angular Momentum onto Circumplanetary Disks


We investigate gas accretion flow onto a circumplanetary disk from a protoplanetary disk in detail by using high-resolution three-dimensional nested-grid hydrodynamic simulations, in order to provide a basis of formation processes of satellites around giant planets. Based on detailed analyses of gas accretion flow, we find that most of gas accretion onto circumplanetary disks occurs nearly vertically toward the disk surface from high altitude, which generates a shock surface at several scale heights of the circumplanetary disk. The gas that has passed through the shock surface moves inward because its specific angular momentum is smaller than that of the local Keplerian rotation, while gas near the midplane in the protoplanetary disk cannot accrete to the circumplanetary disk. Gas near the midplane within the planet’s Hill sphere spirals outward and escapes from the Hill sphere through the two Lagrangian points L and L. We also analyze fluxes of accreting mass and angular momentum in detail and find that the distributions of the fluxes onto the disk surface are well described by power-law functions and that a large fraction of gas accretion occurs at the outer region of the disk, i.e., at about 0.1 times the Hill radius. The nature of power-law functions indicates that, other than the outer edge, there is no specific radius where gas accretion is concentrated. These source functions of mass and angular momentum in the circumplanetary disk would provide us with useful constraints on the structure and evolution of the circumplanetary disk, which is important for satellite formation.

Subject headings:
hydrodynamics — methods: numerical — planets and satellites: formation — protoplanetary disks — shock waves

1. Introduction

Satellite systems around the giant planets in our solar system are commonly seen. They are thought to have formed in circumplanetary disks, which are believed to have existed around giant planets during their gas capturing growing stage.

In earlier works, formation process of satellite systems have been considered based on a minimum mass subnebula (MMSN) model, in which satellites form from a disk that contains sufficient solid mass with solar composition for reproducing the current satellite systems (e.g., Lunine & Stevenson (1982)), as an analog of the minimum mass solar nebula model (Hayashi, 1981). However, it was suggested that the MMSN model has difficulty in reproducing current satellite systems around Jupiter and Saturn (Canup & Ward, 2002). One of severe problems is that the model leads to much higher temperature than that of HO ice sublimation at the current regular satellite region, which means that ice, which is the main component of the satellites, cannot be used as building material of the satellites.

In order to overcome the difficulties of the MMSN-type models which assumes a closed and static disk, alternative models have been developed. Canup & Ward (2002) proposed a model in which an accretion disk with continuous supply of gas and solid is considered as a proto-satellite disk. This model is based on results of hydrodynamic simulations of gas capturing process of giant planets (e.g., Lubow et al., 1999; D’Angelo et al., 2002), which demonstrated that gas accretion from the protoplanetary disk toward the parent planets occurs through a circumplanetary disk. In this model, the surface density and temperature of the circumplanetary disk is kept lower than assumed in the MMSN model, and thus HO ice can be used as a solid building material of satellites. Based on comparison among time scales of different processes, they concluded that formation of the Galilean satellites can be best explained if the circumplanetary disk had one order of magnitude lower gas surface density than the MMSN disk, with slow gas accretion rate corresponding to Jupiter’s growth time scale longer than a few yr (see also Sasaki et al. (2010)). On the other hand, Mosqueira & Estrada (2003) proposed another disk model which consists of two components, i.e., inner MMSN-type massive disk and outer low-density extended disk. The model reproduced, for example, three inner Galilean satellites as well as only partially differentiated Callisto. Although these models seem to reproduce the current satellite systems of the giant planets in our solar system, they needed to assume parameters for the disk structure such as surface density profile as a basis of physical processes of satellite formation.

The structure of a circumplanetary disk is closely related to the gas accretion process of giant planets. Gas flow around proto-giant-planets in protoplanetary disks has been studied using hydrodynamic simulations. Earlier two-dimensional simulations showed that the gas in the Hill sphere rotates in the prograde direction (Miki, 1982; Sekiya et al., 1987; Korycansky & Papaloizou, 1996) and a pair of spiral shocks stands in the circumplanetary disks (e.g., Kley, 1999; Lubow et al., 1999; Tanigawa & Watanabe, 2002; D’Angelo et al., 2002). However, the scale height of the proto-planetary disk for Jupiter-sized planets is comparable to the Hill radius, thus simulations with two-dimension approximation cannot capture the feature of the accretion flow. Recent three-dimensional hydrodynamic simulations revealed that the two-dimensional picture for circumplanetary disks is not appropriate for the flow in the Hill sphere (e.g., D’Angelo et al., 2003; Bate et al., 2003; Machida et al., 2008; Paardekooper & Mellema, 2008; Ayliffe & Bate, 2009b; Coradini et al., 2010). One of the most important features newly found in three-dimensional calculations is that, in the region of circumplanetary disks, the gas accretes nearly vertically downward toward the midplane from high altitude.

Although some studies using hydrodynamic simulations were able to obtain the structures of circumplanetary disks, the direct use of results of hydrodynamic simulations for the structure, such as surface density, is problematic. One of the reasons is that disk-like shear-dominant flow is susceptible to numerical viscosity, which is intrinsic nature of numerical hydrodynamic calculation, no matter if SPH method or mesh based method is used, thus it is difficult to obtain reliable structure of the disk directly from such simulations. In particular, orbital radius of the rotating gas in disks basically changes only slightly with time, thus numerical error tends to accumulate easily in simulation of long-term evolution. Also timescale required for low viscous disks such as circumplanetary disks to reach steady state is much longer than the typical dynamical time of the fluid. While high-resolution hydrodynamic simulations is needed to resolve the structure of circumplanetary disks, their long-term evolution is difficult to follow with such time-consuming simulation. In addition, physical (non-numerical) viscosity in the disk is not well understood, and we need to assume specific viscosity models that are hard to justify.

As an alternative approach, in the present work, we examine gas accretion flow onto circumplanetary disks from proto-planetary disks, in order to determine gas accretion rate as a function of distance from parent planets. Unlike the structure in a rotation disk, the accretion flow onto circumplanetary disks is not susceptible to numerical viscosity. Also, because the circumplanetary disks are located at the downstream of super-sonic accretion flow, the accretion flow itself is hardly affected by the circumplanetary disk structure, which depends on poorly-known effective viscosity. In §2, we describe our settings of hydrodynamic simulations. In §3, results and analyses of the simulations are shown. We discuss implication of our analyses of the accretion flow in §4. We summarize our results in §5.

2. Methods

2.1. Settings and Basic Equations

We consider a situation in which a growing giant planet embedded in a protoplanetary disk has induced the nucleated instability and the gas of the disk accretes dynamically onto the planet. We take local Cartesian coordinates rotating with the planet, -axis corresponding to the radial direction, -axis the revolving direction of the planet, and -axis normal to the disk midplane. The planet is located at the origin.

We adopt the local approximation, in which tidal potential is linearized and curvature of the protoplanetary disk is neglected. This approximation is valid as long as the Hill radius is much smaller than the orbital radius of the planet. The orbit of the planet is assumed to be fixed circular and co-planar with the disk midplane. The gas is assumed to be inviscid and isothermal. Magnetic field and self-gravity of the gas are neglected.

We employ equation of continuity, equation of motion for compressive inviscid gas, and equation of state for isothermal gas to simulate the gas motion. In order to understand physics clearly, we use these equations in a non-dimensional form. We normalize time by the inverse of the Keplerian angular velocity , length by scale height , and mass by unperturbed surface density of the protoplanetary disk gas times , where is the gravitational constant, is the solar mass, is the semi-major axis of the planet, and is the sound speed of the disk gas on the planet orbit. The sound speed is unity in our normalization. Normalized quantities are denoted with tilde on the variables in this paper. The normalized equations can be written as




In the above is distance from the origin, is normalized Hill radius, , is the planet mass, and is a unit vector in the -direction. The third term of the right hand side in Eq. (4) is added so that at the two Lagrangian points L and L (i.e., ). The normalized quantities can be written, for example, as , , , and for velocity, density, pressure, and potential, respectively. Full equation set before the normalization can be found in Machida et al. (2008). In the normalized system, is the only characteristic physical parameter of the system, if we do not consider physical radius of the planet. In this paper, we focus on the case with , which roughly corresponds to Jupiter mass at 5AU.

2.2. Numerical Modeling

We employ a three-dimensional nested-grid hydrodynamic simulation code (e.g., Machida et al., 2005, 2006), which was originally developed to explore star formation process by collapse of molecular cloud core (Matsumoto & Hanawa, 2003). Size of the whole computational domain is (24,24,6). The domain has a symmetry about plane, and the covered region in simulation is , , and . We set 11 levels for the nested-grid system, and the numbers of grids in each level are . The level of the nested grid is denoted by and is the largest grid level. Increment of by one reduces the size of the computational domain to half in all directions, keeping its center at the planet. The finest grid size is thus , which corresponds to about 1/4 of the present Jovian radius at 5.2 AU.

We adopt two types of artificial gravitational weakening, which allows us to avoid non-physical crash of the computation without essentially changing results. One is the widely-used weakening for region adjacent to the planet to avoid singularity of the planet gravity. The modified potential is given by


where is so-called smoothing length and we set , which corresponds to two grid size of the finest grid level (). The other gravitational weakening is also applied for the -component of the tidal force in the high- region. In this region, tidal force in the -direction (toward the mid-plane) is very strong, and causes too large density gradient to be described by the grid size we set in the code. Thus we connect two parabolic curves smoothly at as


where we set . Initial density profile in the -direction is modified accordingly as


for hydrostatic equilibrium under the potential given by Eq.(8). This weakening does not cause any significant effect for the region where we need to observe because the mass above is negligible (0.0063 % of the total mass for the hydrostatic structure).

As boundary conditions, we set the flow at to be the unperturbed one, and the mirror condition is applied at and . As for the boundaries at and (inflow region), we use mixed boundary conditions as follows. When , we set unperturbed condition, which is the same as the initial condition, and change it gradually to periodic boundary condition until based on linear interpolation with respect to time. When , periodic boundary condition is applied to all region at . For the initial condition, we adopt unperturbed density and velocity; , for velocity and , for density, which is identical to that of Machida et al. (2008).

In order to obtain high-resolution results efficiently, calculations are started with only the largest grid level and the number of levels is increased with time. The second largest level () starts at , at , and the higher levels () are started in order from by considering relaxation degree in each level. In this paper, we will show a snapshot of the flow of , when the accretion flow is nearly in the steady state.

We set sink cells around the origin in order to see the pure effect of gas accretion flow without the effect of the planet body, as well as to mimic gas accretion phase onto the planet that follows the nucleated instability (e.g., Mizuno, 1980; Bodenheimer & Pollack, 1986; Ikoma et al., 2000). In the sink cells, gas is removed at a rate that corresponds to . We set .

3. Results

3.1. Path to Circumplanetary Disks

First, we examine the overall structure of accreting gas flow onto a circumplanetary disk. Figure 1 shows streamlines from outside of the Hill sphere (i.e., from a protoplanetary disk) toward the planet. The starting point of a streamline is defined by . Four panels show streamlines starting from four different height; 0.0, 0.5, 1.0, 1.5. Based on the destination of streamlines, we divide the flow into three regions: passing region where gas approaches the planet and passes by it without making U-turn, horseshoe region where gas crosses the planet orbit with U-turn during the encounter, and accreting region where gas is accreted onto the circumplanetary disk. We find that there are no accreting streamlines on the midplane (), while some streamlines starting from off-midplane region accretes to the circumplanetary disk. This implies that the gas near the midplane in proto-planetary disks is harder to accrete to circumplanetary disks and planets. This is confirmed by Figure 2, which shows streamlines starting from three different radii 0.05, 0.1, 0.2, that is, . We can see that gas on the midplane in spirals outward and escapes from the Hill sphere within a short timescale through one of the two Lagrangian points (L or L), and that the outward radial velocity decreases with decreasing distance from the planet (see also Fig. 7).

The vertical structure of the flow is more clearly demonstrated in Figure 3. In this figure, initial positions of streamlines on the plane at are classified into three regions by their destinations. We define the passing region where streamlines reach the boundary at with , horseshoe region where streamlines reach the boundary at with , and accreting region where streamlines end up within a sphere with radius ; we set . We can clearly see that there is no accretion band in the midplane, while off-midplane site has an accretion band with a significant width. Note that the classification of the three regions does not depend on as long as . This implies that, even in the Hill sphere, the gas in the region stays there temporarily and is not necessarily captured by the planet. We will discuss the flow in this region in §4.1.

3.2. Disk Structure and Gas Motion

Figure 4 shows azimuthally averaged density, velocity, and specific angular momentum of the gas with three different scales. We first recognize that the gas clearly forms a disk-like structure; the density distribution is concentrated near the midplane, almost Keplerian rotation is realized in the high density region (judging from the contour lines of specific angular momentum), and the flow velocity in the direction of the region is very low. The density profile in the -direction can be roughly described by hydrostatic equilibrium in the region where and , where is scale height of the circumplanetary disk defined by


From this dependence a flare-up disk is expected, which is seen indeed in Fig. 4. Above the disk surface, large downward velocity, which is almost free fall velocity, is observed; this shows that gas is indeed accreted directly onto the disk surface. The accreting gas forms a shock surface, which stands at . By analyzing streamlines, we confirmed that these gas elements are actually originated from off-midplane gas (mostly from ) in the protoplanetary disk. Also, the contour lines above the shock surface are aligned with the velocity vectors, which means that specific angular momentum does not change very much before the gas elements hit the disk surface in this inner region, where three-body effect is weak enough to be neglected.

Figure 5 shows radial mass flux (where ), as a function of the azimuth and elevation angles at spheres with four different radii , 0.3, 0.1, 0.03. This figure clearly shows that the accretion manner strongly depends on both these angles and is not spherically symmetric at all. For example, the flux at the sphere shows that the mass flux can be both inward and outward directions near the midplane (where , where is elevation angle) depending on azimuth angle , while it is always inward (i.e., negative flux) at high . The two maxima near and 180 on the midplane () in the case of corresponds to the outflow from the Hill sphere through the two Lagrangian points L and L shown in Fig. 2. The flux on the other three spheres also shows that there are two positive maxima and negative minima near the midplane, which implies the formation of a two-arm spiral structure in the circumplanetary disk. The range of the elevation angle that corresponds to outward flux shrinks with decreasing , which indicates the change of disk aspect ratio with .

Although Figure 5 shows that the gas at high elevation angle is always falling radially inward, it is difficult to judge the direction of net mass flux in the disk (corresponding to low ) from these plots. In order to examine the net mass flux, next we calculate azimuthally integrated mass flux at spheres with several radii as a function of (Figure 6(a)) :


We find that the direction of the net radial flux near the midplane ( in the case of ), where the gas moves both inward and outward through the sphere in the region depending on the longitude, is actually outward, while it moves inward for higher elevation angle ( in the case of ), which is obvious in Fig. 5. This basic profile does not change with , but where , the profile becomes sharper, that is, the region where gas moves outward is limited in narrower elevation angle, which corresponds to the dependence of disk aspect ratio on radius. This clear outward motion around the midplane for wide range of shows another evidence of the unbound state of the gas where , which is already shown in the above (Fig. 2). In the cases with and 0.03, there is a narrow band where the mass flux is inward (negative) with peak at and , respectively. Figure 6(b) shows azimuthally averaged radial velocity:


This quantitatively shows that the velocity of the gas elements accreting toward the disk nearly vertically at high elevation angle ( in the case with ) is nearly free-fall velocity (), which is faster than the sound speed (). Near the midplane, radial velocity is very small, which corresponds to the hydrostatic disk region described above (Fig. 4). Between the two regions, there is a transition region, which is under the shock standing at the disk surface, and the gas in the region is not in dynamically equilibrium state as a part of the disk. The range of of the transition region corresponds to the narrow band with inward mass flux described above. Thus, this is a kind of layered accretion and its mechanism will be discussed in § 3.3.

Figure 7 shows radial velocity in the midplane normalized by local Keplerian velocity defined by


We find that radial velocity is very small but takes on positive values for a wide range of radii (). Significant positive values at correspond to the outward motion which is already shown in Fig. 2. Non-negligible positive values of at may reflect non-steady state small-scale structure in the region, or artifact of averaging operation over discrete grids whose size may not be sufficiently small as compared to the distance from the origin. Formation of the circumplanetary disk with nearly Keplerian motion is also demonstrated in Fig. 8, which shows azimuthally averaged specific angular momentum at the midplane. We clearly see that the rotation velocity is nearly Keplerian for a wide range of radii (). Keplerian rotation is realized when the ratio of thermal energy of the gas to potential energy of the planet is much smaller than unity (or pressure force is much weaker than gravity force). Eq. (10) indicates that , which is square root of the ratio, becomes smaller with decreasing , which is consistent with our result.

Figure 9 shows the plots of azimuthally averaged density at the midplane and surface density. We can see that the profiles at can be fitted roughly by power-law functions as and , respectively. In an equilibrium state, we have . The above two functions imply that this relationship is approximately satisfied and the hydrostatic equilibrium is realized. Deviation from the power-law functions for the inner region may be due to insufficient computation time as compared to the time required to reach a steady state2, or due to the effect of sink cells around the center. Here we emphasize again that the distributions of mass and angular momentum of accreting gas that we present in § 3.3 are more important and useful quantities less affected by numerical procedures, as we mentioned in § 1. Note that Machida (2009) introduced an idea of centrifugal barrier to explain the peak of surface density. However, we do not find any physical reason for the existence of such a barrier because, as we will show in §3.3, distributions of accreting mass and angular momentum onto the circumplanetary disk are well described by power-law functions, which do not have typical lengths, such as centrifugal radius.

Figures 10 and 11 show variation of physical quantities along two streamlines starting from different heights that correspond to the passing and accreting regions. Overall variation is shown in Fig. 10, while close-up view near the shock front is shown in Fig. 11. The key quantity in this figure is the Bernoulli integral given by


where we included potential energy. This quantity is conserved along each streamline except shock surfaces.

First we examine a streamline in the midplane plane (Figs. 10(a) and 11(a)). The starting point of the streamline is , which is in the passing region (red region in Figs. 1 and 3); keeps positive values not far from 2, and monotonically decreases with increasing , where is the distance from the starting points of streamlines. With approaching to the planet (corresponding to ), the gas element climbs the tidal potential slope decreasing its kinetic energy, and density slightly decreases. Although these quantities change with , Bernoulli integral , which is the sum of the three quantities, is constant until , which shows there are no shocks along the path until this point. However, there is a shock surface at where kinetic energy decreases and density increases discontinuously, and decreases as a result of shock dissipation. This shock surface stands from the Hill sphere toward both inside and outside of the planet orbit, forming spiral density waves around the central star. After passage of the shock, remains nearly constant, although other quantities change with increasing . This means that the gas element on the streamline underwent only one shock with modest strength.

Next we examine an off-midplane streamline (Figs. 10(b) and 11(b)). The starting point is , which is just above the previous starting point with the same horizontal position, and this corresponds to the accreting region (blue region in Figs. 1 and 3). The quantities change in a similar way to the case in the midplane before , where the gas element hits the shock surface forming around the Hill sphere. We have non-zero values in this case but it is kept almost constant until it encounters the shock, which shows that the flow before the shock is laminar. When the gas element passes through the shock surface at , basic behaviour is the same as the case, but the gas element changes the direction slightly upward (positive direction) at the shock (see the line showing the change in in Fig. 11(b)). This is because the shock surface forms a bow shock and is curved upward, thus the gas element of off-midplane, which passed an oblique shock, gains upward momentum. After that, the gas element falls steeply toward the mid-plane and hits the surface of the circumplanetary disk, forming another shock at . The point of fall is at , which is sufficiently close to the planet to make use of its gravity for acceleration. Thus Mach number of the gas when it hits the disk surface is much larger than unity, and strong shock dissipation is caused. Owing to this strong dissipation, the gas becomes captured within the Hill sphere and merged as a part of disk. Note that the gradual decrease of before arises presumably from collisions between falling gas elements on the way to the surface, especially between two gas elements originated from and in the protoplanetary disk (see also Fig. 12 and the corresponding description in the text). After the passage of the shock forming above the circumplanetary disk (), becomes almost constant and azimuth angle changes with a constant rate, which means that the gas is now rotating around the planet, i.e., the gas element has become a part of the circumplanetary disk. Note that specific angular momentum is not constant along the streamlines at all, although it is sometimes assumed to be a conserved quantity.

Now we consider why the gas in the midplane cannot accrete to the circumplanetary disk. As described above, effective shock dissipation is required for gas elements to become accreted onto the circumplanetary disk, and the planet gravity is essential for the effective acceleration. When gas element consumes potential energy without significant increase of gas density, kinetic energy of the gas increases effectively (see Eq. (14)) and the gas element can have a high Mach number, which leads to strong shock dissipation. However, in order for the gas in the midplane to consume planet gravitational energy, the gas element has to pass through the high density region of the circumplanetary disk. The increase of density absorbs potential energy and prevents gas from effective acceleration. Gas near the midplane is thus difficult to form strong shocks and difficult to dissipate energy effectively, which prevents it from being captured by the planet gravity. Therefore, the gas in the midplane is difficult to be accreted into the circumplanetary disk.

One may think that passing many weak shocks during rotating around the planet in the circumplanetary disk would contribute to inward migration of the rotating gas. However, the decrease in is a third-order small quantity with respect to (where is Mach number), so such weak shocks are not effective for energy dissipation and inward migration. Note also that acceleration by the tidal force in the -direction alone is not sufficiently strong to form strong shocks. This is because the gas is in hydrostatic equilibrium supported by thermal pressure, so that the resultant accelerated velocity should not be much larger than the sound speed, being able to produce only weak shocks. We therefore conclude that the strong shock produced by direct infall onto the surface of the inner part of the circumplanetary disk and associated energy dissipation is essential for the accretion of the gas from protoplanetary disk onto the circumplanetary disk.

3.3. Distribution of Accreting Mass and Angular Momentum onto the Disk

Next, we examine accretion rate of mass and angular momentum as a function of distance from the planet. First we define the surface at which physical quantities immediately prior to the accretion onto the surface of a circumplanetary disk is measured. The surface is axisymmetric about the -axis and the height of the surface from the midplane is denoted by . We define the height of the surface by a linear function of the local scale height of the circumplanetary disk as


where and are constants. We set and so that the surface is above but not far from the shock surface of the circumplanetary disk in the region with .

Next we define mass flux onto the circumplanetary disk by


where and are the density and velocity vector at the surface defined by Eq. (15), is an upward unit vector normal to the surface, , and the negative sign is added so that becomes positive when the gas passes the surface downward. Note that appears in the denominator in Eq. (16) because is defined as a flux at the surface per unit area of the plane, not a flux per unit area of the surface.

Figure 12 shows several quantities on the surface projected on the plane. We find that the density of the accreting gas at the surface has peaks at , while the mass flux has peaks at , which are located closer to the planet than the density maxima. This is because accreting velocity, which is almost free-fall velocity, is faster at the inner region. These two peaks of the two-arm structure correspond to the places where gas elements coming from interior () and exterior () to the planet orbit collide with each other. As for the -component of specific angular momentum around the planet measured in the inertial frame (where the term arises from the rotation of the coordinate system), its value is positive (i.e., prograde) in the whole region. The distribution is roughly axisymmetric with some elongation in the -direction within , and the value increases with increasing radius. As for the angular momentum flux, the above axisymmetric property breaks and has peaks at , which simply corresponds to the two-arm structure observed in density or mass flux. Note that the negative values observed in in the outer region () are due to the fact that the gas passes through the surface from beneath. This is because, in this region , the ratio of the scale height to the radial distance of the circum-planetary disk is not much smaller than unity anymore and is comparable to accordingly. This leads to negative for a gas element whose horizontal velocity is comparable or larger than vertical velocity even when the vertical velocity is negative (downward). The surface of the disk in this region is hard to define, thus the way to measure mass flux described above may not be appropriate in this region.

For a better understanding of the distribution of accreting gas, we next show azimuthally-averaged quantities at the surface defined by Eq. (15). First we examine mass accretion rate. Figure 13 shows the plots of azimuthally-averaged mass flux onto the disk, defined by


Figure 14 shows the plots of the cumulative mass accretion rate, i.e., the mass accreted per unit time through the part of the surface where horizontal distance from the planet is smaller than , i.e.,


We can see from Fig. 13 that is almost constant from all the way to the very center (). Cumulative mass accretion rate is thus roughly proportional to , which is a power-law function and thus does not have any typical length except for the outer end of the downward accretion on the disk surface (at in this case). This distribution shows that the main contribution of mass accretion comes from the outer region () and there is no specific radius where gas accretion is concentrated somewhere in between. Note that dotted lines in Figures 13 and 14 show quantities calculated by Eqs. (17) and (18), respectively, integration being performed over only regions with . Deviation from each solid line can be only seen at , where gas can pass through the surface even upward depending on azimuth angle (see Fig. 12). This is the reason why (solid line) decreases with increasing radius at although it is a cumulative quantity. The treatment of our analysis in this region is a difficult problem and will be discussed later (§4.1).

The distribution of (or ) provides us with useful information but it is not sufficient for an understanding of the mass distribution in a circumplanetary disk. In general, the radial distance from the planet of the location where gas is first accreted on the disk is different from its final orbit after it settles in as a part of the Keplerian rotating disk. What connects the two is angular momentum of the gas. Figure 15 shows azimuthally-averaged specific angular momentum at the surface given by


We can see that is always smaller than the specific angular momentum for Keplerian rotation . This means that the gas does not have enough angular momentum for Keplerian rotation at the radius, which would lead inward migration of the gas after the fall on the disk surface, while the gas near the midplane moves outward, as mentioned above (Fig. 7). This explains the inward stream of the layer just under the shock of the disk surface (Fig. 6). In addition, is nearly proportional to , which is steeper than , and thus the ratio of to decreases with decreasing . This suggests that the gas accreted in the outer region of the disk () tends to keep the position when it rotates as a part of the circumplanetary disk, whereas the gas accreted in the inner region tends to move further inward in order to achieve balance between centrifugal and gravitational forces.

We now have the mass flux and angular momentum in the accretion flow, which allows us to estimate the effective distribution of accreting gas elements, assuming their redistribution to radial distances where their specific angular momentum matches that of the local Keplerian rotation. Let denote the radius where the gas with specific angular momentum is rotating with the Keplerian velocity, i.e.,


Then the effective distribution of azimuthally-averaged mass flux after redistribution of the gas with angular momentum conservation can be written as




The corresponding cumulative mass accretion rate is


The plots of and are shown in Figures 13 and 14. We can see that is roughly proportional to and is to . By comparison with , we find that is more center-concentrated distribution, but still outer region is a dominant source of gas accretion in the sense that . Since angular momentum of a gas element after the shock of the disk surface is not necessarily conserved until the gas reaches the radius where it rotates in the Keplerian velocity, the distributions of and should be regarded as the case of the most center-concentrated limit.

4. Discussion

4.1. Size of Circumplanetary Disks

The size of a circumplanetary disk is closely related to the location of satellite formation. However, it is not easy to define the outer edge of the disk, and there have been several attempts. One natural way would be to define the disk edge based on its density distribution. But the density is monotonically and smoothly decreasing with increasing radius toward the Hill sphere, so it is difficult to define the edge from density distribution. Ayliffe & Bate (2009b), who performed three-dimensional hydrodynamic simulations, suggested a criterion for the disk edge based on angular momentum distribution. Their simulations show that specific angular momentum of the disk gas has a peak at a certain radial location, whereas specific angular momentum of a Keplerian disk increases monotonically. Their Figure 2 suggests that the turnover point is about 1/3 of the Hill radius, and they defined the disk edge by the radial location of this peak. Martin & Lubow (2011) examined periodic orbits of a particle around a planet under the influence of the gravitational force from a central star and the planet, and found that, as the size of the orbit is increased, the orbits start crossing with each other at . They inferred that this corresponds to the location of the disk’s outer edge where tidal torque of central star’s gravity becomes strong, and called it tidal truncation radius (). The above value of is in agreement with the point of turnover of specific angular momentum found in their two-dimensional SPH simulation, as well as in the result of Ayliffe & Bate (2009b). Our Figure 8 shows that the turnover point is , which is also roughly in agreement with .3

We here consider an alternative criterion by examining defined by Eq. (13) to define the position of the outer edge. We observed that azimuthally-averaged radial velocity in the midplane is positive (outward) in almost all the region, although the value is very small for (Fig. 7), and the outward velocity significantly increases at (see also Fig. 2). Since gas at moves outward and escapes from the Hill sphere quickly, the region can be regarded as outside of the circumplanetary disk. We thus define the disk edge as the radial location where starts increasing significantly and has non-negligible positive value; in our case shown in Fig. 7. Note that a two-dimensional SPH simulation shows that negative torque is exerted on the outer disk (Martin & Lubow, 2011), which gives rise to inward flow and is contrary to our results. This might be due to the fact that two-dimensional simulations, which create clearer spiral structure than three-dimensional ones, tend to enhance torque density. In particular, the disk thickness in the outer region (), at which outflow is observed, is very thick (thickness is comparable to radius). Thus, it would be unlikely that our three-dimensional calculation is significantly affected by the negative torque suggested by Martin & Lubow (2011). The distribution of accreting angular momentum also seems to indicate a similar radius for the edge. As we see in Fig. 15, and have different dependence on . If we fit each profile as a single power-law function, they cross each other at in the present case. This roughly agrees with the location of the disk edge defined above based on the significant increase of outward velocity. The gas exterior to this radius has too large angular momentum to achieve Keplerian rotation, thus moves radially outward.

Once angular momentum of accreting gas is obtained, one can calculate mean specific angular momentum of the accretion flow, which is sometimes used in calculating the the so-called centrifugal radius to infer the disk size (e.g., Mosqueira & Estrada, 2003; Ward & Canup, 2010). Mean specific angular momentum of the accretion flow within a radius is given by




We can see from Figure 16 that increases nearly linearly with radius and levels around , where . This corresponds to , which might give us a rough estimation of the disk size. However, it is practically difficult to determine the upper bound of the integral range with respect to , which affects the value of , because the circumplanetary disks is smoothly connected to the proto-planetary disks. Also, unlike the case of a particle, angular momentum of a gas element is not a conserved quantity along a flow in principle, because gas is continuum medium which can transfer angular momentum through waves. This also makes it difficult to determine precisely, and thus as well.

4.2. Picture of Gas Accretion Flow onto Circumplanetary Disks

Figure 17 shows a schematic picture of gas accretion flow onto a circumplanetary disk based on the results obtained in the present work. Gas accretion occurs mostly downward from high altitude with high incident angle (see Fig. 4). The accreting gas is accelerated by the planet gravity to have almost free-fall velocity of the planet (Fig. 6b). The value of angular momentum of the accreting gas normalized by the local Keplerian angular momentum is lower in the inner region of the disk and higher in the outer region (Fig. 15). The falling gas reaches the shock surface formed on the top of the circumplanetary disk.

Gas near the midplane (especially where ) is almost in Keplerian rotation and hydrostatic equilibrium in the -direction. In this region, radial velocity is very small and it is difficult for the gas to accrete inward through the disk midplane, as mentioned above. On the other hand, the gas at in the midplane shows significant outflow, which eventually escapes from the Hill sphere (see Fig. 2). Thus inflow from high altitude and outflow near the midplane co-exist, and they do not interfere with each other. This means that there is a circulation across the Hill sphere and fresh (protoplanetary) gas is always supplied in the region. In addition, the two Lagrangian points L and L, which are often thought to be the most likely points of inflow, are actually the points of outflow, even in the gas accretion stage (Fig. 5). This seems consistent with Klahr & Kley (2006), who also performed three-dimensional hydrodynamic simulation with no explicit physical viscosity, but with lower resolution and with radiation. They showed that accretion mainly occurs via the poles of the planet and no inflow along the equatorial plane, which is quite similar to ours. But they explained that the circulation in their simulation is driven by accretion heating, which we do not consider. Thus it is not clear if the driving mechanisms of the circulation are the same. Note that Ayliffe & Bate (2009b) shows that the vast majority of the mass flows into the Hill sphere near the equator, which seems to be inconsistent with our result. There are several differences in setting between theirs and ours, so it is not easy to judge which factor causes the difference. Disk thickness would be one of the reasons for that, even though our disk thickness is not very different from theirs. Another possibility is gap formation. When a deep gap is formed, the contribution of the vertically accreting gas would become less significant, which reduces the difference. However, here we think that viscosity is the main reason for the difference. When there is explicit viscosity, the circumplanetary disk gas should transfer its angular momentum outward and most of gas would move inward accordingly, but this is not necessarily true in the inviscid limit. Actually, inviscid simulations by Klahr & Kley (2006) showed similar results with ours. Also, inviscid limit might not be bad because circumplanetary disks are likely to be MRI inactive in most cases (Fujii et al., 2011).

Since dust particles tend to settle down toward the midplane, gas accretion flow from high altitude is likely dust-poor gas, which diminishes dust-to-gas ratio in the circumplanetary disk. On the other hand, since gas in the outflow region () rotates significantly slower than Keplerian velocity around the planet, dust particles would migrate inward quickly. Thus the outflow in the midplane is also likely dust-poor gas, which would be the source of solid material in the circumplanetary disk and would enhance dust-to-gas ratio in the disk. Dust-to-gas ratio in the circumplanetary disk is one of the most important factors for satellite formation processes (Canup & Ward, 2006), thus these two filtering effects would become important for satellite formation. In addition, studies on size evolution of solid material in protoplanetary disks, such as Kobayashi & Tanaka (2010), would also be important for the filtering effects. Further studies are needed on this issue.

As we see in Figure 6, there is the layer just under the shock at the disk surface where gas moves inward. This layered accretion is explained by the fact that angular momentum of the accreting gas onto the disk surface is lower than that of the Keplerian rotation at the radius (Fig. 15). The inward stream in the layer might play an important role for the net mass accretion toward the planet in circumplanetary disks. This accretion picture should be explored more quantitatively by further high-resolution simulations, using a code with lower artificial viscosity such as the one with polar coordinates for numerical grids.

4.3. Effects of Gap Formation

One important issue to be addressed is the effect of gap, which is a lower density annulus region near the planet orbit in the protoplanetary disk. We observed the gap as a slightly lower density band formed around in our simulations. However, in the last stage of giant planet formation, a giant planet would become massive enough to create a deep gap, which would be able to truncate its growth. Since the deep gap is associated with steep density gradient at the edge, gap formation may affect accretion flow and circumplanetary disk formation. However, at the edge of such a deep gap, the gas density changes over a radial distance comparable to the disk scale height, while the width of the accretion band onto the circumplanetary disk is much narrower than the scale height (Fig. 3). Therefore, we think that the gap formation would not affect the qualitative feature of the accretion flow. Another possible effect of gap formation is the disk size. As described in §4.1, several mechanisms are proposed to explain the disk size. If the disk size is determined by the tidal effect as suggested by Martin & Lubow (2011), it should not be affected by gap formation, except when the Bondi radius is smaller than the Hill radius. On the other hand, if the disk size is determined by the radially outward velocity of the gas as described in §4.1, lower density around the Hill sphere and thus stronger radial pressure gradient near the disk edge may enhance the outflow from the circumplanetary disk, and the disk size may become smaller. In any case, effects of gap formation should be examined in future works to check the validity of the results shown in this paper.

5. Summary

In order to understand the structure of circumplanetary disks, we performed high-resolution hydrodynamic simulation and analyzed gas accretion flow in detail. We confirmed that gas accretion onto circumplanetary disks occurs in a manner that the gas is accreted from high altitude toward the disk surface downward with large incident angle, which was suggested by previous studies (e.g., D’Angelo et al., 2003; Bate et al., 2003), whereas Ayliffe & Bate (2009b) showed that, in terms of mass flux, the downward accretion flow is not significant, which is inconsistent with our results. We found that the gas that has passed through the shock surface moves inward because its specific angular momentum is smaller than that of Keplerian rotation, whereas gas accretion through the midplane does not occur. While the net gas accretion across the Hill sphere is inward, outflow through near the two Lagrangian points is observed, although the regions around the two Lagrangian points have been, from the point of view of the potential energy, thought as a main accretion channel from protoplanetary disks. This outflow was not observed in previous hydrodynamic simulations for Jupiter-sized planets (e.g., Bate et al., 2003; Ayliffe & Bate, 2009b). Outward radial velocity of the gas near the midplane significantly increases at a radial distance of about 0.2 times the Hill radius from the planet, and the gas can escape from the Hill sphere within a short period of time.

We also obtained the distribution of mass and angular momentum of accreting gas onto the surface of circumplanetary disks. We found that the accretion rates of mass and angular momentum can be well described by power-law functions. This distribution would be useful in the study of satellite formation, for example, a radially one-dimensional viscous-evolution model for circumplanetary disks, such as Ward & Canup (2010) and Martin & Lubow (2011). Recent development of viscous modeling for MRI turbulence in protoplanetary and circumplanetary disks (Okuzumi & Hirose, 2011; Fujii et al., 2011) would also contribute to construct more realistic models for the circumplanetary disks. However, in order to understand satellite formation processes, long term evolution of circumplanetary disks, which is determined by the evolution of accretion rate from protoplanetary disks to circumplanetary disks, is necessary and thus global evolution of protoplanetary disks with embedded giant planets until the complete dissipation of protoplanetary disks needs to be examined. Together with such global models (e.g., Tanigawa & Ikoma, 2007; Ward & Canup, 2010), we will be able to obtain long-term evolution of circumplanetary disks, which would provide better understanding of satellite formation processes.

Our results demonstrate that gas accretion toward and within circumplanetary disks has a complicated vertical structure. The width of the accretion band toward a circumplanetary disk depends on the initial height of gas elements (Fig. 3). Also, after accretion onto the disk, the gas just below the shock surface migrates inward, while the gas near the midplane moves radially outward (Fig. 17). Such vertical heterogeneity of the flow may have significant influence on the dynamical evolution of solid bodies in the circumplanetary disk, which we will examine in our future work. We also need to address the effect of viscosity and local calculation, which would affect the flow quantitatively.

We are grateful to Hidekazu Tanaka, Hiroshi Kobayashi, and Satoshi Okuzumi for valuable comments. T. T. also thanks Sei-ichiro Watanabe and Ryuji Morishima for continuous encouragement to advance this. This work was supported by Center for Planetary Science running under the auspices of the MEXT Global COE Program entitled “Foundation of International Center for Planetary Science”. We are also grateful for the support by JSPS and NASA’s Origins of Solar Systems Program. Numerical calculations were carried out on NEC SX-9 at Center for Computational Astrophysics, CfCA, of National Astronomical Observatory of Japan. A part of the figures were produced by GFD-DENNOU Library.
Figure 1.— Streamlines starting from four different heights (), with and . Interval of the starting points is 0.05 in the -direction. Green, blue, and red curves show streamlines in the horseshoe, accretion, and passing region, respectively.
Figure 2.— Streamlines in the midplane of a circumplanetary disk. Left panel shows a streamline starting from in the midplane . Dashed line shows contour line of , which passes the two Lagrangian points L and L . Colors show potential ; red regions are and blue regions are . Right panel shows a closer view of three streamlines starting from (same as the left panel), 0.1, and 0.05, respectively.
Figure 3.— Classification of starting points of streamlines on the plane at . Streamlines starting from the red region reach the boundary with (passing region), those in green region make U-turn and reach the boundary with (horseshoe region), and those in blue region become trapped in (accretion region).
Figure 4.— Circumplanetary disk structure in the plane with three different spatial scales. Log density is shown with colors, specific angular momentum is shown with contour lines, and gas velocity is expressed with arrows. All quantities are averaged in the azimuthal direction.
Figure 5.— Mass flux on spheres with radii 1.0, 0.3, 0.1, and 0.03. Horizontal axis is azimuth angle and vertical axis is elevation angle on the sphere. Outward flux is defined to be positive. and correspond to the sub-solar and anti-solar points, respectively.

Figure 6.— Left panel shows azimuthally integrated mass flux at the radii 1.0 (solid red line), 0.3 (dashed green line), 0.1 (dot-dashed blue line), 0.03 (dotted purple line) as a function of elevation angle . Right panel shows azimuthally averaged radial velocity at the same radii.
Figure 7.— Azimuthally averaged radial velocity in the midplane normalized by local Keplerian velocity , which corresponds to angle between the velocity vector of the flow and the circular orbit at the point.
Figure 8.— Azimuthally averaged specific angular momentum at the midplane of the circumplanetary disk. Solid line shows specific angular momentum measured in the rotating frame of the Hill coordinate, dashed line is the case when it is measured in the inertial frame, and dotted line shows specific angular momentum of Keplerian rotation around the planet.

Figure 9.— Azimuthally-averaged density at the midplane (left). Azimuthally averaged surface density (right), both as a function of the radial distance in the midplane.
Figure 10.— Quantities along two streamlines starting from two different heights. Horizontal axis is distance along each streamline from the starting point. Starting points for upper and lower panel are and , respectively. Thick solid (red) line shows Bernoulli integral, thick dashed (green) line shows kinetic energy, thick dotted (blue) line shows potential energy, and thick dot-dashed line (purple) shows logarithm of density. Thin dotted, dashed, dot-dashed lines show , , , respectively. The streamline shown in the lower panel is accreted onto the circumplanetary disk, while that in the upper panel is not. Note that the streamline in the case of is truncated at a point where the gas element starts rotating around the planet in the circumplanetary disk.
Figure 11.— Close up view of Fig. 10. In addition to the quantities shown in Fig. 10, azimuth angle divided by (thin dashed blue), specific angular momentum (thin dot-dot-dashed green), and logarithm of (thin solid red) are also shown.
Figure 12.— Distribution of the quantities associated with the accretion flow at the surface defined by . Density , mass flux through the disk surface , specific angular momentum , and angular momentum flux through the surface .
Figure 13.— Azimuthally averaged mass fluxes through the surface as a function of radius. Solid line shows and dashed line shows . Dotted line shows but excluding the region where from the integration Eq. (17). is a flux directly observed at the surface, while is an effective flux assuming redistribution of radial location where specific angular momentum matches that of the local Keplerian rotation. (see Eqs. (17) and (21)).
Figure 14.— Mass accretion rate onto the disk surface within a circle of radius. Solid line shows (Eq. (18)) and dashed line shows (Eq. (23)). Dotted line shows but excluding the region where from the integration of Eq. (18).
Figure 15.— Azimuthally-averaged specific angular momentum of the accretion flow onto circumplanetary disks denoted by in the text (solid line). Dashed line is , which is specific angular momentum for Keplerian rotation around the planet.
Figure 16.— Mean specific angular momentum of the gas accreting onto the part of the circumplanetary disk within the radial distance .
Figure 17.— Schematic picture of flow structure of circumplanetary disks.


  1. slugcomment: Accepted for publication in ApJ
  2. Time required to reach steady state can be estimated as , which is longer than the duration with full level calculation. Here we used (see Fig. 14).
  3. In this comparison, we assume that specific angular momentum in these works were not measured in inertial frame, but in the rotating frame.


  1. Ayliffe, B. A., & Bate, M. R. 2009, MNRAS, 397, 657
  2. Bate, M. R., Lubow, S. H., Ogilvie, G. I., & Miller, K. A. 2003, MNRAS, 341, 213
  3. Bodenheimer, P., & Pollack, J. B. 1986, Icarus, 67, 391
  4. Canup, R. M., & Ward, W. R. 2002, AJ, 124, 3404
  5. Canup, R. M., & Ward, W. R. 2006, Nature, 441, 834
  6. Coradini, A., Magni, G., & Turrini, D. 2010, Space Sci. Rev., 153, 411
  7. D’Angelo, G., Henning, T., & Kley, W. 2002, A&A, 385, 647
  8. D’Angelo, G., Kley, W., & Henning, T. 2003, ApJ, 586, 540
  9. Fujii, Y. I., Okuzumi, S., & Inutsuka, S. 2011, ApJ, 743, 53
  10. Hayashi, C. 1981, Prog. Theor. Phys. Suppl., 70, 35
  11. Ikoma, M., Nakazawa, K., & Emori, H. 2000, ApJ, 537, 1013
  12. Klahr, H. H., & Kley, W. 2006, A&A, 445, 747
  13. Kley, W. 1999, MNRAS, 303, 696
  14. Kobayashi, H., & Tanaka, H. 2010, Icarus, 206, 735
  15. Korycansky, D. G., & Papaloizou, J. C. B. 1996, ApJS, 105, 181
  16. Lubow, S. H., Seibert, M., & Artymowicz, P. 1999, ApJ, 526, 1001
  17. Lunine, J. I., & Stevenson, D. J. 1982, Icarus, 52, 14
  18. Machida, M. N., Matsumoto, T., Tomisaka, K., & Hanawa, T. 2005, MNRAS, 362, 369
  19. Machida, M. N., Matsumoto, T., Hanawa, T. & Tomisaka, K. 2006, ApJ, 645, 1227
  20. Machida, M. N., Kokubo, E., Inutsuka, S., & Matsumoto, T. 2008, ApJ, 685, 1220
  21. Machida, M. N. 2009, MNRAS, 392, 514
  22. Martin, R. G., & Lubow, S. H. 2011, MNRAS, 413, 1447
  23. Matsumoto, T., & Hanawa, T. 2003, ApJ, 595, 913
  24. Miki, S. 1982, Prog. Theor. Phys., 67, 1053
  25. Mizuno, H. 1980, Prog. Theor. Phys., 64, 544
  26. Mosqueira, I., & Estrada, P. R. 2003, Icarus, 163, 198
  27. Okuzumi, S., & Hirose, S. 2011, ApJ, 742, 65
  28. Paardekooper, S.-J., & Mellema, G. 2008, A&A, 478, 245
  29. Sasaki, T., Steward, G. R., & Ida, S. 2010, ApJ, 714, 1052
  30. Sekiya, M., Miyama, S., & Hayashi. C. 1987, EM&P, 39, 1
  31. Tanigawa, T., & Watanabe, S. 2002, ApJ, 580, 506
  32. Tanigawa, T., & Ikoma, M. 2007, ApJ, 667, 557
  33. Ward, W. R., & Canup, R. M. 2010, AJ, 140, 1168
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description