ACCRETION OF JUPITER-MASS PLANETS IN THE LIMIT OF VANISHING VISCOSITY

How exactly giant planets form is still one of the most puzzling questions in today's planetary science. The most popular giant planet formation theory is the core accretion model by Bodenheimer & Pollack (1986). There are three main stages of formation in this model. First, a planetary core is formed and starts attracting the gas within its Bondi radius. When the core reaches 12–16 M_⊕, the gas envelope starts to contract quasi-statically while the accretion rates increase (Pollack et al. 1996). This second stage takes the longest time, of the order of a few million years. The final, runaway gas accretion phase starts when the envelope and core masses are approximately equal. This phase should not stop until the planet has opened a deep gap in the gas disk. However, this happens only when the planet reaches a mass of 5–10 Jupiter masses (Kley 1999; Lubow & D'Angelo 2006). Thus, Jovian-mass planets can double their mass on an order of 10⁴ yr.

Thus, one should expect to observe a dichotomy in the mass distribution of planets; planets should be either smaller than ∼30 Earth masses, i.e., those that did not reach the phase of runaway gas accretion, or larger than a few Jupiter masses, i.e., those that entered and completed the fast runaway growth phase. Planets in between these two mass categories should be extremely rare. This is the converse of what is observed (e.g., Mayor et al. 2011). Thus, there is a need to understand what sets the final mass of a giant planet.

An obvious possibility to stop accretion is that the gas disappears while the planet is still growing. However, the lifetime of gaseous protoplanetary disks is of the order of a few million years (Haisch et al. 2001), which is much longer than the runaway growth timescale (10⁴ yr). It is very unlikely that the disappearance of the disk can happen at the right time to stop the runaway growth of the planet. Another possibility is that a planet cannot accrete more gas than what is delivered to its orbit by viscous accretion, i.e., it cannot grow faster than the star accretion rate. In general, the accretion rate observed in protoplanetary disks is on the order of 10⁻⁸–10⁻⁷ M_☉ yr⁻¹. This would allow the accretion of Jupiter's atmosphere in 10⁴–10⁵ yr, which is too short relative to the lifetime of the disk. If one requires that Jupiter takes a million years to accrete its envelope, then its runaway growth needs to be limited by a stellar accretion rate of 10⁻⁹ M_☉ yr⁻¹. However, at this very low rate, the disk photoevaporates rapidly (i.e., a few 10⁵ yr, see, e.g., Koepferl et al. 2013; Gorti et al. 2009; Szulágyi et al. 2012). Thus a Jupiter mass of gas is unlikely to be accreted by the planet. A very accurate tuning between the viscous accretion rate, the photoevaporation rate, and the runaway growth seems to be needed to allow a planet to grow to Jupiter mass but not beyond this limit. Something must still be missing from the picture.

What we need is a mechanism that slows down runaway growth. So it occurs on a timescale comparable to the disk's lifetime. A possibility is that the circumplanetary disk (CPD) acts as a regulator of gas accretion rate onto the planet. Before the gas is accreted by the planet, it has to pass through the CPD because of angular momentum conservation. The actual accretion rate of the planet then depends on the timescale for angular momentum transport within the CPD. If the CPD has a very low viscosity, then the transport of angular momentum through this disk is inefficient and gas accretes onto the planet at a slow rate. In this situation, the observed mass spectrum of the giant planets is set by the competition between gas accretion and gas dissipation (Rivier et al. 2012).

There are good reasons to think that the viscosity is very low in the CPD (although see Gressel et al. 2013 for an opposite view). First, the planets are thought to be formed in a dead zone of the circumstellar disk, where the viscosity is low (Thommes et al. 2008; Martin & Lubow 2011). Second, because the CPD is shadowed by the circumstellar disk and by the remaining gas in the gap, its irradiation geometry may be unfavorable for ionization (Turner et al. 2010, 2013). Finally, the large orbital frequencies in the CPD make the magnetic Reynolds number too small to derive the magneto-rotational instability (MRI; Fujii et al. 2011; Turner et al. 2013; Fujii et al. 2013).

Motivated by these considerations, in this paper we study the dynamics of the CPD in detail.

Our simulations are similar to those in Machida et al. (2010) and Tanigawa et al. (2012), with one main difference. Instead of using a local shearing sheet approximation, we perform global disk simulations. This better allows us to study the connection between the circumstellar disk and the CPD, i.e., opening of a gap, gas flow through the gap, etc. Moreover, we investigate in more detail the accretion rate of a Jovian-mass planet in the limit of vanishing viscosity. To do this, it is not enough to perform simulations with no prescribed viscosity as in Tanigawa et al. (2012), because every numerical simulation is affected by numerical viscosity. We need to identify the various accretion mechanisms and distinguish between those dependent on viscosity and those independent of viscosity (i.e., polar inflow from the circumstellar disk, loss of angular momentum due to shocks, stellar torque exerted on the CPD, etc.) and evaluate their magnitude.

A well-known crucial issue for simulating gas accretion onto a planet is the choice the equation of state (EOS). Several works have stressed the need to use an adiabatic EOS—possibly complemented by a recipe for radiative cooling—in order to study planet accretion (D'Angelo et al. 2003; Klahr & Kley 2006; Paardekooper & Mellema 2008; Ayliffe & Bate 2009). However, the differences with the isothermal EOS are fundamental for small-mass planets (up to Saturn's mass) but less for Jupiter-mass planets. In the latter case, the flow of gas is mostly dominated by the planet's gravity. CPDs definitely form around Jupiter-mass planets, and the differences between isothermal and adiabatic simulations are mostly limited to the mass of the CPD and its scale height (D'Angelo et al. 2003; Ayliffe & Bate 2009; Gressel et al. 2013). Thus, we prefer to use the isothermal EOS, with temperature dependence on stellar distance (hereafter locally isothermal), for multiple reasons. The first is that we wish to focus our paper on the role of numerical viscosity and on viscosity-independent transport mechanisms within the CPD, which have never been thoroughly discussed before; these considerations should be independent of the EOS assumed. Second, radiative simulations imply additional, badly constrained parameters such as those in the prescription for the opacity laws (e.g., Ayliffe & Bate 2009; Bitsch et al. 2013). We want to focus the discussion on the objectives stated above without distraction. Third, we wish to make direct comparisons particularly concerning the differences between our global disk simulations and shearing-sheet studies (Machida et al. 2010; Tanigawa et al. 2012), and the latter have been done with isothermal EOS. Finally, this paper is the first in a series of future studies; therefore, we wish to begin with the most simple case and build on it incrementally. Nevertheless, for each result that we present, we will state to what extent we expect it to be valid or different in a non-isothermal context.

The paper is structured as follows. In Section 2, we describe the setup of our hydrodynamic simulations. This is followed by the results on the structure of the CPD in Section 3. Then, Section 4 discusses our findings on the accretion mechanism. Section 5 reports discussions and perspectives. Finally, Section 6 summarizes the conclusions of our work.

2.1. Physical Model

2.2. Numerical Model

For the simulations, we used a three-dimensional (3D) nested-grid code, called Jupiter. The Jupiter code solves the Riemann problem at every cell boundary (Toro 2009) to ensure the conservation of mass and the three components of momentum:

$$\begin{equation} \rho _t+\nabla \cdot (p{{\mathbb v}})=0 \end{equation} \tag{ 1 }$$

$$\begin{equation} (\rho v)_t+\nabla \cdot (\rho {{\mathbb v}} \otimes {{\mathbb v}} + p {\mathbb {I}})=0, \end{equation} \tag{ 2 }$$

where ρ is the density, p is the pressure, $\mathbb {v}$ is the velocity vector, and $\mathbb {I}$ indicates the identity matrix. The use of a Riemann-solver makes Jupiter particularly suited to treat shocks, contrary to the van Leer method (van Leer 1977). The Riemann-solvers implemented in the Jupiter code are approximated solvers based on the exact solution: a Two-Shock solver and a Two-Rarefaction solver (de Val-Borro et al. 2006).

The timestep in the simulation is adapted by the code in order to satisfy the Courant–Friedrichs–Lewy condition for all levels of mesh resolution:

$$\begin{equation} C=\Delta t \sum _{i=1}^{3} \frac{v_{x_{i}}}{\Delta x_i} \le 1.0, \end{equation} \tag{ 3 }$$

where C is the Courant number, i represents the number of dimensions, x_i means the spatial variables, and v indicates the velocity. The timestep at a given level can be the same as the timestep on the higher resolution level, or it can be twice that timestep, in which case two iterations are performed on the finer level while one iteration is done on the coarser level. This latter technique is called timestep subcycling. We use an adaptive subcycling procedure, which will be described in a forthcoming publication, in order to obtain the maximum speed up of the code (the highest possible ratio of physical time over wall clock time).

The full viscous stress tensor is implemented in the code in three geometries: Cartesian, cylindrical, and spherical. The spherical implementation, which we use in this work, has been tested thoroughly in a prior work (Fromang et al. 2011).

We employed a system of 8 nested grids, where at level 0 (i.e., in the coarsest grid) the resolution was 628 × 208 × 15 cells for the directions of azimuth, radius, and colatitude, respectively. Each additional grid was added after the gas reached a stationary configuration, and each of them was centered on the planet. The size of the cells in a grid at a given level was 1/2 in each spatial direction of the cell size of the next larger grid. Table 1 contains the number of cells on each level and the grid boundaries. In the finest level, the cell length was 7.82 × 10⁻⁵a, which is 0.113% of the Hill radius of the planet and 87% of the radius of the present day Jupiter, assuming the planet orbits 5 AU away from the star. The cells had the same length in every directions (i.e., they were cubes), and the radial spacing in between them was arithmetic.

Table 1. Number of Cells on Different Grid Levels

Level	No. of Cells in	No. of Cells in	No. of Cells in	Boundaries of the	Boundaries of the	Boundaries of the
	Azimuth	Radius	Colatitude	Levels in Azimuth (rad)	Levels in Radius (a)	Levels in Colatitude (rad)
0	628	208	15	[−π, π]	[0.41, 2.49]	[1.42, π/2]
1	112	112	24	[−0.27735, 0.27735]	[0.72264, 1.27735]	[1.451041, π/2]
2	112	112	40	[−0.138675, 0.138675]	[0.861325, 1.138675]	[1.47137, π/2]
3	112	110	56	[−0.0693375, 0.0693375]	[0.9306625, 1.0693375]	[1.5014588, π/2]
4	112	110	56	[−0.03466875, 0.03466875]	[0.96533125,1.03468875]	[1.5361276, π/2]
5	112	112	56	[−0.017334375, 0.017334375]	[0.98266562, 1.017334375]	[1.553462, π/2]
6	112	112	56	[−0.0086671875, 0.00866719]	[0.99133281, 1.0086671875]	[1.5621291, π/2]
7	112	110	56	[−0.0043335938, 0.00433359]	[0.99566641, 1.0043335938]	[1.566462737, π/2]

Download table as:  ASCIITypeset image

On level 0, we used reflecting boundary conditions except in the azimuthal direction where we used periodic conditions. The communication between the grids at level i and i + 1 (where i = 0, ..., 7) were done through ghost cells with multi-linear interpolation.

To test the effects of numerical viscosity, we also ran a simulation with a twice finer resolution (1256  ×  416  ×  30), which we call hereafter the "high-resolution simulation." Because the simulation is extremely slow at this resolution, we did not start it from time zero but rather from the output of the nominal resolution at 238 orbits; we rebinned the gas on the new grids and then ran the code for an additional 10 orbits.

The planet was not modeled but was treated as a point mass placed in the corner of four cells on the midplane. In order to avoid a singularity, the planetary potential was smoothed as

$$\begin{eqnarray} U_{p} &=& -\frac{G M_p}{\sqrt{x_d^2+y_d^2+z_d^2+{r_s}^2}}, \end{eqnarray} \tag{ 4 }$$

where x_d = x − x_p, y_d = y − y_p, and z_d = z − z_p are the distance-vector components from the planet in Cartesian coordinates. The smoothing length r_s was set equal to the cell size in levels 0–4. From level 5 on, we used 2 cell sizes. Moreover, when introducing levels 6 and 7, we progressively decreased r_s from the value used in the previous level to its desired final value. For example, when introducing level 7, first the same smoothing length was applied as on level 6 (which is equal to 4 cell sizes on level 7), but then it was decreased in time with a sinusoidal function until it reached r_s = 2 cells. This technique was done to allow the gas to adapt to a gradual change of the gravitational potential.

Because of the isothermal character of our fluid equations, a huge amount of mass tends to pile up in the few cells neighboring the point-mass planet. This causes numerical instability. Thus, we applied a density cut: if the volume density reached 1.42 × 10⁵ (in the code's units), then the volume density of that cell was limited to this value (hereafter we refer to this as the "mass cut"). We keep track of the mass removed in this operation, from which we compute the planet accretion rate. However, the mass of the planet that enters in the gravitational potential was not changed. This is because we are interested in the accretion rate of a Jupiter-mass planet and not in the growth of the planet itself.

In this section, we describe our results about the CPD structure: its vertical structure, its radial structure, and the radial flow in the midplane. All the results that we present are from our nominal simulation unless we specify otherwise.

Because the grid at level 0 covers the circumstellar disk globally, the planet can open a gap around its orbit (see Figure 1). This was not the case in the simulations of Machida et al. (2010) and Tanigawa et al. (2012) because of the shearing-sheet approximation they adopted. Figure 1 also shows that the density wave launched by the planet in the circumstellar disk smoothly joins the CPD and spirals into it down to the planet (see also Figure 10).

Zoom In Zoom Out Reset image size

3.1. The Vertical Structure of the Circumplanetary Disk

We start by discussing the vertical structure of the CPD. For this purpose, it is convenient to characterize the CPD based on the z-component (in Cartesian coordinates) of the specific angular momentum with respect to the planet, normalized to the Keplerian value:

$$\begin{equation} L_z=\frac{x_dv_y-y_dv_x+(x_d^2+y_d^2) \Omega }{\sqrt{G M_{P} \sqrt{x_d^2+y_d^2}}}, \end{equation} \tag{ 5 }$$

where v_x and v_y are the velocity-components transformed to Cartesian coordinates in the corotating frame.

Figure 2 represents a vertical slice at azimuth = 0.0 of the L_z distribution in the neighborhood of the planet, which is located at the center of the upper axis. We see that L_z rapidly drops from ∼1 to ∼0.5 at a location where the density shows a clear discontinuity (see Figure 3). Therefore, hereafter we define the CPD as the region where L_z is larger than 0.65 (see the corresponding contour line in Figure 2). Even this value is arbitrary; however, given the steep gradient of L_z near the surface of the disk, changing this threshold would not change significantly the results presented below.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 2 shows that the gas in the midplane and near the midplane is sub-Keplerian (similarly to Tanigawa et al. 2012; Uribe et al. 2013). However, note that near the upper layer of the disk the flow is slightly super-Keplerian, i.e., in the region bounded by the contour line L_z = 1. This is due to the fact that the disk is very flared, so that the radial pressure gradient near its surface is positive. In the viscous simulation, however, this super-Keplerian near-surface layer does not exist. This is due to the higher viscosity that limits the vertical shear in the CPD.

The gas located below the CPD is falling toward the CPD with a large vertical velocity, as indicated by the arrows on Figure 3 (see also Ayliffe & Bate 2009 and Tanigawa et al. 2012). As pointed out in Tanigawa et al. (2012), the sharp vertical boundary of the CPD clearly visible in the L_z and density maps is due to a shock front generated by the vertical influx. As in Tanigawa et al. (2012), we also note from Figure 2 that the vertical inflow hits the CPD with a value of L_z that is much smaller than that characterizing the CPD at the same location. Thus, the vertical inflow slows down the rotation of the CPD, promoting radial infall at the surface of the CPD.

We find that the vertical influx has also a strong influence on the aspect ratio of the CPD. First of all, as a reminder, the pressure scale height of the CPD at hydrostatic equilibrium is H_CPD ≡ c_s/ω, where $\omega =\sqrt{GM_p/d^3}$ is the angular velocity around the planet, and $d=\sqrt{\vphantom{A^A}\smash{{{x_d^2+y_d^2}}}}$ indicates the distance from the planet. The sound speed (c_s = 0.05 r^−1/2) is almost constant in the CPD in our locally isothermal simulation. So, we expect the aspect ratio of the CPD to be $H_{{\rm CPD}}/d= ({c_s}/{\sqrt{GM_p}})\,d^{1/2} \,\approx \, 1.6(d/a)^{1/2}=\, 0.16\,(d/0.01a)^{1/2}$, which is very thick and flared. As we show in Figure 4, the surface of the CPD defined by z_CPD (i.e., the uppermost z-coordinate where L_z ≈ 0.65) is indeed strongly flared, but its aspect ratio has also a strong dependence on the azimuth relative to the planet. In fact, z_CPD/d is changing from ∼20% to >100%. To our knowledge, this wavy surface has not been described yet in the literature. The wavy surface pattern is due to the dynamical pressure of the vertical mass inflow, which is not uniform in planetocentric azimuth (see Figure 5); it is maximal along an axis close to the axis of the spiral arm. Now we take a Lagrangian approach and consider fluid elements in the CPD orbiting on circles centered on the planet; because the pressure due to the vertical inflow has two maxima, they will feel a maximum of pressure twice per orbit. However, it takes a time t_delay ≈ H_CPD/c_s for the gas in the CPD to react to the pressure pulse. In this time, the fluid elements rotate by an angle θ_delay = ω × t_delay, which is the angle between the axis marking the minimum height of the CPD and that marking the maximum pressure. Because ω × t_delay = ω/H_CPD/c_s = 1, this angle is independent of the distance from the planet d. The comparison of Figures 4 and 5 clearly shows an angle of the order of unity (in radians) between the maximum pressure and the minimum CPD height. A toy model is presented in the Appendix about how the pressure of the vertical influx leads to the observed structure of the CPD.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 6 shows the vertical density distribution in the CPD from z = 0 to z = z_CPD at a given radius for various values of the azimuth. The mass is conserved along an orbital period, so the integral of each density curve is the same. On top of the expected equilibrium Gaussian shape, one can note oscillations with two knot points where the density does not change with azimuth. This is reminiscent of stationary waves.

Zoom In Zoom Out Reset image size

Simulations implementing an adiabatic EOS (e.g., Ayliffe & Bate 2009, 2012; Gressel et al. 2013) also find that the vertical inflow is the main feeding mechanism for the CPD. The CPD, however, has a larger scale height as it is hotter and the boundary between the disk and the vertical flow is less sharp than in our isothermal simulations. We will come back to this last, important issue in Section 5.

3.2. Radial Structure of the Circumplanetary Disk

We now move to discuss the orbital motion inside the CPD and the location of its outer radial boundary.

The specific normalized angular momentum L_z on the midplane declines with the distance from the planet (Ayliffe & Bate 2009; Tanigawa et al. 2012), but it does not show a steep gradient like the one in the vertical direction near the surface of the disk. Therefore, previous authors assumed arbitrary limits in L_z, obtaining different radial extensions for the CPD.

We think that it is more meaningful to look at the orbital motion of the fluid elements, defining the CPD as the region where the orbits are basically circular. Quasi-circular orbits may have a small value of L_z if the CPD is strongly sub-Keplerian due to a steep radial pressure gradient, but they are clearly part of a disk.

In order to visualize easily where the orbits of the disk are quasi-circular, we proceeded as follows. First, we calculated the semi-major axis a and the eccentricity e in every cell from the cell's coordinates, the recorded velocities, and the planetary potential; then we plotted the apocenter of the orbit Q = a(1 + e) versus the planetocentric radius d of the cell. If at a given radius every cell, whatever its planetocentric azimuth, appears to be at the apocenter (Q = d), it obviously means the streamline in the disk is circular, although sub-Keplerian. We can see on Figure 7 that this is the case up to ∼0.48 Hill radii. If we use this definition for the radial extent of the CPD, then the disk is a bit wider than the previously recorded radial extensions of ∼0.1–0.3 Hill radii (Tanigawa et al. 2012; Ayliffe & Bate 2009).

Zoom In Zoom Out Reset image size

We remark, however, that the eccentricity of the streamlines in the disk depends on the viscosity. In fact, as we will see in Section 3.3, the streamlines are eccentric if they are shocked at the passage through the wave generated by the stellar tide. The smaller is the effective viscosity—prescribed or numeric—in the CPD, the closer to the planet the wave propagates and shocks. Thus, defining the CPD as the region where streamlines are circular may lead to the uncomfortable situation that the disk may become vanishingly small in the ideal limit of zero viscosity. In fact, in our viscous simulation, circular orbits extend up to ∼0.55 Hill radii, and in the high resolution simulation, which halves the numerical viscosity, they extend only up to ∼0.28 Hill radii. This is an important point that should be kept in mind when analyzing the results of simulations, regardless if conducted with an isothermal or adiabatic EOS.

Alternatively, we may define the radial extent of the CPD as the largest circle from which streamlines wrap around the planet at least once before becoming unbound, in either the forward or backward integration. If we adopt this definition, the radius of our CPD is approximately 1/2 to 3/4 of the Hill radius.

On Figure 8, our CPD's column density (∫ρdz) profile can be seen. We have less massive CPD than Tanigawa et al. (2012), probably because our global disk simulation contained a planetary gap in contrary to the sheering sheet box simulations. Instead, the column density at 0.1 Hill radius in our CPD (∼100 g cm⁻² for Jupiter at 5 AU in a MMSN) is comparable to that in the radiative simulations with reduced opacity of (Ayliffe & Bate 2009) and with the most viscous simulation in D'Angelo et al. (2003).

Zoom In Zoom Out Reset image size

3.3. Flow in the Midplane of the Circumplanetary Disk

There is a debate in the literature about the direction of the radial flow on the midplane of the CPD. Ayliffe & Bate (2009) found inflow in their simulations, while Tanigawa et al. (2012), Klahr & Kley (2006), and Ayliffe & Bate (2012) found outflow.

We find that the direction of the radial flow on the midplane of the CPD depends strongly on viscosity. In our viscous simulation, the flow is outward, as shown by the streamlines plotted in Figure 9. The outflow near the midplane—together with inflow in the upper layers—is indeed typical of a 3D viscous-accretion disk (see Urpin 1984; Siemiginowska 1988; Kley & Lin 1992; Rozyczka et al. 1994; Regev & Gitelman 2002; Takeuchi & Lin 2002).

Zoom In Zoom Out Reset image size

In our nominal simulation, instead, the net flow is inward. This is due to two reasons: (1) the effective viscosity is smaller and (2) the flow suffers more pronounced shocks when crossing the spiral density wave. The latter issue is well visible in Figure 10. The shocks correspond to the points where the streamlines change abruptly direction. Look in particular at the accreting streamline on the figure. When it encounters the wave for the first time, the streamline changes abruptly direction relative to the position of the planet. The streamline now makes a hyperbolic arc around the planet. If unperturbed, it would leave the planet's sphere of influence, but it is shocked again when crossing the outer branch of the density wave at the apex of its trajectory. The shock deviates the motion once again and reduces its angular momentum relative to the planet. The streamline now makes a downward arc around the planet with a large eccentricity, but then it is shocked again and again, every half orbit. Each shock causes a loss in angular momentum, so that the streamline spirals toward the planet.

Zoom In Zoom Out Reset image size

The shocks were also visible in Figure 9, but they were less pronounced. The net flow is the result of the competition between the viscous stress, which pushes the flow outward, and the shocks, which cause angular momentum losses. In the viscous simulation, the former wins; in our nominal simulation, the latter win. The same competition should occur also for CPDs with adiabatic EOS. Shocks are weaker in that case (D'Angelo et al. 2003), but in the limit of zero viscosity they should dominate nevertheless.

It is unclear to us why Tanigawa et al. (2012) found outflow in their simulation, which had no prescribed viscosity as in our nominal case. Possibly, the numerical viscosity in their simulation was higher than in ours, or alternatively, the fact that no gap opened in their simulation changed substantially the local dynamics.

We stress, however, that the discussion about the direction of the midplane flow in the CPD is mostly academic. In fact, even in the case of inflow on the midplane, the midplane flow accounts for only 10% of the total delivered material to the CPD. We derived this percentage through the following procedure. First, we plotted the azimuthally averaged mass flux on circles at different planetocentric radii. The mass flux is increasing with decreasing distance because mass is continuously added to the CPD from the vertical direction. Considering a distance of 0.58 Hill radii from the planet, which corresponds to the largest radius at which all streamlines are accreting (see Figure 10), the mass flux is 10% of the planet's accretion rate (see "mass cut" accretion rate in Section 4 for details). Thus, the remaining 90% of the accretion has to come from the vertical direction.

After having analyzed in detail the dynamics in the vicinity of the planet, we are now ready to discuss the planet's accretion rate.

First, we checked whether we reached a stationary state at the end of the simulation by comparing the radial mass fluxes (averaged over azimuth and integrated vertically) obtained at different output times throughout the circumstellar disk and the CPD. Having concluded positively that a quasi stationary state was reached, we then checked whether the flux of mass onto the planet was a simple consequence of the mass flux toward the star in the circumstellar disk. We found, as discussed more in details in (Morbidelli et al. 2013), that the flux of gas toward the planet is due to the flux of gas into the gap from both of its sides. Thus, the planet accretion rate would not be zero even in an equilibrium disk without any net mass flux to the star. Our disk is indeed very close to an equilibrium disk for α = 0; the flux of gas toward the star in our disk is not significant and therefore not correlated to the accretion of the planet.

In order to measure the accretion rate in the simulation, we measured how much mass was removed through the mass cut. In our nominal simulation, after reaching a stationary state, we found a large accretion rate, namely $\dot{M} = 2 \times 10^{-7}\, {M_*} \Omega$. Again, the results scale linearly with Σ₀ and the relationship derived in Section 2 is $\Sigma _0=\Sigma _{{\rm MMSN}} \sqrt{a_{{\rm jup}}/a} (M_*/M_{\odot })$. Moreover, $\Omega =\sqrt{GM_*/a^3}=\sqrt{{M_*}/{M_{\odot }}({GM_{\odot }}/{(1\,\rm {AU})^3})({(1\,\rm {AU})^3}/{a^3}})$. Plugging in these will lead to $\dot{M} = 5.51 \times 10^{-7}\,{M_{\odot }}\, {\rm {yr^{-1}}} \times (M_*/M_{\odot })^{1/2} \times (\Sigma _0/\Sigma _{{\rm MMSN}}) \times (a/1\,\rm {AU})^{-1}$. If the planet is at 5.2 AU, this corresponds to ∼10⁻⁴ Jupiter masses per year. We argue that this high accretion rate is due to numerical viscosity. In fact, in the high resolution simulation, where the numerical viscosity is halved, the accretion rate is reduced by a factor of two.

Interestingly, in the high resolution simulation, the mass in the CPD is basically the same as in the nominal simulation (see Figure 8). This is because the polar inflow is also reduced by a factor of two. This is at first surprising because in 2D disks at low viscosity the width and depth of a gap is independent of viscosity (Crida et al. 2006). However, 3D gaps behave differently. The detailed analysis of the gas dynamics in a 3D gap will be the object of another paper (Morbidelli et al. 2013). However, in brief, the dynamics of a gap in a 3D disk is characterized by an interesting circulation: the gas flows into the gap from the surface of the circumstellar disk, then precipitates toward the midplane. In doing this, it falls either to the CPD or gets kicked out by the planet and goes back into the circumstellar disk. The flow into the gap at the disk's surface is dominated by the numerical viscosity and therefore changes by a factor of two from the nominal to the high resolution simulation.

This result shows that it is not possible to assess the accretion rate of a planet in the low-viscosity limit just by using simulations with no prescribed viscosity. This has to be kept in mind regardless of the EOS used in the simulations. Instead, we need to identify and quantify the accretion mechanisms that are independent of viscosity. However, as a reminder, our analysis is based on isothermal simulations; if a more realistic EOS is used, the quantitative relevance of each mechanisms may change. The analysis below, therefore, should be regarded as a proof of concept, useful also for future radiative studies and not for its quantitative results (see also Section 5).

Here we envision a scenario in which the circumstellar disk has a layered structure, with a dead zone near the midplane and an active viscous layer near the surface, in agreement with MRI studies (Gammie 1996). We envision also that the CPD is mostly MRI-inactive, in agreement with (Fujii et al. 2011; Turner et al. 2013; Fujii et al. 2013), so we investigate the planet accretion rate in the limit of vanishing viscosity in the CPD.

A first mechanism of accretion, independent of the viscosity in the CPD, is the vertical inflow. We stress that the vertical inflow is sustained by the flux of gas in the active layer of the circumstellar disk (Morbidelli et al. 2013), so that it should exist also if the planet forms in a dead zone and the CPD is inviscid. We have seen in Section 3.1 that the vertical inflow has a specific angular momentum significantly smaller than the CPD. Because the contact of the inflow and the CPD happens through a shock, the inflow subtracts angular momentum from the CPD even in the limit of zero viscosity. Nevertheless, if the specific angular momentum of the incoming gas is larger than that corresponding to an orbit at the surface of the planet, the inflow cannot promote accretion onto the planet. Therefore, the mass accreted by the planet cannot be larger than the mass carried by the inflow of gas with a specific angular momentum smaller than ($j<\sqrt{GM_pR_p}$) (Tanigawa et al. 2012). We call this subset of the vertical inflow a "polar inflow".

To estimate the accretion rate due to the polar inflow, we proceed as follows. We set the radius of the planet to be equal to twice the current radius of Jupiter. This is because the planet at the accretion time was much hotter and therefore its radius was inflated by more or less a factor of two (Guillot et al. 2004). Also, we refer to the viscous simulation. The reason is that, as we said above, the vertical infall is fed by the gas entering into the gap at the surface of the circumstellar disk, and the latter should be MRI active. With these settings, we find an accretion rate of 4 × 10⁻⁹ × M_*Ω = 11.02 × 10⁻⁹ M_☉ yr⁻¹ × (M_*/M_☉)^1/2 × (Σ₀/Σ_MMSN) × (a/1 AU)⁻¹, i.e., 2 × 10⁻⁶ Jupiter masses yr⁻¹ with the usual scalings, and it should scale linearly with α. This estimate is one order of magnitude smaller than in Tanigawa et al. (2012), presumably because of the fact that in our simulations the planet opened a gap.

The second accretion process that does not depend on the viscosity in the CPD is the loss of angular momentum in the CPD due to the torque exerted by the star through the spiral density wave (Martin & Lubow 2011; Rivier et al. 2012). This torque was already considered in Rivier et al. (2012) in their 2D simulations. The authors there assumed that in the inviscid case the torque is deposited only in the very inner part of the CPD. However, as we have seen in Section 3.3, the wave shocks and removes angular momentum also in the outer part of the CPD. The fact that the wave does not seem to shock in the inner part of the CPD is probably an artifact of numerical viscosity, which increases approaching the planet and smears out the density contrasts, consequently erasing the wave and its shock front. Because the simulation does not allow us to resolve where in the CPD the torque is deposited, in order to provide an upper bound of the planet's accretion rate promoted by the stellar torque we adopt the following simple recipe. We integrate the stellar torque from the planet to the radius where it becomes positive, which is basically at the edge of the CPD; then, we estimate the fraction of the CPD mass accreted per unit time as the fraction between the integrated stellar torque and the total angular momentum in the disk.

For both the nominal and the high resolution simulation, we derive that the stellar torque promotes the accretion of 3 × 10⁻³ of the mass of the CPD per planet's orbital period, i.e., 2.5 × 10⁻⁴ of the CPD mass per year, if the planet is at 5.2 AU around a solar-mass star. Because of the fact that this result is independent of the numerical resolution, we are confident of its robustness.

The mass of the CPD in our simulation is only 4 × 10⁻⁴ M_J. However, if the disk could not accrete onto the planet as fast as in our simulation because of the lack of viscosity, the gas would pile up into the disk, increasing the CPD mass. How massive the disk can become cannot be studied using isothermal simulations and will be the object of a future study. In Rivier et al. (2012), it was estimated analytically that the maximum mass of the CPD is ∼10⁻³ M_J; at this mass, its vertical pressure gradient becomes large enough to stop the vertical inflow, so that the mass of the CPD cannot grow further. This estimate is probably valid only at the order of magnitude level. However, even assuming a CPD mass of 0.01 M_J, the stellar torque would imply an accretion rate of only 2.5 × 10⁻⁶ M_J yr⁻¹, i.e., a mass doubling time of 400,000 yr. This timescale is comparable to that of the photoevaporation of the circumstellar disk. If this result is confirmed in future, more realistic studies, it implies that if giant planets form toward the end of the disk's lifetime, the competition between the planet's accretion timescale and the disk removal timescale might explain the wide range of masses observed for giant planets.

In this paper, the simulations were all isothermal. Previous studies showed that for small planets (∼10  M_⊕) the flow near the planet strongly depends on the EOS (Paardekooper & Mellema 2008, Ayliffe & Bate 2009, Nelson & Ruffert 2013), but for Jupiter-mass planets the accretion rate in non-isothermal simulations is close to that in isothermal calculations (Machida et al. 2010).

However, we suspect that the Machida et al. (2010) result is due to the large planet's accretion rate, which is a consequence of the numerical viscosity. At the level of detail at which we explored the local dynamics in this paper, we expect that the EOS would strongly influence the results at a quantitative level. In particular, in the limit of vanishing viscosity, the gas should pile up in the CPD, and an adiabatic EOS with flux limited energy transfer is expected to change significantly the final equilibrium structure of the CPD relative to the isothermal EOS.

The issue of the pileup of material in the CPD is crucial to estimate the planet's accretion rate in the inviscid limit. If the mass of the CPD becomes large, the stellar torque can be sufficient to promote a fast accretion onto the planet. Martin & Lubow (2012) suggested that the disk may become gravitationally unstable, which would cause FU Orionis-like accretion bursts onto the planet. In fact, works comparing the results in isothermal and adiabatic simulations, such as D'Angelo et al. (2003) or Gressel et al. (2013), show that the CPD tends to be less massive if heating and cooling effects are taken into account. However, these works are affected by a large viscosity—prescribed or numeric—which prevents the pileup of mass in the CPD. It should be investigated as to what actually happens in the ideal inviscid limit.

Nevertheless, it may be possible that the CPD becomes so hot and vertically extended that it does not allow the accretion of new material from the vertical inflow. The velocity of the vertical inflow could become sub-sonic; there would be no shock at the surface of the CPD, and the flow could be diverted by the pressure gradient. Indeed, adiabatic simulations like Ayliffe & Bate (2009) and Gressel et al. (2013) show that the boundary between the CPD and the vertical flow is less sharp than in our study, suggesting a weakening of the shock front. If the vertical inflow is diverted before that, the CPD becomes gravitationally unstable, and then a steady state equilibrium can be reached. As in Section 4, the accretion rate onto the planet will depend on the mass of this steady state CPD and the stellar torque, but the quantitative estimate will presumably be different from the one achieved in this paper. We also note that adiabatic simulations (e.g., Gressel et al. 2013) show that the spiral wave launched by the star in the CPD is much less prominent than in the isothermal case, which would reduce the stellar torque.

In the future, it will be necessary to investigate how the results change from the quantitative point of view if a more realistic EOS is used. To study the pileup in the CPD, however, one will still have the problem of numerical viscosity. If the latter promotes the accretion of material from the CPD to the planet, the final mass distribution in the CPD will not be the same as in the ideal inviscid case. Particular care will be needed to address this issue.

In this paper, we studied the dynamics of gas in the vicinity of a Jupiter-mass planet and the properties of the CPD. For this purpose we used the Jupiter code, a 3D nested-grid hydrodynamical code. We performed locally isothermal simulations with two prescribed α viscosities (α = 0.004 and α = 0) and, for α = 0, with two different resolutions.

Our results confirm those of Ayliffe & Bate (2009), Machida et al. (2010), and Tanigawa et al. (2012) concerning the vertical inflow and the CPD vertical structure. We have pointed out, however, that the CPD upper layer is wavy, i.e., the aspect ratio of the CPD changes with planetocentric azimuth, because of the pressure of the inhomogeneous vertical inflow. In a reference frame rotating with the gas around the planet (at a given radius), this pressure exerts a periodic perturbation, leading to the formation of a stationary wave in the CPD vertical structure. We also were able to reduce the viscosity more than previous local box simulations; in our inviscid simulation, the shocks were more pronounced.

We found that CPD is mostly sub-Keplerian, similar to Tanigawa et al. (2012), and Uribe et al. (2013), except in its upper layer, where it can be slightly super-Keplerian because of the significant flaring of the disk. The radial extent of the disk where the streamlines are quasi-circular depends on viscosity and, if α = 0, also on numerical resolution. The smaller the effective viscosity is, the smaller the circular portion of the disk is.

We found that the flow in the CPD midplane is inward if α = 0, in contrast with Tanigawa et al. (2012) and Ayliffe & Bate (2012). In this case, the gas flow in the CPD is crossing the spiral density wave twice in every orbit, and each crossing leads to the loss of angular momentum due to a shock. Thus, the flow spirals down to the planet. Nevertheless, we showed that the radial inflow of mass through the outer boundary of the CPD is only 10% of the gas influx coming from the vertical direction. Instead, in the case of the viscous simulation with α = 0.004, the flow is spiraling outward in the midplane. Therefore, one can conclude that the viscosity determines the directions of the flow in the CPD.

Our simulation resulted in a high planetary accretion rate, namely $\dot{M} = 1\times 10^{-4}$ Jupiter masses per year for a Jupiter-mass planet at 5.2 AU in a MMSN; however, we showed that this high rate is due to numerical viscosity. We identified that the main accretion mechanisms, independent of viscosity, is the torque exerted by the star onto the CPD. We found that the stellar torque promotes the accretion of 2.5 × 10⁻⁴ of the mass of the CPD per year, assuming a planet's orbital period of 12 yr. However, we cannot provide a reliable estimate of the mass of the CPD with our isothermal simulations, particularly in the limit of vanishing viscosity, which could lead to a significant pileup of material in the CPD. An order of magnitude analytic estimate in Rivier et al. (2012) reported a CPD mass of ∼10⁻³ M_J. Even assuming a CPD mass of 0.01 M_J, the stellar torque would lead to an accretion rate of only 2.5 × 10⁻⁶ M_J yr⁻¹. In other words, a Jupiter would build up in 400,000 yr with this accretion rate. This timescale is comparable to the removal timescale of the circumstellar disk gas (e.g., Koepferl et al. 2013; Gorti et al. 2009; Szulágyi et al. 2012).

Although future simulations implementing a realistic, non-isothermal EOS are needed to achieve a reliable quantitative estimate of the planet's accretion rate in the limit of vanishing viscosity, many conceptual results of this paper, particularly those on the role of numerical viscosity and viscosity-independent transport mechanisms in the CPD, should be valid also in a more realistic context.

The main result presented in this paper is encouraging. The similarity between planet accretion and disk removal timescales suggests that if the giant planets form toward the end of the disk's lifetime, the competition between the planet's accretion process and disk's photoevaporation could explain the observed, wide range of giant planet masses.

This work is part of the MOJO project, supported by the French ANR. We are grateful to T. Tanigawa and B. Bitsch for helpful discussions, and we thank S. Jacobson for the language corrections. J. Szulágyi acknowledges support from the Capital Fund Management's J. P. Aguilar Grant. A. Morbidelli also thanks the German Helmholtz Alliance for providing support through its Planetary Evolution and Life program. F. Masset acknowledges UNAM support through the grant PAPIIT IA101113.

Here we present a toy model for the wavy structure of the CPD. Our goal is to schematize with some basic physical considerations the CPD's reaction to the pressure of the vertical inflow, which is leading to this wavy disk surface. This toy model might help one to better understand the process without giving an exact (i.e., complex) physical model, since this is not the goal of this paper.

In the reference frame rotating with the fluid elements in the CPD, the periodic excitation by the vertical inflow's pressure creates a stationary wave in the vertical structure of the CPD. The solution for the acoustic wave equation

$$\begin{equation} \frac{1}{c_s^2}\frac{\partial ^2p}{\partial t^2}=\frac{\partial ^2p}{\partial z^2} \end{equation} \tag{ A1 }$$

for a stationary wave can be written in the following form:

$$\begin{equation} p(z,t)=2p_0 e^{i\nu t}\cos (kz), \end{equation} \tag{ A2 }$$

where t stands for the time, z represents the vertical coordinate, ν indicates the wave frequency, and k means the wavenumber. The term cos (kz) does not involve a phase because the CPD is supposed to be symmetric with respect to the midplane, so that we have ∂p/∂z = 0 at z = 0. Putting this equation back into the wave equation, we get k = ν/c_s = 2π/λ, where c_s indicates the sound speed and λ is the wavelength. Because the pressure exerted by the vertical inflow has a frequency that is twice the planetocentric orbital frequency (ν = 2ω), then at z_max the equation can be written as $p(z_{{\rm max}},t)=p_{z_{{\rm max}}} e^{i2\omega t}$. If we set this equal to Equation (A2), then we get λ = πc_s/ω = πH_CPD.

This is precisely what is seen in Figure 6. The figure shows the vertical profile of the volume density in the CPD for various values of the azimuth, at a distance of d ≈ 0.058 Hill radii from the planet, where H_CPD/d ≈ 0.1. The z coordinate is normalized by H_CPD. The profiles oscillate around the well-known Gaussian hydrostatic equilibrium profile. One can see two knots, where all curves intersect, corresponding to the locations in z where the amplitude of the wave is zero, namely corresponding to cos (kz) = 0: z = (π/4)H_CPD and z = (3π/4)H_CPD. The distance between the two knots is λ/2 = π/2H_CPD. The computation of the cumulated mass along the curves shown in Figure 6 reveals that >97% of the disk mass is below the knot at (3/4)λ = (3/4)πH_CPD ≈ 2.4H_CPD. Thus, the extreme "waviness" of the surface of the disk observed in Figure 4 concerns solely an "atmosphere" of the disk accounting only for <3% of the disk mass.